codekingpro/portable-devtools
114k
1"""
2==============================================================
3Hermite Series, "Physicists" (:mod:`numpy.polynomial.hermite`)
4==============================================================
5
6This module provides a number of objects (mostly functions) useful for
7dealing with Hermite series, including a `Hermite` 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 Hermite
18
19Constants
20---------
21.. autosummary::
22 :toctree: generated/
23
24 hermdomain
25 hermzero
26 hermone
27 hermx
28
29Arithmetic
30----------
31.. autosummary::
32 :toctree: generated/
33
34 hermadd
35 hermsub
36 hermmulx
37 hermmul
38 hermdiv
39 hermpow
40 hermval
41 hermval2d
42 hermval3d
43 hermgrid2d
44 hermgrid3d
45
46Calculus
47--------
48.. autosummary::
49 :toctree: generated/
50
51 hermder
52 hermint
53
54Misc Functions
55--------------
56.. autosummary::
57 :toctree: generated/
58
59 hermfromroots
60 hermroots
61 hermvander
62 hermvander2d
63 hermvander3d
64 hermgauss
65 hermweight
66 hermcompanion
67 hermfit
68 hermtrim
69 hermline
70 herm2poly
71 poly2herm
72
73See also
74--------
75`numpy.polynomial`
76
77"""
78import numpy as np
79
80from . import polyutils as pu
81from ._polybase import ABCPolyBase
82
83__all__ = [
84 'hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline', 'hermadd',
85 'hermsub', 'hermmulx', 'hermmul', 'hermdiv', 'hermpow', 'hermval',
86 'hermder', 'hermint', 'herm2poly', 'poly2herm', 'hermfromroots',
87 'hermvander', 'hermfit', 'hermtrim', 'hermroots', 'Hermite',
88 'hermval2d', 'hermval3d', 'hermgrid2d', 'hermgrid3d', 'hermvander2d',
89 'hermvander3d', 'hermcompanion', 'hermgauss', 'hermweight']
90
91hermtrim = pu.trimcoef
92
93
94def poly2herm(pol):
95 """
96 poly2herm(pol)
97
98 Convert a polynomial to a Hermite series.
99
100 Convert an array representing the coefficients of a polynomial (relative
101 to the "standard" basis) ordered from lowest degree to highest, to an
102 array of the coefficients of the equivalent Hermite series, ordered
103 from lowest to highest degree.
104
105 Parameters
106 ----------
107 pol : array_like
108 1-D array containing the polynomial coefficients
109
110 Returns
111 -------
112 c : ndarray
113 1-D array containing the coefficients of the equivalent Hermite
114 series.
115
116 See Also
117 --------
118 herm2poly
119
120 Notes
121 -----
122 The easy way to do conversions between polynomial basis sets
123 is to use the convert method of a class instance.
124
125 Examples
126 --------
127 >>> from numpy.polynomial.hermite import poly2herm
128 >>> poly2herm(np.arange(4))
129 array([1. , 2.75 , 0.5 , 0.375])
130
131 """
132 [pol] = pu.as_series([pol])
133 deg = len(pol) - 1
134 res = 0
135 for i in range(deg, -1, -1):
136 res = hermadd(hermmulx(res), pol[i])
137 return res
138
139
140def herm2poly(c):
141 """
142 Convert a Hermite series to a polynomial.
143
144 Convert an array representing the coefficients of a Hermite series,
145 ordered from lowest degree to highest, to an array of the coefficients
146 of the equivalent polynomial (relative to the "standard" basis) ordered
147 from lowest to highest degree.
148
149 Parameters
150 ----------
151 c : array_like
152 1-D array containing the Hermite series coefficients, ordered
153 from lowest order term to highest.
154
155 Returns
156 -------
157 pol : ndarray
158 1-D array containing the coefficients of the equivalent polynomial
159 (relative to the "standard" basis) ordered from lowest order term
160 to highest.
161
162 See Also
163 --------
164 poly2herm
165
166 Notes
167 -----
168 The easy way to do conversions between polynomial basis sets
169 is to use the convert method of a class instance.
170
171 Examples
172 --------
173 >>> from numpy.polynomial.hermite import herm2poly
174 >>> herm2poly([ 1. , 2.75 , 0.5 , 0.375])
175 array([0., 1., 2., 3.])
176
177 """
178 from .polynomial import polyadd, polymulx, polysub
179
180 [c] = pu.as_series([c])
181 n = len(c)
182 if n == 1:
183 return c
184 if n == 2:
185 c[1] *= 2
186 return c
187 else:
188 c0 = c[-2]
189 c1 = c[-1]
190 # i is the current degree of c1
191 for i in range(n - 1, 1, -1):
192 tmp = c0
193 c0 = polysub(c[i - 2], c1 * (2 * (i - 1)))
194 c1 = polyadd(tmp, polymulx(c1) * 2)
195 return polyadd(c0, polymulx(c1) * 2)
196
197
198#
199# These are constant arrays are of integer type so as to be compatible
200# with the widest range of other types, such as Decimal.
201#
202
203# Hermite
204hermdomain = np.array([-1., 1.])
205
206# Hermite coefficients representing zero.
207hermzero = np.array([0])
208
209# Hermite coefficients representing one.
210hermone = np.array([1])
211
212# Hermite coefficients representing the identity x.
213hermx = np.array([0, 1 / 2])
214
215
216def hermline(off, scl):
217 """
218 Hermite series whose graph is a straight line.
219
220
221
222 Parameters
223 ----------
224 off, scl : scalars
225 The specified line is given by ``off + scl*x``.
226
227 Returns
228 -------
229 y : ndarray
230 This module's representation of the Hermite series for
231 ``off + scl*x``.
232
233 See Also
234 --------
235 numpy.polynomial.polynomial.polyline
236 numpy.polynomial.chebyshev.chebline
237 numpy.polynomial.legendre.legline
238 numpy.polynomial.laguerre.lagline
239 numpy.polynomial.hermite_e.hermeline
240
241 Examples
242 --------
243 >>> from numpy.polynomial.hermite import hermline, hermval
244 >>> hermval(0,hermline(3, 2))
245 3.0
246 >>> hermval(1,hermline(3, 2))
247 5.0
248
249 """
250 if scl != 0:
251 return np.array([off, scl / 2])
252 else:
253 return np.array([off])
254
255
256def hermfromroots(roots):
257 """
258 Generate a Hermite series with given roots.
259
260 The function returns the coefficients of the polynomial
261
262 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
263
264 in Hermite form, where the :math:`r_n` are the roots specified in `roots`.
265 If a zero has multiplicity n, then it must appear in `roots` n times.
266 For instance, if 2 is a root of multiplicity three and 3 is a root of
267 multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The
268 roots can appear in any order.
269
270 If the returned coefficients are `c`, then
271
272 .. math:: p(x) = c_0 + c_1 * H_1(x) + ... + c_n * H_n(x)
273
274 The coefficient of the last term is not generally 1 for monic
275 polynomials in Hermite form.
276
277 Parameters
278 ----------
279 roots : array_like
280 Sequence containing the roots.
281
282 Returns
283 -------
284 out : ndarray
285 1-D array of coefficients. If all roots are real then `out` is a
286 real array, if some of the roots are complex, then `out` is complex
287 even if all the coefficients in the result are real (see Examples
288 below).
289
290 See Also
291 --------
292 numpy.polynomial.polynomial.polyfromroots
293 numpy.polynomial.legendre.legfromroots
294 numpy.polynomial.laguerre.lagfromroots
295 numpy.polynomial.chebyshev.chebfromroots
296 numpy.polynomial.hermite_e.hermefromroots
297
298 Examples
299 --------
300 >>> from numpy.polynomial.hermite import hermfromroots, hermval
301 >>> coef = hermfromroots((-1, 0, 1))
302 >>> hermval((-1, 0, 1), coef)
303 array([0., 0., 0.])
304 >>> coef = hermfromroots((-1j, 1j))
305 >>> hermval((-1j, 1j), coef)
306 array([0.+0.j, 0.+0.j])
307
308 """
309 return pu._fromroots(hermline, hermmul, roots)
310
311
312def hermadd(c1, c2):
313 """
314 Add one Hermite series to another.
315
316 Returns the sum of two Hermite series `c1` + `c2`. The arguments
317 are sequences of coefficients ordered from lowest order term to
318 highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
319
320 Parameters
321 ----------
322 c1, c2 : array_like
323 1-D arrays of Hermite series coefficients ordered from low to
324 high.
325
326 Returns
327 -------
328 out : ndarray
329 Array representing the Hermite series of their sum.
330
331 See Also
332 --------
333 hermsub, hermmulx, hermmul, hermdiv, hermpow
334
335 Notes
336 -----
337 Unlike multiplication, division, etc., the sum of two Hermite series
338 is a Hermite series (without having to "reproject" the result onto
339 the basis set) so addition, just like that of "standard" polynomials,
340 is simply "component-wise."
341
342 Examples
343 --------
344 >>> from numpy.polynomial.hermite import hermadd
345 >>> hermadd([1, 2, 3], [1, 2, 3, 4])
346 array([2., 4., 6., 4.])
347
348 """
349 return pu._add(c1, c2)
350
351
352def hermsub(c1, c2):
353 """
354 Subtract one Hermite series from another.
355
356 Returns the difference of two Hermite series `c1` - `c2`. The
357 sequences of coefficients are from lowest order term to highest, i.e.,
358 [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
359
360 Parameters
361 ----------
362 c1, c2 : array_like
363 1-D arrays of Hermite series coefficients ordered from low to
364 high.
365
366 Returns
367 -------
368 out : ndarray
369 Of Hermite series coefficients representing their difference.
370
371 See Also
372 --------
373 hermadd, hermmulx, hermmul, hermdiv, hermpow
374
375 Notes
376 -----
377 Unlike multiplication, division, etc., the difference of two Hermite
378 series is a Hermite series (without having to "reproject" the result
379 onto the basis set) so subtraction, just like that of "standard"
380 polynomials, is simply "component-wise."
381
382 Examples
383 --------
384 >>> from numpy.polynomial.hermite import hermsub
385 >>> hermsub([1, 2, 3, 4], [1, 2, 3])
386 array([0., 0., 0., 4.])
387
388 """
389 return pu._sub(c1, c2)
390
391
392def hermmulx(c):
393 """Multiply a Hermite series by x.
394
395 Multiply the Hermite series `c` by x, where x is the independent
396 variable.
397
398
399 Parameters
400 ----------
401 c : array_like
402 1-D array of Hermite series coefficients ordered from low to
403 high.
404
405 Returns
406 -------
407 out : ndarray
408 Array representing the result of the multiplication.
409
410 See Also
411 --------
412 hermadd, hermsub, hermmul, hermdiv, hermpow
413
414 Notes
415 -----
416 The multiplication uses the recursion relationship for Hermite
417 polynomials in the form
418
419 .. math::
420
421 xP_i(x) = (P_{i + 1}(x)/2 + i*P_{i - 1}(x))
422
423 Examples
424 --------
425 >>> from numpy.polynomial.hermite import hermmulx
426 >>> hermmulx([1, 2, 3])
427 array([2. , 6.5, 1. , 1.5])
428
429 """
430 # c is a trimmed copy
431 [c] = pu.as_series([c])
432 # The zero series needs special treatment
433 if len(c) == 1 and c[0] == 0:
434 return c
435
436 prd = np.empty(len(c) + 1, dtype=c.dtype)
437 prd[0] = c[0] * 0
438 prd[1] = c[0] / 2
439 for i in range(1, len(c)):
440 prd[i + 1] = c[i] / 2
441 prd[i - 1] += c[i] * i
442 return prd
443
444
445def hermmul(c1, c2):
446 """
447 Multiply one Hermite series by another.
448
449 Returns the product of two Hermite series `c1` * `c2`. The arguments
450 are sequences of coefficients, from lowest order "term" to highest,
451 e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
452
453 Parameters
454 ----------
455 c1, c2 : array_like
456 1-D arrays of Hermite series coefficients ordered from low to
457 high.
458
459 Returns
460 -------
461 out : ndarray
462 Of Hermite series coefficients representing their product.
463
464 See Also
465 --------
466 hermadd, hermsub, hermmulx, hermdiv, hermpow
467
468 Notes
469 -----
470 In general, the (polynomial) product of two C-series results in terms
471 that are not in the Hermite polynomial basis set. Thus, to express
472 the product as a Hermite series, it is necessary to "reproject" the
473 product onto said basis set, which may produce "unintuitive" (but
474 correct) results; see Examples section below.
475
476 Examples
477 --------
478 >>> from numpy.polynomial.hermite import hermmul
479 >>> hermmul([1, 2, 3], [0, 1, 2])
480 array([52., 29., 52., 7., 6.])
481
482 """
483 # s1, s2 are trimmed copies
484 [c1, c2] = pu.as_series([c1, c2])
485
486 if len(c1) > len(c2):
487 c = c2
488 xs = c1
489 else:
490 c = c1
491 xs = c2
492
493 if len(c) == 1:
494 c0 = c[0] * xs
495 c1 = 0
496 elif len(c) == 2:
497 c0 = c[0] * xs
498 c1 = c[1] * xs
499 else:
500 nd = len(c)
501 c0 = c[-2] * xs
502 c1 = c[-1] * xs
503 for i in range(3, len(c) + 1):
504 tmp = c0
505 nd = nd - 1
506 c0 = hermsub(c[-i] * xs, c1 * (2 * (nd - 1)))
507 c1 = hermadd(tmp, hermmulx(c1) * 2)
508 return hermadd(c0, hermmulx(c1) * 2)
509
510
511def hermdiv(c1, c2):
512 """
513 Divide one Hermite series by another.
514
515 Returns the quotient-with-remainder of two Hermite series
516 `c1` / `c2`. The arguments are sequences of coefficients from lowest
517 order "term" to highest, e.g., [1,2,3] represents the series
518 ``P_0 + 2*P_1 + 3*P_2``.
519
520 Parameters
521 ----------
522 c1, c2 : array_like
523 1-D arrays of Hermite series coefficients ordered from low to
524 high.
525
526 Returns
527 -------
528 [quo, rem] : ndarrays
529 Of Hermite series coefficients representing the quotient and
530 remainder.
531
532 See Also
533 --------
534 hermadd, hermsub, hermmulx, hermmul, hermpow
535
536 Notes
537 -----
538 In general, the (polynomial) division of one Hermite series by another
539 results in quotient and remainder terms that are not in the Hermite
540 polynomial basis set. Thus, to express these results as a Hermite
541 series, it is necessary to "reproject" the results onto the Hermite
542 basis set, which may produce "unintuitive" (but correct) results; see
543 Examples section below.
544
545 Examples
546 --------
547 >>> from numpy.polynomial.hermite import hermdiv
548 >>> hermdiv([ 52., 29., 52., 7., 6.], [0, 1, 2])
549 (array([1., 2., 3.]), array([0.]))
550 >>> hermdiv([ 54., 31., 52., 7., 6.], [0, 1, 2])
551 (array([1., 2., 3.]), array([2., 2.]))
552 >>> hermdiv([ 53., 30., 52., 7., 6.], [0, 1, 2])
553 (array([1., 2., 3.]), array([1., 1.]))
554
555 """
556 return pu._div(hermmul, c1, c2)
557
558
559def hermpow(c, pow, maxpower=16):
560 """Raise a Hermite series to a power.
561
562 Returns the Hermite series `c` raised to the power `pow`. The
563 argument `c` is a sequence of coefficients ordered from low to high.
564 i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.``
565
566 Parameters
567 ----------
568 c : array_like
569 1-D array of Hermite series coefficients ordered from low to
570 high.
571 pow : integer
572 Power to which the series will be raised
573 maxpower : integer, optional
574 Maximum power allowed. This is mainly to limit growth of the series
575 to unmanageable size. Default is 16
576
577 Returns
578 -------
579 coef : ndarray
580 Hermite series of power.
581
582 See Also
583 --------
584 hermadd, hermsub, hermmulx, hermmul, hermdiv
585
586 Examples
587 --------
588 >>> from numpy.polynomial.hermite import hermpow
589 >>> hermpow([1, 2, 3], 2)
590 array([81., 52., 82., 12., 9.])
591
592 """
593 return pu._pow(hermmul, c, pow, maxpower)
594
595
596def hermder(c, m=1, scl=1, axis=0):
597 """
598 Differentiate a Hermite series.
599
600 Returns the Hermite series coefficients `c` differentiated `m` times
601 along `axis`. At each iteration the result is multiplied by `scl` (the
602 scaling factor is for use in a linear change of variable). The argument
603 `c` is an array of coefficients from low to high degree along each
604 axis, e.g., [1,2,3] represents the series ``1*H_0 + 2*H_1 + 3*H_2``
605 while [[1,2],[1,2]] represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) +
606 2*H_0(x)*H_1(y) + 2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is
607 ``y``.
608
609 Parameters
610 ----------
611 c : array_like
612 Array of Hermite series coefficients. If `c` is multidimensional the
613 different axis correspond to different variables with the degree in
614 each axis given by the corresponding index.
615 m : int, optional
616 Number of derivatives taken, must be non-negative. (Default: 1)
617 scl : scalar, optional
618 Each differentiation is multiplied by `scl`. The end result is
619 multiplication by ``scl**m``. This is for use in a linear change of
620 variable. (Default: 1)
621 axis : int, optional
622 Axis over which the derivative is taken. (Default: 0).
623
624 Returns
625 -------
626 der : ndarray
627 Hermite series of the derivative.
628
629 See Also
630 --------
631 hermint
632
633 Notes
634 -----
635 In general, the result of differentiating a Hermite series does not
636 resemble the same operation on a power series. Thus the result of this
637 function may be "unintuitive," albeit correct; see Examples section
638 below.
639
640 Examples
641 --------
642 >>> from numpy.polynomial.hermite import hermder
643 >>> hermder([ 1. , 0.5, 0.5, 0.5])
644 array([1., 2., 3.])
645 >>> hermder([-0.5, 1./2., 1./8., 1./12., 1./16.], m=2)
646 array([1., 2., 3.])
647
648 """
649 c = np.array(c, ndmin=1, copy=True)
650 if c.dtype.char in '?bBhHiIlLqQpP':
651 c = c.astype(np.double)
652 cnt = pu._as_int(m, "the order of derivation")
653 iaxis = pu._as_int(axis, "the axis")
654 if cnt < 0:
655 raise ValueError("The order of derivation must be non-negative")
656 iaxis = np.lib.array_utils.normalize_axis_index(iaxis, c.ndim)
657
658 if cnt == 0:
659 return c
660
661 c = np.moveaxis(c, iaxis, 0)
662 n = len(c)
663 if cnt >= n:
664 c = c[:1] * 0
665 else:
666 for i in range(cnt):
667 n = n - 1
668 c *= scl
669 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
670 for j in range(n, 0, -1):
671 der[j - 1] = (2 * j) * c[j]
672 c = der
673 c = np.moveaxis(c, 0, iaxis)
674 return c
675
676
677def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
678 """
679 Integrate a Hermite series.
680
681 Returns the Hermite series coefficients `c` integrated `m` times from
682 `lbnd` along `axis`. At each iteration the resulting series is
683 **multiplied** by `scl` and an integration constant, `k`, is added.
684 The scaling factor is for use in a linear change of variable. ("Buyer
685 beware": note that, depending on what one is doing, one may want `scl`
686 to be the reciprocal of what one might expect; for more information,
687 see the Notes section below.) The argument `c` is an array of
688 coefficients from low to high degree along each axis, e.g., [1,2,3]
689 represents the series ``H_0 + 2*H_1 + 3*H_2`` while [[1,2],[1,2]]
690 represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) + 2*H_0(x)*H_1(y) +
691 2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
692
693 Parameters
694 ----------
695 c : array_like
696 Array of Hermite series coefficients. If c is multidimensional the
697 different axis correspond to different variables with the degree in
698 each axis given by the corresponding index.
699 m : int, optional
700 Order of integration, must be positive. (Default: 1)
701 k : {[], list, scalar}, optional
702 Integration constant(s). The value of the first integral at
703 ``lbnd`` is the first value in the list, the value of the second
704 integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
705 default), all constants are set to zero. If ``m == 1``, a single
706 scalar can be given instead of a list.
707 lbnd : scalar, optional
708 The lower bound of the integral. (Default: 0)
709 scl : scalar, optional
710 Following each integration the result is *multiplied* by `scl`
711 before the integration constant is added. (Default: 1)
712 axis : int, optional
713 Axis over which the integral is taken. (Default: 0).
714
715 Returns
716 -------
717 S : ndarray
718 Hermite series coefficients of the integral.
719
720 Raises
721 ------
722 ValueError
723 If ``m < 0``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
724 ``np.ndim(scl) != 0``.
725
726 See Also
727 --------
728 hermder
729
730 Notes
731 -----
732 Note that the result of each integration is *multiplied* by `scl`.
733 Why is this important to note? Say one is making a linear change of
734 variable :math:`u = ax + b` in an integral relative to `x`. Then
735 :math:`dx = du/a`, so one will need to set `scl` equal to
736 :math:`1/a` - perhaps not what one would have first thought.
737
738 Also note that, in general, the result of integrating a C-series needs
739 to be "reprojected" onto the C-series basis set. Thus, typically,
740 the result of this function is "unintuitive," albeit correct; see
741 Examples section below.
742
743 Examples
744 --------
745 >>> from numpy.polynomial.hermite import hermint
746 >>> hermint([1,2,3]) # integrate once, value 0 at 0.
747 array([1. , 0.5, 0.5, 0.5])
748 >>> hermint([1,2,3], m=2) # integrate twice, value & deriv 0 at 0
749 array([-0.5 , 0.5 , 0.125 , 0.08333333, 0.0625 ]) # may vary
750 >>> hermint([1,2,3], k=1) # integrate once, value 1 at 0.
751 array([2. , 0.5, 0.5, 0.5])
752 >>> hermint([1,2,3], lbnd=-1) # integrate once, value 0 at -1
753 array([-2. , 0.5, 0.5, 0.5])
754 >>> hermint([1,2,3], m=2, k=[1,2], lbnd=-1)
755 array([ 1.66666667, -0.5 , 0.125 , 0.08333333, 0.0625 ]) # may vary
756
757 """
758 c = np.array(c, ndmin=1, copy=True)
759 if c.dtype.char in '?bBhHiIlLqQpP':
760 c = c.astype(np.double)
761 if not np.iterable(k):
762 k = [k]
763 cnt = pu._as_int(m, "the order of integration")
764 iaxis = pu._as_int(axis, "the axis")
765 if cnt < 0:
766 raise ValueError("The order of integration must be non-negative")
767 if len(k) > cnt:
768 raise ValueError("Too many integration constants")
769 if np.ndim(lbnd) != 0:
770 raise ValueError("lbnd must be a scalar.")
771 if np.ndim(scl) != 0:
772 raise ValueError("scl must be a scalar.")
773 iaxis = np.lib.array_utils.normalize_axis_index(iaxis, c.ndim)
774
775 if cnt == 0:
776 return c
777
778 c = np.moveaxis(c, iaxis, 0)
779 k = list(k) + [0] * (cnt - len(k))
780 for i in range(cnt):
781 n = len(c)
782 c *= scl
783 if n == 1 and np.all(c[0] == 0):
784 c[0] += k[i]
785 else:
786 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
787 tmp[0] = c[0] * 0
788 tmp[1] = c[0] / 2
789 for j in range(1, n):
790 tmp[j + 1] = c[j] / (2 * (j + 1))
791 tmp[0] += k[i] - hermval(lbnd, tmp)
792 c = tmp
793 c = np.moveaxis(c, 0, iaxis)
794 return c
795
796
797def hermval(x, c, tensor=True):
798 """
799 Evaluate a Hermite series at points x.
800
801 If `c` is of length ``n + 1``, this function returns the value:
802
803 .. math:: p(x) = c_0 * H_0(x) + c_1 * H_1(x) + ... + c_n * H_n(x)
804
805 The parameter `x` is converted to an array only if it is a tuple or a
806 list, otherwise it is treated as a scalar. In either case, either `x`
807 or its elements must support multiplication and addition both with
808 themselves and with the elements of `c`.
809
810 If `c` is a 1-D array, then ``p(x)`` will have the same shape as `x`. If
811 `c` is multidimensional, then the shape of the result depends on the
812 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
813 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
814 scalars have shape (,).
815
816 Trailing zeros in the coefficients will be used in the evaluation, so
817 they should be avoided if efficiency is a concern.
818
819 Parameters
820 ----------
821 x : array_like, compatible object
822 If `x` is a list or tuple, it is converted to an ndarray, otherwise
823 it is left unchanged and treated as a scalar. In either case, `x`
824 or its elements must support addition and multiplication with
825 themselves and with the elements of `c`.
826 c : array_like
827 Array of coefficients ordered so that the coefficients for terms of
828 degree n are contained in c[n]. If `c` is multidimensional the
829 remaining indices enumerate multiple polynomials. In the two
830 dimensional case the coefficients may be thought of as stored in
831 the columns of `c`.
832 tensor : boolean, optional
833 If True, the shape of the coefficient array is extended with ones
834 on the right, one for each dimension of `x`. Scalars have dimension 0
835 for this action. The result is that every column of coefficients in
836 `c` is evaluated for every element of `x`. If False, `x` is broadcast
837 over the columns of `c` for the evaluation. This keyword is useful
838 when `c` is multidimensional. The default value is True.
839
840 Returns
841 -------
842 values : ndarray, algebra_like
843 The shape of the return value is described above.
844
845 See Also
846 --------
847 hermval2d, hermgrid2d, hermval3d, hermgrid3d
848
849 Notes
850 -----
851 The evaluation uses Clenshaw recursion, aka synthetic division.
852
853 Examples
854 --------
855 >>> from numpy.polynomial.hermite import hermval
856 >>> coef = [1,2,3]
857 >>> hermval(1, coef)
858 11.0
859 >>> hermval([[1,2],[3,4]], coef)
860 array([[ 11., 51.],
861 [115., 203.]])
862
863 """
864 c = np.array(c, ndmin=1, copy=None)
865 if c.dtype.char in '?bBhHiIlLqQpP':
866 c = c.astype(np.double)
867 if isinstance(x, (tuple, list)):
868 x = np.asarray(x)
869 if isinstance(x, np.ndarray) and tensor:
870 c = c.reshape(c.shape + (1,) * x.ndim)
871
872 x2 = x * 2
873 if len(c) == 1:
874 c0 = c[0]
875 c1 = 0
876 elif len(c) == 2:
877 c0 = c[0]
878 c1 = c[1]
879 else:
880 nd = len(c)
881 c0 = c[-2]
882 c1 = c[-1]
883 for i in range(3, len(c) + 1):
884 tmp = c0
885 nd = nd - 1
886 c0 = c[-i] - c1 * (2 * (nd - 1))
887 c1 = tmp + c1 * x2
888 return c0 + c1 * x2
889
890
891def hermval2d(x, y, c):
892 """
893 Evaluate a 2-D Hermite series at points (x, y).
894
895 This function returns the values:
896
897 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * H_i(x) * H_j(y)
898
899 The parameters `x` and `y` are converted to arrays only if they are
900 tuples or a lists, otherwise they are treated as a scalars and they
901 must have the same shape after conversion. In either case, either `x`
902 and `y` or their elements must support multiplication and addition both
903 with themselves and with the elements of `c`.
904
905 If `c` is a 1-D array a one is implicitly appended to its shape to make
906 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
907
908 Parameters
909 ----------
910 x, y : array_like, compatible objects
911 The two dimensional series is evaluated at the points ``(x, y)``,
912 where `x` and `y` must have the same shape. If `x` or `y` is a list
913 or tuple, it is first converted to an ndarray, otherwise it is left
914 unchanged and if it isn't an ndarray it is treated as a scalar.
915 c : array_like
916 Array of coefficients ordered so that the coefficient of the term
917 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
918 dimension greater than two the remaining indices enumerate multiple
919 sets of coefficients.
920
921 Returns
922 -------
923 values : ndarray, compatible object
924 The values of the two dimensional polynomial at points formed with
925 pairs of corresponding values from `x` and `y`.
926
927 See Also
928 --------
929 hermval, hermgrid2d, hermval3d, hermgrid3d
930
931 Examples
932 --------
933 >>> from numpy.polynomial.hermite import hermval2d
934 >>> x = [1, 2]
935 >>> y = [4, 5]
936 >>> c = [[1, 2, 3], [4, 5, 6]]
937 >>> hermval2d(x, y, c)
938 array([1035., 2883.])
939
940 """
941 return pu._valnd(hermval, c, x, y)
942
943
944def hermgrid2d(x, y, c):
945 """
946 Evaluate a 2-D Hermite series on the Cartesian product of x and y.
947
948 This function returns the values:
949
950 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * H_i(a) * H_j(b)
951
952 where the points ``(a, b)`` consist of all pairs formed by taking
953 `a` from `x` and `b` from `y`. The resulting points form a grid with
954 `x` in the first dimension and `y` in the second.
955
956 The parameters `x` and `y` are converted to arrays only if they are
957 tuples or a lists, otherwise they are treated as a scalars. In either
958 case, either `x` and `y` or their elements must support multiplication
959 and addition both with themselves and with the elements of `c`.
960
961 If `c` has fewer than two dimensions, ones are implicitly appended to
962 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
963 x.shape.
964
965 Parameters
966 ----------
967 x, y : array_like, compatible objects
968 The two dimensional series is evaluated at the points in the
969 Cartesian product of `x` and `y`. If `x` or `y` is a list or
970 tuple, it is first converted to an ndarray, otherwise it is left
971 unchanged and, if it isn't an ndarray, it is treated as a scalar.
972 c : array_like
973 Array of coefficients ordered so that the coefficients for terms of
974 degree i,j are contained in ``c[i,j]``. If `c` has dimension
975 greater than two the remaining indices enumerate multiple sets of
976 coefficients.
977
978 Returns
979 -------
980 values : ndarray, compatible object
981 The values of the two dimensional polynomial at points in the Cartesian
982 product of `x` and `y`.
983
984 See Also
985 --------
986 hermval, hermval2d, hermval3d, hermgrid3d
987
988 Examples
989 --------
990 >>> from numpy.polynomial.hermite import hermgrid2d
991 >>> x = [1, 2, 3]
992 >>> y = [4, 5]
993 >>> c = [[1, 2, 3], [4, 5, 6]]
994 >>> hermgrid2d(x, y, c)
995 array([[1035., 1599.],
996 [1867., 2883.],
997 [2699., 4167.]])
998
999 """
1000 return pu._gridnd(hermval, c, x, y)
1001
1002
1003def hermval3d(x, y, z, c):
1004 """
1005 Evaluate a 3-D Hermite series at points (x, y, z).
1006
1007 This function returns the values:
1008
1009 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * H_i(x) * H_j(y) * H_k(z)
1010
1011 The parameters `x`, `y`, and `z` are converted to arrays only if
1012 they are tuples or a lists, otherwise they are treated as a scalars and
1013 they must have the same shape after conversion. In either case, either
1014 `x`, `y`, and `z` or their elements must support multiplication and
1015 addition both with themselves and with the elements of `c`.
1016
1017 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
1018 shape to make it 3-D. The shape of the result will be c.shape[3:] +
1019 x.shape.
1020
1021 Parameters
1022 ----------
1023 x, y, z : array_like, compatible object
1024 The three dimensional series is evaluated at the points
1025 ``(x, y, z)``, where `x`, `y`, and `z` must have the same shape. If
1026 any of `x`, `y`, or `z` is a list or tuple, it is first converted
1027 to an ndarray, otherwise it is left unchanged and if it isn't an
1028 ndarray it is treated as a scalar.
1029 c : array_like
1030 Array of coefficients ordered so that the coefficient of the term of
1031 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
1032 greater than 3 the remaining indices enumerate multiple sets of
1033 coefficients.
1034
1035 Returns
1036 -------
1037 values : ndarray, compatible object
1038 The values of the multidimensional polynomial on points formed with
1039 triples of corresponding values from `x`, `y`, and `z`.
1040
1041 See Also
1042 --------
1043 hermval, hermval2d, hermgrid2d, hermgrid3d
1044
1045 Examples
1046 --------
1047 >>> from numpy.polynomial.hermite import hermval3d
1048 >>> x = [1, 2]
1049 >>> y = [4, 5]
1050 >>> z = [6, 7]
1051 >>> c = [[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]]
1052 >>> hermval3d(x, y, z, c)
1053 array([ 40077., 120131.])
1054
1055 """
1056 return pu._valnd(hermval, c, x, y, z)
1057
1058
1059def hermgrid3d(x, y, z, c):
1060 """
1061 Evaluate a 3-D Hermite series on the Cartesian product of x, y, and z.
1062
1063 This function returns the values:
1064
1065 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * H_i(a) * H_j(b) * H_k(c)
1066
1067 where the points ``(a, b, c)`` consist of all triples formed by taking
1068 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1069 a grid with `x` in the first dimension, `y` in the second, and `z` in
1070 the third.
1071
1072 The parameters `x`, `y`, and `z` are converted to arrays only if they
1073 are tuples or a lists, otherwise they are treated as a scalars. In
1074 either case, either `x`, `y`, and `z` or their elements must support
1075 multiplication and addition both with themselves and with the elements
1076 of `c`.
1077
1078 If `c` has fewer than three dimensions, ones are implicitly appended to
1079 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1080 x.shape + y.shape + z.shape.
1081
1082 Parameters
1083 ----------
1084 x, y, z : array_like, compatible objects
1085 The three dimensional series is evaluated at the points in the
1086 Cartesian product of `x`, `y`, and `z`. If `x`, `y`, or `z` is a
1087 list or tuple, it is first converted to an ndarray, otherwise it is
1088 left unchanged and, if it isn't an ndarray, it is treated as a
1089 scalar.
1090 c : array_like
1091 Array of coefficients ordered so that the coefficients for terms of
1092 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1093 greater than two the remaining indices enumerate multiple sets of
1094 coefficients.
1095
1096 Returns
1097 -------
1098 values : ndarray, compatible object
1099 The values of the two dimensional polynomial at points in the Cartesian
1100 product of `x` and `y`.
1101
1102 See Also
1103 --------
1104 hermval, hermval2d, hermgrid2d, hermval3d
1105
1106 Examples
1107 --------
1108 >>> from numpy.polynomial.hermite import hermgrid3d
1109 >>> x = [1, 2]
1110 >>> y = [4, 5]
1111 >>> z = [6, 7]
1112 >>> c = [[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]]
1113 >>> hermgrid3d(x, y, z, c)
1114 array([[[ 40077., 54117.],
1115 [ 49293., 66561.]],
1116 [[ 72375., 97719.],
1117 [ 88975., 120131.]]])
1118
1119 """
1120 return pu._gridnd(hermval, c, x, y, z)
1121
1122
1123def hermvander(x, deg):
1124 """Pseudo-Vandermonde matrix of given degree.
1125
1126 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
1127 `x`. The pseudo-Vandermonde matrix is defined by
1128
1129 .. math:: V[..., i] = H_i(x),
1130
1131 where ``0 <= i <= deg``. The leading indices of `V` index the elements of
1132 `x` and the last index is the degree of the Hermite polynomial.
1133
1134 If `c` is a 1-D array of coefficients of length ``n + 1`` and `V` is the
1135 array ``V = hermvander(x, n)``, then ``np.dot(V, c)`` and
1136 ``hermval(x, c)`` are the same up to roundoff. This equivalence is
1137 useful both for least squares fitting and for the evaluation of a large
1138 number of Hermite series of the same degree and sample points.
1139
1140 Parameters
1141 ----------
1142 x : array_like
1143 Array of points. The dtype is converted to float64 or complex128
1144 depending on whether any of the elements are complex. If `x` is
1145 scalar it is converted to a 1-D array.
1146 deg : int
1147 Degree of the resulting matrix.
1148
1149 Returns
1150 -------
1151 vander : ndarray
1152 The pseudo-Vandermonde matrix. The shape of the returned matrix is
1153 ``x.shape + (deg + 1,)``, where The last index is the degree of the
1154 corresponding Hermite polynomial. The dtype will be the same as
1155 the converted `x`.
1156
1157 Examples
1158 --------
1159 >>> import numpy as np
1160 >>> from numpy.polynomial.hermite import hermvander
1161 >>> x = np.array([-1, 0, 1])
1162 >>> hermvander(x, 3)
1163 array([[ 1., -2., 2., 4.],
1164 [ 1., 0., -2., -0.],
1165 [ 1., 2., 2., -4.]])
1166
1167 """
1168 ideg = pu._as_int(deg, "deg")
1169 if ideg < 0:
1170 raise ValueError("deg must be non-negative")
1171
1172 x = np.array(x, copy=None, ndmin=1) + 0.0
1173 dims = (ideg + 1,) + x.shape
1174 dtyp = x.dtype
1175 v = np.empty(dims, dtype=dtyp)
1176 v[0] = x * 0 + 1
1177 if ideg > 0:
1178 x2 = x * 2
1179 v[1] = x2
1180 for i in range(2, ideg + 1):
1181 v[i] = (v[i - 1] * x2 - v[i - 2] * (2 * (i - 1)))
1182 return np.moveaxis(v, 0, -1)
1183
1184
1185def hermvander2d(x, y, deg):
1186 """Pseudo-Vandermonde matrix of given degrees.
1187
1188 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1189 points ``(x, y)``. The pseudo-Vandermonde matrix is defined by
1190
1191 .. math:: V[..., (deg[1] + 1)*i + j] = H_i(x) * H_j(y),
1192
1193 where ``0 <= i <= deg[0]`` and ``0 <= j <= deg[1]``. The leading indices of
1194 `V` index the points ``(x, y)`` and the last index encodes the degrees of
1195 the Hermite polynomials.
1196
1197 If ``V = hermvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1198 correspond to the elements of a 2-D coefficient array `c` of shape
1199 (xdeg + 1, ydeg + 1) in the order
1200