Add the orthogonal_polynomials module with the Legendre functions.

diff --git a/mjo/all.py b/mjo/all.py
index 83ecf4559dd628a53db9c2253fa16e14d325a99d..2a2e73becd6bd4e391ef981ff8ca6edfc072d88d 100644 (file)
--- a/mjo/all.py
+++ b/mjo/all.py
@@ -5,6 +5,7 @@ in his script. Instead, he can just `from mjo.all import *`.
 
 from interpolation import *
 from misc import *
+from orthogonal_polynomials import *
 from plot import *
 from symbol_sequence import *
 from symbolic import *

diff --git a/mjo/orthogonal_polynomials.py b/mjo/orthogonal_polynomials.py
new file mode 100644 (file)
index 0000000..c56d3d0
--- /dev/null
+++ b/mjo/orthogonal_polynomials.py
@@ -0,0 +1,123 @@
+from sage.all import *
+from sage.symbolic.expression import is_Expression
+
+def standard_legendre_p(n, x):
+    """
+    Returns the ``n``th Legendre polynomial of the first kind over the
+    interval [-1, 1] with respect to ``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.
+
+    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 should agree with Maxima for all `n`::
+
+        sage: eq = lambda k: bool(standard_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(x=10)
+        1
+
+    """
+    if not n in ZZ:
+        raise TypeError('n must be a natural number')
+
+    if n < 0:
+        raise ValueError('n must be nonnegative')
+
+    # Ensures that 1/(2**n) is not integer division.
+    n = ZZ(n)
+    dn = 1/(2**n)
+
+    def c(m):
+        return binomial(n,m)*binomial(n, n-m)
+
+    def g(m):
+        return ( ((x - 1)**(n-m)) * (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) ])
+
+    # If `x` is a symbolic expression, we want to return a symbolic
+    # expression (even if that expression is e.g. `1`).
+    if is_Expression(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