From cf113d36bb50229eeb063e50a41d26f7201511ed Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sat, 17 Nov 2012 20:15:26 -0500 Subject: [PATCH] Combine the standard/general Legendre polynomial functions. --- mjo/orthogonal_polynomials.py | 117 ++++++++++++++-------------------- 1 file changed, 49 insertions(+), 68 deletions(-) diff --git a/mjo/orthogonal_polynomials.py b/mjo/orthogonal_polynomials.py index c56d3d0..3f6b227 100644 --- 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 -- 2.43.2