Mercurial > hg > vamp-build-and-test
comparison DEPENDENCIES/mingw32/Python27/Lib/site-packages/numpy/polynomial/legendre.py @ 87:2a2c65a20a8b
Add Python libs and headers
author | Chris Cannam |
---|---|
date | Wed, 25 Feb 2015 14:05:22 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
86:413a9d26189e | 87:2a2c65a20a8b |
---|---|
1 """ | |
2 Legendre Series (:mod: `numpy.polynomial.legendre`) | |
3 =================================================== | |
4 | |
5 .. currentmodule:: numpy.polynomial.polynomial | |
6 | |
7 This module provides a number of objects (mostly functions) useful for | |
8 dealing with Legendre series, including a `Legendre` class that | |
9 encapsulates the usual arithmetic operations. (General information | |
10 on how this module represents and works with such polynomials is in the | |
11 docstring for its "parent" sub-package, `numpy.polynomial`). | |
12 | |
13 Constants | |
14 --------- | |
15 | |
16 .. autosummary:: | |
17 :toctree: generated/ | |
18 | |
19 legdomain Legendre series default domain, [-1,1]. | |
20 legzero Legendre series that evaluates identically to 0. | |
21 legone Legendre series that evaluates identically to 1. | |
22 legx Legendre series for the identity map, ``f(x) = x``. | |
23 | |
24 Arithmetic | |
25 ---------- | |
26 | |
27 .. autosummary:: | |
28 :toctree: generated/ | |
29 | |
30 legmulx multiply a Legendre series in P_i(x) by x. | |
31 legadd add two Legendre series. | |
32 legsub subtract one Legendre series from another. | |
33 legmul multiply two Legendre series. | |
34 legdiv divide one Legendre series by another. | |
35 legpow raise a Legendre series to an positive integer power | |
36 legval evaluate a Legendre series at given points. | |
37 legval2d evaluate a 2D Legendre series at given points. | |
38 legval3d evaluate a 3D Legendre series at given points. | |
39 leggrid2d evaluate a 2D Legendre series on a Cartesian product. | |
40 leggrid3d evaluate a 3D Legendre series on a Cartesian product. | |
41 | |
42 Calculus | |
43 -------- | |
44 | |
45 .. autosummary:: | |
46 :toctree: generated/ | |
47 | |
48 legder differentiate a Legendre series. | |
49 legint integrate a Legendre series. | |
50 | |
51 Misc Functions | |
52 -------------- | |
53 | |
54 .. autosummary:: | |
55 :toctree: generated/ | |
56 | |
57 legfromroots create a Legendre series with specified roots. | |
58 legroots find the roots of a Legendre series. | |
59 legvander Vandermonde-like matrix for Legendre polynomials. | |
60 legvander2d Vandermonde-like matrix for 2D power series. | |
61 legvander3d Vandermonde-like matrix for 3D power series. | |
62 leggauss Gauss-Legendre quadrature, points and weights. | |
63 legweight Legendre weight function. | |
64 legcompanion symmetrized companion matrix in Legendre form. | |
65 legfit least-squares fit returning a Legendre series. | |
66 legtrim trim leading coefficients from a Legendre series. | |
67 legline Legendre series representing given straight line. | |
68 leg2poly convert a Legendre series to a polynomial. | |
69 poly2leg convert a polynomial to a Legendre series. | |
70 | |
71 Classes | |
72 ------- | |
73 Legendre A Legendre series class. | |
74 | |
75 See also | |
76 -------- | |
77 numpy.polynomial.polynomial | |
78 numpy.polynomial.chebyshev | |
79 numpy.polynomial.laguerre | |
80 numpy.polynomial.hermite | |
81 numpy.polynomial.hermite_e | |
82 | |
83 """ | |
84 from __future__ import division, absolute_import, print_function | |
85 | |
86 import warnings | |
87 import numpy as np | |
88 import numpy.linalg as la | |
89 | |
90 from . import polyutils as pu | |
91 from ._polybase import ABCPolyBase | |
92 | |
93 __all__ = [ | |
94 'legzero', 'legone', 'legx', 'legdomain', 'legline', 'legadd', | |
95 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow', 'legval', 'legder', | |
96 'legint', 'leg2poly', 'poly2leg', 'legfromroots', 'legvander', | |
97 'legfit', 'legtrim', 'legroots', 'Legendre', 'legval2d', 'legval3d', | |
98 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d', 'legcompanion', | |
99 'leggauss', 'legweight'] | |
100 | |
101 legtrim = pu.trimcoef | |
102 | |
103 | |
104 def poly2leg(pol): | |
105 """ | |
106 Convert a polynomial to a Legendre series. | |
107 | |
108 Convert an array representing the coefficients of a polynomial (relative | |
109 to the "standard" basis) ordered from lowest degree to highest, to an | |
110 array of the coefficients of the equivalent Legendre series, ordered | |
111 from lowest to highest degree. | |
112 | |
113 Parameters | |
114 ---------- | |
115 pol : array_like | |
116 1-D array containing the polynomial coefficients | |
117 | |
118 Returns | |
119 ------- | |
120 c : ndarray | |
121 1-D array containing the coefficients of the equivalent Legendre | |
122 series. | |
123 | |
124 See Also | |
125 -------- | |
126 leg2poly | |
127 | |
128 Notes | |
129 ----- | |
130 The easy way to do conversions between polynomial basis sets | |
131 is to use the convert method of a class instance. | |
132 | |
133 Examples | |
134 -------- | |
135 >>> from numpy import polynomial as P | |
136 >>> p = P.Polynomial(np.arange(4)) | |
137 >>> p | |
138 Polynomial([ 0., 1., 2., 3.], [-1., 1.]) | |
139 >>> c = P.Legendre(P.poly2leg(p.coef)) | |
140 >>> c | |
141 Legendre([ 1. , 3.25, 1. , 0.75], [-1., 1.]) | |
142 | |
143 """ | |
144 [pol] = pu.as_series([pol]) | |
145 deg = len(pol) - 1 | |
146 res = 0 | |
147 for i in range(deg, -1, -1): | |
148 res = legadd(legmulx(res), pol[i]) | |
149 return res | |
150 | |
151 | |
152 def leg2poly(c): | |
153 """ | |
154 Convert a Legendre series to a polynomial. | |
155 | |
156 Convert an array representing the coefficients of a Legendre series, | |
157 ordered from lowest degree to highest, to an array of the coefficients | |
158 of the equivalent polynomial (relative to the "standard" basis) ordered | |
159 from lowest to highest degree. | |
160 | |
161 Parameters | |
162 ---------- | |
163 c : array_like | |
164 1-D array containing the Legendre series coefficients, ordered | |
165 from lowest order term to highest. | |
166 | |
167 Returns | |
168 ------- | |
169 pol : ndarray | |
170 1-D array containing the coefficients of the equivalent polynomial | |
171 (relative to the "standard" basis) ordered from lowest order term | |
172 to highest. | |
173 | |
174 See Also | |
175 -------- | |
176 poly2leg | |
177 | |
178 Notes | |
179 ----- | |
180 The easy way to do conversions between polynomial basis sets | |
181 is to use the convert method of a class instance. | |
182 | |
183 Examples | |
184 -------- | |
185 >>> c = P.Legendre(range(4)) | |
186 >>> c | |
187 Legendre([ 0., 1., 2., 3.], [-1., 1.]) | |
188 >>> p = c.convert(kind=P.Polynomial) | |
189 >>> p | |
190 Polynomial([-1. , -3.5, 3. , 7.5], [-1., 1.]) | |
191 >>> P.leg2poly(range(4)) | |
192 array([-1. , -3.5, 3. , 7.5]) | |
193 | |
194 | |
195 """ | |
196 from .polynomial import polyadd, polysub, polymulx | |
197 | |
198 [c] = pu.as_series([c]) | |
199 n = len(c) | |
200 if n < 3: | |
201 return c | |
202 else: | |
203 c0 = c[-2] | |
204 c1 = c[-1] | |
205 # i is the current degree of c1 | |
206 for i in range(n - 1, 1, -1): | |
207 tmp = c0 | |
208 c0 = polysub(c[i - 2], (c1*(i - 1))/i) | |
209 c1 = polyadd(tmp, (polymulx(c1)*(2*i - 1))/i) | |
210 return polyadd(c0, polymulx(c1)) | |
211 | |
212 # | |
213 # These are constant arrays are of integer type so as to be compatible | |
214 # with the widest range of other types, such as Decimal. | |
215 # | |
216 | |
217 # Legendre | |
218 legdomain = np.array([-1, 1]) | |
219 | |
220 # Legendre coefficients representing zero. | |
221 legzero = np.array([0]) | |
222 | |
223 # Legendre coefficients representing one. | |
224 legone = np.array([1]) | |
225 | |
226 # Legendre coefficients representing the identity x. | |
227 legx = np.array([0, 1]) | |
228 | |
229 | |
230 def legline(off, scl): | |
231 """ | |
232 Legendre series whose graph is a straight line. | |
233 | |
234 | |
235 | |
236 Parameters | |
237 ---------- | |
238 off, scl : scalars | |
239 The specified line is given by ``off + scl*x``. | |
240 | |
241 Returns | |
242 ------- | |
243 y : ndarray | |
244 This module's representation of the Legendre series for | |
245 ``off + scl*x``. | |
246 | |
247 See Also | |
248 -------- | |
249 polyline, chebline | |
250 | |
251 Examples | |
252 -------- | |
253 >>> import numpy.polynomial.legendre as L | |
254 >>> L.legline(3,2) | |
255 array([3, 2]) | |
256 >>> L.legval(-3, L.legline(3,2)) # should be -3 | |
257 -3.0 | |
258 | |
259 """ | |
260 if scl != 0: | |
261 return np.array([off, scl]) | |
262 else: | |
263 return np.array([off]) | |
264 | |
265 | |
266 def legfromroots(roots): | |
267 """ | |
268 Generate a Legendre series with given roots. | |
269 | |
270 The function returns the coefficients of the polynomial | |
271 | |
272 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), | |
273 | |
274 in Legendre form, where the `r_n` are the roots specified in `roots`. | |
275 If a zero has multiplicity n, then it must appear in `roots` n times. | |
276 For instance, if 2 is a root of multiplicity three and 3 is a root of | |
277 multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The | |
278 roots can appear in any order. | |
279 | |
280 If the returned coefficients are `c`, then | |
281 | |
282 .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x) | |
283 | |
284 The coefficient of the last term is not generally 1 for monic | |
285 polynomials in Legendre form. | |
286 | |
287 Parameters | |
288 ---------- | |
289 roots : array_like | |
290 Sequence containing the roots. | |
291 | |
292 Returns | |
293 ------- | |
294 out : ndarray | |
295 1-D array of coefficients. If all roots are real then `out` is a | |
296 real array, if some of the roots are complex, then `out` is complex | |
297 even if all the coefficients in the result are real (see Examples | |
298 below). | |
299 | |
300 See Also | |
301 -------- | |
302 polyfromroots, chebfromroots, lagfromroots, hermfromroots, | |
303 hermefromroots. | |
304 | |
305 Examples | |
306 -------- | |
307 >>> import numpy.polynomial.legendre as L | |
308 >>> L.legfromroots((-1,0,1)) # x^3 - x relative to the standard basis | |
309 array([ 0. , -0.4, 0. , 0.4]) | |
310 >>> j = complex(0,1) | |
311 >>> L.legfromroots((-j,j)) # x^2 + 1 relative to the standard basis | |
312 array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) | |
313 | |
314 """ | |
315 if len(roots) == 0: | |
316 return np.ones(1) | |
317 else: | |
318 [roots] = pu.as_series([roots], trim=False) | |
319 roots.sort() | |
320 p = [legline(-r, 1) for r in roots] | |
321 n = len(p) | |
322 while n > 1: | |
323 m, r = divmod(n, 2) | |
324 tmp = [legmul(p[i], p[i+m]) for i in range(m)] | |
325 if r: | |
326 tmp[0] = legmul(tmp[0], p[-1]) | |
327 p = tmp | |
328 n = m | |
329 return p[0] | |
330 | |
331 | |
332 def legadd(c1, c2): | |
333 """ | |
334 Add one Legendre series to another. | |
335 | |
336 Returns the sum of two Legendre series `c1` + `c2`. The arguments | |
337 are sequences of coefficients ordered from lowest order term to | |
338 highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. | |
339 | |
340 Parameters | |
341 ---------- | |
342 c1, c2 : array_like | |
343 1-D arrays of Legendre series coefficients ordered from low to | |
344 high. | |
345 | |
346 Returns | |
347 ------- | |
348 out : ndarray | |
349 Array representing the Legendre series of their sum. | |
350 | |
351 See Also | |
352 -------- | |
353 legsub, legmul, legdiv, legpow | |
354 | |
355 Notes | |
356 ----- | |
357 Unlike multiplication, division, etc., the sum of two Legendre series | |
358 is a Legendre series (without having to "reproject" the result onto | |
359 the basis set) so addition, just like that of "standard" polynomials, | |
360 is simply "component-wise." | |
361 | |
362 Examples | |
363 -------- | |
364 >>> from numpy.polynomial import legendre as L | |
365 >>> c1 = (1,2,3) | |
366 >>> c2 = (3,2,1) | |
367 >>> L.legadd(c1,c2) | |
368 array([ 4., 4., 4.]) | |
369 | |
370 """ | |
371 # c1, c2 are trimmed copies | |
372 [c1, c2] = pu.as_series([c1, c2]) | |
373 if len(c1) > len(c2): | |
374 c1[:c2.size] += c2 | |
375 ret = c1 | |
376 else: | |
377 c2[:c1.size] += c1 | |
378 ret = c2 | |
379 return pu.trimseq(ret) | |
380 | |
381 | |
382 def legsub(c1, c2): | |
383 """ | |
384 Subtract one Legendre series from another. | |
385 | |
386 Returns the difference of two Legendre series `c1` - `c2`. The | |
387 sequences of coefficients are from lowest order term to highest, i.e., | |
388 [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. | |
389 | |
390 Parameters | |
391 ---------- | |
392 c1, c2 : array_like | |
393 1-D arrays of Legendre series coefficients ordered from low to | |
394 high. | |
395 | |
396 Returns | |
397 ------- | |
398 out : ndarray | |
399 Of Legendre series coefficients representing their difference. | |
400 | |
401 See Also | |
402 -------- | |
403 legadd, legmul, legdiv, legpow | |
404 | |
405 Notes | |
406 ----- | |
407 Unlike multiplication, division, etc., the difference of two Legendre | |
408 series is a Legendre series (without having to "reproject" the result | |
409 onto the basis set) so subtraction, just like that of "standard" | |
410 polynomials, is simply "component-wise." | |
411 | |
412 Examples | |
413 -------- | |
414 >>> from numpy.polynomial import legendre as L | |
415 >>> c1 = (1,2,3) | |
416 >>> c2 = (3,2,1) | |
417 >>> L.legsub(c1,c2) | |
418 array([-2., 0., 2.]) | |
419 >>> L.legsub(c2,c1) # -C.legsub(c1,c2) | |
420 array([ 2., 0., -2.]) | |
421 | |
422 """ | |
423 # c1, c2 are trimmed copies | |
424 [c1, c2] = pu.as_series([c1, c2]) | |
425 if len(c1) > len(c2): | |
426 c1[:c2.size] -= c2 | |
427 ret = c1 | |
428 else: | |
429 c2 = -c2 | |
430 c2[:c1.size] += c1 | |
431 ret = c2 | |
432 return pu.trimseq(ret) | |
433 | |
434 | |
435 def legmulx(c): | |
436 """Multiply a Legendre series by x. | |
437 | |
438 Multiply the Legendre series `c` by x, where x is the independent | |
439 variable. | |
440 | |
441 | |
442 Parameters | |
443 ---------- | |
444 c : array_like | |
445 1-D array of Legendre series coefficients ordered from low to | |
446 high. | |
447 | |
448 Returns | |
449 ------- | |
450 out : ndarray | |
451 Array representing the result of the multiplication. | |
452 | |
453 Notes | |
454 ----- | |
455 The multiplication uses the recursion relationship for Legendre | |
456 polynomials in the form | |
457 | |
458 .. math:: | |
459 | |
460 xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) | |
461 | |
462 """ | |
463 # c is a trimmed copy | |
464 [c] = pu.as_series([c]) | |
465 # The zero series needs special treatment | |
466 if len(c) == 1 and c[0] == 0: | |
467 return c | |
468 | |
469 prd = np.empty(len(c) + 1, dtype=c.dtype) | |
470 prd[0] = c[0]*0 | |
471 prd[1] = c[0] | |
472 for i in range(1, len(c)): | |
473 j = i + 1 | |
474 k = i - 1 | |
475 s = i + j | |
476 prd[j] = (c[i]*j)/s | |
477 prd[k] += (c[i]*i)/s | |
478 return prd | |
479 | |
480 | |
481 def legmul(c1, c2): | |
482 """ | |
483 Multiply one Legendre series by another. | |
484 | |
485 Returns the product of two Legendre series `c1` * `c2`. The arguments | |
486 are sequences of coefficients, from lowest order "term" to highest, | |
487 e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. | |
488 | |
489 Parameters | |
490 ---------- | |
491 c1, c2 : array_like | |
492 1-D arrays of Legendre series coefficients ordered from low to | |
493 high. | |
494 | |
495 Returns | |
496 ------- | |
497 out : ndarray | |
498 Of Legendre series coefficients representing their product. | |
499 | |
500 See Also | |
501 -------- | |
502 legadd, legsub, legdiv, legpow | |
503 | |
504 Notes | |
505 ----- | |
506 In general, the (polynomial) product of two C-series results in terms | |
507 that are not in the Legendre polynomial basis set. Thus, to express | |
508 the product as a Legendre series, it is necessary to "reproject" the | |
509 product onto said basis set, which may produce "unintuitive" (but | |
510 correct) results; see Examples section below. | |
511 | |
512 Examples | |
513 -------- | |
514 >>> from numpy.polynomial import legendre as L | |
515 >>> c1 = (1,2,3) | |
516 >>> c2 = (3,2) | |
517 >>> P.legmul(c1,c2) # multiplication requires "reprojection" | |
518 array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) | |
519 | |
520 """ | |
521 # s1, s2 are trimmed copies | |
522 [c1, c2] = pu.as_series([c1, c2]) | |
523 | |
524 if len(c1) > len(c2): | |
525 c = c2 | |
526 xs = c1 | |
527 else: | |
528 c = c1 | |
529 xs = c2 | |
530 | |
531 if len(c) == 1: | |
532 c0 = c[0]*xs | |
533 c1 = 0 | |
534 elif len(c) == 2: | |
535 c0 = c[0]*xs | |
536 c1 = c[1]*xs | |
537 else: | |
538 nd = len(c) | |
539 c0 = c[-2]*xs | |
540 c1 = c[-1]*xs | |
541 for i in range(3, len(c) + 1): | |
542 tmp = c0 | |
543 nd = nd - 1 | |
544 c0 = legsub(c[-i]*xs, (c1*(nd - 1))/nd) | |
545 c1 = legadd(tmp, (legmulx(c1)*(2*nd - 1))/nd) | |
546 return legadd(c0, legmulx(c1)) | |
547 | |
548 | |
549 def legdiv(c1, c2): | |
550 """ | |
551 Divide one Legendre series by another. | |
552 | |
553 Returns the quotient-with-remainder of two Legendre series | |
554 `c1` / `c2`. The arguments are sequences of coefficients from lowest | |
555 order "term" to highest, e.g., [1,2,3] represents the series | |
556 ``P_0 + 2*P_1 + 3*P_2``. | |
557 | |
558 Parameters | |
559 ---------- | |
560 c1, c2 : array_like | |
561 1-D arrays of Legendre series coefficients ordered from low to | |
562 high. | |
563 | |
564 Returns | |
565 ------- | |
566 quo, rem : ndarrays | |
567 Of Legendre series coefficients representing the quotient and | |
568 remainder. | |
569 | |
570 See Also | |
571 -------- | |
572 legadd, legsub, legmul, legpow | |
573 | |
574 Notes | |
575 ----- | |
576 In general, the (polynomial) division of one Legendre series by another | |
577 results in quotient and remainder terms that are not in the Legendre | |
578 polynomial basis set. Thus, to express these results as a Legendre | |
579 series, it is necessary to "reproject" the results onto the Legendre | |
580 basis set, which may produce "unintuitive" (but correct) results; see | |
581 Examples section below. | |
582 | |
583 Examples | |
584 -------- | |
585 >>> from numpy.polynomial import legendre as L | |
586 >>> c1 = (1,2,3) | |
587 >>> c2 = (3,2,1) | |
588 >>> L.legdiv(c1,c2) # quotient "intuitive," remainder not | |
589 (array([ 3.]), array([-8., -4.])) | |
590 >>> c2 = (0,1,2,3) | |
591 >>> L.legdiv(c2,c1) # neither "intuitive" | |
592 (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) | |
593 | |
594 """ | |
595 # c1, c2 are trimmed copies | |
596 [c1, c2] = pu.as_series([c1, c2]) | |
597 if c2[-1] == 0: | |
598 raise ZeroDivisionError() | |
599 | |
600 lc1 = len(c1) | |
601 lc2 = len(c2) | |
602 if lc1 < lc2: | |
603 return c1[:1]*0, c1 | |
604 elif lc2 == 1: | |
605 return c1/c2[-1], c1[:1]*0 | |
606 else: | |
607 quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype) | |
608 rem = c1 | |
609 for i in range(lc1 - lc2, - 1, -1): | |
610 p = legmul([0]*i + [1], c2) | |
611 q = rem[-1]/p[-1] | |
612 rem = rem[:-1] - q*p[:-1] | |
613 quo[i] = q | |
614 return quo, pu.trimseq(rem) | |
615 | |
616 | |
617 def legpow(c, pow, maxpower=16): | |
618 """Raise a Legendre series to a power. | |
619 | |
620 Returns the Legendre series `c` raised to the power `pow`. The | |
621 arguement `c` is a sequence of coefficients ordered from low to high. | |
622 i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` | |
623 | |
624 Parameters | |
625 ---------- | |
626 c : array_like | |
627 1-D array of Legendre series coefficients ordered from low to | |
628 high. | |
629 pow : integer | |
630 Power to which the series will be raised | |
631 maxpower : integer, optional | |
632 Maximum power allowed. This is mainly to limit growth of the series | |
633 to unmanageable size. Default is 16 | |
634 | |
635 Returns | |
636 ------- | |
637 coef : ndarray | |
638 Legendre series of power. | |
639 | |
640 See Also | |
641 -------- | |
642 legadd, legsub, legmul, legdiv | |
643 | |
644 Examples | |
645 -------- | |
646 | |
647 """ | |
648 # c is a trimmed copy | |
649 [c] = pu.as_series([c]) | |
650 power = int(pow) | |
651 if power != pow or power < 0: | |
652 raise ValueError("Power must be a non-negative integer.") | |
653 elif maxpower is not None and power > maxpower: | |
654 raise ValueError("Power is too large") | |
655 elif power == 0: | |
656 return np.array([1], dtype=c.dtype) | |
657 elif power == 1: | |
658 return c | |
659 else: | |
660 # This can be made more efficient by using powers of two | |
661 # in the usual way. | |
662 prd = c | |
663 for i in range(2, power + 1): | |
664 prd = legmul(prd, c) | |
665 return prd | |
666 | |
667 | |
668 def legder(c, m=1, scl=1, axis=0): | |
669 """ | |
670 Differentiate a Legendre series. | |
671 | |
672 Returns the Legendre series coefficients `c` differentiated `m` times | |
673 along `axis`. At each iteration the result is multiplied by `scl` (the | |
674 scaling factor is for use in a linear change of variable). The argument | |
675 `c` is an array of coefficients from low to high degree along each | |
676 axis, e.g., [1,2,3] represents the series ``1*L_0 + 2*L_1 + 3*L_2`` | |
677 while [[1,2],[1,2]] represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + | |
678 2*L_0(x)*L_1(y) + 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is | |
679 ``y``. | |
680 | |
681 Parameters | |
682 ---------- | |
683 c : array_like | |
684 Array of Legendre series coefficients. If c is multidimensional the | |
685 different axis correspond to different variables with the degree in | |
686 each axis given by the corresponding index. | |
687 m : int, optional | |
688 Number of derivatives taken, must be non-negative. (Default: 1) | |
689 scl : scalar, optional | |
690 Each differentiation is multiplied by `scl`. The end result is | |
691 multiplication by ``scl**m``. This is for use in a linear change of | |
692 variable. (Default: 1) | |
693 axis : int, optional | |
694 Axis over which the derivative is taken. (Default: 0). | |
695 | |
696 .. versionadded:: 1.7.0 | |
697 | |
698 Returns | |
699 ------- | |
700 der : ndarray | |
701 Legendre series of the derivative. | |
702 | |
703 See Also | |
704 -------- | |
705 legint | |
706 | |
707 Notes | |
708 ----- | |
709 In general, the result of differentiating a Legendre series does not | |
710 resemble the same operation on a power series. Thus the result of this | |
711 function may be "unintuitive," albeit correct; see Examples section | |
712 below. | |
713 | |
714 Examples | |
715 -------- | |
716 >>> from numpy.polynomial import legendre as L | |
717 >>> c = (1,2,3,4) | |
718 >>> L.legder(c) | |
719 array([ 6., 9., 20.]) | |
720 >>> L.legder(c, 3) | |
721 array([ 60.]) | |
722 >>> L.legder(c, scl=-1) | |
723 array([ -6., -9., -20.]) | |
724 >>> L.legder(c, 2,-1) | |
725 array([ 9., 60.]) | |
726 | |
727 """ | |
728 c = np.array(c, ndmin=1, copy=1) | |
729 if c.dtype.char in '?bBhHiIlLqQpP': | |
730 c = c.astype(np.double) | |
731 cnt, iaxis = [int(t) for t in [m, axis]] | |
732 | |
733 if cnt != m: | |
734 raise ValueError("The order of derivation must be integer") | |
735 if cnt < 0: | |
736 raise ValueError("The order of derivation must be non-negative") | |
737 if iaxis != axis: | |
738 raise ValueError("The axis must be integer") | |
739 if not -c.ndim <= iaxis < c.ndim: | |
740 raise ValueError("The axis is out of range") | |
741 if iaxis < 0: | |
742 iaxis += c.ndim | |
743 | |
744 if cnt == 0: | |
745 return c | |
746 | |
747 c = np.rollaxis(c, iaxis) | |
748 n = len(c) | |
749 if cnt >= n: | |
750 c = c[:1]*0 | |
751 else: | |
752 for i in range(cnt): | |
753 n = n - 1 | |
754 c *= scl | |
755 der = np.empty((n,) + c.shape[1:], dtype=c.dtype) | |
756 for j in range(n, 2, -1): | |
757 der[j - 1] = (2*j - 1)*c[j] | |
758 c[j - 2] += c[j] | |
759 if n > 1: | |
760 der[1] = 3*c[2] | |
761 der[0] = c[1] | |
762 c = der | |
763 c = np.rollaxis(c, 0, iaxis + 1) | |
764 return c | |
765 | |
766 | |
767 def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0): | |
768 """ | |
769 Integrate a Legendre series. | |
770 | |
771 Returns the Legendre series coefficients `c` integrated `m` times from | |
772 `lbnd` along `axis`. At each iteration the resulting series is | |
773 **multiplied** by `scl` and an integration constant, `k`, is added. | |
774 The scaling factor is for use in a linear change of variable. ("Buyer | |
775 beware": note that, depending on what one is doing, one may want `scl` | |
776 to be the reciprocal of what one might expect; for more information, | |
777 see the Notes section below.) The argument `c` is an array of | |
778 coefficients from low to high degree along each axis, e.g., [1,2,3] | |
779 represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]] | |
780 represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) + | |
781 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``. | |
782 | |
783 Parameters | |
784 ---------- | |
785 c : array_like | |
786 Array of Legendre series coefficients. If c is multidimensional the | |
787 different axis correspond to different variables with the degree in | |
788 each axis given by the corresponding index. | |
789 m : int, optional | |
790 Order of integration, must be positive. (Default: 1) | |
791 k : {[], list, scalar}, optional | |
792 Integration constant(s). The value of the first integral at | |
793 ``lbnd`` is the first value in the list, the value of the second | |
794 integral at ``lbnd`` is the second value, etc. If ``k == []`` (the | |
795 default), all constants are set to zero. If ``m == 1``, a single | |
796 scalar can be given instead of a list. | |
797 lbnd : scalar, optional | |
798 The lower bound of the integral. (Default: 0) | |
799 scl : scalar, optional | |
800 Following each integration the result is *multiplied* by `scl` | |
801 before the integration constant is added. (Default: 1) | |
802 axis : int, optional | |
803 Axis over which the integral is taken. (Default: 0). | |
804 | |
805 .. versionadded:: 1.7.0 | |
806 | |
807 Returns | |
808 ------- | |
809 S : ndarray | |
810 Legendre series coefficient array of the integral. | |
811 | |
812 Raises | |
813 ------ | |
814 ValueError | |
815 If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or | |
816 ``np.isscalar(scl) == False``. | |
817 | |
818 See Also | |
819 -------- | |
820 legder | |
821 | |
822 Notes | |
823 ----- | |
824 Note that the result of each integration is *multiplied* by `scl`. | |
825 Why is this important to note? Say one is making a linear change of | |
826 variable :math:`u = ax + b` in an integral relative to `x`. Then | |
827 .. math::`dx = du/a`, so one will need to set `scl` equal to | |
828 :math:`1/a` - perhaps not what one would have first thought. | |
829 | |
830 Also note that, in general, the result of integrating a C-series needs | |
831 to be "reprojected" onto the C-series basis set. Thus, typically, | |
832 the result of this function is "unintuitive," albeit correct; see | |
833 Examples section below. | |
834 | |
835 Examples | |
836 -------- | |
837 >>> from numpy.polynomial import legendre as L | |
838 >>> c = (1,2,3) | |
839 >>> L.legint(c) | |
840 array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) | |
841 >>> L.legint(c, 3) | |
842 array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, | |
843 -1.73472348e-18, 1.90476190e-02, 9.52380952e-03]) | |
844 >>> L.legint(c, k=3) | |
845 array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) | |
846 >>> L.legint(c, lbnd=-2) | |
847 array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) | |
848 >>> L.legint(c, scl=2) | |
849 array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) | |
850 | |
851 """ | |
852 c = np.array(c, ndmin=1, copy=1) | |
853 if c.dtype.char in '?bBhHiIlLqQpP': | |
854 c = c.astype(np.double) | |
855 if not np.iterable(k): | |
856 k = [k] | |
857 cnt, iaxis = [int(t) for t in [m, axis]] | |
858 | |
859 if cnt != m: | |
860 raise ValueError("The order of integration must be integer") | |
861 if cnt < 0: | |
862 raise ValueError("The order of integration must be non-negative") | |
863 if len(k) > cnt: | |
864 raise ValueError("Too many integration constants") | |
865 if iaxis != axis: | |
866 raise ValueError("The axis must be integer") | |
867 if not -c.ndim <= iaxis < c.ndim: | |
868 raise ValueError("The axis is out of range") | |
869 if iaxis < 0: | |
870 iaxis += c.ndim | |
871 | |
872 if cnt == 0: | |
873 return c | |
874 | |
875 c = np.rollaxis(c, iaxis) | |
876 k = list(k) + [0]*(cnt - len(k)) | |
877 for i in range(cnt): | |
878 n = len(c) | |
879 c *= scl | |
880 if n == 1 and np.all(c[0] == 0): | |
881 c[0] += k[i] | |
882 else: | |
883 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype) | |
884 tmp[0] = c[0]*0 | |
885 tmp[1] = c[0] | |
886 if n > 1: | |
887 tmp[2] = c[1]/3 | |
888 for j in range(2, n): | |
889 t = c[j]/(2*j + 1) | |
890 tmp[j + 1] = t | |
891 tmp[j - 1] -= t | |
892 tmp[0] += k[i] - legval(lbnd, tmp) | |
893 c = tmp | |
894 c = np.rollaxis(c, 0, iaxis + 1) | |
895 return c | |
896 | |
897 | |
898 def legval(x, c, tensor=True): | |
899 """ | |
900 Evaluate a Legendre series at points x. | |
901 | |
902 If `c` is of length `n + 1`, this function returns the value: | |
903 | |
904 .. math:: p(x) = c_0 * L_0(x) + c_1 * L_1(x) + ... + c_n * L_n(x) | |
905 | |
906 The parameter `x` is converted to an array only if it is a tuple or a | |
907 list, otherwise it is treated as a scalar. In either case, either `x` | |
908 or its elements must support multiplication and addition both with | |
909 themselves and with the elements of `c`. | |
910 | |
911 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If | |
912 `c` is multidimensional, then the shape of the result depends on the | |
913 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + | |
914 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that | |
915 scalars have shape (,). | |
916 | |
917 Trailing zeros in the coefficients will be used in the evaluation, so | |
918 they should be avoided if efficiency is a concern. | |
919 | |
920 Parameters | |
921 ---------- | |
922 x : array_like, compatible object | |
923 If `x` is a list or tuple, it is converted to an ndarray, otherwise | |
924 it is left unchanged and treated as a scalar. In either case, `x` | |
925 or its elements must support addition and multiplication with | |
926 with themselves and with the elements of `c`. | |
927 c : array_like | |
928 Array of coefficients ordered so that the coefficients for terms of | |
929 degree n are contained in c[n]. If `c` is multidimensional the | |
930 remaining indices enumerate multiple polynomials. In the two | |
931 dimensional case the coefficients may be thought of as stored in | |
932 the columns of `c`. | |
933 tensor : boolean, optional | |
934 If True, the shape of the coefficient array is extended with ones | |
935 on the right, one for each dimension of `x`. Scalars have dimension 0 | |
936 for this action. The result is that every column of coefficients in | |
937 `c` is evaluated for every element of `x`. If False, `x` is broadcast | |
938 over the columns of `c` for the evaluation. This keyword is useful | |
939 when `c` is multidimensional. The default value is True. | |
940 | |
941 .. versionadded:: 1.7.0 | |
942 | |
943 Returns | |
944 ------- | |
945 values : ndarray, algebra_like | |
946 The shape of the return value is described above. | |
947 | |
948 See Also | |
949 -------- | |
950 legval2d, leggrid2d, legval3d, leggrid3d | |
951 | |
952 Notes | |
953 ----- | |
954 The evaluation uses Clenshaw recursion, aka synthetic division. | |
955 | |
956 Examples | |
957 -------- | |
958 | |
959 """ | |
960 c = np.array(c, ndmin=1, copy=0) | |
961 if c.dtype.char in '?bBhHiIlLqQpP': | |
962 c = c.astype(np.double) | |
963 if isinstance(x, (tuple, list)): | |
964 x = np.asarray(x) | |
965 if isinstance(x, np.ndarray) and tensor: | |
966 c = c.reshape(c.shape + (1,)*x.ndim) | |
967 | |
968 if len(c) == 1: | |
969 c0 = c[0] | |
970 c1 = 0 | |
971 elif len(c) == 2: | |
972 c0 = c[0] | |
973 c1 = c[1] | |
974 else: | |
975 nd = len(c) | |
976 c0 = c[-2] | |
977 c1 = c[-1] | |
978 for i in range(3, len(c) + 1): | |
979 tmp = c0 | |
980 nd = nd - 1 | |
981 c0 = c[-i] - (c1*(nd - 1))/nd | |
982 c1 = tmp + (c1*x*(2*nd - 1))/nd | |
983 return c0 + c1*x | |
984 | |
985 | |
986 def legval2d(x, y, c): | |
987 """ | |
988 Evaluate a 2-D Legendre series at points (x, y). | |
989 | |
990 This function returns the values: | |
991 | |
992 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * L_i(x) * L_j(y) | |
993 | |
994 The parameters `x` and `y` are converted to arrays only if they are | |
995 tuples or a lists, otherwise they are treated as a scalars and they | |
996 must have the same shape after conversion. In either case, either `x` | |
997 and `y` or their elements must support multiplication and addition both | |
998 with themselves and with the elements of `c`. | |
999 | |
1000 If `c` is a 1-D array a one is implicitly appended to its shape to make | |
1001 it 2-D. The shape of the result will be c.shape[2:] + x.shape. | |
1002 | |
1003 Parameters | |
1004 ---------- | |
1005 x, y : array_like, compatible objects | |
1006 The two dimensional series is evaluated at the points `(x, y)`, | |
1007 where `x` and `y` must have the same shape. If `x` or `y` is a list | |
1008 or tuple, it is first converted to an ndarray, otherwise it is left | |
1009 unchanged and if it isn't an ndarray it is treated as a scalar. | |
1010 c : array_like | |
1011 Array of coefficients ordered so that the coefficient of the term | |
1012 of multi-degree i,j is contained in ``c[i,j]``. If `c` has | |
1013 dimension greater than two the remaining indices enumerate multiple | |
1014 sets of coefficients. | |
1015 | |
1016 Returns | |
1017 ------- | |
1018 values : ndarray, compatible object | |
1019 The values of the two dimensional Legendre series at points formed | |
1020 from pairs of corresponding values from `x` and `y`. | |
1021 | |
1022 See Also | |
1023 -------- | |
1024 legval, leggrid2d, legval3d, leggrid3d | |
1025 | |
1026 Notes | |
1027 ----- | |
1028 | |
1029 .. versionadded::1.7.0 | |
1030 | |
1031 """ | |
1032 try: | |
1033 x, y = np.array((x, y), copy=0) | |
1034 except: | |
1035 raise ValueError('x, y are incompatible') | |
1036 | |
1037 c = legval(x, c) | |
1038 c = legval(y, c, tensor=False) | |
1039 return c | |
1040 | |
1041 | |
1042 def leggrid2d(x, y, c): | |
1043 """ | |
1044 Evaluate a 2-D Legendre series on the Cartesian product of x and y. | |
1045 | |
1046 This function returns the values: | |
1047 | |
1048 .. math:: p(a,b) = \sum_{i,j} c_{i,j} * L_i(a) * L_j(b) | |
1049 | |
1050 where the points `(a, b)` consist of all pairs formed by taking | |
1051 `a` from `x` and `b` from `y`. The resulting points form a grid with | |
1052 `x` in the first dimension and `y` in the second. | |
1053 | |
1054 The parameters `x` and `y` are converted to arrays only if they are | |
1055 tuples or a lists, otherwise they are treated as a scalars. In either | |
1056 case, either `x` and `y` or their elements must support multiplication | |
1057 and addition both with themselves and with the elements of `c`. | |
1058 | |
1059 If `c` has fewer than two dimensions, ones are implicitly appended to | |
1060 its shape to make it 2-D. The shape of the result will be c.shape[2:] + | |
1061 x.shape + y.shape. | |
1062 | |
1063 Parameters | |
1064 ---------- | |
1065 x, y : array_like, compatible objects | |
1066 The two dimensional series is evaluated at the points in the | |
1067 Cartesian product of `x` and `y`. If `x` or `y` is a list or | |
1068 tuple, it is first converted to an ndarray, otherwise it is left | |
1069 unchanged and, if it isn't an ndarray, it is treated as a scalar. | |
1070 c : array_like | |
1071 Array of coefficients ordered so that the coefficient of the term of | |
1072 multi-degree i,j is contained in `c[i,j]`. If `c` has dimension | |
1073 greater than two the remaining indices enumerate multiple sets of | |
1074 coefficients. | |
1075 | |
1076 Returns | |
1077 ------- | |
1078 values : ndarray, compatible object | |
1079 The values of the two dimensional Chebyshev series at points in the | |
1080 Cartesian product of `x` and `y`. | |
1081 | |
1082 See Also | |
1083 -------- | |
1084 legval, legval2d, legval3d, leggrid3d | |
1085 | |
1086 Notes | |
1087 ----- | |
1088 | |
1089 .. versionadded::1.7.0 | |
1090 | |
1091 """ | |
1092 c = legval(x, c) | |
1093 c = legval(y, c) | |
1094 return c | |
1095 | |
1096 | |
1097 def legval3d(x, y, z, c): | |
1098 """ | |
1099 Evaluate a 3-D Legendre series at points (x, y, z). | |
1100 | |
1101 This function returns the values: | |
1102 | |
1103 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * L_i(x) * L_j(y) * L_k(z) | |
1104 | |
1105 The parameters `x`, `y`, and `z` are converted to arrays only if | |
1106 they are tuples or a lists, otherwise they are treated as a scalars and | |
1107 they must have the same shape after conversion. In either case, either | |
1108 `x`, `y`, and `z` or their elements must support multiplication and | |
1109 addition both with themselves and with the elements of `c`. | |
1110 | |
1111 If `c` has fewer than 3 dimensions, ones are implicitly appended to its | |
1112 shape to make it 3-D. The shape of the result will be c.shape[3:] + | |
1113 x.shape. | |
1114 | |
1115 Parameters | |
1116 ---------- | |
1117 x, y, z : array_like, compatible object | |
1118 The three dimensional series is evaluated at the points | |
1119 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If | |
1120 any of `x`, `y`, or `z` is a list or tuple, it is first converted | |
1121 to an ndarray, otherwise it is left unchanged and if it isn't an | |
1122 ndarray it is treated as a scalar. | |
1123 c : array_like | |
1124 Array of coefficients ordered so that the coefficient of the term of | |
1125 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension | |
1126 greater than 3 the remaining indices enumerate multiple sets of | |
1127 coefficients. | |
1128 | |
1129 Returns | |
1130 ------- | |
1131 values : ndarray, compatible object | |
1132 The values of the multidimensional polynomial on points formed with | |
1133 triples of corresponding values from `x`, `y`, and `z`. | |
1134 | |
1135 See Also | |
1136 -------- | |
1137 legval, legval2d, leggrid2d, leggrid3d | |
1138 | |
1139 Notes | |
1140 ----- | |
1141 | |
1142 .. versionadded::1.7.0 | |
1143 | |
1144 """ | |
1145 try: | |
1146 x, y, z = np.array((x, y, z), copy=0) | |
1147 except: | |
1148 raise ValueError('x, y, z are incompatible') | |
1149 | |
1150 c = legval(x, c) | |
1151 c = legval(y, c, tensor=False) | |
1152 c = legval(z, c, tensor=False) | |
1153 return c | |
1154 | |
1155 | |
1156 def leggrid3d(x, y, z, c): | |
1157 """ | |
1158 Evaluate a 3-D Legendre series on the Cartesian product of x, y, and z. | |
1159 | |
1160 This function returns the values: | |
1161 | |
1162 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * L_i(a) * L_j(b) * L_k(c) | |
1163 | |
1164 where the points `(a, b, c)` consist of all triples formed by taking | |
1165 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form | |
1166 a grid with `x` in the first dimension, `y` in the second, and `z` in | |
1167 the third. | |
1168 | |
1169 The parameters `x`, `y`, and `z` are converted to arrays only if they | |
1170 are tuples or a lists, otherwise they are treated as a scalars. In | |
1171 either case, either `x`, `y`, and `z` or their elements must support | |
1172 multiplication and addition both with themselves and with the elements | |
1173 of `c`. | |
1174 | |
1175 If `c` has fewer than three dimensions, ones are implicitly appended to | |
1176 its shape to make it 3-D. The shape of the result will be c.shape[3:] + | |
1177 x.shape + y.shape + z.shape. | |
1178 | |
1179 Parameters | |
1180 ---------- | |
1181 x, y, z : array_like, compatible objects | |
1182 The three dimensional series is evaluated at the points in the | |
1183 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a | |
1184 list or tuple, it is first converted to an ndarray, otherwise it is | |
1185 left unchanged and, if it isn't an ndarray, it is treated as a | |
1186 scalar. | |
1187 c : array_like | |
1188 Array of coefficients ordered so that the coefficients for terms of | |
1189 degree i,j are contained in ``c[i,j]``. If `c` has dimension | |
1190 greater than two the remaining indices enumerate multiple sets of | |
1191 coefficients. | |
1192 | |
1193 Returns | |
1194 ------- | |
1195 values : ndarray, compatible object | |
1196 The values of the two dimensional polynomial at points in the Cartesian | |
1197 product of `x` and `y`. | |
1198 | |
1199 See Also | |
1200 -------- | |
1201 legval, legval2d, leggrid2d, legval3d | |
1202 | |
1203 Notes | |
1204 ----- | |
1205 | |
1206 .. versionadded::1.7.0 | |
1207 | |
1208 """ | |
1209 c = legval(x, c) | |
1210 c = legval(y, c) | |
1211 c = legval(z, c) | |
1212 return c | |
1213 | |
1214 | |
1215 def legvander(x, deg): | |
1216 """Pseudo-Vandermonde matrix of given degree. | |
1217 | |
1218 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points | |
1219 `x`. The pseudo-Vandermonde matrix is defined by | |
1220 | |
1221 .. math:: V[..., i] = L_i(x) | |
1222 | |
1223 where `0 <= i <= deg`. The leading indices of `V` index the elements of | |
1224 `x` and the last index is the degree of the Legendre polynomial. | |
1225 | |
1226 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the | |
1227 array ``V = legvander(x, n)``, then ``np.dot(V, c)`` and | |
1228 ``legval(x, c)`` are the same up to roundoff. This equivalence is | |
1229 useful both for least squares fitting and for the evaluation of a large | |
1230 number of Legendre series of the same degree and sample points. | |
1231 | |
1232 Parameters | |
1233 ---------- | |
1234 x : array_like | |
1235 Array of points. The dtype is converted to float64 or complex128 | |
1236 depending on whether any of the elements are complex. If `x` is | |
1237 scalar it is converted to a 1-D array. | |
1238 deg : int | |
1239 Degree of the resulting matrix. | |
1240 | |
1241 Returns | |
1242 ------- | |
1243 vander : ndarray | |
1244 The pseudo-Vandermonde matrix. The shape of the returned matrix is | |
1245 ``x.shape + (deg + 1,)``, where The last index is the degree of the | |
1246 corresponding Legendre polynomial. The dtype will be the same as | |
1247 the converted `x`. | |
1248 | |
1249 """ | |
1250 ideg = int(deg) | |
1251 if ideg != deg: | |
1252 raise ValueError("deg must be integer") | |
1253 if ideg < 0: | |
1254 raise ValueError("deg must be non-negative") | |
1255 | |
1256 x = np.array(x, copy=0, ndmin=1) + 0.0 | |
1257 dims = (ideg + 1,) + x.shape | |
1258 dtyp = x.dtype | |
1259 v = np.empty(dims, dtype=dtyp) | |
1260 # Use forward recursion to generate the entries. This is not as accurate | |
1261 # as reverse recursion in this application but it is more efficient. | |
1262 v[0] = x*0 + 1 | |
1263 if ideg > 0: | |
1264 v[1] = x | |
1265 for i in range(2, ideg + 1): | |
1266 v[i] = (v[i-1]*x*(2*i - 1) - v[i-2]*(i - 1))/i | |
1267 return np.rollaxis(v, 0, v.ndim) | |
1268 | |
1269 | |
1270 def legvander2d(x, y, deg): | |
1271 """Pseudo-Vandermonde matrix of given degrees. | |
1272 | |
1273 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample | |
1274 points `(x, y)`. The pseudo-Vandermonde matrix is defined by | |
1275 | |
1276 .. math:: V[..., deg[1]*i + j] = L_i(x) * L_j(y), | |
1277 | |
1278 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of | |
1279 `V` index the points `(x, y)` and the last index encodes the degrees of | |
1280 the Legendre polynomials. | |
1281 | |
1282 If ``V = legvander2d(x, y, [xdeg, ydeg])``, then the columns of `V` | |
1283 correspond to the elements of a 2-D coefficient array `c` of shape | |
1284 (xdeg + 1, ydeg + 1) in the order | |
1285 | |
1286 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ... | |
1287 | |
1288 and ``np.dot(V, c.flat)`` and ``legval2d(x, y, c)`` will be the same | |
1289 up to roundoff. This equivalence is useful both for least squares | |
1290 fitting and for the evaluation of a large number of 2-D Legendre | |
1291 series of the same degrees and sample points. | |
1292 | |
1293 Parameters | |
1294 ---------- | |
1295 x, y : array_like | |
1296 Arrays of point coordinates, all of the same shape. The dtypes | |
1297 will be converted to either float64 or complex128 depending on | |
1298 whether any of the elements are complex. Scalars are converted to | |
1299 1-D arrays. | |
1300 deg : list of ints | |
1301 List of maximum degrees of the form [x_deg, y_deg]. | |
1302 | |
1303 Returns | |
1304 ------- | |
1305 vander2d : ndarray | |
1306 The shape of the returned matrix is ``x.shape + (order,)``, where | |
1307 :math:`order = (deg[0]+1)*(deg([1]+1)`. The dtype will be the same | |
1308 as the converted `x` and `y`. | |
1309 | |
1310 See Also | |
1311 -------- | |
1312 legvander, legvander3d. legval2d, legval3d | |
1313 | |
1314 Notes | |
1315 ----- | |
1316 | |
1317 .. versionadded::1.7.0 | |
1318 | |
1319 """ | |
1320 ideg = [int(d) for d in deg] | |
1321 is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)] | |
1322 if is_valid != [1, 1]: | |
1323 raise ValueError("degrees must be non-negative integers") | |
1324 degx, degy = ideg | |
1325 x, y = np.array((x, y), copy=0) + 0.0 | |
1326 | |
1327 vx = legvander(x, degx) | |
1328 vy = legvander(y, degy) | |
1329 v = vx[..., None]*vy[..., None,:] | |
1330 return v.reshape(v.shape[:-2] + (-1,)) | |
1331 | |
1332 | |
1333 def legvander3d(x, y, z, deg): | |
1334 """Pseudo-Vandermonde matrix of given degrees. | |
1335 | |
1336 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample | |
1337 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`, | |
1338 then The pseudo-Vandermonde matrix is defined by | |
1339 | |
1340 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = L_i(x)*L_j(y)*L_k(z), | |
1341 | |
1342 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading | |
1343 indices of `V` index the points `(x, y, z)` and the last index encodes | |
1344 the degrees of the Legendre polynomials. | |
1345 | |
1346 If ``V = legvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns | |
1347 of `V` correspond to the elements of a 3-D coefficient array `c` of | |
1348 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order | |
1349 | |
1350 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},... | |
1351 | |
1352 and ``np.dot(V, c.flat)`` and ``legval3d(x, y, z, c)`` will be the | |
1353 same up to roundoff. This equivalence is useful both for least squares | |
1354 fitting and for the evaluation of a large number of 3-D Legendre | |
1355 series of the same degrees and sample points. | |
1356 | |
1357 Parameters | |
1358 ---------- | |
1359 x, y, z : array_like | |
1360 Arrays of point coordinates, all of the same shape. The dtypes will | |
1361 be converted to either float64 or complex128 depending on whether | |
1362 any of the elements are complex. Scalars are converted to 1-D | |
1363 arrays. | |
1364 deg : list of ints | |
1365 List of maximum degrees of the form [x_deg, y_deg, z_deg]. | |
1366 | |
1367 Returns | |
1368 ------- | |
1369 vander3d : ndarray | |
1370 The shape of the returned matrix is ``x.shape + (order,)``, where | |
1371 :math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will | |
1372 be the same as the converted `x`, `y`, and `z`. | |
1373 | |
1374 See Also | |
1375 -------- | |
1376 legvander, legvander3d. legval2d, legval3d | |
1377 | |
1378 Notes | |
1379 ----- | |
1380 | |
1381 .. versionadded::1.7.0 | |
1382 | |
1383 """ | |
1384 ideg = [int(d) for d in deg] | |
1385 is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)] | |
1386 if is_valid != [1, 1, 1]: | |
1387 raise ValueError("degrees must be non-negative integers") | |
1388 degx, degy, degz = ideg | |
1389 x, y, z = np.array((x, y, z), copy=0) + 0.0 | |
1390 | |
1391 vx = legvander(x, degx) | |
1392 vy = legvander(y, degy) | |
1393 vz = legvander(z, degz) | |
1394 v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:] | |
1395 return v.reshape(v.shape[:-3] + (-1,)) | |
1396 | |
1397 | |
1398 def legfit(x, y, deg, rcond=None, full=False, w=None): | |
1399 """ | |
1400 Least squares fit of Legendre series to data. | |
1401 | |
1402 Return the coefficients of a Legendre series of degree `deg` that is the | |
1403 least squares fit to the data values `y` given at points `x`. If `y` is | |
1404 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple | |
1405 fits are done, one for each column of `y`, and the resulting | |
1406 coefficients are stored in the corresponding columns of a 2-D return. | |
1407 The fitted polynomial(s) are in the form | |
1408 | |
1409 .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x), | |
1410 | |
1411 where `n` is `deg`. | |
1412 | |
1413 Parameters | |
1414 ---------- | |
1415 x : array_like, shape (M,) | |
1416 x-coordinates of the M sample points ``(x[i], y[i])``. | |
1417 y : array_like, shape (M,) or (M, K) | |
1418 y-coordinates of the sample points. Several data sets of sample | |
1419 points sharing the same x-coordinates can be fitted at once by | |
1420 passing in a 2D-array that contains one dataset per column. | |
1421 deg : int | |
1422 Degree of the fitting polynomial | |
1423 rcond : float, optional | |
1424 Relative condition number of the fit. Singular values smaller than | |
1425 this relative to the largest singular value will be ignored. The | |
1426 default value is len(x)*eps, where eps is the relative precision of | |
1427 the float type, about 2e-16 in most cases. | |
1428 full : bool, optional | |
1429 Switch determining nature of return value. When it is False (the | |
1430 default) just the coefficients are returned, when True diagnostic | |
1431 information from the singular value decomposition is also returned. | |
1432 w : array_like, shape (`M`,), optional | |
1433 Weights. If not None, the contribution of each point | |
1434 ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the | |
1435 weights are chosen so that the errors of the products ``w[i]*y[i]`` | |
1436 all have the same variance. The default value is None. | |
1437 | |
1438 .. versionadded:: 1.5.0 | |
1439 | |
1440 Returns | |
1441 ------- | |
1442 coef : ndarray, shape (M,) or (M, K) | |
1443 Legendre coefficients ordered from low to high. If `y` was 2-D, | |
1444 the coefficients for the data in column k of `y` are in column | |
1445 `k`. | |
1446 | |
1447 [residuals, rank, singular_values, rcond] : list | |
1448 These values are only returned if `full` = True | |
1449 | |
1450 resid -- sum of squared residuals of the least squares fit | |
1451 rank -- the numerical rank of the scaled Vandermonde matrix | |
1452 sv -- singular values of the scaled Vandermonde matrix | |
1453 rcond -- value of `rcond`. | |
1454 | |
1455 For more details, see `linalg.lstsq`. | |
1456 | |
1457 Warns | |
1458 ----- | |
1459 RankWarning | |
1460 The rank of the coefficient matrix in the least-squares fit is | |
1461 deficient. The warning is only raised if `full` = False. The | |
1462 warnings can be turned off by | |
1463 | |
1464 >>> import warnings | |
1465 >>> warnings.simplefilter('ignore', RankWarning) | |
1466 | |
1467 See Also | |
1468 -------- | |
1469 chebfit, polyfit, lagfit, hermfit, hermefit | |
1470 legval : Evaluates a Legendre series. | |
1471 legvander : Vandermonde matrix of Legendre series. | |
1472 legweight : Legendre weight function (= 1). | |
1473 linalg.lstsq : Computes a least-squares fit from the matrix. | |
1474 scipy.interpolate.UnivariateSpline : Computes spline fits. | |
1475 | |
1476 Notes | |
1477 ----- | |
1478 The solution is the coefficients of the Legendre series `p` that | |
1479 minimizes the sum of the weighted squared errors | |
1480 | |
1481 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2, | |
1482 | |
1483 where :math:`w_j` are the weights. This problem is solved by setting up | |
1484 as the (typically) overdetermined matrix equation | |
1485 | |
1486 .. math:: V(x) * c = w * y, | |
1487 | |
1488 where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the | |
1489 coefficients to be solved for, `w` are the weights, and `y` are the | |
1490 observed values. This equation is then solved using the singular value | |
1491 decomposition of `V`. | |
1492 | |
1493 If some of the singular values of `V` are so small that they are | |
1494 neglected, then a `RankWarning` will be issued. This means that the | |
1495 coefficient values may be poorly determined. Using a lower order fit | |
1496 will usually get rid of the warning. The `rcond` parameter can also be | |
1497 set to a value smaller than its default, but the resulting fit may be | |
1498 spurious and have large contributions from roundoff error. | |
1499 | |
1500 Fits using Legendre series are usually better conditioned than fits | |
1501 using power series, but much can depend on the distribution of the | |
1502 sample points and the smoothness of the data. If the quality of the fit | |
1503 is inadequate splines may be a good alternative. | |
1504 | |
1505 References | |
1506 ---------- | |
1507 .. [1] Wikipedia, "Curve fitting", | |
1508 http://en.wikipedia.org/wiki/Curve_fitting | |
1509 | |
1510 Examples | |
1511 -------- | |
1512 | |
1513 """ | |
1514 order = int(deg) + 1 | |
1515 x = np.asarray(x) + 0.0 | |
1516 y = np.asarray(y) + 0.0 | |
1517 | |
1518 # check arguments. | |
1519 if deg < 0: | |
1520 raise ValueError("expected deg >= 0") | |
1521 if x.ndim != 1: | |
1522 raise TypeError("expected 1D vector for x") | |
1523 if x.size == 0: | |
1524 raise TypeError("expected non-empty vector for x") | |
1525 if y.ndim < 1 or y.ndim > 2: | |
1526 raise TypeError("expected 1D or 2D array for y") | |
1527 if len(x) != len(y): | |
1528 raise TypeError("expected x and y to have same length") | |
1529 | |
1530 # set up the least squares matrices in transposed form | |
1531 lhs = legvander(x, deg).T | |
1532 rhs = y.T | |
1533 if w is not None: | |
1534 w = np.asarray(w) + 0.0 | |
1535 if w.ndim != 1: | |
1536 raise TypeError("expected 1D vector for w") | |
1537 if len(x) != len(w): | |
1538 raise TypeError("expected x and w to have same length") | |
1539 # apply weights. Don't use inplace operations as they | |
1540 # can cause problems with NA. | |
1541 lhs = lhs * w | |
1542 rhs = rhs * w | |
1543 | |
1544 # set rcond | |
1545 if rcond is None: | |
1546 rcond = len(x)*np.finfo(x.dtype).eps | |
1547 | |
1548 # Determine the norms of the design matrix columns. | |
1549 if issubclass(lhs.dtype.type, np.complexfloating): | |
1550 scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1)) | |
1551 else: | |
1552 scl = np.sqrt(np.square(lhs).sum(1)) | |
1553 scl[scl == 0] = 1 | |
1554 | |
1555 # Solve the least squares problem. | |
1556 c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) | |
1557 c = (c.T/scl).T | |
1558 | |
1559 # warn on rank reduction | |
1560 if rank != order and not full: | |
1561 msg = "The fit may be poorly conditioned" | |
1562 warnings.warn(msg, pu.RankWarning) | |
1563 | |
1564 if full: | |
1565 return c, [resids, rank, s, rcond] | |
1566 else: | |
1567 return c | |
1568 | |
1569 | |
1570 def legcompanion(c): | |
1571 """Return the scaled companion matrix of c. | |
1572 | |
1573 The basis polynomials are scaled so that the companion matrix is | |
1574 symmetric when `c` is an Legendre basis polynomial. This provides | |
1575 better eigenvalue estimates than the unscaled case and for basis | |
1576 polynomials the eigenvalues are guaranteed to be real if | |
1577 `numpy.linalg.eigvalsh` is used to obtain them. | |
1578 | |
1579 Parameters | |
1580 ---------- | |
1581 c : array_like | |
1582 1-D array of Legendre series coefficients ordered from low to high | |
1583 degree. | |
1584 | |
1585 Returns | |
1586 ------- | |
1587 mat : ndarray | |
1588 Scaled companion matrix of dimensions (deg, deg). | |
1589 | |
1590 Notes | |
1591 ----- | |
1592 | |
1593 .. versionadded::1.7.0 | |
1594 | |
1595 """ | |
1596 # c is a trimmed copy | |
1597 [c] = pu.as_series([c]) | |
1598 if len(c) < 2: | |
1599 raise ValueError('Series must have maximum degree of at least 1.') | |
1600 if len(c) == 2: | |
1601 return np.array([[-c[0]/c[1]]]) | |
1602 | |
1603 n = len(c) - 1 | |
1604 mat = np.zeros((n, n), dtype=c.dtype) | |
1605 scl = 1./np.sqrt(2*np.arange(n) + 1) | |
1606 top = mat.reshape(-1)[1::n+1] | |
1607 bot = mat.reshape(-1)[n::n+1] | |
1608 top[...] = np.arange(1, n)*scl[:n-1]*scl[1:n] | |
1609 bot[...] = top | |
1610 mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*(n/(2*n - 1)) | |
1611 return mat | |
1612 | |
1613 | |
1614 def legroots(c): | |
1615 """ | |
1616 Compute the roots of a Legendre series. | |
1617 | |
1618 Return the roots (a.k.a. "zeros") of the polynomial | |
1619 | |
1620 .. math:: p(x) = \\sum_i c[i] * L_i(x). | |
1621 | |
1622 Parameters | |
1623 ---------- | |
1624 c : 1-D array_like | |
1625 1-D array of coefficients. | |
1626 | |
1627 Returns | |
1628 ------- | |
1629 out : ndarray | |
1630 Array of the roots of the series. If all the roots are real, | |
1631 then `out` is also real, otherwise it is complex. | |
1632 | |
1633 See Also | |
1634 -------- | |
1635 polyroots, chebroots, lagroots, hermroots, hermeroots | |
1636 | |
1637 Notes | |
1638 ----- | |
1639 The root estimates are obtained as the eigenvalues of the companion | |
1640 matrix, Roots far from the origin of the complex plane may have large | |
1641 errors due to the numerical instability of the series for such values. | |
1642 Roots with multiplicity greater than 1 will also show larger errors as | |
1643 the value of the series near such points is relatively insensitive to | |
1644 errors in the roots. Isolated roots near the origin can be improved by | |
1645 a few iterations of Newton's method. | |
1646 | |
1647 The Legendre series basis polynomials aren't powers of ``x`` so the | |
1648 results of this function may seem unintuitive. | |
1649 | |
1650 Examples | |
1651 -------- | |
1652 >>> import numpy.polynomial.legendre as leg | |
1653 >>> leg.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0, all real roots | |
1654 array([-0.85099543, -0.11407192, 0.51506735]) | |
1655 | |
1656 """ | |
1657 # c is a trimmed copy | |
1658 [c] = pu.as_series([c]) | |
1659 if len(c) < 2: | |
1660 return np.array([], dtype=c.dtype) | |
1661 if len(c) == 2: | |
1662 return np.array([-c[0]/c[1]]) | |
1663 | |
1664 m = legcompanion(c) | |
1665 r = la.eigvals(m) | |
1666 r.sort() | |
1667 return r | |
1668 | |
1669 | |
1670 def leggauss(deg): | |
1671 """ | |
1672 Gauss-Legendre quadrature. | |
1673 | |
1674 Computes the sample points and weights for Gauss-Legendre quadrature. | |
1675 These sample points and weights will correctly integrate polynomials of | |
1676 degree :math:`2*deg - 1` or less over the interval :math:`[-1, 1]` with | |
1677 the weight function :math:`f(x) = 1`. | |
1678 | |
1679 Parameters | |
1680 ---------- | |
1681 deg : int | |
1682 Number of sample points and weights. It must be >= 1. | |
1683 | |
1684 Returns | |
1685 ------- | |
1686 x : ndarray | |
1687 1-D ndarray containing the sample points. | |
1688 y : ndarray | |
1689 1-D ndarray containing the weights. | |
1690 | |
1691 Notes | |
1692 ----- | |
1693 | |
1694 .. versionadded::1.7.0 | |
1695 | |
1696 The results have only been tested up to degree 100, higher degrees may | |
1697 be problematic. The weights are determined by using the fact that | |
1698 | |
1699 .. math:: w_k = c / (L'_n(x_k) * L_{n-1}(x_k)) | |
1700 | |
1701 where :math:`c` is a constant independent of :math:`k` and :math:`x_k` | |
1702 is the k'th root of :math:`L_n`, and then scaling the results to get | |
1703 the right value when integrating 1. | |
1704 | |
1705 """ | |
1706 ideg = int(deg) | |
1707 if ideg != deg or ideg < 1: | |
1708 raise ValueError("deg must be a non-negative integer") | |
1709 | |
1710 # first approximation of roots. We use the fact that the companion | |
1711 # matrix is symmetric in this case in order to obtain better zeros. | |
1712 c = np.array([0]*deg + [1]) | |
1713 m = legcompanion(c) | |
1714 x = la.eigvals(m) | |
1715 x.sort() | |
1716 | |
1717 # improve roots by one application of Newton | |
1718 dy = legval(x, c) | |
1719 df = legval(x, legder(c)) | |
1720 x -= dy/df | |
1721 | |
1722 # compute the weights. We scale the factor to avoid possible numerical | |
1723 # overflow. | |
1724 fm = legval(x, c[1:]) | |
1725 fm /= np.abs(fm).max() | |
1726 df /= np.abs(df).max() | |
1727 w = 1/(fm * df) | |
1728 | |
1729 # for Legendre we can also symmetrize | |
1730 w = (w + w[::-1])/2 | |
1731 x = (x - x[::-1])/2 | |
1732 | |
1733 # scale w to get the right value | |
1734 w *= 2. / w.sum() | |
1735 | |
1736 return x, w | |
1737 | |
1738 | |
1739 def legweight(x): | |
1740 """ | |
1741 Weight function of the Legendre polynomials. | |
1742 | |
1743 The weight function is :math:`1` and the interval of integration is | |
1744 :math:`[-1, 1]`. The Legendre polynomials are orthogonal, but not | |
1745 normalized, with respect to this weight function. | |
1746 | |
1747 Parameters | |
1748 ---------- | |
1749 x : array_like | |
1750 Values at which the weight function will be computed. | |
1751 | |
1752 Returns | |
1753 ------- | |
1754 w : ndarray | |
1755 The weight function at `x`. | |
1756 | |
1757 Notes | |
1758 ----- | |
1759 | |
1760 .. versionadded::1.7.0 | |
1761 | |
1762 """ | |
1763 w = x*0.0 + 1.0 | |
1764 return w | |
1765 | |
1766 # | |
1767 # Legendre series class | |
1768 # | |
1769 | |
1770 class Legendre(ABCPolyBase): | |
1771 """A Legendre series class. | |
1772 | |
1773 The Legendre class provides the standard Python numerical methods | |
1774 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the | |
1775 attributes and methods listed in the `ABCPolyBase` documentation. | |
1776 | |
1777 Parameters | |
1778 ---------- | |
1779 coef : array_like | |
1780 Legendre coefficients in order of increasing degree, i.e., | |
1781 ``(1, 2, 3)`` gives ``1*P_0(x) + 2*P_1(x) + 3*P_2(x)``. | |
1782 domain : (2,) array_like, optional | |
1783 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped | |
1784 to the interval ``[window[0], window[1]]`` by shifting and scaling. | |
1785 The default value is [-1, 1]. | |
1786 window : (2,) array_like, optional | |
1787 Window, see `domain` for its use. The default value is [-1, 1]. | |
1788 | |
1789 .. versionadded:: 1.6.0 | |
1790 | |
1791 """ | |
1792 # Virtual Functions | |
1793 _add = staticmethod(legadd) | |
1794 _sub = staticmethod(legsub) | |
1795 _mul = staticmethod(legmul) | |
1796 _div = staticmethod(legdiv) | |
1797 _pow = staticmethod(legpow) | |
1798 _val = staticmethod(legval) | |
1799 _int = staticmethod(legint) | |
1800 _der = staticmethod(legder) | |
1801 _fit = staticmethod(legfit) | |
1802 _line = staticmethod(legline) | |
1803 _roots = staticmethod(legroots) | |
1804 _fromroots = staticmethod(legfromroots) | |
1805 | |
1806 # Virtual properties | |
1807 nickname = 'leg' | |
1808 domain = np.array(legdomain) | |
1809 window = np.array(legdomain) |