1"""
2===================================================================
3HermiteE Series, "Probabilists" (:mod:`numpy.polynomial.hermite_e`)
4===================================================================
5
6This module provides a number of objects (mostly functions) useful for
7dealing with Hermite_e series, including a `HermiteE` 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 HermiteE
18
19Constants
20---------
21.. autosummary::
22 :toctree: generated/
23
24 hermedomain
25 hermezero
26 hermeone
27 hermex
28
29Arithmetic
30----------
31.. autosummary::
32 :toctree: generated/
33
34 hermeadd
35 hermesub
36 hermemulx
37 hermemul
38 hermediv
39 hermepow
40 hermeval
41 hermeval2d
42 hermeval3d
43 hermegrid2d
44 hermegrid3d
45
46Calculus
47--------
48.. autosummary::
49 :toctree: generated/
50
51 hermeder
52 hermeint
53
54Misc Functions
55--------------
56.. autosummary::
57 :toctree: generated/
58
59 hermefromroots
60 hermeroots
61 hermevander
62 hermevander2d
63 hermevander3d
64 hermegauss
65 hermeweight
66 hermecompanion
67 hermefit
68 hermetrim
69 hermeline
70 herme2poly
71 poly2herme
72
73See also
74--------
75`numpy.polynomial`
76
77"""
78import numpy as np
79import numpy.linalg as la
80from numpy.core.multiarray import normalize_axis_index
81
82from . import polyutils as pu
83from ._polybase import ABCPolyBase
84
85__all__ = [
86 'hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline',
87 'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv',
88 'hermepow', 'hermeval', 'hermeder', 'hermeint', 'herme2poly',
89 'poly2herme', 'hermefromroots', 'hermevander', 'hermefit', 'hermetrim',
90 'hermeroots', 'HermiteE', 'hermeval2d', 'hermeval3d', 'hermegrid2d',
91 'hermegrid3d', 'hermevander2d', 'hermevander3d', 'hermecompanion',
92 'hermegauss', 'hermeweight']
93
94hermetrim = pu.trimcoef
95
96
97def poly2herme(pol):
98 """
99 poly2herme(pol)
100
101 Convert a polynomial to a Hermite series.
102
103 Convert an array representing the coefficients of a polynomial (relative
104 to the "standard" basis) ordered from lowest degree to highest, to an
105 array of the coefficients of the equivalent Hermite series, ordered
106 from lowest to highest degree.
107
108 Parameters
109 ----------
110 pol : array_like
111 1-D array containing the polynomial coefficients
112
113 Returns
114 -------
115 c : ndarray
116 1-D array containing the coefficients of the equivalent Hermite
117 series.
118
119 See Also
120 --------
121 herme2poly
122
123 Notes
124 -----
125 The easy way to do conversions between polynomial basis sets
126 is to use the convert method of a class instance.
127
128 Examples
129 --------
130 >>> from numpy.polynomial.hermite_e import poly2herme
131 >>> poly2herme(np.arange(4))
132 array([ 2., 10., 2., 3.])
133
134 """
135 [pol] = pu.as_series([pol])
136 deg = len(pol) - 1
137 res = 0
138 for i in range(deg, -1, -1):
139 res = hermeadd(hermemulx(res), pol[i])
140 return res
141
142
143def herme2poly(c):
144 """
145 Convert a Hermite series to a polynomial.
146
147 Convert an array representing the coefficients of a Hermite series,
148 ordered from lowest degree to highest, to an array of the coefficients
149 of the equivalent polynomial (relative to the "standard" basis) ordered
150 from lowest to highest degree.
151
152 Parameters
153 ----------
154 c : array_like
155 1-D array containing the Hermite series coefficients, ordered
156 from lowest order term to highest.
157
158 Returns
159 -------
160 pol : ndarray
161 1-D array containing the coefficients of the equivalent polynomial
162 (relative to the "standard" basis) ordered from lowest order term
163 to highest.
164
165 See Also
166 --------
167 poly2herme
168
169 Notes
170 -----
171 The easy way to do conversions between polynomial basis sets
172 is to use the convert method of a class instance.
173
174 Examples
175 --------
176 >>> from numpy.polynomial.hermite_e import herme2poly
177 >>> herme2poly([ 2., 10., 2., 3.])
178 array([0., 1., 2., 3.])
179
180 """
181 from .polynomial import polyadd, polysub, polymulx
182
183 [c] = pu.as_series([c])
184 n = len(c)
185 if n == 1:
186 return c
187 if n == 2:
188 return c
189 else:
190 c0 = c[-2]
191 c1 = c[-1]
192 # i is the current degree of c1
193 for i in range(n - 1, 1, -1):
194 tmp = c0
195 c0 = polysub(c[i - 2], c1*(i - 1))
196 c1 = polyadd(tmp, polymulx(c1))
197 return polyadd(c0, polymulx(c1))
198
199#
200# These are constant arrays are of integer type so as to be compatible
201# with the widest range of other types, such as Decimal.
202#
203
204# Hermite
205hermedomain = np.array([-1, 1])
206
207# Hermite coefficients representing zero.
208hermezero = np.array([0])
209
210# Hermite coefficients representing one.
211hermeone = np.array([1])
212
213# Hermite coefficients representing the identity x.
214hermex = np.array([0, 1])
215
216
217def hermeline(off, scl):
218 """
219 Hermite series whose graph is a straight line.
220
221 Parameters
222 ----------
223 off, scl : scalars
224 The specified line is given by ``off + scl*x``.
225
226 Returns
227 -------
228 y : ndarray
229 This module's representation of the Hermite series for
230 ``off + scl*x``.
231
232 See Also
233 --------
234 numpy.polynomial.polynomial.polyline
235 numpy.polynomial.chebyshev.chebline
236 numpy.polynomial.legendre.legline
237 numpy.polynomial.laguerre.lagline
238 numpy.polynomial.hermite.hermline
239
240 Examples
241 --------
242 >>> from numpy.polynomial.hermite_e import hermeline
243 >>> from numpy.polynomial.hermite_e import hermeline, hermeval
244 >>> hermeval(0,hermeline(3, 2))
245 3.0
246 >>> hermeval(1,hermeline(3, 2))
247 5.0
248
249 """
250 if scl != 0:
251 return np.array([off, scl])
252 else:
253 return np.array([off])
254
255
256def hermefromroots(roots):
257 """
258 Generate a HermiteE 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 HermiteE form, where the `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 * He_1(x) + ... + c_n * He_n(x)
273
274 The coefficient of the last term is not generally 1 for monic
275 polynomials in HermiteE 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.hermite.hermfromroots
296 numpy.polynomial.chebyshev.chebfromroots
297
298 Examples
299 --------
300 >>> from numpy.polynomial.hermite_e import hermefromroots, hermeval
301 >>> coef = hermefromroots((-1, 0, 1))
302 >>> hermeval((-1, 0, 1), coef)
303 array([0., 0., 0.])
304 >>> coef = hermefromroots((-1j, 1j))
305 >>> hermeval((-1j, 1j), coef)
306 array([0.+0.j, 0.+0.j])
307
308 """
309 return pu._fromroots(hermeline, hermemul, roots)
310
311
312def hermeadd(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 hermesub, hermemulx, hermemul, hermediv, hermepow
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_e import hermeadd
345 >>> hermeadd([1, 2, 3], [1, 2, 3, 4])
346 array([2., 4., 6., 4.])
347
348 """
349 return pu._add(c1, c2)
350
351
352def hermesub(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 hermeadd, hermemulx, hermemul, hermediv, hermepow
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_e import hermesub
385 >>> hermesub([1, 2, 3, 4], [1, 2, 3])
386 array([0., 0., 0., 4.])
387
388 """
389 return pu._sub(c1, c2)
390
391
392def hermemulx(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 Notes
411 -----
412 The multiplication uses the recursion relationship for Hermite
413 polynomials in the form
414
415 .. math::
416
417 xP_i(x) = (P_{i + 1}(x) + iP_{i - 1}(x)))
418
419 Examples
420 --------
421 >>> from numpy.polynomial.hermite_e import hermemulx
422 >>> hermemulx([1, 2, 3])
423 array([2., 7., 2., 3.])
424
425 """
426 # c is a trimmed copy
427 [c] = pu.as_series([c])
428 # The zero series needs special treatment
429 if len(c) == 1 and c[0] == 0:
430 return c
431
432 prd = np.empty(len(c) + 1, dtype=c.dtype)
433 prd[0] = c[0]*0
434 prd[1] = c[0]
435 for i in range(1, len(c)):
436 prd[i + 1] = c[i]
437 prd[i - 1] += c[i]*i
438 return prd
439
440
441def hermemul(c1, c2):
442 """
443 Multiply one Hermite series by another.
444
445 Returns the product of two Hermite series `c1` * `c2`. The arguments
446 are sequences of coefficients, from lowest order "term" to highest,
447 e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
448
449 Parameters
450 ----------
451 c1, c2 : array_like
452 1-D arrays of Hermite series coefficients ordered from low to
453 high.
454
455 Returns
456 -------
457 out : ndarray
458 Of Hermite series coefficients representing their product.
459
460 See Also
461 --------
462 hermeadd, hermesub, hermemulx, hermediv, hermepow
463
464 Notes
465 -----
466 In general, the (polynomial) product of two C-series results in terms
467 that are not in the Hermite polynomial basis set. Thus, to express
468 the product as a Hermite series, it is necessary to "reproject" the
469 product onto said basis set, which may produce "unintuitive" (but
470 correct) results; see Examples section below.
471
472 Examples
473 --------
474 >>> from numpy.polynomial.hermite_e import hermemul
475 >>> hermemul([1, 2, 3], [0, 1, 2])
476 array([14., 15., 28., 7., 6.])
477
478 """
479 # s1, s2 are trimmed copies
480 [c1, c2] = pu.as_series([c1, c2])
481
482 if len(c1) > len(c2):
483 c = c2
484 xs = c1
485 else:
486 c = c1
487 xs = c2
488
489 if len(c) == 1:
490 c0 = c[0]*xs
491 c1 = 0
492 elif len(c) == 2:
493 c0 = c[0]*xs
494 c1 = c[1]*xs
495 else:
496 nd = len(c)
497 c0 = c[-2]*xs
498 c1 = c[-1]*xs
499 for i in range(3, len(c) + 1):
500 tmp = c0
501 nd = nd - 1
502 c0 = hermesub(c[-i]*xs, c1*(nd - 1))
503 c1 = hermeadd(tmp, hermemulx(c1))
504 return hermeadd(c0, hermemulx(c1))
505
506
507def hermediv(c1, c2):
508 """
509 Divide one Hermite series by another.
510
511 Returns the quotient-with-remainder of two Hermite series
512 `c1` / `c2`. The arguments are sequences of coefficients from lowest
513 order "term" to highest, e.g., [1,2,3] represents the series
514 ``P_0 + 2*P_1 + 3*P_2``.
515
516 Parameters
517 ----------
518 c1, c2 : array_like
519 1-D arrays of Hermite series coefficients ordered from low to
520 high.
521
522 Returns
523 -------
524 [quo, rem] : ndarrays
525 Of Hermite series coefficients representing the quotient and
526 remainder.
527
528 See Also
529 --------
530 hermeadd, hermesub, hermemulx, hermemul, hermepow
531
532 Notes
533 -----
534 In general, the (polynomial) division of one Hermite series by another
535 results in quotient and remainder terms that are not in the Hermite
536 polynomial basis set. Thus, to express these results as a Hermite
537 series, it is necessary to "reproject" the results onto the Hermite
538 basis set, which may produce "unintuitive" (but correct) results; see
539 Examples section below.
540
541 Examples
542 --------
543 >>> from numpy.polynomial.hermite_e import hermediv
544 >>> hermediv([ 14., 15., 28., 7., 6.], [0, 1, 2])
545 (array([1., 2., 3.]), array([0.]))
546 >>> hermediv([ 15., 17., 28., 7., 6.], [0, 1, 2])
547 (array([1., 2., 3.]), array([1., 2.]))
548
549 """
550 return pu._div(hermemul, c1, c2)
551
552
553def hermepow(c, pow, maxpower=16):
554 """Raise a Hermite series to a power.
555
556 Returns the Hermite series `c` raised to the power `pow`. The
557 argument `c` is a sequence of coefficients ordered from low to high.
558 i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.``
559
560 Parameters
561 ----------
562 c : array_like
563 1-D array of Hermite series coefficients ordered from low to
564 high.
565 pow : integer
566 Power to which the series will be raised
567 maxpower : integer, optional
568 Maximum power allowed. This is mainly to limit growth of the series
569 to unmanageable size. Default is 16
570
571 Returns
572 -------
573 coef : ndarray
574 Hermite series of power.
575
576 See Also
577 --------
578 hermeadd, hermesub, hermemulx, hermemul, hermediv
579
580 Examples
581 --------
582 >>> from numpy.polynomial.hermite_e import hermepow
583 >>> hermepow([1, 2, 3], 2)
584 array([23., 28., 46., 12., 9.])
585
586 """
587 return pu._pow(hermemul, c, pow, maxpower)
588
589
590def hermeder(c, m=1, scl=1, axis=0):
591 """
592 Differentiate a Hermite_e series.
593
594 Returns the series coefficients `c` differentiated `m` times along
595 `axis`. At each iteration the result is multiplied by `scl` (the
596 scaling factor is for use in a linear change of variable). The argument
597 `c` is an array of coefficients from low to high degree along each
598 axis, e.g., [1,2,3] represents the series ``1*He_0 + 2*He_1 + 3*He_2``
599 while [[1,2],[1,2]] represents ``1*He_0(x)*He_0(y) + 1*He_1(x)*He_0(y)
600 + 2*He_0(x)*He_1(y) + 2*He_1(x)*He_1(y)`` if axis=0 is ``x`` and axis=1
601 is ``y``.
602
603 Parameters
604 ----------
605 c : array_like
606 Array of Hermite_e series coefficients. If `c` is multidimensional
607 the different axis correspond to different variables with the
608 degree in each axis given by the corresponding index.
609 m : int, optional
610 Number of derivatives taken, must be non-negative. (Default: 1)
611 scl : scalar, optional
612 Each differentiation is multiplied by `scl`. The end result is
613 multiplication by ``scl**m``. This is for use in a linear change of
614 variable. (Default: 1)
615 axis : int, optional
616 Axis over which the derivative is taken. (Default: 0).
617
618 .. versionadded:: 1.7.0
619
620 Returns
621 -------
622 der : ndarray
623 Hermite series of the derivative.
624
625 See Also
626 --------
627 hermeint
628
629 Notes
630 -----
631 In general, the result of differentiating a Hermite series does not
632 resemble the same operation on a power series. Thus the result of this
633 function may be "unintuitive," albeit correct; see Examples section
634 below.
635
636 Examples
637 --------
638 >>> from numpy.polynomial.hermite_e import hermeder
639 >>> hermeder([ 1., 1., 1., 1.])
640 array([1., 2., 3.])
641 >>> hermeder([-0.25, 1., 1./2., 1./3., 1./4 ], m=2)
642 array([1., 2., 3.])
643
644 """
645 c = np.array(c, ndmin=1, copy=True)
646 if c.dtype.char in '?bBhHiIlLqQpP':
647 c = c.astype(np.double)
648 cnt = pu._deprecate_as_int(m, "the order of derivation")
649 iaxis = pu._deprecate_as_int(axis, "the axis")
650 if cnt < 0:
651 raise ValueError("The order of derivation must be non-negative")
652 iaxis = normalize_axis_index(iaxis, c.ndim)
653
654 if cnt == 0:
655 return c
656
657 c = np.moveaxis(c, iaxis, 0)
658 n = len(c)
659 if cnt >= n:
660 return c[:1]*0
661 else:
662 for i in range(cnt):
663 n = n - 1
664 c *= scl
665 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
666 for j in range(n, 0, -1):
667 der[j - 1] = j*c[j]
668 c = der
669 c = np.moveaxis(c, 0, iaxis)
670 return c
671
672
673def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
674 """
675 Integrate a Hermite_e series.
676
677 Returns the Hermite_e series coefficients `c` integrated `m` times from
678 `lbnd` along `axis`. At each iteration the resulting series is
679 **multiplied** by `scl` and an integration constant, `k`, is added.
680 The scaling factor is for use in a linear change of variable. ("Buyer
681 beware": note that, depending on what one is doing, one may want `scl`
682 to be the reciprocal of what one might expect; for more information,
683 see the Notes section below.) The argument `c` is an array of
684 coefficients from low to high degree along each axis, e.g., [1,2,3]
685 represents the series ``H_0 + 2*H_1 + 3*H_2`` while [[1,2],[1,2]]
686 represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) + 2*H_0(x)*H_1(y) +
687 2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
688
689 Parameters
690 ----------
691 c : array_like
692 Array of Hermite_e series coefficients. If c is multidimensional
693 the different axis correspond to different variables with the
694 degree in each axis given by the corresponding index.
695 m : int, optional
696 Order of integration, must be positive. (Default: 1)
697 k : {[], list, scalar}, optional
698 Integration constant(s). The value of the first integral at
699 ``lbnd`` is the first value in the list, the value of the second
700 integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
701 default), all constants are set to zero. If ``m == 1``, a single
702 scalar can be given instead of a list.
703 lbnd : scalar, optional
704 The lower bound of the integral. (Default: 0)
705 scl : scalar, optional
706 Following each integration the result is *multiplied* by `scl`
707 before the integration constant is added. (Default: 1)
708 axis : int, optional
709 Axis over which the integral is taken. (Default: 0).
710
711 .. versionadded:: 1.7.0
712
713 Returns
714 -------
715 S : ndarray
716 Hermite_e series coefficients of the integral.
717
718 Raises
719 ------
720 ValueError
721 If ``m < 0``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
722 ``np.ndim(scl) != 0``.
723
724 See Also
725 --------
726 hermeder
727
728 Notes
729 -----
730 Note that the result of each integration is *multiplied* by `scl`.
731 Why is this important to note? Say one is making a linear change of
732 variable :math:`u = ax + b` in an integral relative to `x`. Then
733 :math:`dx = du/a`, so one will need to set `scl` equal to
734 :math:`1/a` - perhaps not what one would have first thought.
735
736 Also note that, in general, the result of integrating a C-series needs
737 to be "reprojected" onto the C-series basis set. Thus, typically,
738 the result of this function is "unintuitive," albeit correct; see
739 Examples section below.
740
741 Examples
742 --------
743 >>> from numpy.polynomial.hermite_e import hermeint
744 >>> hermeint([1, 2, 3]) # integrate once, value 0 at 0.
745 array([1., 1., 1., 1.])
746 >>> hermeint([1, 2, 3], m=2) # integrate twice, value & deriv 0 at 0
747 array([-0.25 , 1. , 0.5 , 0.33333333, 0.25 ]) # may vary
748 >>> hermeint([1, 2, 3], k=1) # integrate once, value 1 at 0.
749 array([2., 1., 1., 1.])
750 >>> hermeint([1, 2, 3], lbnd=-1) # integrate once, value 0 at -1
751 array([-1., 1., 1., 1.])
752 >>> hermeint([1, 2, 3], m=2, k=[1, 2], lbnd=-1)
753 array([ 1.83333333, 0. , 0.5 , 0.33333333, 0.25 ]) # may vary
754
755 """
756 c = np.array(c, ndmin=1, copy=True)
757 if c.dtype.char in '?bBhHiIlLqQpP':
758 c = c.astype(np.double)
759 if not np.iterable(k):
760 k = [k]
761 cnt = pu._deprecate_as_int(m, "the order of integration")
762 iaxis = pu._deprecate_as_int(axis, "the axis")
763 if cnt < 0:
764 raise ValueError("The order of integration must be non-negative")
765 if len(k) > cnt:
766 raise ValueError("Too many integration constants")
767 if np.ndim(lbnd) != 0:
768 raise ValueError("lbnd must be a scalar.")
769 if np.ndim(scl) != 0:
770 raise ValueError("scl must be a scalar.")
771 iaxis = normalize_axis_index(iaxis, c.ndim)
772
773 if cnt == 0:
774 return c
775
776 c = np.moveaxis(c, iaxis, 0)
777 k = list(k) + [0]*(cnt - len(k))
778 for i in range(cnt):
779 n = len(c)
780 c *= scl
781 if n == 1 and np.all(c[0] == 0):
782 c[0] += k[i]
783 else:
784 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
785 tmp[0] = c[0]*0
786 tmp[1] = c[0]
787 for j in range(1, n):
788 tmp[j + 1] = c[j]/(j + 1)
789 tmp[0] += k[i] - hermeval(lbnd, tmp)
790 c = tmp
791 c = np.moveaxis(c, 0, iaxis)
792 return c
793
794
795def hermeval(x, c, tensor=True):
796 """
797 Evaluate an HermiteE series at points x.
798
799 If `c` is of length `n + 1`, this function returns the value:
800
801 .. math:: p(x) = c_0 * He_0(x) + c_1 * He_1(x) + ... + c_n * He_n(x)
802
803 The parameter `x` is converted to an array only if it is a tuple or a
804 list, otherwise it is treated as a scalar. In either case, either `x`
805 or its elements must support multiplication and addition both with
806 themselves and with the elements of `c`.
807
808 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
809 `c` is multidimensional, then the shape of the result depends on the
810 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
811 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
812 scalars have shape (,).
813
814 Trailing zeros in the coefficients will be used in the evaluation, so
815 they should be avoided if efficiency is a concern.
816
817 Parameters
818 ----------
819 x : array_like, compatible object
820 If `x` is a list or tuple, it is converted to an ndarray, otherwise
821 it is left unchanged and treated as a scalar. In either case, `x`
822 or its elements must support addition and multiplication with
823 with themselves and with the elements of `c`.
824 c : array_like
825 Array of coefficients ordered so that the coefficients for terms of
826 degree n are contained in c[n]. If `c` is multidimensional the
827 remaining indices enumerate multiple polynomials. In the two
828 dimensional case the coefficients may be thought of as stored in
829 the columns of `c`.
830 tensor : boolean, optional
831 If True, the shape of the coefficient array is extended with ones
832 on the right, one for each dimension of `x`. Scalars have dimension 0
833 for this action. The result is that every column of coefficients in
834 `c` is evaluated for every element of `x`. If False, `x` is broadcast
835 over the columns of `c` for the evaluation. This keyword is useful
836 when `c` is multidimensional. The default value is True.
837
838 .. versionadded:: 1.7.0
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 hermeval2d, hermegrid2d, hermeval3d, hermegrid3d
848
849 Notes
850 -----
851 The evaluation uses Clenshaw recursion, aka synthetic division.
852
853 Examples
854 --------
855 >>> from numpy.polynomial.hermite_e import hermeval
856 >>> coef = [1,2,3]
857 >>> hermeval(1, coef)
858 3.0
859 >>> hermeval([[1,2],[3,4]], coef)
860 array([[ 3., 14.],
861 [31., 54.]])
862
863 """
864 c = np.array(c, ndmin=1, copy=False)
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 if len(c) == 1:
873 c0 = c[0]
874 c1 = 0
875 elif len(c) == 2:
876 c0 = c[0]
877 c1 = c[1]
878 else:
879 nd = len(c)
880 c0 = c[-2]
881 c1 = c[-1]
882 for i in range(3, len(c) + 1):
883 tmp = c0
884 nd = nd - 1
885 c0 = c[-i] - c1*(nd - 1)
886 c1 = tmp + c1*x
887 return c0 + c1*x
888
889
890def hermeval2d(x, y, c):
891 """
892 Evaluate a 2-D HermiteE series at points (x, y).
893
894 This function returns the values:
895
896 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * He_i(x) * He_j(y)
897
898 The parameters `x` and `y` are converted to arrays only if they are
899 tuples or a lists, otherwise they are treated as a scalars and they
900 must have the same shape after conversion. In either case, either `x`
901 and `y` or their elements must support multiplication and addition both
902 with themselves and with the elements of `c`.
903
904 If `c` is a 1-D array a one is implicitly appended to its shape to make
905 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
906
907 Parameters
908 ----------
909 x, y : array_like, compatible objects
910 The two dimensional series is evaluated at the points `(x, y)`,
911 where `x` and `y` must have the same shape. If `x` or `y` is a list
912 or tuple, it is first converted to an ndarray, otherwise it is left
913 unchanged and if it isn't an ndarray it is treated as a scalar.
914 c : array_like
915 Array of coefficients ordered so that the coefficient of the term
916 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
917 dimension greater than two the remaining indices enumerate multiple
918 sets of coefficients.
919
920 Returns
921 -------
922 values : ndarray, compatible object
923 The values of the two dimensional polynomial at points formed with
924 pairs of corresponding values from `x` and `y`.
925
926 See Also
927 --------
928 hermeval, hermegrid2d, hermeval3d, hermegrid3d
929
930 Notes
931 -----
932
933 .. versionadded:: 1.7.0
934
935 """
936 return pu._valnd(hermeval, c, x, y)
937
938
939def hermegrid2d(x, y, c):
940 """
941 Evaluate a 2-D HermiteE series on the Cartesian product of x and y.
942
943 This function returns the values:
944
945 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * H_i(a) * H_j(b)
946
947 where the points `(a, b)` consist of all pairs formed by taking
948 `a` from `x` and `b` from `y`. The resulting points form a grid with
949 `x` in the first dimension and `y` in the second.
950
951 The parameters `x` and `y` are converted to arrays only if they are
952 tuples or a lists, otherwise they are treated as a scalars. In either
953 case, either `x` and `y` or their elements must support multiplication
954 and addition both with themselves and with the elements of `c`.
955
956 If `c` has fewer than two dimensions, ones are implicitly appended to
957 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
958 x.shape.
959
960 Parameters
961 ----------
962 x, y : array_like, compatible objects
963 The two dimensional series is evaluated at the points in the
964 Cartesian product of `x` and `y`. If `x` or `y` is a list or
965 tuple, it is first converted to an ndarray, otherwise it is left
966 unchanged and, if it isn't an ndarray, it is treated as a scalar.
967 c : array_like
968 Array of coefficients ordered so that the coefficients for terms of
969 degree i,j are contained in ``c[i,j]``. If `c` has dimension
970 greater than two the remaining indices enumerate multiple sets of
971 coefficients.
972
973 Returns
974 -------
975 values : ndarray, compatible object
976 The values of the two dimensional polynomial at points in the Cartesian
977 product of `x` and `y`.
978
979 See Also
980 --------
981 hermeval, hermeval2d, hermeval3d, hermegrid3d
982
983 Notes
984 -----
985
986 .. versionadded:: 1.7.0
987
988 """
989 return pu._gridnd(hermeval, c, x, y)
990
991
992def hermeval3d(x, y, z, c):
993 """
994 Evaluate a 3-D Hermite_e series at points (x, y, z).
995
996 This function returns the values:
997
998 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * He_i(x) * He_j(y) * He_k(z)
999
1000 The parameters `x`, `y`, and `z` are converted to arrays only if
1001 they are tuples or a lists, otherwise they are treated as a scalars and
1002 they must have the same shape after conversion. In either case, either
1003 `x`, `y`, and `z` or their elements must support multiplication and
1004 addition both with themselves and with the elements of `c`.
1005
1006 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
1007 shape to make it 3-D. The shape of the result will be c.shape[3:] +
1008 x.shape.
1009
1010 Parameters
1011 ----------
1012 x, y, z : array_like, compatible object
1013 The three dimensional series is evaluated at the points
1014 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
1015 any of `x`, `y`, or `z` is a list or tuple, it is first converted
1016 to an ndarray, otherwise it is left unchanged and if it isn't an
1017 ndarray it is treated as a scalar.
1018 c : array_like
1019 Array of coefficients ordered so that the coefficient of the term of
1020 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
1021 greater than 3 the remaining indices enumerate multiple sets of
1022 coefficients.
1023
1024 Returns
1025 -------
1026 values : ndarray, compatible object
1027 The values of the multidimensional polynomial on points formed with
1028 triples of corresponding values from `x`, `y`, and `z`.
1029
1030 See Also
1031 --------
1032 hermeval, hermeval2d, hermegrid2d, hermegrid3d
1033
1034 Notes
1035 -----
1036
1037 .. versionadded:: 1.7.0
1038
1039 """
1040 return pu._valnd(hermeval, c, x, y, z)
1041
1042
1043def hermegrid3d(x, y, z, c):
1044 """
1045 Evaluate a 3-D HermiteE series on the Cartesian product of x, y, and z.
1046
1047 This function returns the values:
1048
1049 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * He_i(a) * He_j(b) * He_k(c)
1050
1051 where the points `(a, b, c)` consist of all triples formed by taking
1052 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1053 a grid with `x` in the first dimension, `y` in the second, and `z` in
1054 the third.
1055
1056 The parameters `x`, `y`, and `z` are converted to arrays only if they
1057 are tuples or a lists, otherwise they are treated as a scalars. In
1058 either case, either `x`, `y`, and `z` or their elements must support
1059 multiplication and addition both with themselves and with the elements
1060 of `c`.
1061
1062 If `c` has fewer than three dimensions, ones are implicitly appended to
1063 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1064 x.shape + y.shape + z.shape.
1065
1066 Parameters
1067 ----------
1068 x, y, z : array_like, compatible objects
1069 The three dimensional series is evaluated at the points in the
1070 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
1071 list or tuple, it is first converted to an ndarray, otherwise it is
1072 left unchanged and, if it isn't an ndarray, it is treated as a
1073 scalar.
1074 c : array_like
1075 Array of coefficients ordered so that the coefficients for terms of
1076 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1077 greater than two the remaining indices enumerate multiple sets of
1078 coefficients.
1079
1080 Returns
1081 -------
1082 values : ndarray, compatible object
1083 The values of the two dimensional polynomial at points in the Cartesian
1084 product of `x` and `y`.
1085
1086 See Also
1087 --------
1088 hermeval, hermeval2d, hermegrid2d, hermeval3d
1089
1090 Notes
1091 -----
1092
1093 .. versionadded:: 1.7.0
1094
1095 """
1096 return pu._gridnd(hermeval, c, x, y, z)
1097
1098
1099def hermevander(x, deg):
1100 """Pseudo-Vandermonde matrix of given degree.
1101
1102 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
1103 `x`. The pseudo-Vandermonde matrix is defined by
1104
1105 .. math:: V[..., i] = He_i(x),
1106
1107 where `0 <= i <= deg`. The leading indices of `V` index the elements of
1108 `x` and the last index is the degree of the HermiteE polynomial.
1109
1110 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
1111 array ``V = hermevander(x, n)``, then ``np.dot(V, c)`` and
1112 ``hermeval(x, c)`` are the same up to roundoff. This equivalence is
1113 useful both for least squares fitting and for the evaluation of a large
1114 number of HermiteE series of the same degree and sample points.
1115
1116 Parameters
1117 ----------
1118 x : array_like
1119 Array of points. The dtype is converted to float64 or complex128
1120 depending on whether any of the elements are complex. If `x` is
1121 scalar it is converted to a 1-D array.
1122 deg : int
1123 Degree of the resulting matrix.
1124
1125 Returns
1126 -------
1127 vander : ndarray
1128 The pseudo-Vandermonde matrix. The shape of the returned matrix is
1129 ``x.shape + (deg + 1,)``, where The last index is the degree of the
1130 corresponding HermiteE polynomial. The dtype will be the same as
1131 the converted `x`.
1132
1133 Examples
1134 --------
1135 >>> from numpy.polynomial.hermite_e import hermevander
1136 >>> x = np.array([-1, 0, 1])
1137 >>> hermevander(x, 3)
1138 array([[ 1., -1., 0., 2.],
1139 [ 1., 0., -1., -0.],
1140 [ 1., 1., 0., -2.]])
1141
1142 """
1143 ideg = pu._deprecate_as_int(deg, "deg")
1144 if ideg < 0:
1145 raise ValueError("deg must be non-negative")
1146
1147 x = np.array(x, copy=False, ndmin=1) + 0.0
1148 dims = (ideg + 1,) + x.shape
1149 dtyp = x.dtype
1150 v = np.empty(dims, dtype=dtyp)
1151 v[0] = x*0 + 1
1152 if ideg > 0:
1153 v[1] = x
1154 for i in range(2, ideg + 1):
1155 v[i] = (v[i-1]*x - v[i-2]*(i - 1))
1156 return np.moveaxis(v, 0, -1)
1157
1158
1159def hermevander2d(x, y, deg):
1160 """Pseudo-Vandermonde matrix of given degrees.
1161
1162 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1163 points `(x, y)`. The pseudo-Vandermonde matrix is defined by
1164
1165 .. math:: V[..., (deg[1] + 1)*i + j] = He_i(x) * He_j(y),
1166
1167 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
1168 `V` index the points `(x, y)` and the last index encodes the degrees of
1169 the HermiteE polynomials.
1170
1171 If ``V = hermevander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1172 correspond to the elements of a 2-D coefficient array `c` of shape
1173 (xdeg + 1, ydeg + 1) in the order
1174
1175 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
1176
1177 and ``np.dot(V, c.flat)`` and ``hermeval2d(x, y, c)`` will be the same
1178 up to roundoff. This equivalence is useful both for least squares
1179 fitting and for the evaluation of a large number of 2-D HermiteE
1180 series of the same degrees and sample points.
1181
1182 Parameters
1183 ----------
1184 x, y : array_like
1185 Arrays of point coordinates, all of the same shape. The dtypes
1186 will be converted to either float64 or complex128 depending on
1187 whether any of the elements are complex. Scalars are converted to
1188 1-D arrays.
1189 deg : list of ints
1190 List of maximum degrees of the form [x_deg, y_deg].
1191
1192 Returns
1193 -------
1194 vander2d : ndarray
1195 The shape of the returned matrix is ``x.shape + (order,)``, where
1196 :math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same
1197 as the converted `x` and `y`.
1198
1199 See Also
1200 --------
1201 hermevander, hermevander3d, hermeval2d, hermeval3d
1202
1203 Notes
1204 -----
1205
1206 .. versionadded:: 1.7.0
1207
1208 """
1209 return pu._vander_nd_flat((hermevander, hermevander), (x, y), deg)
1210
1211
1212def hermevander3d(x, y, z, deg):
1213 """Pseudo-Vandermonde matrix of given degrees.
1214
1215 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1216 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
1217 then Hehe pseudo-Vandermonde matrix is defined by
1218
1219 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = He_i(x)*He_j(y)*He_k(z),
1220
1221 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
1222 indices of `V` index the points `(x, y, z)` and the last index encodes
1223 the degrees of the HermiteE polynomials.
1224
1225 If ``V = hermevander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
1226 of `V` correspond to the elements of a 3-D coefficient array `c` of
1227 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
1228
1229 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
1230
1231 and ``np.dot(V, c.flat)`` and ``hermeval3d(x, y, z, c)`` will be the
1232 same up to roundoff. This equivalence is useful both for least squares
1233 fitting and for the evaluation of a large number of 3-D HermiteE
1234 series of the same degrees and sample points.
1235
1236 Parameters
1237 ----------
1238 x, y, z : array_like
1239 Arrays of point coordinates, all of the same shape. The dtypes will
1240 be converted to either float64 or complex128 depending on whether
1241 any of the elements are complex. Scalars are converted to 1-D
1242 arrays.
1243 deg : list of ints
1244 List of maximum degrees of the form [x_deg, y_deg, z_deg].
1245
1246 Returns
1247 -------
1248 vander3d : ndarray
1249 The shape of the returned matrix is ``x.shape + (order,)``, where
1250 :math:`order = (deg[0]+1)*(deg[1]+1)*(deg[2]+1)`. The dtype will
1251 be the same as the converted `x`, `y`, and `z`.
1252
1253 See Also
1254 --------
1255 hermevander, hermevander3d, hermeval2d, hermeval3d
1256
1257 Notes
1258 -----
1259
1260 .. versionadded:: 1.7.0
1261
1262 """
1263 return pu._vander_nd_flat((hermevander, hermevander, hermevander), (x, y, z), deg)
1264
1265
1266def hermefit(x, y, deg, rcond=None, full=False, w=None):
1267 """
1268 Least squares fit of Hermite series to data.
1269
1270 Return the coefficients of a HermiteE series of degree `deg` that is
1271 the least squares fit to the data values `y` given at points `x`. If
1272 `y` is 1-D the returned coefficients will also be 1-D. If `y` is 2-D
1273 multiple fits are done, one for each column of `y`, and the resulting
1274 coefficients are stored in the corresponding columns of a 2-D return.
1275 The fitted polynomial(s) are in the form
1276
1277 .. math:: p(x) = c_0 + c_1 * He_1(x) + ... + c_n * He_n(x),
1278
1279 where `n` is `deg`.
1280
1281 Parameters
1282 ----------
1283 x : array_like, shape (M,)
1284 x-coordinates of the M sample points ``(x[i], y[i])``.
1285 y : array_like, shape (M,) or (M, K)
1286 y-coordinates of the sample points. Several data sets of sample
1287 points sharing the same x-coordinates can be fitted at once by
1288 passing in a 2D-array that contains one dataset per column.
1289 deg : int or 1-D array_like
1290 Degree(s) of the fitting polynomials. If `deg` is a single integer
1291 all terms up to and including the `deg`'th term are included in the
1292 fit. For NumPy versions >= 1.11.0 a list of integers specifying the
1293 degrees of the terms to include may be used instead.
1294 rcond : float, optional
1295 Relative condition number of the fit. Singular values smaller than
1296 this relative to the largest singular value will be ignored. The
1297 default value is len(x)*eps, where eps is the relative precision of
1298 the float type, about 2e-16 in most cases.
1299 full : bool, optional
1300 Switch determining nature of return value. When it is False (the
1301 default) just the coefficients are returned, when True diagnostic
1302 information from the singular value decomposition is also returned.
1303 w : array_like, shape (`M`,), optional
1304 Weights. If not None, the weight ``w[i]`` applies to the unsquared
1305 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
1306 chosen so that the errors of the products ``w[i]*y[i]`` all have the
1307 same variance. When using inverse-variance weighting, use
1308 ``w[i] = 1/sigma(y[i])``. The default value is None.
1309
1310 Returns
1311 -------
1312 coef : ndarray, shape (M,) or (M, K)
1313 Hermite coefficients ordered from low to high. If `y` was 2-D,
1314 the coefficients for the data in column k of `y` are in column
1315 `k`.
1316
1317 [residuals, rank, singular_values, rcond] : list
1318 These values are only returned if ``full == True``
1319
1320 - residuals -- sum of squared residuals of the least squares fit
1321 - rank -- the numerical rank of the scaled Vandermonde matrix
1322 - singular_values -- singular values of the scaled Vandermonde matrix
1323 - rcond -- value of `rcond`.
1324
1325 For more details, see `numpy.linalg.lstsq`.
1326
1327 Warns
1328 -----
1329 RankWarning
1330 The rank of the coefficient matrix in the least-squares fit is
1331 deficient. The warning is only raised if ``full = False``. The
1332 warnings can be turned off by
1333
1334 >>> import warnings
1335 >>> warnings.simplefilter('ignore', np.RankWarning)
1336
1337 See Also
1338 --------
1339 numpy.polynomial.chebyshev.chebfit
1340 numpy.polynomial.legendre.legfit
1341 numpy.polynomial.polynomial.polyfit
1342 numpy.polynomial.hermite.hermfit
1343 numpy.polynomial.laguerre.lagfit
1344 hermeval : Evaluates a Hermite series.
1345 hermevander : pseudo Vandermonde matrix of Hermite series.
1346 hermeweight : HermiteE weight function.
1347 numpy.linalg.lstsq : Computes a least-squares fit from the matrix.
1348 scipy.interpolate.UnivariateSpline : Computes spline fits.
1349
1350 Notes
1351 -----
1352 The solution is the coefficients of the HermiteE series `p` that
1353 minimizes the sum of the weighted squared errors
1354
1355 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
1356
1357 where the :math:`w_j` are the weights. This problem is solved by
1358 setting up the (typically) overdetermined matrix equation
1359
1360 .. math:: V(x) * c = w * y,
1361
1362 where `V` is the pseudo Vandermonde matrix of `x`, the elements of `c`
1363 are the coefficients to be solved for, and the elements of `y` are the
1364 observed values. This equation is then solved using the singular value
1365 decomposition of `V`.
1366
1367 If some of the singular values of `V` are so small that they are
1368 neglected, then a `RankWarning` will be issued. This means that the
1369 coefficient values may be poorly determined. Using a lower order fit
1370 will usually get rid of the warning. The `rcond` parameter can also be
1371 set to a value smaller than its default, but the resulting fit may be
1372 spurious and have large contributions from roundoff error.
1373
1374 Fits using HermiteE series are probably most useful when the data can
1375 be approximated by ``sqrt(w(x)) * p(x)``, where `w(x)` is the HermiteE
1376 weight. In that case the weight ``sqrt(w(x[i]))`` should be used
1377 together with data values ``y[i]/sqrt(w(x[i]))``. The weight function is
1378 available as `hermeweight`.
1379
1380 References
1381 ----------
1382 .. [1] Wikipedia, "Curve fitting",
1383 https://en.wikipedia.org/wiki/Curve_fitting
1384
1385 Examples
1386 --------
1387 >>> from numpy.polynomial.hermite_e import hermefit, hermeval
1388 >>> x = np.linspace(-10, 10)
1389 >>> np.random.seed(123)
1390 >>> err = np.random.randn(len(x))/10
1391 >>> y = hermeval(x, [1, 2, 3]) + err
1392 >>> hermefit(x, y, 2)
1393 array([ 1.01690445, 1.99951418, 2.99948696]) # may vary
1394
1395 """
1396 return pu._fit(hermevander, x, y, deg, rcond, full, w)
1397
1398
1399def hermecompanion(c):
1400 """
1401 Return the scaled companion matrix of c.
1402
1403 The basis polynomials are scaled so that the companion matrix is
1404 symmetric when `c` is an HermiteE basis polynomial. This provides
1405 better eigenvalue estimates than the unscaled case and for basis
1406 polynomials the eigenvalues are guaranteed to be real if
1407 `numpy.linalg.eigvalsh` is used to obtain them.
1408
1409 Parameters
1410 ----------
1411 c : array_like
1412 1-D array of HermiteE series coefficients ordered from low to high
1413 degree.
1414
1415 Returns
1416 -------
1417 mat : ndarray
1418 Scaled companion matrix of dimensions (deg, deg).
1419
1420 Notes
1421 -----
1422
1423 .. versionadded:: 1.7.0
1424
1425 """
1426 # c is a trimmed copy
1427 [c] = pu.as_series([c])
1428 if len(c) < 2:
1429 raise ValueError('Series must have maximum degree of at least 1.')
1430 if len(c) == 2:
1431 return np.array([[-c[0]/c[1]]])
1432
1433 n = len(c) - 1
1434 mat = np.zeros((n, n), dtype=c.dtype)
1435 scl = np.hstack((1., 1./np.sqrt(np.arange(n - 1, 0, -1))))
1436 scl = np.multiply.accumulate(scl)[::-1]
1437 top = mat.reshape(-1)[1::n+1]
1438 bot = mat.reshape(-1)[n::n+1]
1439 top[...] = np.sqrt(np.arange(1, n))
1440 bot[...] = top
1441 mat[:, -1] -= scl*c[:-1]/c[-1]
1442 return mat
1443
1444
1445def hermeroots(c):
1446 """
1447 Compute the roots of a HermiteE series.
1448
1449 Return the roots (a.k.a. "zeros") of the polynomial
1450
1451 .. math:: p(x) = \\sum_i c[i] * He_i(x).
1452
1453 Parameters
1454 ----------
1455 c : 1-D array_like
1456 1-D array of coefficients.
1457
1458 Returns
1459 -------
1460 out : ndarray
1461 Array of the roots of the series. If all the roots are real,
1462 then `out` is also real, otherwise it is complex.
1463
1464 See Also
1465 --------
1466 numpy.polynomial.polynomial.polyroots
1467 numpy.polynomial.legendre.legroots
1468 numpy.polynomial.laguerre.lagroots
1469 numpy.polynomial.hermite.hermroots
1470 numpy.polynomial.chebyshev.chebroots
1471
1472 Notes
1473 -----
1474 The root estimates are obtained as the eigenvalues of the companion
1475 matrix, Roots far from the origin of the complex plane may have large
1476 errors due to the numerical instability of the series for such
1477 values. Roots with multiplicity greater than 1 will also show larger
1478 errors as the value of the series near such points is relatively
1479 insensitive to errors in the roots. Isolated roots near the origin can
1480 be improved by a few iterations of Newton's method.
1481
1482 The HermiteE series basis polynomials aren't powers of `x` so the
1483 results of this function may seem unintuitive.
1484
1485 Examples
1486 --------
1487 >>> from numpy.polynomial.hermite_e import hermeroots, hermefromroots
1488 >>> coef = hermefromroots([-1, 0, 1])
1489 >>> coef
1490 array([0., 2., 0., 1.])
1491 >>> hermeroots(coef)
1492 array([-1., 0., 1.]) # may vary
1493
1494 """
1495 # c is a trimmed copy
1496 [c] = pu.as_series([c])
1497 if len(c) <= 1:
1498 return np.array([], dtype=c.dtype)
1499 if len(c) == 2:
1500 return np.array([-c[0]/c[1]])
1501
1502 # rotated companion matrix reduces error
1503 m = hermecompanion(c)[::-1,::-1]
1504 r = la.eigvals(m)
1505 r.sort()
1506 return r
1507
1508
1509def _normed_hermite_e_n(x, n):
1510 """
1511 Evaluate a normalized HermiteE polynomial.
1512
1513 Compute the value of the normalized HermiteE polynomial of degree ``n``
1514 at the points ``x``.
1515
1516
1517 Parameters
1518 ----------
1519 x : ndarray of double.
1520 Points at which to evaluate the function
1521 n : int
1522 Degree of the normalized HermiteE function to be evaluated.
1523
1524 Returns
1525 -------
1526 values : ndarray
1527 The shape of the return value is described above.
1528
1529 Notes
1530 -----
1531 .. versionadded:: 1.10.0
1532
1533 This function is needed for finding the Gauss points and integration
1534 weights for high degrees. The values of the standard HermiteE functions
1535 overflow when n >= 207.
1536
1537 """
1538 if n == 0:
1539 return np.full(x.shape, 1/np.sqrt(np.sqrt(2*np.pi)))
1540
1541 c0 = 0.
1542 c1 = 1./np.sqrt(np.sqrt(2*np.pi))
1543 nd = float(n)
1544 for i in range(n - 1):
1545 tmp = c0
1546 c0 = -c1*np.sqrt((nd - 1.)/nd)
1547 c1 = tmp + c1*x*np.sqrt(1./nd)
1548 nd = nd - 1.0
1549 return c0 + c1*x
1550
1551
1552def hermegauss(deg):
1553 """
1554 Gauss-HermiteE quadrature.
1555
1556 Computes the sample points and weights for Gauss-HermiteE quadrature.
1557 These sample points and weights will correctly integrate polynomials of
1558 degree :math:`2*deg - 1` or less over the interval :math:`[-\\inf, \\inf]`
1559 with the weight function :math:`f(x) = \\exp(-x^2/2)`.
1560
1561 Parameters
1562 ----------
1563 deg : int
1564 Number of sample points and weights. It must be >= 1.
1565
1566 Returns
1567 -------
1568 x : ndarray
1569 1-D ndarray containing the sample points.
1570 y : ndarray
1571 1-D ndarray containing the weights.
1572
1573 Notes
1574 -----
1575
1576 .. versionadded:: 1.7.0
1577
1578 The results have only been tested up to degree 100, higher degrees may
1579 be problematic. The weights are determined by using the fact that
1580
1581 .. math:: w_k = c / (He'_n(x_k) * He_{n-1}(x_k))
1582
1583 where :math:`c` is a constant independent of :math:`k` and :math:`x_k`
1584 is the k'th root of :math:`He_n`, and then scaling the results to get
1585 the right value when integrating 1.
1586
1587 """
1588 ideg = pu._deprecate_as_int(deg, "deg")
1589 if ideg <= 0:
1590 raise ValueError("deg must be a positive integer")
1591
1592 # first approximation of roots. We use the fact that the companion
1593 # matrix is symmetric in this case in order to obtain better zeros.
1594 c = np.array([0]*deg + [1])
1595 m = hermecompanion(c)
1596 x = la.eigvalsh(m)
1597
1598 # improve roots by one application of Newton
1599 dy = _normed_hermite_e_n(x, ideg)
1600 df = _normed_hermite_e_n(x, ideg - 1) * np.sqrt(ideg)
1601 x -= dy/df
1602
1603 # compute the weights. We scale the factor to avoid possible numerical
1604 # overflow.
1605 fm = _normed_hermite_e_n(x, ideg - 1)
1606 fm /= np.abs(fm).max()
1607 w = 1/(fm * fm)
1608
1609 # for Hermite_e we can also symmetrize
1610 w = (w + w[::-1])/2
1611 x = (x - x[::-1])/2
1612
1613 # scale w to get the right value
1614 w *= np.sqrt(2*np.pi) / w.sum()
1615
1616 return x, w
1617
1618
1619def hermeweight(x):
1620 """Weight function of the Hermite_e polynomials.
1621
1622 The weight function is :math:`\\exp(-x^2/2)` and the interval of
1623 integration is :math:`[-\\inf, \\inf]`. the HermiteE polynomials are
1624 orthogonal, but not normalized, with respect to this weight function.
1625
1626 Parameters
1627 ----------
1628 x : array_like
1629 Values at which the weight function will be computed.
1630
1631 Returns
1632 -------
1633 w : ndarray
1634 The weight function at `x`.
1635
1636 Notes
1637 -----
1638
1639 .. versionadded:: 1.7.0
1640
1641 """
1642 w = np.exp(-.5*x**2)
1643 return w
1644
1645
1646#
1647# HermiteE series class
1648#
1649
1650class HermiteE(ABCPolyBase):
1651 """An HermiteE series class.
1652
1653 The HermiteE class provides the standard Python numerical methods
1654 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
1655 attributes and methods listed in the `ABCPolyBase` documentation.
1656
1657 Parameters
1658 ----------
1659 coef : array_like
1660 HermiteE coefficients in order of increasing degree, i.e,
1661 ``(1, 2, 3)`` gives ``1*He_0(x) + 2*He_1(X) + 3*He_2(x)``.
1662 domain : (2,) array_like, optional
1663 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
1664 to the interval ``[window[0], window[1]]`` by shifting and scaling.
1665 The default value is [-1, 1].
1666 window : (2,) array_like, optional
1667 Window, see `domain` for its use. The default value is [-1, 1].
1668
1669 .. versionadded:: 1.6.0
1670
1671 """
1672 # Virtual Functions
1673 _add = staticmethod(hermeadd)
1674 _sub = staticmethod(hermesub)
1675 _mul = staticmethod(hermemul)
1676 _div = staticmethod(hermediv)
1677 _pow = staticmethod(hermepow)
1678 _val = staticmethod(hermeval)
1679 _int = staticmethod(hermeint)
1680 _der = staticmethod(hermeder)
1681 _fit = staticmethod(hermefit)
1682 _line = staticmethod(hermeline)
1683 _roots = staticmethod(hermeroots)
1684 _fromroots = staticmethod(hermefromroots)
1685
1686 # Virtual properties
1687 domain = np.array(hermedomain)
1688 window = np.array(hermedomain)
1689 basis_name = 'He'