laguerre.py 50 KB

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