legendre.py 50 KB

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