mjo/orthogonal_polynomials.py

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