]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/orthogonal_polynomials.py
hurwitz.py: add conjugate() method for matrices.
[sage.d.git] / mjo / orthogonal_polynomials.py
1 from sage.all import *
2
3 def legendre_p(n, x, a = -1, b = 1):
4 r"""
5 Returns the ``n``th Legendre polynomial of the first kind over the
6 interval [a, b] with respect to ``x``.
7
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`.
12
13 INPUT:
14
15 * ``n`` -- The index of the polynomial.
16
17 * ``x`` -- Either the variable to use as the independent
18 variable in the polynomial, or a point at which to evaluate
19 the polynomial.
20
21 * ``a`` -- The "left" endpoint of the interval. Must be a real number.
22
23 * ``b`` -- The "right" endpoint of the interval. Must be a real number.
24
25 OUTPUT:
26
27 If ``x`` is a variable, a polynomial (symbolic expression) will be
28 returned. Otherwise, the value of the ``n``th polynomial at ``x``
29 will be returned.
30
31 SETUP::
32
33 sage: from mjo.orthogonal_polynomials import legendre_p
34
35 EXAMPLES:
36
37 Create the standard Legendre polynomials in `x`::
38
39 sage: legendre_p(0,x)
40 1
41 sage: legendre_p(1,x)
42 x
43
44 Reuse the variable from a polynomial ring::
45 sage: P.<t> = QQ[]
46 sage: legendre_p(2,t)
47 3/2*t^2 - 1/2
48
49 If ``x`` is a real number, the result should be as well::
50
51 sage: legendre_p(3, 1.1)
52 1.67750000000000
53
54 Similarly for complex numbers::
55
56 sage: legendre_p(3, 1 + I)
57 7/2*I - 13/2
58
59 Even matrices work::
60
61 sage: legendre_p(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
62 [-179 242]
63 [-484 547]
64
65 And finite field elements::
66
67 sage: legendre_p(3, GF(11)(5))
68 8
69
70 Solve a simple least squares problem over `[-\pi, \pi]`::
71
72 sage: a = -pi
73 sage: b = pi
74 sage: def inner_product(v1, v2):
75 ....: return integrate(v1*v2, x, a, b)
76 sage: def norm(v):
77 ....: return sqrt(inner_product(v,v))
78 sage: def project(basis, v):
79 ....: return sum( inner_product(v, b)*b/norm(b)**2
80 ....: for b in basis)
81 sage: f = sin(x)
82 sage: legendre_basis = [ legendre_p(k, x, a, b) for k in range(4) ]
83 sage: proj = project(legendre_basis, f)
84 sage: proj.simplify_trig()
85 5/2*(7*(pi^2 - 15)*x^3 - 3*(pi^4 - 21*pi^2)*x)/pi^6
86
87 TESTS:
88
89 We should agree with Maxima for all `n`::
90
91 sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x))
92 sage: all( eq(k) for k in range(20) ) # long time
93 True
94
95 We can evaluate the result of the zeroth polynomial::
96
97 sage: f = legendre_p(0,x)
98 sage: f(x=10)
99 1
100
101 We should have |P(a)| = |P(b)| = 1 for all a,b::
102
103 sage: a = RR.random_element()
104 sage: b = RR.random_element()
105 sage: k = ZZ.random_element(20)
106 sage: P = legendre_p(k, x, a, b)
107 sage: abs(P(x=a)) # abs tol 1e-12
108 1
109 sage: abs(P(x=b)) # abs tol 1e-12
110 1
111
112 Two different polynomials should be orthogonal with respect to the
113 inner product over `[a,b]`. Note that this test can fail if QQ is
114 replaced with RR, since integrate() can return NaN::
115
116 sage: a = QQ.random_element()
117 sage: b = QQ.random_element()
118 sage: j = ZZ.random_element(20)
119 sage: k = j + 1
120 sage: Pj = legendre_p(j, x, a, b)
121 sage: Pk = legendre_p(k, x, a, b)
122 sage: integrate(Pj*Pk, x, a, b) # abs tol 1e-12
123 0
124
125 The first few polynomials shifted to [0,1] are known to be::
126
127 sage: p0 = 1
128 sage: p1 = 2*x - 1
129 sage: p2 = 6*x^2 - 6*x + 1
130 sage: p3 = 20*x^3 - 30*x^2 + 12*x - 1
131 sage: bool(legendre_p(0, x, 0, 1) == p0)
132 True
133 sage: bool(legendre_p(1, x, 0, 1) == p1)
134 True
135 sage: bool(legendre_p(2, x, 0, 1) == p2)
136 True
137 sage: bool(legendre_p(3, x, 0, 1) == p3)
138 True
139
140 The zeroth polynomial is an element of the ring that we're working
141 in::
142
143 sage: legendre_p(0, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
144 [1 0]
145 [0 1]
146
147 """
148 if not n in ZZ:
149 raise TypeError('n must be a natural number')
150
151 if n < 0:
152 raise ValueError('n must be nonnegative')
153
154 if not (a in RR and b in RR):
155 raise TypeError('both `a` and `b` must be real numbers')
156
157 if n == 0:
158 # Easy base case, save time. Attempt to return a value in the
159 # same field/ring as `x`.
160 return x.parent()(1)
161
162 # Even though we know a,b are real we use a different ring. We
163 # prefer ZZ so that we can support division of finite field
164 # elements by (a-b). Eventually this should be supported for QQ as
165 # well, although it does not work at the moment. The preference of
166 # SR over RR is to return something attractive when e.g. a=pi.
167 if a in ZZ:
168 a = ZZ(a)
169 else:
170 a = SR(a)
171
172 if b in ZZ:
173 b = ZZ(b)
174 else:
175 b = SR(b)
176
177 # Ensure that (2**n) is an element of ZZ. This is used later --
178 # we can divide finite field elements by integers but we can't
179 # multiply them by rationals at the moment.
180 n = ZZ(n)
181
182 def phi(t):
183 # This is an affine map from [a,b] into [-1,1] and so
184 # preserves orthogonality.
185 return (a + b - 2*t)/(a - b)
186
187 def c(m):
188 return binomial(n,m)*binomial(n, n-m)
189
190 def g(m):
191 # As given in A&S, but with `x` replaced by `phi(x)`.
192 return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m )
193
194 # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
195 # Also massaged to support finite field elements.
196 P = sum( c(m)*g(m) for m in range(n+1) )/(2**n)
197
198 return P