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     """
  70     if not n in ZZ:
  71         raise TypeError('n must be a natural number')
  72
  73     if n < 0:
  74         raise ValueError('n must be nonnegative')
  75
  76     if not (a in RR and b in RR):
  77         raise TypeError('both `a` and `b` must be real numbers')
  78
  79     a = RR(a)
  80     b = RR(b)
  81     n = ZZ(n) # Ensure that 1/(2**n) is not integer division.
  82     dn = 1/(2**n)
  83
  84     def phi(t):
  85         # This is an affine map from [a,b] into [-1,1] and so preserves
  86         # orthogonality.
  87         return (2 / (b-a))*t + 1 - (2*b)/(b-a)
  88
  89     def c(m):
  90         return binomial(n,m)*binomial(n, n-m)
  91
  92     def g(m):
  93         # As given in A&S, but with `x` replaced by `phi(x)`.
  94         return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m )
  95
  96     # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
  97     P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ])
  98
  99     # If `x` is a symbolic expression, we want to return a symbolic
 100     # expression (even if that expression is e.g. `1`).
 101     if is_Expression(x):
 102         P = SR(P)
 103
 104     return P