from sage.all import *
def legendre_p(n, x, a = -1, b = 1):
- """
+ r"""
Returns the ``n``th Legendre polynomial of the first kind over the
interval [a, b] with respect to ``x``.
returned. Otherwise, the value of the ``n``th polynomial at ``x``
will be returned.
+ SETUP::
+
+ sage: from mjo.orthogonal_polynomials import legendre_p
+
EXAMPLES:
Create the standard Legendre polynomials in `x`::
Reuse the variable from a polynomial ring::
sage: P.<t> = QQ[]
- sage: legendre_p(2,t).simplify_rational()
+ sage: legendre_p(2,t)
3/2*t^2 - 1/2
If ``x`` is a real number, the result should be as well::
[-179 242]
[-484 547]
+ And finite field elements::
+
+ sage: legendre_p(3, GF(11)(5))
+ 8
+
+ Solve a simple least squares problem over `[-\pi, \pi]`::
+
+ 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
+
TESTS:
We should agree with Maxima for all `n`::
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
+ sage: all( eq(k) for k in range(20) ) # long time
True
We can evaluate the result of the zeroth polynomial::
# same field/ring as `x`.
return x.parent()(1)
- # Even though we know a,b are real we use the symbolic ring. This
- # lets us return pretty expressions where possible.
- a = SR(a)
- b = SR(b)
- n = ZZ(n) # Ensure that 1/(2**n) is not integer division.
- dn = 1/(2**n)
+ # 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)
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)
+ # 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)
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) ])
+ # Also massaged to support finite field elements.
+ P = sum( c(m)*g(m) for m in range(n+1) )/(2**n)
return P