hermite.py 50 KB

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