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)