mjo/orthogonal_polynomials.py

   1 from sage.all import *
   2
   3 def legendre_p(n, x, a = -1, b = 1):
   4     """
   5     Returns the ``n``th Legendre polynomial of the first kind over the
   6     interval [a, b] with respect to ``x``.
   7
   8     When [a,b] is not [-1,1], we scale the standard Legendre
   9     polynomial (which is defined over [-1,1]) via an affine map. The
  10     resulting polynomials are still orthogonal and possess the
  11     property that `P(a) = P(b) = 1`.
  12
  13     INPUT:
  14
  15       * ``n`` -- The index of the polynomial.
  16
  17       * ``x`` -- Either the variable to use as the independent
  18         variable in the polynomial, or a point at which to evaluate
  19         the polynomial.
  20
  21       * ``a`` -- The "left" endpoint of the interval. Must be a real number.
  22
  23       * ``b`` -- The "right" endpoint of the interval. Must be a real number.
  24
  25     OUTPUT:
  26
  27     If ``x`` is a variable, a polynomial (symbolic expression) will be
  28     returned. Otherwise, the value of the ``n``th polynomial at ``x``
  29     will be returned.
  30
  31     EXAMPLES:
  32
  33     Create the standard Legendre polynomials in `x`::
  34
  35         sage: legendre_p(0,x)
  36         1
  37         sage: legendre_p(1,x)
  38         x
  39
  40     Reuse the variable from a polynomial ring::
  41         sage: P.<t> = QQ[]
  42         sage: legendre_p(2,t).simplify_rational()
  43         3/2*t^2 - 1/2
  44
  45     If ``x`` is a real number, the result should be as well::
  46
  47         sage: legendre_p(3, 1.1)
  48         1.67750000000000
  49
  50     Similarly for complex numbers::
  51
  52         sage: legendre_p(3, 1 + I)
  53         7/2*I - 13/2
  54
  55     Even matrices work::
  56
  57         sage: legendre_p(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
  58         [-179  242]
  59         [-484  547]
  60
  61     TESTS:
  62
  63     We should agree with Maxima for all `n`::
  64
  65         sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x))
  66         sage: all([eq(k) for k in range(0,20) ]) # long time
  67         True
  68
  69     We can evaluate the result of the zeroth polynomial::
  70
  71         sage: f = legendre_p(0,x)
  72         sage: f(x=10)
  73         1
  74
  75     We should have |P(a)| = |P(b)| = 1 for all a,b::
  76
  77         sage: a = RR.random_element()
  78         sage: b = RR.random_element()
  79         sage: k = ZZ.random_element(20)
  80         sage: P = legendre_p(k, x, a, b)
  81         sage: abs(P(x=a)) # abs tol 1e-12
  82         1
  83         sage: abs(P(x=b)) # abs tol 1e-12
  84         1
  85
  86     Two different polynomials should be orthogonal with respect to the
  87     inner product over `[a,b]`. Note that this test can fail if QQ is
  88     replaced with RR, since integrate() can return NaN::
  89
  90         sage: a = QQ.random_element()
  91         sage: b = QQ.random_element()
  92         sage: j = ZZ.random_element(20)
  93         sage: k = j + 1
  94         sage: Pj = legendre_p(j, x, a, b)
  95         sage: Pk = legendre_p(k, x, a, b)
  96         sage: integrate(Pj*Pk, x, a, b) # abs tol 1e-12
  97         0
  98
  99     The first few polynomials shifted to [0,1] are known to be::
 100
 101         sage: p0 = 1
 102         sage: p1 = 2*x - 1
 103         sage: p2 = 6*x^2 - 6*x + 1
 104         sage: p3 = 20*x^3 - 30*x^2 + 12*x - 1
 105         sage: bool(legendre_p(0, x, 0, 1) == p0)
 106         True
 107         sage: bool(legendre_p(1, x, 0, 1) == p1)
 108         True
 109         sage: bool(legendre_p(2, x, 0, 1) == p2)
 110         True
 111         sage: bool(legendre_p(3, x, 0, 1) == p3)
 112         True
 113
 114     The zeroth polynomial is an element of the ring that we're working
 115     in::
 116
 117         sage: legendre_p(0, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
 118         [1 0]
 119         [0 1]
 120
 121     """
 122     if not n in ZZ:
 123         raise TypeError('n must be a natural number')
 124
 125     if n < 0:
 126         raise ValueError('n must be nonnegative')
 127
 128     if not (a in RR and b in RR):
 129         raise TypeError('both `a` and `b` must be real numbers')
 130
 131     if n == 0:
 132         # Easy base case, save time. Attempt to return a value in the
 133         # same field/ring as `x`.
 134         return x.parent()(1)
 135
 136     # Even though we know a,b are real we use the symbolic ring. This
 137     # lets us return pretty expressions where possible.
 138     a = SR(a)
 139     b = SR(b)
 140     n = ZZ(n) # Ensure that 1/(2**n) is not integer division.
 141     dn = 1/(2**n)
 142
 143     def phi(t):
 144         # This is an affine map from [a,b] into [-1,1] and so preserves
 145         # orthogonality.
 146         return (2 / (b-a))*t + 1 - (2*b)/(b-a)
 147
 148     def c(m):
 149         return binomial(n,m)*binomial(n, n-m)
 150
 151     def g(m):
 152         # As given in A&S, but with `x` replaced by `phi(x)`.
 153         return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m )
 154
 155     # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
 156     P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ])
 157
 158     return P