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