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:
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.<t> = 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
projects / sage.d.git / blobdiff