]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/orthogonal_polynomials.py
3 def legendre_p(n
, x
, a
= -1, b
= 1):
5 Returns the ``n``th Legendre polynomial of the first kind over the
6 interval [a, b] with respect to ``x``.
8 When [a,b] is not [-1,1], we scale the standard Legendre
9 polynomial (which is defined over [-1,1]) via an affine map. The
10 resulting polynomials are still orthogonal and possess the
11 property that `P(a) = P(b) = 1`.
15 * ``n`` -- The index of the polynomial.
17 * ``x`` -- Either the variable to use as the independent
18 variable in the polynomial, or a point at which to evaluate
21 * ``a`` -- The "left" endpoint of the interval. Must be a real number.
23 * ``b`` -- The "right" endpoint of the interval. Must be a real number.
27 If ``x`` is a variable, a polynomial (symbolic expression) will be
28 returned. Otherwise, the value of the ``n``th polynomial at ``x``
33 sage: from mjo.orthogonal_polynomials import legendre_p
37 Create the standard Legendre polynomials in `x`::
44 Reuse the variable from a polynomial ring::
49 If ``x`` is a real number, the result should be as well::
51 sage: legendre_p(3, 1.1)
54 Similarly for complex numbers::
56 sage: legendre_p(3, 1 + I)
61 sage: legendre_p(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
65 And finite field elements::
67 sage: legendre_p(3, GF(11)(5))
70 Solve a simple least squares problem over `[-\pi, \pi]`::
74 sage: def inner_product(v1, v2):
75 ... return integrate(v1*v2, x, a, b)
78 ... return sqrt(inner_product(v,v))
80 sage: def project(basis, v):
81 ... return sum([ inner_product(v, b)*b/norm(b)**2
85 sage: legendre_basis = [ legendre_p(k, x, a, b) for k in range(0,4) ]
86 sage: proj = project(legendre_basis, f)
87 sage: proj.simplify_trig()
88 5/2*(7*(pi^2 - 15)*x^3 - 3*(pi^4 - 21*pi^2)*x)/pi^6
92 We should agree with Maxima for all `n`::
94 sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x))
95 sage: all([eq(k) for k in range(0,20) ]) # long time
98 We can evaluate the result of the zeroth polynomial::
100 sage: f = legendre_p(0,x)
104 We should have |P(a)| = |P(b)| = 1 for all a,b::
106 sage: a = RR.random_element()
107 sage: b = RR.random_element()
108 sage: k = ZZ.random_element(20)
109 sage: P = legendre_p(k, x, a, b)
110 sage: abs(P(x=a)) # abs tol 1e-12
112 sage: abs(P(x=b)) # abs tol 1e-12
115 Two different polynomials should be orthogonal with respect to the
116 inner product over `[a,b]`. Note that this test can fail if QQ is
117 replaced with RR, since integrate() can return NaN::
119 sage: a = QQ.random_element()
120 sage: b = QQ.random_element()
121 sage: j = ZZ.random_element(20)
123 sage: Pj = legendre_p(j, x, a, b)
124 sage: Pk = legendre_p(k, x, a, b)
125 sage: integrate(Pj*Pk, x, a, b) # abs tol 1e-12
128 The first few polynomials shifted to [0,1] are known to be::
132 sage: p2 = 6*x^2 - 6*x + 1
133 sage: p3 = 20*x^3 - 30*x^2 + 12*x - 1
134 sage: bool(legendre_p(0, x, 0, 1) == p0)
136 sage: bool(legendre_p(1, x, 0, 1) == p1)
138 sage: bool(legendre_p(2, x, 0, 1) == p2)
140 sage: bool(legendre_p(3, x, 0, 1) == p3)
143 The zeroth polynomial is an element of the ring that we're working
146 sage: legendre_p(0, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
152 raise TypeError('n must be a natural number')
155 raise ValueError('n must be nonnegative')
157 if not (a
in RR
and b
in RR
):
158 raise TypeError('both `a` and `b` must be real numbers')
161 # Easy base case, save time. Attempt to return a value in the
162 # same field/ring as `x`.
165 # Even though we know a,b are real we use a different ring. We
166 # prefer ZZ so that we can support division of finite field
167 # elements by (a-b). Eventually this should be supported for QQ as
168 # well, although it does not work at the moment. The preference of
169 # SR over RR is to return something attractive when e.g. a=pi.
180 # Ensure that (2**n) is an element of ZZ. This is used later --
181 # we can divide finite field elements by integers but we can't
182 # multiply them by rationals at the moment.
186 # This is an affine map from [a,b] into [-1,1] and so
187 # preserves orthogonality.
188 return (a
+ b
- 2*t
)/(a
- b
)
191 return binomial(n
,m
)*binomial(n
, n
-m
)
194 # As given in A&S, but with `x` replaced by `phi(x)`.
195 return ( ((phi(x
) - 1)**(n
-m
)) * (phi(x
) + 1)**m
)
197 # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
198 # Also massaged to support finite field elements.
199 P
= sum([ c(m
)*g(m
) for m
in range(0,n
+1) ])/(2**n
)