laguerre.py 51 KB

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