laguerre.py 49 KB

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