mjo/orthogonal_polynomials.py

   1 from sage.all import *
   2 from sage.symbolic.expression import is_Expression
   3
   4 def legendre_p(n, x, a = -1, b = 1):
   5     """
   6     Returns the ``n``th Legendre polynomial of the first kind over the
   7     interval [a, b] with respect to ``x``.
   8
   9     When [a,b] is not [-1,1], we scale the standard Legendre
  10     polynomial (which is defined over [-1,1]) via an affine map. The
  11     resulting polynomials are still orthogonal and possess the
  12     property that `P(a) = P(b) = 1`.
  13
  14     INPUT:
  15
  16       * ``n`` -- The index of the polynomial.
  17
  18       * ``x`` -- Either the variable to use as the independent
  19         variable in the polynomial, or a point at which to evaluate
  20         the polynomial.
  21
  22       * ``a`` -- The "left" endpoint of the interval. Must be a real number.
  23
  24       * ``b`` -- The "right" endpoint of the interval. Must be a real number.
  25
  26     OUTPUT:
  27
  28     If ``x`` is a variable, a polynomial (symbolic expression) will be
  29     returned. Otherwise, the value of the ``n``th polynomial at ``x``
  30     will be returned.
  31
  32     TESTS:
  33
  34     We should agree with Maxima for all `n`::
  35
  36         sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x))
  37         sage: all([eq(k) for k in range(0,20) ]) # long time
  38         True
  39
  40     We can evaluate the result of the zeroth polynomial::
  41
  42         sage: f = legendre_p(0,x)
  43         sage: f(x=10)
  44         1
  45
  46     We should have |P(a)| = |P(b)| = 1 for all a,b::
  47
  48         sage: a = RR.random_element()
  49         sage: b = RR.random_element()
  50         sage: k = ZZ.random_element(20)
  51         sage: P = legendre_p(k, x, a, b)
  52         sage: abs(P(x=a)) # abs tol 1e-12
  53         1
  54         sage: abs(P(x=b)) # abs tol 1e-12
  55         1
  56
  57     Two different polynomials should be orthogonal with respect to the
  58     inner product over [a,b]::
  59
  60         sage: a = RR.random_element()
  61         sage: b = RR.random_element()
  62         sage: j = ZZ.random_element(20)
  63         sage: k = j + 1
  64         sage: Pj = legendre_p(j, x, a, b)
  65         sage: Pk = legendre_p(k, x, a, b)
  66         sage: integrate(Pj*Pk, x, a, b) # abs tol 1e-12
  67         0
  68
  69     The first few polynomials shifted to [0,1] are known to be::
  70
  71         sage: p0 = 1
  72         sage: p1 = 2*x - 1
  73         sage: p2 = 6*x^2 - 6*x + 1
  74         sage: p3 = 20*x^3 - 30*x^2 + 12*x - 1
  75         sage: bool(legendre_p(0, x, 0, 1) == p0)
  76         True
  77         sage: bool(legendre_p(1, x, 0, 1) == p1)
  78         True
  79         sage: bool(legendre_p(2, x, 0, 1) == p2)
  80         True
  81         sage: bool(legendre_p(3, x, 0, 1) == p3)
  82         True
  83
  84     """
  85     if not n in ZZ:
  86         raise TypeError('n must be a natural number')
  87
  88     if n < 0:
  89         raise ValueError('n must be nonnegative')
  90
  91     if not (a in RR and b in RR):
  92         raise TypeError('both `a` and `b` must be real numbers')
  93
  94     a = RR(a)
  95     b = RR(b)
  96     n = ZZ(n) # Ensure that 1/(2**n) is not integer division.
  97     dn = 1/(2**n)
  98
  99     def phi(t):
 100         # This is an affine map from [a,b] into [-1,1] and so preserves
 101         # orthogonality.
 102         return (2 / (b-a))*t + 1 - (2*b)/(b-a)
 103
 104     def c(m):
 105         return binomial(n,m)*binomial(n, n-m)
 106
 107     def g(m):
 108         # As given in A&S, but with `x` replaced by `phi(x)`.
 109         return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m )
 110
 111     # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
 112     P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ])
 113
 114     # If `x` is a symbolic expression, we want to return a symbolic
 115     # expression (even if that expression is e.g. `1`).
 116     if is_Expression(x):
 117         P = SR(P)
 118
 119     return P