codekingpro/portable-devtools
114k
1"""
2==================================================
3Legendre Series (:mod:`numpy.polynomial.legendre`)
4==================================================
5
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`).
11
12Classes
13-------
14.. autosummary::
15 :toctree: generated/
16
17 Legendre
18
19Constants
20---------
21
22.. autosummary::
23 :toctree: generated/
24
25 legdomain
26 legzero
27 legone
28 legx
29
30Arithmetic
31----------
32
33.. autosummary::
34 :toctree: generated/
35
36 legadd
37 legsub
38 legmulx
39 legmul
40 legdiv
41 legpow
42 legval
43 legval2d
44 legval3d
45 leggrid2d
46 leggrid3d
47
48Calculus
49--------
50
51.. autosummary::
52 :toctree: generated/
53
54 legder
55 legint
56
57Misc Functions
58--------------
59
60.. autosummary::
61 :toctree: generated/
62
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
76
77See also
78--------
79numpy.polynomial
80
81"""
82import numpy as np
83
84from . import polyutils as pu
85from ._polybase import ABCPolyBase
86
87__all__ = [
88 'legzero', 'legone', 'legx', 'legdomain', 'legline', 'legadd',
89 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow', 'legval', 'legder',
90 'legint', 'leg2poly', 'poly2leg', 'legfromroots', 'legvander',
91 'legfit', 'legtrim', 'legroots', 'Legendre', 'legval2d', 'legval3d',
92 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d', 'legcompanion',
93 'leggauss', 'legweight']
94
95legtrim = pu.trimcoef
96
97
98def poly2leg(pol):
99 """
100 Convert a polynomial to a Legendre series.
101
102 Convert an array representing the coefficients of a polynomial (relative
103 to the "standard" basis) ordered from lowest degree to highest, to an
104 array of the coefficients of the equivalent Legendre series, ordered
105 from lowest to highest degree.
106
107 Parameters
108 ----------
109 pol : array_like
110 1-D array containing the polynomial coefficients
111
112 Returns
113 -------
114 c : ndarray
115 1-D array containing the coefficients of the equivalent Legendre
116 series.
117
118 See Also
119 --------
120 leg2poly
121
122 Notes
123 -----
124 The easy way to do conversions between polynomial basis sets
125 is to use the convert method of a class instance.
126
127 Examples
128 --------
129 >>> import numpy as np
130 >>> from numpy import polynomial as P
131 >>> p = P.Polynomial(np.arange(4))
132 >>> p
133 Polynomial([0., 1., 2., 3.], domain=[-1., 1.], window=[-1., 1.], ...
134 >>> c = P.Legendre(P.legendre.poly2leg(p.coef))
135 >>> c
136 Legendre([ 1. , 3.25, 1. , 0.75], domain=[-1, 1], window=[-1, 1]) # may vary
137
138 """
139 [pol] = pu.as_series([pol])
140 deg = len(pol) - 1
141 res = 0
142 for i in range(deg, -1, -1):
143 res = legadd(legmulx(res), pol[i])
144 return res
145
146
147def leg2poly(c):
148 """
149 Convert a Legendre series to a polynomial.
150
151 Convert an array representing the coefficients of a Legendre series,
152 ordered from lowest degree to highest, to an array of the coefficients
153 of the equivalent polynomial (relative to the "standard" basis) ordered
154 from lowest to highest degree.
155
156 Parameters
157 ----------
158 c : array_like
159 1-D array containing the Legendre series coefficients, ordered
160 from lowest order term to highest.
161
162 Returns
163 -------
164 pol : ndarray
165 1-D array containing the coefficients of the equivalent polynomial
166 (relative to the "standard" basis) ordered from lowest order term
167 to highest.
168
169 See Also
170 --------
171 poly2leg
172
173 Notes
174 -----
175 The easy way to do conversions between polynomial basis sets
176 is to use the convert method of a class instance.
177
178 Examples
179 --------
180 >>> from numpy import polynomial as P
181 >>> c = P.Legendre(range(4))
182 >>> c
183 Legendre([0., 1., 2., 3.], domain=[-1., 1.], window=[-1., 1.], symbol='x')
184 >>> p = c.convert(kind=P.Polynomial)
185 >>> p
186 Polynomial([-1. , -3.5, 3. , 7.5], domain=[-1., 1.], window=[-1., ...
187 >>> P.legendre.leg2poly(range(4))
188 array([-1. , -3.5, 3. , 7.5])
189
190
191 """
192 from .polynomial import polyadd, polymulx, polysub
193
194 [c] = pu.as_series([c])
195 n = len(c)
196 if n < 3:
197 return c
198 else:
199 c0 = c[-2]
200 c1 = c[-1]
201 # i is the current degree of c1
202 for i in range(n - 1, 1, -1):
203 tmp = c0
204 c0 = polysub(c[i - 2], (c1 * (i - 1)) / i)
205 c1 = polyadd(tmp, (polymulx(c1) * (2 * i - 1)) / i)
206 return polyadd(c0, polymulx(c1))
207
208
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#
213
214# Legendre
215legdomain = np.array([-1., 1.])
216
217# Legendre coefficients representing zero.
218legzero = np.array([0])
219
220# Legendre coefficients representing one.
221legone = np.array([1])
222
223# Legendre coefficients representing the identity x.
224legx = np.array([0, 1])
225
226
227def legline(off, scl):
228 """
229 Legendre series whose graph is a straight line.
230
231
232
233 Parameters
234 ----------
235 off, scl : scalars
236 The specified line is given by ``off + scl*x``.
237
238 Returns
239 -------
240 y : ndarray
241 This module's representation of the Legendre series for
242 ``off + scl*x``.
243
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
251
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
259
260 """
261 if scl != 0:
262 return np.array([off, scl])
263 else:
264 return np.array([off])
265
266
267def legfromroots(roots):
268 """
269 Generate a Legendre series with given roots.
270
271 The function returns the coefficients of the polynomial
272
273 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
274
275 in Legendre form, where the :math:`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.
280
281 If the returned coefficients are `c`, then
282
283 .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x)
284
285 The coefficient of the last term is not generally 1 for monic
286 polynomials in Legendre form.
287
288 Parameters
289 ----------
290 roots : array_like
291 Sequence containing the roots.
292
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).
300
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
308
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
317
318 """
319 return pu._fromroots(legline, legmul, roots)
320
321
322def legadd(c1, c2):
323 """
324 Add one Legendre series to another.
325
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``.
329
330 Parameters
331 ----------
332 c1, c2 : array_like
333 1-D arrays of Legendre series coefficients ordered from low to
334 high.
335
336 Returns
337 -------
338 out : ndarray
339 Array representing the Legendre series of their sum.
340
341 See Also
342 --------
343 legsub, legmulx, legmul, legdiv, legpow
344
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."
351
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.])
359
360 """
361 return pu._add(c1, c2)
362
363
364def legsub(c1, c2):
365 """
366 Subtract one Legendre series from another.
367
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``.
371
372 Parameters
373 ----------
374 c1, c2 : array_like
375 1-D arrays of Legendre series coefficients ordered from low to
376 high.
377
378 Returns
379 -------
380 out : ndarray
381 Of Legendre series coefficients representing their difference.
382
383 See Also
384 --------
385 legadd, legmulx, legmul, legdiv, legpow
386
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."
393
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.])
403
404 """
405 return pu._sub(c1, c2)
406
407
408def legmulx(c):
409 """Multiply a Legendre series by x.
410
411 Multiply the Legendre series `c` by x, where x is the independent
412 variable.
413
414
415 Parameters
416 ----------
417 c : array_like
418 1-D array of Legendre series coefficients ordered from low to
419 high.
420
421 Returns
422 -------
423 out : ndarray
424 Array representing the result of the multiplication.
425
426 See Also
427 --------
428 legadd, legsub, legmul, legdiv, legpow
429
430 Notes
431 -----
432 The multiplication uses the recursion relationship for Legendre
433 polynomials in the form
434
435 .. math::
436
437 xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1)
438
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
444
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
451
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
462
463
464def legmul(c1, c2):
465 """
466 Multiply one Legendre series by another.
467
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``.
471
472 Parameters
473 ----------
474 c1, c2 : array_like
475 1-D arrays of Legendre series coefficients ordered from low to
476 high.
477
478 Returns
479 -------
480 out : ndarray
481 Of Legendre series coefficients representing their product.
482
483 See Also
484 --------
485 legadd, legsub, legmulx, legdiv, legpow
486
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.
494
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
502
503 """
504 # s1, s2 are trimmed copies
505 [c1, c2] = pu.as_series([c1, c2])
506
507 if len(c1) > len(c2):
508 c = c2
509 xs = c1
510 else:
511 c = c1
512 xs = c2
513
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))
530
531
532def legdiv(c1, c2):
533 """
534 Divide one Legendre series by another.
535
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``.
540
541 Parameters
542 ----------
543 c1, c2 : array_like
544 1-D arrays of Legendre series coefficients ordered from low to
545 high.
546
547 Returns
548 -------
549 quo, rem : ndarrays
550 Of Legendre series coefficients representing the quotient and
551 remainder.
552
553 See Also
554 --------
555 legadd, legsub, legmulx, legmul, legpow
556
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.
565
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
576
577 """
578 return pu._div(legmul, c1, c2)
579
580
581def legpow(c, pow, maxpower=16):
582 """Raise a Legendre series to a power.
583
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.``
587
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
598
599 Returns
600 -------
601 coef : ndarray
602 Legendre series of power.
603
604 See Also
605 --------
606 legadd, legsub, legmulx, legmul, legdiv
607
608 """
609 return pu._pow(legmul, c, pow, maxpower)
610
611
612def legder(c, m=1, scl=1, axis=0):
613 """
614 Differentiate a Legendre series.
615
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``.
624
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).
639
640 Returns
641 -------
642 der : ndarray
643 Legendre series of the derivative.
644
645 See Also
646 --------
647 legint
648
649 Notes
650 -----
651 In general, the result of differentiating a Legendre series does not
652 resemble the same operation on a power series. Thus the result of this
653 function may be "unintuitive," albeit correct; see Examples section
654 below.
655
656 Examples
657 --------
658 >>> from numpy.polynomial import legendre as L
659 >>> c = (1,2,3,4)
660 >>> L.legder(c)
661 array([ 6., 9., 20.])
662 >>> L.legder(c, 3)
663 array([60.])
664 >>> L.legder(c, scl=-1)
665 array([ -6., -9., -20.])
666 >>> L.legder(c, 2,-1)
667 array([ 9., 60.])
668
669 """
670 c = np.array(c, ndmin=1, copy=True)
671 if c.dtype.char in '?bBhHiIlLqQpP':
672 c = c.astype(np.double)
673 cnt = pu._as_int(m, "the order of derivation")
674 iaxis = pu._as_int(axis, "the axis")
675 if cnt < 0:
676 raise ValueError("The order of derivation must be non-negative")
677 iaxis = np.lib.array_utils.normalize_axis_index(iaxis, c.ndim)
678
679 if cnt == 0:
680 return c
681
682 c = np.moveaxis(c, iaxis, 0)
683 n = len(c)
684 if cnt >= n:
685 c = c[:1] * 0
686 else:
687 for i in range(cnt):
688 n = n - 1
689 c *= scl
690 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
691 for j in range(n, 2, -1):
692 der[j - 1] = (2 * j - 1) * c[j]
693 c[j - 2] += c[j]
694 if n > 1:
695 der[1] = 3 * c[2]
696 der[0] = c[1]
697 c = der
698 c = np.moveaxis(c, 0, iaxis)
699 return c
700
701
702def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
703 """
704 Integrate a Legendre series.
705
706 Returns the Legendre series coefficients `c` integrated `m` times from
707 `lbnd` along `axis`. At each iteration the resulting series is
708 **multiplied** by `scl` and an integration constant, `k`, is added.
709 The scaling factor is for use in a linear change of variable. ("Buyer
710 beware": note that, depending on what one is doing, one may want `scl`
711 to be the reciprocal of what one might expect; for more information,
712 see the Notes section below.) The argument `c` is an array of
713 coefficients from low to high degree along each axis, e.g., [1,2,3]
714 represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]]
715 represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) +
716 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
717
718 Parameters
719 ----------
720 c : array_like
721 Array of Legendre series coefficients. If c is multidimensional the
722 different axis correspond to different variables with the degree in
723 each axis given by the corresponding index.
724 m : int, optional
725 Order of integration, must be positive. (Default: 1)
726 k : {[], list, scalar}, optional
727 Integration constant(s). The value of the first integral at
728 ``lbnd`` is the first value in the list, the value of the second
729 integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
730 default), all constants are set to zero. If ``m == 1``, a single
731 scalar can be given instead of a list.
732 lbnd : scalar, optional
733 The lower bound of the integral. (Default: 0)
734 scl : scalar, optional
735 Following each integration the result is *multiplied* by `scl`
736 before the integration constant is added. (Default: 1)
737 axis : int, optional
738 Axis over which the integral is taken. (Default: 0).
739
740 Returns
741 -------
742 S : ndarray
743 Legendre series coefficient array of the integral.
744
745 Raises
746 ------
747 ValueError
748 If ``m < 0``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
749 ``np.ndim(scl) != 0``.
750
751 See Also
752 --------
753 legder
754
755 Notes
756 -----
757 Note that the result of each integration is *multiplied* by `scl`.
758 Why is this important to note? Say one is making a linear change of
759 variable :math:`u = ax + b` in an integral relative to `x`. Then
760 :math:`dx = du/a`, so one will need to set `scl` equal to
761 :math:`1/a` - perhaps not what one would have first thought.
762
763 Also note that, in general, the result of integrating a C-series needs
764 to be "reprojected" onto the C-series basis set. Thus, typically,
765 the result of this function is "unintuitive," albeit correct; see
766 Examples section below.
767
768 Examples
769 --------
770 >>> from numpy.polynomial import legendre as L
771 >>> c = (1,2,3)
772 >>> L.legint(c)
773 array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
774 >>> L.legint(c, 3)
775 array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, # may vary
776 -1.73472348e-18, 1.90476190e-02, 9.52380952e-03])
777 >>> L.legint(c, k=3)
778 array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
779 >>> L.legint(c, lbnd=-2)
780 array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
781 >>> L.legint(c, scl=2)
782 array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) # may vary
783
784 """
785 c = np.array(c, ndmin=1, copy=True)
786 if c.dtype.char in '?bBhHiIlLqQpP':
787 c = c.astype(np.double)
788 if not np.iterable(k):
789 k = [k]
790 cnt = pu._as_int(m, "the order of integration")
791 iaxis = pu._as_int(axis, "the axis")
792 if cnt < 0:
793 raise ValueError("The order of integration must be non-negative")
794 if len(k) > cnt:
795 raise ValueError("Too many integration constants")
796 if np.ndim(lbnd) != 0:
797 raise ValueError("lbnd must be a scalar.")
798 if np.ndim(scl) != 0:
799 raise ValueError("scl must be a scalar.")
800 iaxis = np.lib.array_utils.normalize_axis_index(iaxis, c.ndim)
801
802 if cnt == 0:
803 return c
804
805 c = np.moveaxis(c, iaxis, 0)
806 k = list(k) + [0] * (cnt - len(k))
807 for i in range(cnt):
808 n = len(c)
809 c *= scl
810 if n == 1 and np.all(c[0] == 0):
811 c[0] += k[i]
812 else:
813 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
814 tmp[0] = c[0] * 0
815 tmp[1] = c[0]
816 if n > 1:
817 tmp[2] = c[1] / 3
818 for j in range(2, n):
819 t = c[j] / (2 * j + 1)
820 tmp[j + 1] = t
821 tmp[j - 1] -= t
822 tmp[0] += k[i] - legval(lbnd, tmp)
823 c = tmp
824 c = np.moveaxis(c, 0, iaxis)
825 return c
826
827
828def legval(x, c, tensor=True):
829 """
830 Evaluate a Legendre series at points x.
831
832 If `c` is of length ``n + 1``, this function returns the value:
833
834 .. math:: p(x) = c_0 * L_0(x) + c_1 * L_1(x) + ... + c_n * L_n(x)
835
836 The parameter `x` is converted to an array only if it is a tuple or a
837 list, otherwise it is treated as a scalar. In either case, either `x`
838 or its elements must support multiplication and addition both with
839 themselves and with the elements of `c`.
840
841 If `c` is a 1-D array, then ``p(x)`` will have the same shape as `x`. If
842 `c` is multidimensional, then the shape of the result depends on the
843 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
844 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
845 scalars have shape (,).
846
847 Trailing zeros in the coefficients will be used in the evaluation, so
848 they should be avoided if efficiency is a concern.
849
850 Parameters
851 ----------
852 x : array_like, compatible object
853 If `x` is a list or tuple, it is converted to an ndarray, otherwise
854 it is left unchanged and treated as a scalar. In either case, `x`
855 or its elements must support addition and multiplication with
856 themselves and with the elements of `c`.
857 c : array_like
858 Array of coefficients ordered so that the coefficients for terms of
859 degree n are contained in c[n]. If `c` is multidimensional the
860 remaining indices enumerate multiple polynomials. In the two
861 dimensional case the coefficients may be thought of as stored in
862 the columns of `c`.
863 tensor : boolean, optional
864 If True, the shape of the coefficient array is extended with ones
865 on the right, one for each dimension of `x`. Scalars have dimension 0
866 for this action. The result is that every column of coefficients in
867 `c` is evaluated for every element of `x`. If False, `x` is broadcast
868 over the columns of `c` for the evaluation. This keyword is useful
869 when `c` is multidimensional. The default value is True.
870
871 Returns
872 -------
873 values : ndarray, algebra_like
874 The shape of the return value is described above.
875
876 See Also
877 --------
878 legval2d, leggrid2d, legval3d, leggrid3d
879
880 Notes
881 -----
882 The evaluation uses Clenshaw recursion, aka synthetic division.
883
884 """
885 c = np.array(c, ndmin=1, copy=None)
886 if c.dtype.char in '?bBhHiIlLqQpP':
887 c = c.astype(np.double)
888 if isinstance(x, (tuple, list)):
889 x = np.asarray(x)
890 if isinstance(x, np.ndarray) and tensor:
891 c = c.reshape(c.shape + (1,) * x.ndim)
892
893 if len(c) == 1:
894 c0 = c[0]
895 c1 = 0
896 elif len(c) == 2:
897 c0 = c[0]
898 c1 = c[1]
899 else:
900 nd = len(c)
901 c0 = c[-2]
902 c1 = c[-1]
903 for i in range(3, len(c) + 1):
904 tmp = c0
905 nd = nd - 1
906 c0 = c[-i] - c1 * ((nd - 1) / nd)
907 c1 = tmp + c1 * x * ((2 * nd - 1) / nd)
908 return c0 + c1 * x
909
910
911def legval2d(x, y, c):
912 """
913 Evaluate a 2-D Legendre series at points (x, y).
914
915 This function returns the values:
916
917 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * L_i(x) * L_j(y)
918
919 The parameters `x` and `y` are converted to arrays only if they are
920 tuples or a lists, otherwise they are treated as a scalars and they
921 must have the same shape after conversion. In either case, either `x`
922 and `y` or their elements must support multiplication and addition both
923 with themselves and with the elements of `c`.
924
925 If `c` is a 1-D array a one is implicitly appended to its shape to make
926 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
927
928 Parameters
929 ----------
930 x, y : array_like, compatible objects
931 The two dimensional series is evaluated at the points ``(x, y)``,
932 where `x` and `y` must have the same shape. If `x` or `y` is a list
933 or tuple, it is first converted to an ndarray, otherwise it is left
934 unchanged and if it isn't an ndarray it is treated as a scalar.
935 c : array_like
936 Array of coefficients ordered so that the coefficient of the term
937 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
938 dimension greater than two the remaining indices enumerate multiple
939 sets of coefficients.
940
941 Returns
942 -------
943 values : ndarray, compatible object
944 The values of the two dimensional Legendre series at points formed
945 from pairs of corresponding values from `x` and `y`.
946
947 See Also
948 --------
949 legval, leggrid2d, legval3d, leggrid3d
950 """
951 return pu._valnd(legval, c, x, y)
952
953
954def leggrid2d(x, y, c):
955 """
956 Evaluate a 2-D Legendre series on the Cartesian product of x and y.
957
958 This function returns the values:
959
960 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * L_i(a) * L_j(b)
961
962 where the points ``(a, b)`` consist of all pairs formed by taking
963 `a` from `x` and `b` from `y`. The resulting points form a grid with
964 `x` in the first dimension and `y` in the second.
965
966 The parameters `x` and `y` are converted to arrays only if they are
967 tuples or a lists, otherwise they are treated as a scalars. In either
968 case, either `x` and `y` or their elements must support multiplication
969 and addition both with themselves and with the elements of `c`.
970
971 If `c` has fewer than two dimensions, ones are implicitly appended to
972 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
973 x.shape + y.shape.
974
975 Parameters
976 ----------
977 x, y : array_like, compatible objects
978 The two dimensional series is evaluated at the points in the
979 Cartesian product of `x` and `y`. If `x` or `y` is a list or
980 tuple, it is first converted to an ndarray, otherwise it is left
981 unchanged and, if it isn't an ndarray, it is treated as a scalar.
982 c : array_like
983 Array of coefficients ordered so that the coefficient of the term of
984 multi-degree i,j is contained in ``c[i,j]``. If `c` has dimension
985 greater than two the remaining indices enumerate multiple sets of
986 coefficients.
987
988 Returns
989 -------
990 values : ndarray, compatible object
991 The values of the two dimensional Chebyshev series at points in the
992 Cartesian product of `x` and `y`.
993
994 See Also
995 --------
996 legval, legval2d, legval3d, leggrid3d
997 """
998 return pu._gridnd(legval, c, x, y)
999
1000
1001def legval3d(x, y, z, c):
1002 """
1003 Evaluate a 3-D Legendre series at points (x, y, z).
1004
1005 This function returns the values:
1006
1007 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * L_i(x) * L_j(y) * L_k(z)
1008
1009 The parameters `x`, `y`, and `z` are converted to arrays only if
1010 they are tuples or a lists, otherwise they are treated as a scalars and
1011 they must have the same shape after conversion. In either case, either
1012 `x`, `y`, and `z` or their elements must support multiplication and
1013 addition both with themselves and with the elements of `c`.
1014
1015 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
1016 shape to make it 3-D. The shape of the result will be c.shape[3:] +
1017 x.shape.
1018
1019 Parameters
1020 ----------
1021 x, y, z : array_like, compatible object
1022 The three dimensional series is evaluated at the points
1023 ``(x, y, z)``, where `x`, `y`, and `z` must have the same shape. If
1024 any of `x`, `y`, or `z` is a list or tuple, it is first converted
1025 to an ndarray, otherwise it is left unchanged and if it isn't an
1026 ndarray it is treated as a scalar.
1027 c : array_like
1028 Array of coefficients ordered so that the coefficient of the term of
1029 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
1030 greater than 3 the remaining indices enumerate multiple sets of
1031 coefficients.
1032
1033 Returns
1034 -------
1035 values : ndarray, compatible object
1036 The values of the multidimensional polynomial on points formed with
1037 triples of corresponding values from `x`, `y`, and `z`.
1038
1039 See Also
1040 --------
1041 legval, legval2d, leggrid2d, leggrid3d
1042 """
1043 return pu._valnd(legval, c, x, y, z)
1044
1045
1046def leggrid3d(x, y, z, c):
1047 """
1048 Evaluate a 3-D Legendre series on the Cartesian product of x, y, and z.
1049
1050 This function returns the values:
1051
1052 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * L_i(a) * L_j(b) * L_k(c)
1053
1054 where the points ``(a, b, c)`` consist of all triples formed by taking
1055 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1056 a grid with `x` in the first dimension, `y` in the second, and `z` in
1057 the third.
1058
1059 The parameters `x`, `y`, and `z` are converted to arrays only if they
1060 are tuples or a lists, otherwise they are treated as a scalars. In
1061 either case, either `x`, `y`, and `z` or their elements must support
1062 multiplication and addition both with themselves and with the elements
1063 of `c`.
1064
1065 If `c` has fewer than three dimensions, ones are implicitly appended to
1066 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1067 x.shape + y.shape + z.shape.
1068
1069 Parameters
1070 ----------
1071 x, y, z : array_like, compatible objects
1072 The three dimensional series is evaluated at the points in the
1073 Cartesian product of `x`, `y`, and `z`. If `x`, `y`, or `z` is a
1074 list or tuple, it is first converted to an ndarray, otherwise it is
1075 left unchanged and, if it isn't an ndarray, it is treated as a
1076 scalar.
1077 c : array_like
1078 Array of coefficients ordered so that the coefficients for terms of
1079 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1080 greater than two the remaining indices enumerate multiple sets of
1081 coefficients.
1082
1083 Returns
1084 -------
1085 values : ndarray, compatible object
1086 The values of the two dimensional polynomial at points in the Cartesian
1087 product of `x` and `y`.
1088
1089 See Also
1090 --------
1091 legval, legval2d, leggrid2d, legval3d
1092 """
1093 return pu._gridnd(legval, c, x, y, z)
1094
1095
1096def legvander(x, deg):
1097 """Pseudo-Vandermonde matrix of given degree.
1098
1099 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
1100 `x`. The pseudo-Vandermonde matrix is defined by
1101
1102 .. math:: V[..., i] = L_i(x)
1103
1104 where ``0 <= i <= deg``. The leading indices of `V` index the elements of
1105 `x` and the last index is the degree of the Legendre polynomial.
1106
1107 If `c` is a 1-D array of coefficients of length ``n + 1`` and `V` is the
1108 array ``V = legvander(x, n)``, then ``np.dot(V, c)`` and
1109 ``legval(x, c)`` are the same up to roundoff. This equivalence is
1110 useful both for least squares fitting and for the evaluation of a large
1111 number of Legendre series of the same degree and sample points.
1112
1113 Parameters
1114 ----------
1115 x : array_like
1116 Array of points. The dtype is converted to float64 or complex128
1117 depending on whether any of the elements are complex. If `x` is
1118 scalar it is converted to a 1-D array.
1119 deg : int
1120 Degree of the resulting matrix.
1121
1122 Returns
1123 -------
1124 vander : ndarray
1125 The pseudo-Vandermonde matrix. The shape of the returned matrix is
1126 ``x.shape + (deg + 1,)``, where The last index is the degree of the
1127 corresponding Legendre polynomial. The dtype will be the same as
1128 the converted `x`.
1129
1130 """
1131 ideg = pu._as_int(deg, "deg")
1132 if ideg < 0:
1133 raise ValueError("deg must be non-negative")
1134
1135 x = np.array(x, copy=None, ndmin=1) + 0.0
1136 dims = (ideg + 1,) + x.shape
1137 dtyp = x.dtype
1138 v = np.empty(dims, dtype=dtyp)
1139 # Use forward recursion to generate the entries. This is not as accurate
1140 # as reverse recursion in this application but it is more efficient.
1141 v[0] = x * 0 + 1
1142 if ideg > 0:
1143 v[1] = x
1144 for i in range(2, ideg + 1):
1145 v[i] = (v[i - 1] * x * (2 * i - 1) - v[i - 2] * (i - 1)) / i
1146 return np.moveaxis(v, 0, -1)
1147
1148
1149def legvander2d(x, y, deg):
1150 """Pseudo-Vandermonde matrix of given degrees.
1151
1152 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1153 points ``(x, y)``. The pseudo-Vandermonde matrix is defined by
1154
1155 .. math:: V[..., (deg[1] + 1)*i + j] = L_i(x) * L_j(y),
1156
1157 where ``0 <= i <= deg[0]`` and ``0 <= j <= deg[1]``. The leading indices of
1158 `V` index the points ``(x, y)`` and the last index encodes the degrees of
1159 the Legendre polynomials.
1160
1161 If ``V = legvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1162 correspond to the elements of a 2-D coefficient array `c` of shape
1163 (xdeg + 1, ydeg + 1) in the order
1164
1165 .. math:: c_{00}, c_{01}, c_{02}, ... , c_{10}, c_{11}, c_{12}, ...
1166
1167 and ``np.dot(V, c.flat)`` and ``legval2d(x, y, c)`` will be the same
1168 up to roundoff. This equivalence is useful both for least squares
1169 fitting and for the evaluation of a large number of 2-D Legendre
1170 series of the same degrees and sample points.
1171
1172 Parameters
1173 ----------
1174 x, y : array_like
1175 Arrays of point coordinates, all of the same shape. The dtypes
1176 will be converted to either float64 or complex128 depending on
1177 whether any of the elements are complex. Scalars are converted to
1178 1-D arrays.
1179 deg : list of ints
1180 List of maximum degrees of the form [x_deg, y_deg].
1181
1182 Returns
1183 -------
1184 vander2d : ndarray
1185 The shape of the returned matrix is ``x.shape + (order,)``, where
1186 :math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same
1187 as the converted `x` and `y`.
1188
1189 See Also
1190 --------
1191 legvander, legvander3d, legval2d, legval3d
1192 """
1193 return pu._vander_nd_flat((legvander, legvander), (x, y), deg)
1194
1195
1196def legvander3d(x, y, z, deg):
1197 """Pseudo-Vandermonde matrix of given degrees.
1198
1199 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1200 points ``(x, y, z)``. If `l`, `m`, `n` are the given degrees in `x`, `y`, `z`,