laguerre.py 50 KB

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