X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Forthogonal_polynomials.py;h=15c695611fb3e287c357075f45e55eeddf550f7f;hb=a339e89c225bb46379332ecb8b0c50b918d34ac6;hp=c56d3d081af98f9459d4b15650b13f9777d571ad;hpb=351c057951b2d059d4faa754caf63bd2486638c6;p=sage.d.git diff --git a/mjo/orthogonal_polynomials.py b/mjo/orthogonal_polynomials.py index c56d3d0..15c6956 100644 --- a/mjo/orthogonal_polynomials.py +++ b/mjo/orthogonal_polynomials.py @@ -1,10 +1,14 @@ 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): + r""" 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,110 +18,181 @@ 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 returned. Otherwise, the value of the ``n``th polynomial at ``x`` will be returned. - TESTS: + SETUP:: - We should agree with Maxima for all `n`:: + sage: from mjo.orthogonal_polynomials import legendre_p - 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 + EXAMPLES: - We can evaluate the result of the zeroth polynomial:: + Create the standard Legendre polynomials in `x`:: - sage: f = standard_legendre_p(0,x) - sage: f(x=10) + sage: legendre_p(0,x) 1 + sage: legendre_p(1,x) + x - """ - 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) + Reuse the variable from a polynomial ring:: + sage: P. = QQ[] + sage: legendre_p(2,t) + 3/2*t^2 - 1/2 - def c(m): - return binomial(n,m)*binomial(n, n-m) + If ``x`` is a real number, the result should be as well:: - def g(m): - return ( ((x - 1)**(n-m)) * (x + 1)**m ) + sage: legendre_p(3, 1.1) + 1.67750000000000 - # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0. - P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ]) + Similarly for complex numbers:: - # 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) + sage: legendre_p(3, 1 + I) + 7/2*I - 13/2 - return P + Even matrices work:: + sage: legendre_p(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7])) + [-179 242] + [-484 547] -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``. + And finite field elements:: - INPUT: + sage: legendre_p(3, GF(11)(5)) + 8 - * ``n`` -- The index of the polynomial. + Solve a simple least squares problem over `[-\pi, \pi]`:: - * ``x`` -- Either the variable to use as the independent - variable in the polynomial, or a point at which to evaluate - the polynomial. + sage: a = -pi + sage: b = pi + sage: def inner_product(v1, v2): + ....: return integrate(v1*v2, x, a, b) + sage: def norm(v): + ....: return sqrt(inner_product(v,v)) + sage: def project(basis, v): + ....: return sum( inner_product(v, b)*b/norm(b)**2 + ....: for b in basis) + sage: f = sin(x) + sage: legendre_basis = [ legendre_p(k, x, a, b) for k in range(4) ] + sage: proj = project(legendre_basis, f) + sage: proj.simplify_trig() + 5/2*(7*(pi^2 - 15)*x^3 - 3*(pi^4 - 21*pi^2)*x)/pi^6 - * ``a`` -- The "left" endpoint of the interval. Must be a real number. + TESTS: - * ``b`` -- The "right" endpoint of the interval. Must be a real number. + We should agree with Maxima for all `n`:: - OUTPUT: + sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x)) + sage: all( eq(k) for k in range(20) ) # long time + True - 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. + We can evaluate the result of the zeroth polynomial:: - TESTS: + sage: f = legendre_p(0,x) + sage: f(x=10) + 1 - We agree with ``standard_legendre_p`` when `a=b=1`:: + We should have |P(a)| = |P(b)| = 1 for all a,b:: - 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 + 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 - We should have |P(a)| = |P(b)| = 1 for all a,b. + Two different polynomials should be orthogonal with respect to the + inner product over `[a,b]`. Note that this test can fail if QQ is + replaced with RR, since integrate() can return NaN:: 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) + 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 + + The first few polynomials shifted to [0,1] are known to be:: + + sage: p0 = 1 + sage: p1 = 2*x - 1 + sage: p2 = 6*x^2 - 6*x + 1 + sage: p3 = 20*x^3 - 30*x^2 + 12*x - 1 + sage: bool(legendre_p(0, x, 0, 1) == p0) + True + sage: bool(legendre_p(1, x, 0, 1) == p1) True - sage: bool(abs(P(x=b)) == 1) + sage: bool(legendre_p(2, x, 0, 1) == p2) True + sage: bool(legendre_p(3, x, 0, 1) == p3) + True + + The zeroth polynomial is an element of the ring that we're working + in:: + + sage: legendre_p(0, MatrixSpace(ZZ, 2)([1, 2, -4, 7])) + [1 0] + [0 1] """ + if not n in ZZ: + raise TypeError('n must be a natural number') + + if n < 0: + raise ValueError('n must be nonnegative') + if not (a in RR and b in RR): - raise TypeError('both `a` and `b` must be a real numbers') + raise TypeError('both `a` and `b` must be real numbers') + + if n == 0: + # Easy base case, save time. Attempt to return a value in the + # same field/ring as `x`. + return x.parent()(1) + + # Even though we know a,b are real we use a different ring. We + # prefer ZZ so that we can support division of finite field + # elements by (a-b). Eventually this should be supported for QQ as + # well, although it does not work at the moment. The preference of + # SR over RR is to return something attractive when e.g. a=pi. + if a in ZZ: + a = ZZ(a) + else: + a = SR(a) + + if b in ZZ: + b = ZZ(b) + else: + b = SR(b) + + # Ensure that (2**n) is an element of ZZ. This is used later -- + # we can divide finite field elements by integers but we can't + # multiply them by rationals at the moment. + n = ZZ(n) - a = RR(a) - b = RR(b) - t = SR.symbol('t') - P = standard_legendre_p(n, t) + def phi(t): + # This is an affine map from [a,b] into [-1,1] and so + # preserves orthogonality. + return (a + b - 2*t)/(a - b) + + def c(m): + return binomial(n,m)*binomial(n, n-m) - # 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) + def g(m): + # As given in A&S, but with `x` replaced by `phi(x)`. + return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m ) - P_tilde = P(t = phi(x)) + # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0. + # Also massaged to support finite field elements. + P = sum( c(m)*g(m) for m in range(n+1) )/(2**n) - return P_tilde + return P