From 351c057951b2d059d4faa754caf63bd2486638c6 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sat, 17 Nov 2012 15:53:14 -0500 Subject: [PATCH] Add the orthogonal_polynomials module with the Legendre functions. --- mjo/all.py | 1 + mjo/orthogonal_polynomials.py | 123 ++++++++++++++++++++++++++++++++++ 2 files changed, 124 insertions(+) create mode 100644 mjo/orthogonal_polynomials.py diff --git a/mjo/all.py b/mjo/all.py index 83ecf45..2a2e73b 100644 --- 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 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 -- 2.44.2