Combine the standard/general Legendre polynomial functions.

diff --git a/mjo/orthogonal_polynomials.py b/mjo/orthogonal_polynomials.py
index c56d3d081af98f9459d4b15650b13f9777d571ad..3f6b227fb5b660f41de61f10ad2f2edb9542f1e8 100644 (file)
--- a/mjo/orthogonal_polynomials.py
+++ b/mjo/orthogonal_polynomials.py
@@ -1,10 +1,15 @@
 from sage.all import *
 from sage.symbolic.expression import is_Expression
 
-def standard_legendre_p(n, x):
+def legendre_p(n, x, a = -1, b = 1):
     """
     Returns the ``n``th Legendre polynomial of the first kind over the
-    interval [-1, 1] with respect to ``x``.
+    interval [a, b] with respect to ``x``.
+
+    When [a,b] is not [-1,1], we scale the standard Legendre
+    polynomial (which is defined over [-1,1]) via an affine map. The
+    resulting polynomials are still orthogonal and possess the
+    property that `P(a) = P(b) = 1`.
 
     INPUT:
 
@@ -14,6 +19,10 @@ def standard_legendre_p(n, x):
         variable in the polynomial, or a point at which to evaluate
         the polynomial.
 
+      * ``a`` -- The "left" endpoint of the interval. Must be a real number.
+
+      * ``b`` -- The "right" endpoint of the interval. Must be a real number.
+
     OUTPUT:
 
     If ``x`` is a variable, a polynomial (symbolic expression) will be
@@ -24,16 +33,39 @@ def standard_legendre_p(n, x):
 
     We should agree with Maxima for all `n`::
 
-        sage: eq = lambda k: bool(standard_legendre_p(k,x) == legendre_P(k,x))
+        sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x))
         sage: all([eq(k) for k in range(0,20) ]) # long time
         True
 
     We can evaluate the result of the zeroth polynomial::
 
-        sage: f = standard_legendre_p(0,x)
+        sage: f = legendre_p(0,x)
         sage: f(x=10)
         1
 
+    We should have |P(a)| = |P(b)| = 1 for all a,b::
+
+        sage: a = RR.random_element()
+        sage: b = RR.random_element()
+        sage: k = ZZ.random_element(20)
+        sage: P = legendre_p(k, x, a, b)
+        sage: abs(P(x=a)) # abs tol 1e-12
+        1
+        sage: abs(P(x=b)) # abs tol 1e-12
+        1
+
+    Two different polynomials should be orthogonal with respect to the
+    inner product over [a,b]::
+
+        sage: a = RR.random_element()
+        sage: b = RR.random_element()
+        sage: j = ZZ.random_element(20)
+        sage: k = j + 1
+        sage: Pj = legendre_p(j, x, a, b)
+        sage: Pk = legendre_p(k, x, a, b)
+        sage: integrate(Pj*Pk, x, a, b) # abs tol 1e-12
+        0
+
     """
     if not n in ZZ:
         raise TypeError('n must be a natural number')
@@ -41,15 +73,25 @@ def standard_legendre_p(n, x):
     if n < 0:
         raise ValueError('n must be nonnegative')
 
-    # Ensures that 1/(2**n) is not integer division.
-    n = ZZ(n)
+    if not (a in RR and b in RR):
+        raise TypeError('both `a` and `b` must be real numbers')
+
+    a = RR(a)
+    b = RR(b)
+    n = ZZ(n) # Ensure that 1/(2**n) is not integer division.
     dn = 1/(2**n)
 
+    def phi(t):
+        # This is an affine map from [a,b] into [-1,1] and so preserves
+        # orthogonality.
+        return (2 / (b-a))*t + 1 - (2*b)/(b-a)
+
     def c(m):
         return binomial(n,m)*binomial(n, n-m)
 
     def g(m):
-        return ( ((x - 1)**(n-m)) * (x + 1)**m )
+        # As given in A&S, but with `x` replaced by `phi(x)`.
+        return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m )
 
     # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
     P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ])
@@ -60,64 +102,3 @@ def standard_legendre_p(n, x):
         P = SR(P)
 
     return P
-
-
-def legendre_p(n, x, a=-1, b=1):
-    """
-    Return the ``n``th Legendre polynomial over the interval `[a,b]`
-    with respect to the variable ``x``.
-
-    INPUT:
-
-      * ``n`` -- The index of the polynomial.
-
-      * ``x`` -- Either the variable to use as the independent
-        variable in the polynomial, or a point at which to evaluate
-        the polynomial.
-
-      * ``a`` -- The "left" endpoint of the interval. Must be a real number.
-
-      * ``b`` -- The "right" endpoint of the interval. Must be a real number.
-
-    OUTPUT:
-
-    If ``x`` is a variable, a polynomial (symbolic expression) will be
-    returned. Otherwise, the value of the ``n``th polynomial at ``x``
-    will be returned.
-
-    TESTS:
-
-    We agree with ``standard_legendre_p`` when `a=b=1`::
-
-        sage: eq = lambda k: bool(legendre_p(k,x) == standard_legendre_p(k,x))
-        sage: all([ eq(k) for k in range(0, 20) ]) # long time
-        True
-
-    We should have |P(a)| = |P(b)| = 1 for all a,b.
-
-        sage: a = QQ.random_element()
-        sage: b = QQ.random_element()
-        sage: k = ZZ.random_element(20)
-        sage: P = legendre_p(k, x, a, b)
-        sage: bool(abs(P(x=a)) == 1)
-        True
-        sage: bool(abs(P(x=b)) == 1)
-        True
-
-    """
-    if not (a in RR and b in RR):
-        raise TypeError('both `a` and `b` must be a real numbers')
-
-    a = RR(a)
-    b = RR(b)
-    t = SR.symbol('t')
-    P = standard_legendre_p(n, t)
-
-    # This is an affine map from [a,b] into [-1,1] and so preserves
-    # orthogonality. If we define this with 'def' instead of a lambda,
-    # Python segfaults as we evaluate P.
-    phi = lambda y: (2 / (b-a))*y + 1 - (2*b)/(b-a)
-
-    P_tilde = P(t = phi(x))
-
-    return P_tilde