Coverage for /var/srv/projects/api.amasfac.comuna18.com/tmp/venv/lib/python3.9/site-packages/numpy/polynomial/legendre.py: 16%
261 statements
« prev ^ index » next coverage.py v6.4.4, created at 2023-07-17 14:22 -0600
1"""
2==================================================
3Legendre Series (:mod:`numpy.polynomial.legendre`)
4==================================================
6This module provides a number of objects (mostly functions) useful for
7dealing with Legendre series, including a `Legendre` class that
8encapsulates the usual arithmetic operations. (General information
9on how this module represents and works with such polynomials is in the
10docstring for its "parent" sub-package, `numpy.polynomial`).
12Classes
13-------
14.. autosummary::
15 :toctree: generated/
17 Legendre
19Constants
20---------
22.. autosummary::
23 :toctree: generated/
25 legdomain
26 legzero
27 legone
28 legx
30Arithmetic
31----------
33.. autosummary::
34 :toctree: generated/
36 legadd
37 legsub
38 legmulx
39 legmul
40 legdiv
41 legpow
42 legval
43 legval2d
44 legval3d
45 leggrid2d
46 leggrid3d
48Calculus
49--------
51.. autosummary::
52 :toctree: generated/
54 legder
55 legint
57Misc Functions
58--------------
60.. autosummary::
61 :toctree: generated/
63 legfromroots
64 legroots
65 legvander
66 legvander2d
67 legvander3d
68 leggauss
69 legweight
70 legcompanion
71 legfit
72 legtrim
73 legline
74 leg2poly
75 poly2leg
77See also
78--------
79numpy.polynomial
81"""
82import numpy as np
83import numpy.linalg as la
84from numpy.core.multiarray import normalize_axis_index
86from . import polyutils as pu
87from ._polybase import ABCPolyBase
89__all__ = [
90 'legzero', 'legone', 'legx', 'legdomain', 'legline', 'legadd',
91 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow', 'legval', 'legder',
92 'legint', 'leg2poly', 'poly2leg', 'legfromroots', 'legvander',
93 'legfit', 'legtrim', 'legroots', 'Legendre', 'legval2d', 'legval3d',
94 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d', 'legcompanion',
95 'leggauss', 'legweight']
97legtrim = pu.trimcoef
100def poly2leg(pol):
101 """
102 Convert a polynomial to a Legendre series.
104 Convert an array representing the coefficients of a polynomial (relative
105 to the "standard" basis) ordered from lowest degree to highest, to an
106 array of the coefficients of the equivalent Legendre series, ordered
107 from lowest to highest degree.
109 Parameters
110 ----------
111 pol : array_like
112 1-D array containing the polynomial coefficients
114 Returns
115 -------
116 c : ndarray
117 1-D array containing the coefficients of the equivalent Legendre
118 series.
120 See Also
121 --------
122 leg2poly
124 Notes
125 -----
126 The easy way to do conversions between polynomial basis sets
127 is to use the convert method of a class instance.
129 Examples
130 --------
131 >>> from numpy import polynomial as P
132 >>> p = P.Polynomial(np.arange(4))
133 >>> p
134 Polynomial([0., 1., 2., 3.], domain=[-1, 1], window=[-1, 1])
135 >>> c = P.Legendre(P.legendre.poly2leg(p.coef))
136 >>> c
137 Legendre([ 1. , 3.25, 1. , 0.75], domain=[-1, 1], window=[-1, 1]) # may vary
139 """
140 [pol] = pu.as_series([pol])
141 deg = len(pol) - 1
142 res = 0
143 for i in range(deg, -1, -1):
144 res = legadd(legmulx(res), pol[i])
145 return res
148def leg2poly(c):
149 """
150 Convert a Legendre series to a polynomial.
152 Convert an array representing the coefficients of a Legendre series,
153 ordered from lowest degree to highest, to an array of the coefficients
154 of the equivalent polynomial (relative to the "standard" basis) ordered
155 from lowest to highest degree.
157 Parameters
158 ----------
159 c : array_like
160 1-D array containing the Legendre series coefficients, ordered
161 from lowest order term to highest.
163 Returns
164 -------
165 pol : ndarray
166 1-D array containing the coefficients of the equivalent polynomial
167 (relative to the "standard" basis) ordered from lowest order term
168 to highest.
170 See Also
171 --------
172 poly2leg
174 Notes
175 -----
176 The easy way to do conversions between polynomial basis sets
177 is to use the convert method of a class instance.
179 Examples
180 --------
181 >>> from numpy import polynomial as P
182 >>> c = P.Legendre(range(4))
183 >>> c
184 Legendre([0., 1., 2., 3.], domain=[-1, 1], window=[-1, 1])
185 >>> p = c.convert(kind=P.Polynomial)
186 >>> p
187 Polynomial([-1. , -3.5, 3. , 7.5], domain=[-1., 1.], window=[-1., 1.])
188 >>> P.legendre.leg2poly(range(4))
189 array([-1. , -3.5, 3. , 7.5])
192 """
193 from .polynomial import polyadd, polysub, polymulx
195 [c] = pu.as_series([c])
196 n = len(c)
197 if n < 3:
198 return c
199 else:
200 c0 = c[-2]
201 c1 = c[-1]
202 # i is the current degree of c1
203 for i in range(n - 1, 1, -1):
204 tmp = c0
205 c0 = polysub(c[i - 2], (c1*(i - 1))/i)
206 c1 = polyadd(tmp, (polymulx(c1)*(2*i - 1))/i)
207 return polyadd(c0, polymulx(c1))
209#
210# These are constant arrays are of integer type so as to be compatible
211# with the widest range of other types, such as Decimal.
212#
214# Legendre
215legdomain = np.array([-1, 1])
217# Legendre coefficients representing zero.
218legzero = np.array([0])
220# Legendre coefficients representing one.
221legone = np.array([1])
223# Legendre coefficients representing the identity x.
224legx = np.array([0, 1])
227def legline(off, scl):
228 """
229 Legendre series whose graph is a straight line.
233 Parameters
234 ----------
235 off, scl : scalars
236 The specified line is given by ``off + scl*x``.
238 Returns
239 -------
240 y : ndarray
241 This module's representation of the Legendre series for
242 ``off + scl*x``.
244 See Also
245 --------
246 numpy.polynomial.polynomial.polyline
247 numpy.polynomial.chebyshev.chebline
248 numpy.polynomial.laguerre.lagline
249 numpy.polynomial.hermite.hermline
250 numpy.polynomial.hermite_e.hermeline
252 Examples
253 --------
254 >>> import numpy.polynomial.legendre as L
255 >>> L.legline(3,2)
256 array([3, 2])
257 >>> L.legval(-3, L.legline(3,2)) # should be -3
258 -3.0
260 """
261 if scl != 0:
262 return np.array([off, scl])
263 else:
264 return np.array([off])
267def legfromroots(roots):
268 """
269 Generate a Legendre series with given roots.
271 The function returns the coefficients of the polynomial
273 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
275 in Legendre form, where the `r_n` are the roots specified in `roots`.
276 If a zero has multiplicity n, then it must appear in `roots` n times.
277 For instance, if 2 is a root of multiplicity three and 3 is a root of
278 multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The
279 roots can appear in any order.
281 If the returned coefficients are `c`, then
283 .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x)
285 The coefficient of the last term is not generally 1 for monic
286 polynomials in Legendre form.
288 Parameters
289 ----------
290 roots : array_like
291 Sequence containing the roots.
293 Returns
294 -------
295 out : ndarray
296 1-D array of coefficients. If all roots are real then `out` is a
297 real array, if some of the roots are complex, then `out` is complex
298 even if all the coefficients in the result are real (see Examples
299 below).
301 See Also
302 --------
303 numpy.polynomial.polynomial.polyfromroots
304 numpy.polynomial.chebyshev.chebfromroots
305 numpy.polynomial.laguerre.lagfromroots
306 numpy.polynomial.hermite.hermfromroots
307 numpy.polynomial.hermite_e.hermefromroots
309 Examples
310 --------
311 >>> import numpy.polynomial.legendre as L
312 >>> L.legfromroots((-1,0,1)) # x^3 - x relative to the standard basis
313 array([ 0. , -0.4, 0. , 0.4])
314 >>> j = complex(0,1)
315 >>> L.legfromroots((-j,j)) # x^2 + 1 relative to the standard basis
316 array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) # may vary
318 """
319 return pu._fromroots(legline, legmul, roots)
322def legadd(c1, c2):
323 """
324 Add one Legendre series to another.
326 Returns the sum of two Legendre series `c1` + `c2`. The arguments
327 are sequences of coefficients ordered from lowest order term to
328 highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
330 Parameters
331 ----------
332 c1, c2 : array_like
333 1-D arrays of Legendre series coefficients ordered from low to
334 high.
336 Returns
337 -------
338 out : ndarray
339 Array representing the Legendre series of their sum.
341 See Also
342 --------
343 legsub, legmulx, legmul, legdiv, legpow
345 Notes
346 -----
347 Unlike multiplication, division, etc., the sum of two Legendre series
348 is a Legendre series (without having to "reproject" the result onto
349 the basis set) so addition, just like that of "standard" polynomials,
350 is simply "component-wise."
352 Examples
353 --------
354 >>> from numpy.polynomial import legendre as L
355 >>> c1 = (1,2,3)
356 >>> c2 = (3,2,1)
357 >>> L.legadd(c1,c2)
358 array([4., 4., 4.])
360 """
361 return pu._add(c1, c2)
364def legsub(c1, c2):
365 """
366 Subtract one Legendre series from another.
368 Returns the difference of two Legendre series `c1` - `c2`. The
369 sequences of coefficients are from lowest order term to highest, i.e.,
370 [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
372 Parameters
373 ----------
374 c1, c2 : array_like
375 1-D arrays of Legendre series coefficients ordered from low to
376 high.
378 Returns
379 -------
380 out : ndarray
381 Of Legendre series coefficients representing their difference.
383 See Also
384 --------
385 legadd, legmulx, legmul, legdiv, legpow
387 Notes
388 -----
389 Unlike multiplication, division, etc., the difference of two Legendre
390 series is a Legendre series (without having to "reproject" the result
391 onto the basis set) so subtraction, just like that of "standard"
392 polynomials, is simply "component-wise."
394 Examples
395 --------
396 >>> from numpy.polynomial import legendre as L
397 >>> c1 = (1,2,3)
398 >>> c2 = (3,2,1)
399 >>> L.legsub(c1,c2)
400 array([-2., 0., 2.])
401 >>> L.legsub(c2,c1) # -C.legsub(c1,c2)
402 array([ 2., 0., -2.])
404 """
405 return pu._sub(c1, c2)
408def legmulx(c):
409 """Multiply a Legendre series by x.
411 Multiply the Legendre series `c` by x, where x is the independent
412 variable.
415 Parameters
416 ----------
417 c : array_like
418 1-D array of Legendre series coefficients ordered from low to
419 high.
421 Returns
422 -------
423 out : ndarray
424 Array representing the result of the multiplication.
426 See Also
427 --------
428 legadd, legmul, legdiv, legpow
430 Notes
431 -----
432 The multiplication uses the recursion relationship for Legendre
433 polynomials in the form
435 .. math::
437 xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1)
439 Examples
440 --------
441 >>> from numpy.polynomial import legendre as L
442 >>> L.legmulx([1,2,3])
443 array([ 0.66666667, 2.2, 1.33333333, 1.8]) # may vary
445 """
446 # c is a trimmed copy
447 [c] = pu.as_series([c])
448 # The zero series needs special treatment
449 if len(c) == 1 and c[0] == 0:
450 return c
452 prd = np.empty(len(c) + 1, dtype=c.dtype)
453 prd[0] = c[0]*0
454 prd[1] = c[0]
455 for i in range(1, len(c)):
456 j = i + 1
457 k = i - 1
458 s = i + j
459 prd[j] = (c[i]*j)/s
460 prd[k] += (c[i]*i)/s
461 return prd
464def legmul(c1, c2):
465 """
466 Multiply one Legendre series by another.
468 Returns the product of two Legendre series `c1` * `c2`. The arguments
469 are sequences of coefficients, from lowest order "term" to highest,
470 e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
472 Parameters
473 ----------
474 c1, c2 : array_like
475 1-D arrays of Legendre series coefficients ordered from low to
476 high.
478 Returns
479 -------
480 out : ndarray
481 Of Legendre series coefficients representing their product.
483 See Also
484 --------
485 legadd, legsub, legmulx, legdiv, legpow
487 Notes
488 -----
489 In general, the (polynomial) product of two C-series results in terms
490 that are not in the Legendre polynomial basis set. Thus, to express
491 the product as a Legendre series, it is necessary to "reproject" the
492 product onto said basis set, which may produce "unintuitive" (but
493 correct) results; see Examples section below.
495 Examples
496 --------
497 >>> from numpy.polynomial import legendre as L
498 >>> c1 = (1,2,3)
499 >>> c2 = (3,2)
500 >>> L.legmul(c1,c2) # multiplication requires "reprojection"
501 array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) # may vary
503 """
504 # s1, s2 are trimmed copies
505 [c1, c2] = pu.as_series([c1, c2])
507 if len(c1) > len(c2):
508 c = c2
509 xs = c1
510 else:
511 c = c1
512 xs = c2
514 if len(c) == 1:
515 c0 = c[0]*xs
516 c1 = 0
517 elif len(c) == 2:
518 c0 = c[0]*xs
519 c1 = c[1]*xs
520 else:
521 nd = len(c)
522 c0 = c[-2]*xs
523 c1 = c[-1]*xs
524 for i in range(3, len(c) + 1):
525 tmp = c0
526 nd = nd - 1
527 c0 = legsub(c[-i]*xs, (c1*(nd - 1))/nd)
528 c1 = legadd(tmp, (legmulx(c1)*(2*nd - 1))/nd)
529 return legadd(c0, legmulx(c1))
532def legdiv(c1, c2):
533 """
534 Divide one Legendre series by another.
536 Returns the quotient-with-remainder of two Legendre series
537 `c1` / `c2`. The arguments are sequences of coefficients from lowest
538 order "term" to highest, e.g., [1,2,3] represents the series
539 ``P_0 + 2*P_1 + 3*P_2``.
541 Parameters
542 ----------
543 c1, c2 : array_like
544 1-D arrays of Legendre series coefficients ordered from low to
545 high.
547 Returns
548 -------
549 quo, rem : ndarrays
550 Of Legendre series coefficients representing the quotient and
551 remainder.
553 See Also
554 --------
555 legadd, legsub, legmulx, legmul, legpow
557 Notes
558 -----
559 In general, the (polynomial) division of one Legendre series by another
560 results in quotient and remainder terms that are not in the Legendre
561 polynomial basis set. Thus, to express these results as a Legendre
562 series, it is necessary to "reproject" the results onto the Legendre
563 basis set, which may produce "unintuitive" (but correct) results; see
564 Examples section below.
566 Examples
567 --------
568 >>> from numpy.polynomial import legendre as L
569 >>> c1 = (1,2,3)
570 >>> c2 = (3,2,1)
571 >>> L.legdiv(c1,c2) # quotient "intuitive," remainder not
572 (array([3.]), array([-8., -4.]))
573 >>> c2 = (0,1,2,3)
574 >>> L.legdiv(c2,c1) # neither "intuitive"
575 (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) # may vary
577 """
578 return pu._div(legmul, c1, c2)
581def legpow(c, pow, maxpower=16):
582 """Raise a Legendre series to a power.
584 Returns the Legendre series `c` raised to the power `pow`. The
585 argument `c` is a sequence of coefficients ordered from low to high.
586 i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.``
588 Parameters
589 ----------
590 c : array_like
591 1-D array of Legendre series coefficients ordered from low to
592 high.
593 pow : integer
594 Power to which the series will be raised
595 maxpower : integer, optional
596 Maximum power allowed. This is mainly to limit growth of the series
597 to unmanageable size. Default is 16
599 Returns
600 -------
601 coef : ndarray
602 Legendre series of power.
604 See Also
605 --------
606 legadd, legsub, legmulx, legmul, legdiv
608 """
609 return pu._pow(legmul, c, pow, maxpower)
612def legder(c, m=1, scl=1, axis=0):
613 """
614 Differentiate a Legendre series.
616 Returns the Legendre series coefficients `c` differentiated `m` times
617 along `axis`. At each iteration the result is multiplied by `scl` (the
618 scaling factor is for use in a linear change of variable). The argument
619 `c` is an array of coefficients from low to high degree along each
620 axis, e.g., [1,2,3] represents the series ``1*L_0 + 2*L_1 + 3*L_2``
621 while [[1,2],[1,2]] represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) +
622 2*L_0(x)*L_1(y) + 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is
623 ``y``.
625 Parameters
626 ----------
627 c : array_like
628 Array of Legendre series coefficients. If c is multidimensional the
629 different axis correspond to different variables with the degree in
630 each axis given by the corresponding index.
631 m : int, optional
632 Number of derivatives taken, must be non-negative. (Default: 1)
633 scl : scalar, optional
634 Each differentiation is multiplied by `scl`. The end result is
635 multiplication by ``scl**m``. This is for use in a linear change of
636 variable. (Default: 1)
637 axis : int, optional
638 Axis over which the derivative is taken. (Default: 0).
640 .. versionadded:: 1.7.0
642 Returns
643 -------
644 der : ndarray
645 Legendre series of the derivative.
647 See Also
648 --------
649 legint
651 Notes
652 -----
653 In general, the result of differentiating a Legendre series does not
654 resemble the same operation on a power series. Thus the result of this
655 function may be "unintuitive," albeit correct; see Examples section
656 below.
658 Examples
659 --------
660 >>> from numpy.polynomial import legendre as L
661 >>> c = (1,2,3,4)
662 >>> L.legder(c)
663 array([ 6., 9., 20.])
664 >>> L.legder(c, 3)
665 array([60.])
666 >>> L.legder(c, scl=-1)
667 array([ -6., -9., -20.])
668 >>> L.legder(c, 2,-1)
669 array([ 9., 60.])
671 """
672 c = np.array(c, ndmin=1, copy=True)
673 if c.dtype.char in '?bBhHiIlLqQpP':
674 c = c.astype(np.double)
675 cnt = pu._deprecate_as_int(m, "the order of derivation")
676 iaxis = pu._deprecate_as_int(axis, "the axis")
677 if cnt < 0:
678 raise ValueError("The order of derivation must be non-negative")
679 iaxis = normalize_axis_index(iaxis, c.ndim)
681 if cnt == 0:
682 return c
684 c = np.moveaxis(c, iaxis, 0)
685 n = len(c)
686 if cnt >= n:
687 c = c[:1]*0
688 else:
689 for i in range(cnt):
690 n = n - 1
691 c *= scl
692 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
693 for j in range(n, 2, -1):
694 der[j - 1] = (2*j - 1)*c[j]
695 c[j - 2] += c[j]
696 if n > 1:
697 der[1] = 3*c[2]
698 der[0] = c[1]
699 c = der
700 c = np.moveaxis(c, 0, iaxis)
701 return c
704def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
705 """
706 Integrate a Legendre series.
708 Returns the Legendre series coefficients `c` integrated `m` times from
709 `lbnd` along `axis`. At each iteration the resulting series is
710 **multiplied** by `scl` and an integration constant, `k`, is added.
711 The scaling factor is for use in a linear change of variable. ("Buyer
712 beware": note that, depending on what one is doing, one may want `scl`
713 to be the reciprocal of what one might expect; for more information,
714 see the Notes section below.) The argument `c` is an array of
715 coefficients from low to high degree along each axis, e.g., [1,2,3]
716 represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]]
717 represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) +
718 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
720 Parameters
721 ----------
722 c : array_like
723 Array of Legendre series coefficients. If c is multidimensional the
724 different axis correspond to different variables with the degree in
725 each axis given by the corresponding index.
726 m : int, optional
727 Order of integration, must be positive. (Default: 1)
728 k : {[], list, scalar}, optional
729 Integration constant(s). The value of the first integral at
730 ``lbnd`` is the first value in the list, the value of the second
731 integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
732 default), all constants are set to zero. If ``m == 1``, a single
733 scalar can be given instead of a list.
734 lbnd : scalar, optional
735 The lower bound of the integral. (Default: 0)
736 scl : scalar, optional
737 Following each integration the result is *multiplied* by `scl`
738 before the integration constant is added. (Default: 1)
739 axis : int, optional
740 Axis over which the integral is taken. (Default: 0).
742 .. versionadded:: 1.7.0
744 Returns
745 -------
746 S : ndarray
747 Legendre series coefficient array of the integral.
749 Raises
750 ------
751 ValueError
752 If ``m < 0``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
753 ``np.ndim(scl) != 0``.
755 See Also
756 --------
757 legder
759 Notes
760 -----
761 Note that the result of each integration is *multiplied* by `scl`.
762 Why is this important to note? Say one is making a linear change of
763 variable :math:`u = ax + b` in an integral relative to `x`. Then
764 :math:`dx = du/a`, so one will need to set `scl` equal to
765 :math:`1/a` - perhaps not what one would have first thought.
767 Also note that, in general, the result of integrating a C-series needs
768 to be "reprojected" onto the C-series basis set. Thus, typically,
769 the result of this function is "unintuitive," albeit correct; see
770 Examples section below.
772 Examples
773 --------
774 >>> from numpy.polynomial import legendre as L
775 >>> c = (1,2,3)
776 >>> L.legint(c)
777 array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
778 >>> L.legint(c, 3)
779 array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, # may vary
780 -1.73472348e-18, 1.90476190e-02, 9.52380952e-03])
781 >>> L.legint(c, k=3)
782 array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
783 >>> L.legint(c, lbnd=-2)
784 array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
785 >>> L.legint(c, scl=2)
786 array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) # may vary
788 """
789 c = np.array(c, ndmin=1, copy=True)
790 if c.dtype.char in '?bBhHiIlLqQpP':
791 c = c.astype(np.double)
792 if not np.iterable(k):
793 k = [k]
794 cnt = pu._deprecate_as_int(m, "the order of integration")
795 iaxis = pu._deprecate_as_int(axis, "the axis")
796 if cnt < 0:
797 raise ValueError("The order of integration must be non-negative")
798 if len(k) > cnt:
799 raise ValueError("Too many integration constants")
800 if np.ndim(lbnd) != 0:
801 raise ValueError("lbnd must be a scalar.")
802 if np.ndim(scl) != 0:
803 raise ValueError("scl must be a scalar.")
804 iaxis = normalize_axis_index(iaxis, c.ndim)
806 if cnt == 0:
807 return c
809 c = np.moveaxis(c, iaxis, 0)
810 k = list(k) + [0]*(cnt - len(k))
811 for i in range(cnt):
812 n = len(c)
813 c *= scl
814 if n == 1 and np.all(c[0] == 0):
815 c[0] += k[i]
816 else:
817 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
818 tmp[0] = c[0]*0
819 tmp[1] = c[0]
820 if n > 1:
821 tmp[2] = c[1]/3
822 for j in range(2, n):
823 t = c[j]/(2*j + 1)
824 tmp[j + 1] = t
825 tmp[j - 1] -= t
826 tmp[0] += k[i] - legval(lbnd, tmp)
827 c = tmp
828 c = np.moveaxis(c, 0, iaxis)
829 return c
832def legval(x, c, tensor=True):
833 """
834 Evaluate a Legendre series at points x.
836 If `c` is of length `n + 1`, this function returns the value:
838 .. math:: p(x) = c_0 * L_0(x) + c_1 * L_1(x) + ... + c_n * L_n(x)
840 The parameter `x` is converted to an array only if it is a tuple or a
841 list, otherwise it is treated as a scalar. In either case, either `x`
842 or its elements must support multiplication and addition both with
843 themselves and with the elements of `c`.
845 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
846 `c` is multidimensional, then the shape of the result depends on the
847 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
848 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
849 scalars have shape (,).
851 Trailing zeros in the coefficients will be used in the evaluation, so
852 they should be avoided if efficiency is a concern.
854 Parameters
855 ----------
856 x : array_like, compatible object
857 If `x` is a list or tuple, it is converted to an ndarray, otherwise
858 it is left unchanged and treated as a scalar. In either case, `x`
859 or its elements must support addition and multiplication with
860 themselves and with the elements of `c`.
861 c : array_like
862 Array of coefficients ordered so that the coefficients for terms of
863 degree n are contained in c[n]. If `c` is multidimensional the
864 remaining indices enumerate multiple polynomials. In the two
865 dimensional case the coefficients may be thought of as stored in
866 the columns of `c`.
867 tensor : boolean, optional
868 If True, the shape of the coefficient array is extended with ones
869 on the right, one for each dimension of `x`. Scalars have dimension 0
870 for this action. The result is that every column of coefficients in
871 `c` is evaluated for every element of `x`. If False, `x` is broadcast
872 over the columns of `c` for the evaluation. This keyword is useful
873 when `c` is multidimensional. The default value is True.
875 .. versionadded:: 1.7.0
877 Returns
878 -------
879 values : ndarray, algebra_like
880 The shape of the return value is described above.
882 See Also
883 --------
884 legval2d, leggrid2d, legval3d, leggrid3d
886 Notes
887 -----
888 The evaluation uses Clenshaw recursion, aka synthetic division.
890 """
891 c = np.array(c, ndmin=1, copy=False)
892 if c.dtype.char in '?bBhHiIlLqQpP':
893 c = c.astype(np.double)
894 if isinstance(x, (tuple, list)):
895 x = np.asarray(x)
896 if isinstance(x, np.ndarray) and tensor:
897 c = c.reshape(c.shape + (1,)*x.ndim)
899 if len(c) == 1:
900 c0 = c[0]
901 c1 = 0
902 elif len(c) == 2:
903 c0 = c[0]
904 c1 = c[1]
905 else:
906 nd = len(c)
907 c0 = c[-2]
908 c1 = c[-1]
909 for i in range(3, len(c) + 1):
910 tmp = c0
911 nd = nd - 1
912 c0 = c[-i] - (c1*(nd - 1))/nd
913 c1 = tmp + (c1*x*(2*nd - 1))/nd
914 return c0 + c1*x
917def legval2d(x, y, c):
918 """
919 Evaluate a 2-D Legendre series at points (x, y).
921 This function returns the values:
923 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * L_i(x) * L_j(y)
925 The parameters `x` and `y` are converted to arrays only if they are
926 tuples or a lists, otherwise they are treated as a scalars and they
927 must have the same shape after conversion. In either case, either `x`
928 and `y` or their elements must support multiplication and addition both
929 with themselves and with the elements of `c`.
931 If `c` is a 1-D array a one is implicitly appended to its shape to make
932 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
934 Parameters
935 ----------
936 x, y : array_like, compatible objects
937 The two dimensional series is evaluated at the points `(x, y)`,
938 where `x` and `y` must have the same shape. If `x` or `y` is a list
939 or tuple, it is first converted to an ndarray, otherwise it is left
940 unchanged and if it isn't an ndarray it is treated as a scalar.
941 c : array_like
942 Array of coefficients ordered so that the coefficient of the term
943 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
944 dimension greater than two the remaining indices enumerate multiple
945 sets of coefficients.
947 Returns
948 -------
949 values : ndarray, compatible object
950 The values of the two dimensional Legendre series at points formed
951 from pairs of corresponding values from `x` and `y`.
953 See Also
954 --------
955 legval, leggrid2d, legval3d, leggrid3d
957 Notes
958 -----
960 .. versionadded:: 1.7.0
962 """
963 return pu._valnd(legval, c, x, y)
966def leggrid2d(x, y, c):
967 """
968 Evaluate a 2-D Legendre series on the Cartesian product of x and y.
970 This function returns the values:
972 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * L_i(a) * L_j(b)
974 where the points `(a, b)` consist of all pairs formed by taking
975 `a` from `x` and `b` from `y`. The resulting points form a grid with
976 `x` in the first dimension and `y` in the second.
978 The parameters `x` and `y` are converted to arrays only if they are
979 tuples or a lists, otherwise they are treated as a scalars. In either
980 case, either `x` and `y` or their elements must support multiplication
981 and addition both with themselves and with the elements of `c`.
983 If `c` has fewer than two dimensions, ones are implicitly appended to
984 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
985 x.shape + y.shape.
987 Parameters
988 ----------
989 x, y : array_like, compatible objects
990 The two dimensional series is evaluated at the points in the
991 Cartesian product of `x` and `y`. If `x` or `y` is a list or
992 tuple, it is first converted to an ndarray, otherwise it is left
993 unchanged and, if it isn't an ndarray, it is treated as a scalar.
994 c : array_like
995 Array of coefficients ordered so that the coefficient of the term of
996 multi-degree i,j is contained in `c[i,j]`. If `c` has dimension
997 greater than two the remaining indices enumerate multiple sets of
998 coefficients.
1000 Returns
1001 -------
1002 values : ndarray, compatible object
1003 The values of the two dimensional Chebyshev series at points in the
1004 Cartesian product of `x` and `y`.
1006 See Also
1007 --------
1008 legval, legval2d, legval3d, leggrid3d
1010 Notes
1011 -----
1013 .. versionadded:: 1.7.0
1015 """
1016 return pu._gridnd(legval, c, x, y)
1019def legval3d(x, y, z, c):
1020 """
1021 Evaluate a 3-D Legendre series at points (x, y, z).
1023 This function returns the values:
1025 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * L_i(x) * L_j(y) * L_k(z)
1027 The parameters `x`, `y`, and `z` are converted to arrays only if
1028 they are tuples or a lists, otherwise they are treated as a scalars and
1029 they must have the same shape after conversion. In either case, either
1030 `x`, `y`, and `z` or their elements must support multiplication and
1031 addition both with themselves and with the elements of `c`.
1033 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
1034 shape to make it 3-D. The shape of the result will be c.shape[3:] +
1035 x.shape.
1037 Parameters
1038 ----------
1039 x, y, z : array_like, compatible object
1040 The three dimensional series is evaluated at the points
1041 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
1042 any of `x`, `y`, or `z` is a list or tuple, it is first converted
1043 to an ndarray, otherwise it is left unchanged and if it isn't an
1044 ndarray it is treated as a scalar.
1045 c : array_like
1046 Array of coefficients ordered so that the coefficient of the term of
1047 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
1048 greater than 3 the remaining indices enumerate multiple sets of
1049 coefficients.
1051 Returns
1052 -------
1053 values : ndarray, compatible object
1054 The values of the multidimensional polynomial on points formed with
1055 triples of corresponding values from `x`, `y`, and `z`.
1057 See Also
1058 --------
1059 legval, legval2d, leggrid2d, leggrid3d
1061 Notes
1062 -----
1064 .. versionadded:: 1.7.0
1066 """
1067 return pu._valnd(legval, c, x, y, z)
1070def leggrid3d(x, y, z, c):
1071 """
1072 Evaluate a 3-D Legendre series on the Cartesian product of x, y, and z.
1074 This function returns the values:
1076 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * L_i(a) * L_j(b) * L_k(c)
1078 where the points `(a, b, c)` consist of all triples formed by taking
1079 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1080 a grid with `x` in the first dimension, `y` in the second, and `z` in
1081 the third.
1083 The parameters `x`, `y`, and `z` are converted to arrays only if they
1084 are tuples or a lists, otherwise they are treated as a scalars. In
1085 either case, either `x`, `y`, and `z` or their elements must support
1086 multiplication and addition both with themselves and with the elements
1087 of `c`.
1089 If `c` has fewer than three dimensions, ones are implicitly appended to
1090 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1091 x.shape + y.shape + z.shape.
1093 Parameters
1094 ----------
1095 x, y, z : array_like, compatible objects
1096 The three dimensional series is evaluated at the points in the
1097 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
1098 list or tuple, it is first converted to an ndarray, otherwise it is
1099 left unchanged and, if it isn't an ndarray, it is treated as a
1100 scalar.
1101 c : array_like
1102 Array of coefficients ordered so that the coefficients for terms of
1103 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1104 greater than two the remaining indices enumerate multiple sets of
1105 coefficients.
1107 Returns
1108 -------
1109 values : ndarray, compatible object
1110 The values of the two dimensional polynomial at points in the Cartesian
1111 product of `x` and `y`.
1113 See Also
1114 --------
1115 legval, legval2d, leggrid2d, legval3d
1117 Notes
1118 -----
1120 .. versionadded:: 1.7.0
1122 """
1123 return pu._gridnd(legval, c, x, y, z)
1126def legvander(x, deg):
1127 """Pseudo-Vandermonde matrix of given degree.
1129 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
1130 `x`. The pseudo-Vandermonde matrix is defined by
1132 .. math:: V[..., i] = L_i(x)
1134 where `0 <= i <= deg`. The leading indices of `V` index the elements of
1135 `x` and the last index is the degree of the Legendre polynomial.
1137 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
1138 array ``V = legvander(x, n)``, then ``np.dot(V, c)`` and
1139 ``legval(x, c)`` are the same up to roundoff. This equivalence is
1140 useful both for least squares fitting and for the evaluation of a large
1141 number of Legendre series of the same degree and sample points.
1143 Parameters
1144 ----------
1145 x : array_like
1146 Array of points. The dtype is converted to float64 or complex128
1147 depending on whether any of the elements are complex. If `x` is
1148 scalar it is converted to a 1-D array.
1149 deg : int
1150 Degree of the resulting matrix.
1152 Returns
1153 -------
1154 vander : ndarray
1155 The pseudo-Vandermonde matrix. The shape of the returned matrix is
1156 ``x.shape + (deg + 1,)``, where The last index is the degree of the
1157 corresponding Legendre polynomial. The dtype will be the same as
1158 the converted `x`.
1160 """
1161 ideg = pu._deprecate_as_int(deg, "deg")
1162 if ideg < 0:
1163 raise ValueError("deg must be non-negative")
1165 x = np.array(x, copy=False, ndmin=1) + 0.0
1166 dims = (ideg + 1,) + x.shape
1167 dtyp = x.dtype
1168 v = np.empty(dims, dtype=dtyp)
1169 # Use forward recursion to generate the entries. This is not as accurate
1170 # as reverse recursion in this application but it is more efficient.
1171 v[0] = x*0 + 1
1172 if ideg > 0:
1173 v[1] = x
1174 for i in range(2, ideg + 1):
1175 v[i] = (v[i-1]*x*(2*i - 1) - v[i-2]*(i - 1))/i
1176 return np.moveaxis(v, 0, -1)
1179def legvander2d(x, y, deg):
1180 """Pseudo-Vandermonde matrix of given degrees.
1182 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1183 points `(x, y)`. The pseudo-Vandermonde matrix is defined by
1185 .. math:: V[..., (deg[1] + 1)*i + j] = L_i(x) * L_j(y),
1187 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
1188 `V` index the points `(x, y)` and the last index encodes the degrees of
1189 the Legendre polynomials.
1191 If ``V = legvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1192 correspond to the elements of a 2-D coefficient array `c` of shape
1193 (xdeg + 1, ydeg + 1) in the order
1195 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
1197 and ``np.dot(V, c.flat)`` and ``legval2d(x, y, c)`` will be the same
1198 up to roundoff. This equivalence is useful both for least squares
1199 fitting and for the evaluation of a large number of 2-D Legendre
1200 series of the same degrees and sample points.
1202 Parameters
1203 ----------
1204 x, y : array_like
1205 Arrays of point coordinates, all of the same shape. The dtypes
1206 will be converted to either float64 or complex128 depending on
1207 whether any of the elements are complex. Scalars are converted to
1208 1-D arrays.
1209 deg : list of ints
1210 List of maximum degrees of the form [x_deg, y_deg].
1212 Returns
1213 -------
1214 vander2d : ndarray
1215 The shape of the returned matrix is ``x.shape + (order,)``, where
1216 :math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same
1217 as the converted `x` and `y`.
1219 See Also
1220 --------
1221 legvander, legvander3d, legval2d, legval3d
1223 Notes
1224 -----
1226 .. versionadded:: 1.7.0
1228 """
1229 return pu._vander_nd_flat((legvander, legvander), (x, y), deg)
1232def legvander3d(x, y, z, deg):
1233 """Pseudo-Vandermonde matrix of given degrees.
1235 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1236 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
1237 then The pseudo-Vandermonde matrix is defined by
1239 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = L_i(x)*L_j(y)*L_k(z),
1241 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
1242 indices of `V` index the points `(x, y, z)` and the last index encodes
1243 the degrees of the Legendre polynomials.
1245 If ``V = legvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
1246 of `V` correspond to the elements of a 3-D coefficient array `c` of
1247 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
1249 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
1251 and ``np.dot(V, c.flat)`` and ``legval3d(x, y, z, c)`` will be the
1252 same up to roundoff. This equivalence is useful both for least squares
1253 fitting and for the evaluation of a large number of 3-D Legendre
1254 series of the same degrees and sample points.
1256 Parameters
1257 ----------
1258 x, y, z : array_like
1259 Arrays of point coordinates, all of the same shape. The dtypes will
1260 be converted to either float64 or complex128 depending on whether
1261 any of the elements are complex. Scalars are converted to 1-D
1262 arrays.
1263 deg : list of ints
1264 List of maximum degrees of the form [x_deg, y_deg, z_deg].
1266 Returns
1267 -------
1268 vander3d : ndarray
1269 The shape of the returned matrix is ``x.shape + (order,)``, where
1270 :math:`order = (deg[0]+1)*(deg[1]+1)*(deg[2]+1)`. The dtype will
1271 be the same as the converted `x`, `y`, and `z`.
1273 See Also
1274 --------
1275 legvander, legvander3d, legval2d, legval3d
1277 Notes
1278 -----
1280 .. versionadded:: 1.7.0
1282 """
1283 return pu._vander_nd_flat((legvander, legvander, legvander), (x, y, z), deg)
1286def legfit(x, y, deg, rcond=None, full=False, w=None):
1287 """
1288 Least squares fit of Legendre series to data.
1290 Return the coefficients of a Legendre series of degree `deg` that is the
1291 least squares fit to the data values `y` given at points `x`. If `y` is
1292 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
1293 fits are done, one for each column of `y`, and the resulting
1294 coefficients are stored in the corresponding columns of a 2-D return.
1295 The fitted polynomial(s) are in the form
1297 .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x),
1299 where `n` is `deg`.
1301 Parameters
1302 ----------
1303 x : array_like, shape (M,)
1304 x-coordinates of the M sample points ``(x[i], y[i])``.
1305 y : array_like, shape (M,) or (M, K)
1306 y-coordinates of the sample points. Several data sets of sample
1307 points sharing the same x-coordinates can be fitted at once by
1308 passing in a 2D-array that contains one dataset per column.
1309 deg : int or 1-D array_like
1310 Degree(s) of the fitting polynomials. If `deg` is a single integer
1311 all terms up to and including the `deg`'th term are included in the
1312 fit. For NumPy versions >= 1.11.0 a list of integers specifying the
1313 degrees of the terms to include may be used instead.
1314 rcond : float, optional
1315 Relative condition number of the fit. Singular values smaller than
1316 this relative to the largest singular value will be ignored. The
1317 default value is len(x)*eps, where eps is the relative precision of
1318 the float type, about 2e-16 in most cases.
1319 full : bool, optional
1320 Switch determining nature of return value. When it is False (the
1321 default) just the coefficients are returned, when True diagnostic
1322 information from the singular value decomposition is also returned.
1323 w : array_like, shape (`M`,), optional
1324 Weights. If not None, the weight ``w[i]`` applies to the unsquared
1325 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
1326 chosen so that the errors of the products ``w[i]*y[i]`` all have the
1327 same variance. When using inverse-variance weighting, use
1328 ``w[i] = 1/sigma(y[i])``. The default value is None.
1330 .. versionadded:: 1.5.0
1332 Returns
1333 -------
1334 coef : ndarray, shape (M,) or (M, K)
1335 Legendre coefficients ordered from low to high. If `y` was
1336 2-D, the coefficients for the data in column k of `y` are in
1337 column `k`. If `deg` is specified as a list, coefficients for
1338 terms not included in the fit are set equal to zero in the
1339 returned `coef`.
1341 [residuals, rank, singular_values, rcond] : list
1342 These values are only returned if ``full == True``
1344 - residuals -- sum of squared residuals of the least squares fit
1345 - rank -- the numerical rank of the scaled Vandermonde matrix
1346 - singular_values -- singular values of the scaled Vandermonde matrix
1347 - rcond -- value of `rcond`.
1349 For more details, see `numpy.linalg.lstsq`.
1351 Warns
1352 -----
1353 RankWarning
1354 The rank of the coefficient matrix in the least-squares fit is
1355 deficient. The warning is only raised if ``full == False``. The
1356 warnings can be turned off by
1358 >>> import warnings
1359 >>> warnings.simplefilter('ignore', np.RankWarning)
1361 See Also
1362 --------
1363 numpy.polynomial.polynomial.polyfit
1364 numpy.polynomial.chebyshev.chebfit
1365 numpy.polynomial.laguerre.lagfit
1366 numpy.polynomial.hermite.hermfit
1367 numpy.polynomial.hermite_e.hermefit
1368 legval : Evaluates a Legendre series.
1369 legvander : Vandermonde matrix of Legendre series.
1370 legweight : Legendre weight function (= 1).
1371 numpy.linalg.lstsq : Computes a least-squares fit from the matrix.
1372 scipy.interpolate.UnivariateSpline : Computes spline fits.
1374 Notes
1375 -----
1376 The solution is the coefficients of the Legendre series `p` that
1377 minimizes the sum of the weighted squared errors
1379 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
1381 where :math:`w_j` are the weights. This problem is solved by setting up
1382 as the (typically) overdetermined matrix equation
1384 .. math:: V(x) * c = w * y,
1386 where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
1387 coefficients to be solved for, `w` are the weights, and `y` are the
1388 observed values. This equation is then solved using the singular value
1389 decomposition of `V`.
1391 If some of the singular values of `V` are so small that they are
1392 neglected, then a `RankWarning` will be issued. This means that the
1393 coefficient values may be poorly determined. Using a lower order fit
1394 will usually get rid of the warning. The `rcond` parameter can also be
1395 set to a value smaller than its default, but the resulting fit may be
1396 spurious and have large contributions from roundoff error.
1398 Fits using Legendre series are usually better conditioned than fits
1399 using power series, but much can depend on the distribution of the
1400 sample points and the smoothness of the data. If the quality of the fit
1401 is inadequate splines may be a good alternative.
1403 References
1404 ----------
1405 .. [1] Wikipedia, "Curve fitting",
1406 https://en.wikipedia.org/wiki/Curve_fitting
1408 Examples
1409 --------
1411 """
1412 return pu._fit(legvander, x, y, deg, rcond, full, w)
1415def legcompanion(c):
1416 """Return the scaled companion matrix of c.
1418 The basis polynomials are scaled so that the companion matrix is
1419 symmetric when `c` is an Legendre basis polynomial. This provides
1420 better eigenvalue estimates than the unscaled case and for basis
1421 polynomials the eigenvalues are guaranteed to be real if
1422 `numpy.linalg.eigvalsh` is used to obtain them.
1424 Parameters
1425 ----------
1426 c : array_like
1427 1-D array of Legendre series coefficients ordered from low to high
1428 degree.
1430 Returns
1431 -------
1432 mat : ndarray
1433 Scaled companion matrix of dimensions (deg, deg).
1435 Notes
1436 -----
1438 .. versionadded:: 1.7.0
1440 """
1441 # c is a trimmed copy
1442 [c] = pu.as_series([c])
1443 if len(c) < 2:
1444 raise ValueError('Series must have maximum degree of at least 1.')
1445 if len(c) == 2:
1446 return np.array([[-c[0]/c[1]]])
1448 n = len(c) - 1
1449 mat = np.zeros((n, n), dtype=c.dtype)
1450 scl = 1./np.sqrt(2*np.arange(n) + 1)
1451 top = mat.reshape(-1)[1::n+1]
1452 bot = mat.reshape(-1)[n::n+1]
1453 top[...] = np.arange(1, n)*scl[:n-1]*scl[1:n]
1454 bot[...] = top
1455 mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*(n/(2*n - 1))
1456 return mat
1459def legroots(c):
1460 """
1461 Compute the roots of a Legendre series.
1463 Return the roots (a.k.a. "zeros") of the polynomial
1465 .. math:: p(x) = \\sum_i c[i] * L_i(x).
1467 Parameters
1468 ----------
1469 c : 1-D array_like
1470 1-D array of coefficients.
1472 Returns
1473 -------
1474 out : ndarray
1475 Array of the roots of the series. If all the roots are real,
1476 then `out` is also real, otherwise it is complex.
1478 See Also
1479 --------
1480 numpy.polynomial.polynomial.polyroots
1481 numpy.polynomial.chebyshev.chebroots
1482 numpy.polynomial.laguerre.lagroots
1483 numpy.polynomial.hermite.hermroots
1484 numpy.polynomial.hermite_e.hermeroots
1486 Notes
1487 -----
1488 The root estimates are obtained as the eigenvalues of the companion
1489 matrix, Roots far from the origin of the complex plane may have large
1490 errors due to the numerical instability of the series for such values.
1491 Roots with multiplicity greater than 1 will also show larger errors as
1492 the value of the series near such points is relatively insensitive to
1493 errors in the roots. Isolated roots near the origin can be improved by
1494 a few iterations of Newton's method.
1496 The Legendre series basis polynomials aren't powers of ``x`` so the
1497 results of this function may seem unintuitive.
1499 Examples
1500 --------
1501 >>> import numpy.polynomial.legendre as leg
1502 >>> leg.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0, all real roots
1503 array([-0.85099543, -0.11407192, 0.51506735]) # may vary
1505 """
1506 # c is a trimmed copy
1507 [c] = pu.as_series([c])
1508 if len(c) < 2:
1509 return np.array([], dtype=c.dtype)
1510 if len(c) == 2:
1511 return np.array([-c[0]/c[1]])
1513 # rotated companion matrix reduces error
1514 m = legcompanion(c)[::-1,::-1]
1515 r = la.eigvals(m)
1516 r.sort()
1517 return r
1520def leggauss(deg):
1521 """
1522 Gauss-Legendre quadrature.
1524 Computes the sample points and weights for Gauss-Legendre quadrature.
1525 These sample points and weights will correctly integrate polynomials of
1526 degree :math:`2*deg - 1` or less over the interval :math:`[-1, 1]` with
1527 the weight function :math:`f(x) = 1`.
1529 Parameters
1530 ----------
1531 deg : int
1532 Number of sample points and weights. It must be >= 1.
1534 Returns
1535 -------
1536 x : ndarray
1537 1-D ndarray containing the sample points.
1538 y : ndarray
1539 1-D ndarray containing the weights.
1541 Notes
1542 -----
1544 .. versionadded:: 1.7.0
1546 The results have only been tested up to degree 100, higher degrees may
1547 be problematic. The weights are determined by using the fact that
1549 .. math:: w_k = c / (L'_n(x_k) * L_{n-1}(x_k))
1551 where :math:`c` is a constant independent of :math:`k` and :math:`x_k`
1552 is the k'th root of :math:`L_n`, and then scaling the results to get
1553 the right value when integrating 1.
1555 """
1556 ideg = pu._deprecate_as_int(deg, "deg")
1557 if ideg <= 0:
1558 raise ValueError("deg must be a positive integer")
1560 # first approximation of roots. We use the fact that the companion
1561 # matrix is symmetric in this case in order to obtain better zeros.
1562 c = np.array([0]*deg + [1])
1563 m = legcompanion(c)
1564 x = la.eigvalsh(m)
1566 # improve roots by one application of Newton
1567 dy = legval(x, c)
1568 df = legval(x, legder(c))
1569 x -= dy/df
1571 # compute the weights. We scale the factor to avoid possible numerical
1572 # overflow.
1573 fm = legval(x, c[1:])
1574 fm /= np.abs(fm).max()
1575 df /= np.abs(df).max()
1576 w = 1/(fm * df)
1578 # for Legendre we can also symmetrize
1579 w = (w + w[::-1])/2
1580 x = (x - x[::-1])/2
1582 # scale w to get the right value
1583 w *= 2. / w.sum()
1585 return x, w
1588def legweight(x):
1589 """
1590 Weight function of the Legendre polynomials.
1592 The weight function is :math:`1` and the interval of integration is
1593 :math:`[-1, 1]`. The Legendre polynomials are orthogonal, but not
1594 normalized, with respect to this weight function.
1596 Parameters
1597 ----------
1598 x : array_like
1599 Values at which the weight function will be computed.
1601 Returns
1602 -------
1603 w : ndarray
1604 The weight function at `x`.
1606 Notes
1607 -----
1609 .. versionadded:: 1.7.0
1611 """
1612 w = x*0.0 + 1.0
1613 return w
1615#
1616# Legendre series class
1617#
1619class Legendre(ABCPolyBase):
1620 """A Legendre series class.
1622 The Legendre class provides the standard Python numerical methods
1623 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
1624 attributes and methods listed in the `ABCPolyBase` documentation.
1626 Parameters
1627 ----------
1628 coef : array_like
1629 Legendre coefficients in order of increasing degree, i.e.,
1630 ``(1, 2, 3)`` gives ``1*P_0(x) + 2*P_1(x) + 3*P_2(x)``.
1631 domain : (2,) array_like, optional
1632 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
1633 to the interval ``[window[0], window[1]]`` by shifting and scaling.
1634 The default value is [-1, 1].
1635 window : (2,) array_like, optional
1636 Window, see `domain` for its use. The default value is [-1, 1].
1638 .. versionadded:: 1.6.0
1640 """
1641 # Virtual Functions
1642 _add = staticmethod(legadd)
1643 _sub = staticmethod(legsub)
1644 _mul = staticmethod(legmul)
1645 _div = staticmethod(legdiv)
1646 _pow = staticmethod(legpow)
1647 _val = staticmethod(legval)
1648 _int = staticmethod(legint)
1649 _der = staticmethod(legder)
1650 _fit = staticmethod(legfit)
1651 _line = staticmethod(legline)
1652 _roots = staticmethod(legroots)
1653 _fromroots = staticmethod(legfromroots)
1655 # Virtual properties
1656 domain = np.array(legdomain)
1657 window = np.array(legdomain)
1658 basis_name = 'P'