legendre.py 50 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664
  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.legendre.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, 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. """
  484. return pu._pow(legmul, c, pow, maxpower)
  485. def legder(c, m=1, scl=1, axis=0):
  486. """
  487. Differentiate a Legendre series.
  488. Returns the Legendre series coefficients `c` differentiated `m` times
  489. along `axis`. At each iteration the result is multiplied by `scl` (the
  490. scaling factor is for use in a linear change of variable). The argument
  491. `c` is an array of coefficients from low to high degree along each
  492. axis, e.g., [1,2,3] represents the series ``1*L_0 + 2*L_1 + 3*L_2``
  493. while [[1,2],[1,2]] represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) +
  494. 2*L_0(x)*L_1(y) + 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is
  495. ``y``.
  496. Parameters
  497. ----------
  498. c : array_like
  499. Array of Legendre series coefficients. If c is multidimensional the
  500. different axis correspond to different variables with the degree in
  501. each axis given by the corresponding index.
  502. m : int, optional
  503. Number of derivatives taken, must be non-negative. (Default: 1)
  504. scl : scalar, optional
  505. Each differentiation is multiplied by `scl`. The end result is
  506. multiplication by ``scl**m``. This is for use in a linear change of
  507. variable. (Default: 1)
  508. axis : int, optional
  509. Axis over which the derivative is taken. (Default: 0).
  510. .. versionadded:: 1.7.0
  511. Returns
  512. -------
  513. der : ndarray
  514. Legendre series of the derivative.
  515. See Also
  516. --------
  517. legint
  518. Notes
  519. -----
  520. In general, the result of differentiating a Legendre series does not
  521. resemble the same operation on a power series. Thus the result of this
  522. function may be "unintuitive," albeit correct; see Examples section
  523. below.
  524. Examples
  525. --------
  526. >>> from numpy.polynomial import legendre as L
  527. >>> c = (1,2,3,4)
  528. >>> L.legder(c)
  529. array([ 6., 9., 20.])
  530. >>> L.legder(c, 3)
  531. array([60.])
  532. >>> L.legder(c, scl=-1)
  533. array([ -6., -9., -20.])
  534. >>> L.legder(c, 2,-1)
  535. array([ 9., 60.])
  536. """
  537. c = np.array(c, ndmin=1, copy=True)
  538. if c.dtype.char in '?bBhHiIlLqQpP':
  539. c = c.astype(np.double)
  540. cnt = pu._deprecate_as_int(m, "the order of derivation")
  541. iaxis = pu._deprecate_as_int(axis, "the axis")
  542. if cnt < 0:
  543. raise ValueError("The order of derivation must be non-negative")
  544. iaxis = normalize_axis_index(iaxis, c.ndim)
  545. if cnt == 0:
  546. return c
  547. c = np.moveaxis(c, iaxis, 0)
  548. n = len(c)
  549. if cnt >= n:
  550. c = c[:1]*0
  551. else:
  552. for i in range(cnt):
  553. n = n - 1
  554. c *= scl
  555. der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
  556. for j in range(n, 2, -1):
  557. der[j - 1] = (2*j - 1)*c[j]
  558. c[j - 2] += c[j]
  559. if n > 1:
  560. der[1] = 3*c[2]
  561. der[0] = c[1]
  562. c = der
  563. c = np.moveaxis(c, 0, iaxis)
  564. return c
  565. def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
  566. """
  567. Integrate a Legendre series.
  568. Returns the Legendre series coefficients `c` integrated `m` times from
  569. `lbnd` along `axis`. At each iteration the resulting series is
  570. **multiplied** by `scl` and an integration constant, `k`, is added.
  571. The scaling factor is for use in a linear change of variable. ("Buyer
  572. beware": note that, depending on what one is doing, one may want `scl`
  573. to be the reciprocal of what one might expect; for more information,
  574. see the Notes section below.) The argument `c` is an array of
  575. coefficients from low to high degree along each axis, e.g., [1,2,3]
  576. represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]]
  577. represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) +
  578. 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
  579. Parameters
  580. ----------
  581. c : array_like
  582. Array of Legendre series coefficients. If c is multidimensional the
  583. different axis correspond to different variables with the degree in
  584. each axis given by the corresponding index.
  585. m : int, optional
  586. Order of integration, must be positive. (Default: 1)
  587. k : {[], list, scalar}, optional
  588. Integration constant(s). The value of the first integral at
  589. ``lbnd`` is the first value in the list, the value of the second
  590. integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
  591. default), all constants are set to zero. If ``m == 1``, a single
  592. scalar can be given instead of a list.
  593. lbnd : scalar, optional
  594. The lower bound of the integral. (Default: 0)
  595. scl : scalar, optional
  596. Following each integration the result is *multiplied* by `scl`
  597. before the integration constant is added. (Default: 1)
  598. axis : int, optional
  599. Axis over which the integral is taken. (Default: 0).
  600. .. versionadded:: 1.7.0
  601. Returns
  602. -------
  603. S : ndarray
  604. Legendre series coefficient array of the integral.
  605. Raises
  606. ------
  607. ValueError
  608. If ``m < 0``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
  609. ``np.ndim(scl) != 0``.
  610. See Also
  611. --------
  612. legder
  613. Notes
  614. -----
  615. Note that the result of each integration is *multiplied* by `scl`.
  616. Why is this important to note? Say one is making a linear change of
  617. variable :math:`u = ax + b` in an integral relative to `x`. Then
  618. :math:`dx = du/a`, so one will need to set `scl` equal to
  619. :math:`1/a` - perhaps not what one would have first thought.
  620. Also note that, in general, the result of integrating a C-series needs
  621. to be "reprojected" onto the C-series basis set. Thus, typically,
  622. the result of this function is "unintuitive," albeit correct; see
  623. Examples section below.
  624. Examples
  625. --------
  626. >>> from numpy.polynomial import legendre as L
  627. >>> c = (1,2,3)
  628. >>> L.legint(c)
  629. array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
  630. >>> L.legint(c, 3)
  631. array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, # may vary
  632. -1.73472348e-18, 1.90476190e-02, 9.52380952e-03])
  633. >>> L.legint(c, k=3)
  634. array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
  635. >>> L.legint(c, lbnd=-2)
  636. array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary
  637. >>> L.legint(c, scl=2)
  638. array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) # may vary
  639. """
  640. c = np.array(c, ndmin=1, copy=True)
  641. if c.dtype.char in '?bBhHiIlLqQpP':
  642. c = c.astype(np.double)
  643. if not np.iterable(k):
  644. k = [k]
  645. cnt = pu._deprecate_as_int(m, "the order of integration")
  646. iaxis = pu._deprecate_as_int(axis, "the axis")
  647. if cnt < 0:
  648. raise ValueError("The order of integration must be non-negative")
  649. if len(k) > cnt:
  650. raise ValueError("Too many integration constants")
  651. if np.ndim(lbnd) != 0:
  652. raise ValueError("lbnd must be a scalar.")
  653. if np.ndim(scl) != 0:
  654. raise ValueError("scl must be a scalar.")
  655. iaxis = normalize_axis_index(iaxis, c.ndim)
  656. if cnt == 0:
  657. return c
  658. c = np.moveaxis(c, iaxis, 0)
  659. k = list(k) + [0]*(cnt - len(k))
  660. for i in range(cnt):
  661. n = len(c)
  662. c *= scl
  663. if n == 1 and np.all(c[0] == 0):
  664. c[0] += k[i]
  665. else:
  666. tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
  667. tmp[0] = c[0]*0
  668. tmp[1] = c[0]
  669. if n > 1:
  670. tmp[2] = c[1]/3
  671. for j in range(2, n):
  672. t = c[j]/(2*j + 1)
  673. tmp[j + 1] = t
  674. tmp[j - 1] -= t
  675. tmp[0] += k[i] - legval(lbnd, tmp)
  676. c = tmp
  677. c = np.moveaxis(c, 0, iaxis)
  678. return c
  679. def legval(x, c, tensor=True):
  680. """
  681. Evaluate a Legendre series at points x.
  682. If `c` is of length `n + 1`, this function returns the value:
  683. .. math:: p(x) = c_0 * L_0(x) + c_1 * L_1(x) + ... + c_n * L_n(x)
  684. The parameter `x` is converted to an array only if it is a tuple or a
  685. list, otherwise it is treated as a scalar. In either case, either `x`
  686. or its elements must support multiplication and addition both with
  687. themselves and with the elements of `c`.
  688. If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
  689. `c` is multidimensional, then the shape of the result depends on the
  690. value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
  691. x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
  692. scalars have shape (,).
  693. Trailing zeros in the coefficients will be used in the evaluation, so
  694. they should be avoided if efficiency is a concern.
  695. Parameters
  696. ----------
  697. x : array_like, compatible object
  698. If `x` is a list or tuple, it is converted to an ndarray, otherwise
  699. it is left unchanged and treated as a scalar. In either case, `x`
  700. or its elements must support addition and multiplication with
  701. themselves and with the elements of `c`.
  702. c : array_like
  703. Array of coefficients ordered so that the coefficients for terms of
  704. degree n are contained in c[n]. If `c` is multidimensional the
  705. remaining indices enumerate multiple polynomials. In the two
  706. dimensional case the coefficients may be thought of as stored in
  707. the columns of `c`.
  708. tensor : boolean, optional
  709. If True, the shape of the coefficient array is extended with ones
  710. on the right, one for each dimension of `x`. Scalars have dimension 0
  711. for this action. The result is that every column of coefficients in
  712. `c` is evaluated for every element of `x`. If False, `x` is broadcast
  713. over the columns of `c` for the evaluation. This keyword is useful
  714. when `c` is multidimensional. The default value is True.
  715. .. versionadded:: 1.7.0
  716. Returns
  717. -------
  718. values : ndarray, algebra_like
  719. The shape of the return value is described above.
  720. See Also
  721. --------
  722. legval2d, leggrid2d, legval3d, leggrid3d
  723. Notes
  724. -----
  725. The evaluation uses Clenshaw recursion, aka synthetic division.
  726. """
  727. c = np.array(c, ndmin=1, copy=False)
  728. if c.dtype.char in '?bBhHiIlLqQpP':
  729. c = c.astype(np.double)
  730. if isinstance(x, (tuple, list)):
  731. x = np.asarray(x)
  732. if isinstance(x, np.ndarray) and tensor:
  733. c = c.reshape(c.shape + (1,)*x.ndim)
  734. if len(c) == 1:
  735. c0 = c[0]
  736. c1 = 0
  737. elif len(c) == 2:
  738. c0 = c[0]
  739. c1 = c[1]
  740. else:
  741. nd = len(c)
  742. c0 = c[-2]
  743. c1 = c[-1]
  744. for i in range(3, len(c) + 1):
  745. tmp = c0
  746. nd = nd - 1
  747. c0 = c[-i] - (c1*(nd - 1))/nd
  748. c1 = tmp + (c1*x*(2*nd - 1))/nd
  749. return c0 + c1*x
  750. def legval2d(x, y, c):
  751. """
  752. Evaluate a 2-D Legendre series at points (x, y).
  753. This function returns the values:
  754. .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * L_i(x) * L_j(y)
  755. The parameters `x` and `y` are converted to arrays only if they are
  756. tuples or a lists, otherwise they are treated as a scalars and they
  757. must have the same shape after conversion. In either case, either `x`
  758. and `y` or their elements must support multiplication and addition both
  759. with themselves and with the elements of `c`.
  760. If `c` is a 1-D array a one is implicitly appended to its shape to make
  761. it 2-D. The shape of the result will be c.shape[2:] + x.shape.
  762. Parameters
  763. ----------
  764. x, y : array_like, compatible objects
  765. The two dimensional series is evaluated at the points `(x, y)`,
  766. where `x` and `y` must have the same shape. If `x` or `y` is a list
  767. or tuple, it is first converted to an ndarray, otherwise it is left
  768. unchanged and if it isn't an ndarray it is treated as a scalar.
  769. c : array_like
  770. Array of coefficients ordered so that the coefficient of the term
  771. of multi-degree i,j is contained in ``c[i,j]``. If `c` has
  772. dimension greater than two the remaining indices enumerate multiple
  773. sets of coefficients.
  774. Returns
  775. -------
  776. values : ndarray, compatible object
  777. The values of the two dimensional Legendre series at points formed
  778. from pairs of corresponding values from `x` and `y`.
  779. See Also
  780. --------
  781. legval, leggrid2d, legval3d, leggrid3d
  782. Notes
  783. -----
  784. .. versionadded:: 1.7.0
  785. """
  786. return pu._valnd(legval, c, x, y)
  787. def leggrid2d(x, y, c):
  788. """
  789. Evaluate a 2-D Legendre series on the Cartesian product of x and y.
  790. This function returns the values:
  791. .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * L_i(a) * L_j(b)
  792. where the points `(a, b)` consist of all pairs formed by taking
  793. `a` from `x` and `b` from `y`. The resulting points form a grid with
  794. `x` in the first dimension and `y` in the second.
  795. The parameters `x` and `y` are converted to arrays only if they are
  796. tuples or a lists, otherwise they are treated as a scalars. In either
  797. case, either `x` and `y` or their elements must support multiplication
  798. and addition both with themselves and with the elements of `c`.
  799. If `c` has fewer than two dimensions, ones are implicitly appended to
  800. its shape to make it 2-D. The shape of the result will be c.shape[2:] +
  801. x.shape + y.shape.
  802. Parameters
  803. ----------
  804. x, y : array_like, compatible objects
  805. The two dimensional series is evaluated at the points in the
  806. Cartesian product of `x` and `y`. If `x` or `y` is a list or
  807. tuple, it is first converted to an ndarray, otherwise it is left
  808. unchanged and, if it isn't an ndarray, it is treated as a scalar.
  809. c : array_like
  810. Array of coefficients ordered so that the coefficient of the term of
  811. multi-degree i,j is contained in `c[i,j]`. If `c` has dimension
  812. greater than two the remaining indices enumerate multiple sets of
  813. coefficients.
  814. Returns
  815. -------
  816. values : ndarray, compatible object
  817. The values of the two dimensional Chebyshev series at points in the
  818. Cartesian product of `x` and `y`.
  819. See Also
  820. --------
  821. legval, legval2d, legval3d, leggrid3d
  822. Notes
  823. -----
  824. .. versionadded:: 1.7.0
  825. """
  826. return pu._gridnd(legval, c, x, y)
  827. def legval3d(x, y, z, c):
  828. """
  829. Evaluate a 3-D Legendre series at points (x, y, z).
  830. This function returns the values:
  831. .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * L_i(x) * L_j(y) * L_k(z)
  832. The parameters `x`, `y`, and `z` are converted to arrays only if
  833. they are tuples or a lists, otherwise they are treated as a scalars and
  834. they must have the same shape after conversion. In either case, either
  835. `x`, `y`, and `z` or their elements must support multiplication and
  836. addition both with themselves and with the elements of `c`.
  837. If `c` has fewer than 3 dimensions, ones are implicitly appended to its
  838. shape to make it 3-D. The shape of the result will be c.shape[3:] +
  839. x.shape.
  840. Parameters
  841. ----------
  842. x, y, z : array_like, compatible object
  843. The three dimensional series is evaluated at the points
  844. `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
  845. any of `x`, `y`, or `z` is a list or tuple, it is first converted
  846. to an ndarray, otherwise it is left unchanged and if it isn't an
  847. ndarray it is treated as a scalar.
  848. c : array_like
  849. Array of coefficients ordered so that the coefficient of the term of
  850. multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
  851. greater than 3 the remaining indices enumerate multiple sets of
  852. coefficients.
  853. Returns
  854. -------
  855. values : ndarray, compatible object
  856. The values of the multidimensional polynomial on points formed with
  857. triples of corresponding values from `x`, `y`, and `z`.
  858. See Also
  859. --------
  860. legval, legval2d, leggrid2d, leggrid3d
  861. Notes
  862. -----
  863. .. versionadded:: 1.7.0
  864. """
  865. return pu._valnd(legval, c, x, y, z)
  866. def leggrid3d(x, y, z, c):
  867. """
  868. Evaluate a 3-D Legendre series on the Cartesian product of x, y, and z.
  869. This function returns the values:
  870. .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * L_i(a) * L_j(b) * L_k(c)
  871. where the points `(a, b, c)` consist of all triples formed by taking
  872. `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
  873. a grid with `x` in the first dimension, `y` in the second, and `z` in
  874. the third.
  875. The parameters `x`, `y`, and `z` are converted to arrays only if they
  876. are tuples or a lists, otherwise they are treated as a scalars. In
  877. either case, either `x`, `y`, and `z` or their elements must support
  878. multiplication and addition both with themselves and with the elements
  879. of `c`.
  880. If `c` has fewer than three dimensions, ones are implicitly appended to
  881. its shape to make it 3-D. The shape of the result will be c.shape[3:] +
  882. x.shape + y.shape + z.shape.
  883. Parameters
  884. ----------
  885. x, y, z : array_like, compatible objects
  886. The three dimensional series is evaluated at the points in the
  887. Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
  888. list or tuple, it is first converted to an ndarray, otherwise it is
  889. left unchanged and, if it isn't an ndarray, it is treated as a
  890. scalar.
  891. c : array_like
  892. Array of coefficients ordered so that the coefficients for terms of
  893. degree i,j are contained in ``c[i,j]``. If `c` has dimension
  894. greater than two the remaining indices enumerate multiple sets of
  895. coefficients.
  896. Returns
  897. -------
  898. values : ndarray, compatible object
  899. The values of the two dimensional polynomial at points in the Cartesian
  900. product of `x` and `y`.
  901. See Also
  902. --------
  903. legval, legval2d, leggrid2d, legval3d
  904. Notes
  905. -----
  906. .. versionadded:: 1.7.0
  907. """
  908. return pu._gridnd(legval, c, x, y, z)
  909. def legvander(x, deg):
  910. """Pseudo-Vandermonde matrix of given degree.
  911. Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
  912. `x`. The pseudo-Vandermonde matrix is defined by
  913. .. math:: V[..., i] = L_i(x)
  914. where `0 <= i <= deg`. The leading indices of `V` index the elements of
  915. `x` and the last index is the degree of the Legendre polynomial.
  916. If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
  917. array ``V = legvander(x, n)``, then ``np.dot(V, c)`` and
  918. ``legval(x, c)`` are the same up to roundoff. This equivalence is
  919. useful both for least squares fitting and for the evaluation of a large
  920. number of Legendre series of the same degree and sample points.
  921. Parameters
  922. ----------
  923. x : array_like
  924. Array of points. The dtype is converted to float64 or complex128
  925. depending on whether any of the elements are complex. If `x` is
  926. scalar it is converted to a 1-D array.
  927. deg : int
  928. Degree of the resulting matrix.
  929. Returns
  930. -------
  931. vander : ndarray
  932. The pseudo-Vandermonde matrix. The shape of the returned matrix is
  933. ``x.shape + (deg + 1,)``, where The last index is the degree of the
  934. corresponding Legendre polynomial. The dtype will be the same as
  935. the converted `x`.
  936. """
  937. ideg = pu._deprecate_as_int(deg, "deg")
  938. if ideg < 0:
  939. raise ValueError("deg must be non-negative")
  940. x = np.array(x, copy=False, ndmin=1) + 0.0
  941. dims = (ideg + 1,) + x.shape
  942. dtyp = x.dtype
  943. v = np.empty(dims, dtype=dtyp)
  944. # Use forward recursion to generate the entries. This is not as accurate
  945. # as reverse recursion in this application but it is more efficient.
  946. v[0] = x*0 + 1
  947. if ideg > 0:
  948. v[1] = x
  949. for i in range(2, ideg + 1):
  950. v[i] = (v[i-1]*x*(2*i - 1) - v[i-2]*(i - 1))/i
  951. return np.moveaxis(v, 0, -1)
  952. def legvander2d(x, y, deg):
  953. """Pseudo-Vandermonde matrix of given degrees.
  954. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
  955. points `(x, y)`. The pseudo-Vandermonde matrix is defined by
  956. .. math:: V[..., (deg[1] + 1)*i + j] = L_i(x) * L_j(y),
  957. where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
  958. `V` index the points `(x, y)` and the last index encodes the degrees of
  959. the Legendre polynomials.
  960. If ``V = legvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
  961. correspond to the elements of a 2-D coefficient array `c` of shape
  962. (xdeg + 1, ydeg + 1) in the order
  963. .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
  964. and ``np.dot(V, c.flat)`` and ``legval2d(x, y, c)`` will be the same
  965. up to roundoff. This equivalence is useful both for least squares
  966. fitting and for the evaluation of a large number of 2-D Legendre
  967. series of the same degrees and sample points.
  968. Parameters
  969. ----------
  970. x, y : array_like
  971. Arrays of point coordinates, all of the same shape. The dtypes
  972. will be converted to either float64 or complex128 depending on
  973. whether any of the elements are complex. Scalars are converted to
  974. 1-D arrays.
  975. deg : list of ints
  976. List of maximum degrees of the form [x_deg, y_deg].
  977. Returns
  978. -------
  979. vander2d : ndarray
  980. The shape of the returned matrix is ``x.shape + (order,)``, where
  981. :math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same
  982. as the converted `x` and `y`.
  983. See Also
  984. --------
  985. legvander, legvander3d, legval2d, legval3d
  986. Notes
  987. -----
  988. .. versionadded:: 1.7.0
  989. """
  990. return pu._vander_nd_flat((legvander, legvander), (x, y), deg)
  991. def legvander3d(x, y, z, deg):
  992. """Pseudo-Vandermonde matrix of given degrees.
  993. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
  994. points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
  995. then The pseudo-Vandermonde matrix is defined by
  996. .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = L_i(x)*L_j(y)*L_k(z),
  997. where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
  998. indices of `V` index the points `(x, y, z)` and the last index encodes
  999. the degrees of the Legendre polynomials.
  1000. If ``V = legvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
  1001. of `V` correspond to the elements of a 3-D coefficient array `c` of
  1002. shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
  1003. .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
  1004. and ``np.dot(V, c.flat)`` and ``legval3d(x, y, z, c)`` will be the
  1005. same up to roundoff. This equivalence is useful both for least squares
  1006. fitting and for the evaluation of a large number of 3-D Legendre
  1007. series of the same degrees and sample points.
  1008. Parameters
  1009. ----------
  1010. x, y, z : array_like
  1011. Arrays of point coordinates, all of the same shape. The dtypes will
  1012. be converted to either float64 or complex128 depending on whether
  1013. any of the elements are complex. Scalars are converted to 1-D
  1014. arrays.
  1015. deg : list of ints
  1016. List of maximum degrees of the form [x_deg, y_deg, z_deg].
  1017. Returns
  1018. -------
  1019. vander3d : ndarray
  1020. The shape of the returned matrix is ``x.shape + (order,)``, where
  1021. :math:`order = (deg[0]+1)*(deg[1]+1)*(deg[2]+1)`. The dtype will
  1022. be the same as the converted `x`, `y`, and `z`.
  1023. See Also
  1024. --------
  1025. legvander, legvander3d, legval2d, legval3d
  1026. Notes
  1027. -----
  1028. .. versionadded:: 1.7.0
  1029. """
  1030. return pu._vander_nd_flat((legvander, legvander, legvander), (x, y, z), deg)
  1031. def legfit(x, y, deg, rcond=None, full=False, w=None):
  1032. """
  1033. Least squares fit of Legendre series to data.
  1034. Return the coefficients of a Legendre series of degree `deg` that is the
  1035. least squares fit to the data values `y` given at points `x`. If `y` is
  1036. 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
  1037. fits are done, one for each column of `y`, and the resulting
  1038. coefficients are stored in the corresponding columns of a 2-D return.
  1039. The fitted polynomial(s) are in the form
  1040. .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x),
  1041. where `n` is `deg`.
  1042. Parameters
  1043. ----------
  1044. x : array_like, shape (M,)
  1045. x-coordinates of the M sample points ``(x[i], y[i])``.
  1046. y : array_like, shape (M,) or (M, K)
  1047. y-coordinates of the sample points. Several data sets of sample
  1048. points sharing the same x-coordinates can be fitted at once by
  1049. passing in a 2D-array that contains one dataset per column.
  1050. deg : int or 1-D array_like
  1051. Degree(s) of the fitting polynomials. If `deg` is a single integer
  1052. all terms up to and including the `deg`'th term are included in the
  1053. fit. For NumPy versions >= 1.11.0 a list of integers specifying the
  1054. degrees of the terms to include may be used instead.
  1055. rcond : float, optional
  1056. Relative condition number of the fit. Singular values smaller than
  1057. this relative to the largest singular value will be ignored. The
  1058. default value is len(x)*eps, where eps is the relative precision of
  1059. the float type, about 2e-16 in most cases.
  1060. full : bool, optional
  1061. Switch determining nature of return value. When it is False (the
  1062. default) just the coefficients are returned, when True diagnostic
  1063. information from the singular value decomposition is also returned.
  1064. w : array_like, shape (`M`,), optional
  1065. Weights. If not None, the weight ``w[i]`` applies to the unsquared
  1066. residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
  1067. chosen so that the errors of the products ``w[i]*y[i]`` all have the
  1068. same variance. When using inverse-variance weighting, use
  1069. ``w[i] = 1/sigma(y[i])``. The default value is None.
  1070. .. versionadded:: 1.5.0
  1071. Returns
  1072. -------
  1073. coef : ndarray, shape (M,) or (M, K)
  1074. Legendre coefficients ordered from low to high. If `y` was
  1075. 2-D, the coefficients for the data in column k of `y` are in
  1076. column `k`. If `deg` is specified as a list, coefficients for
  1077. terms not included in the fit are set equal to zero in the
  1078. returned `coef`.
  1079. [residuals, rank, singular_values, rcond] : list
  1080. These values are only returned if ``full == True``
  1081. - residuals -- sum of squared residuals of the least squares fit
  1082. - rank -- the numerical rank of the scaled Vandermonde matrix
  1083. - singular_values -- singular values of the scaled Vandermonde matrix
  1084. - rcond -- value of `rcond`.
  1085. For more details, see `numpy.linalg.lstsq`.
  1086. Warns
  1087. -----
  1088. RankWarning
  1089. The rank of the coefficient matrix in the least-squares fit is
  1090. deficient. The warning is only raised if ``full == False``. The
  1091. warnings can be turned off by
  1092. >>> import warnings
  1093. >>> warnings.simplefilter('ignore', np.RankWarning)
  1094. See Also
  1095. --------
  1096. numpy.polynomial.polynomial.polyfit
  1097. numpy.polynomial.chebyshev.chebfit
  1098. numpy.polynomial.laguerre.lagfit
  1099. numpy.polynomial.hermite.hermfit
  1100. numpy.polynomial.hermite_e.hermefit
  1101. legval : Evaluates a Legendre series.
  1102. legvander : Vandermonde matrix of Legendre series.
  1103. legweight : Legendre weight function (= 1).
  1104. numpy.linalg.lstsq : Computes a least-squares fit from the matrix.
  1105. scipy.interpolate.UnivariateSpline : Computes spline fits.
  1106. Notes
  1107. -----
  1108. The solution is the coefficients of the Legendre series `p` that
  1109. minimizes the sum of the weighted squared errors
  1110. .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
  1111. where :math:`w_j` are the weights. This problem is solved by setting up
  1112. as the (typically) overdetermined matrix equation
  1113. .. math:: V(x) * c = w * y,
  1114. where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
  1115. coefficients to be solved for, `w` are the weights, and `y` are the
  1116. observed values. This equation is then solved using the singular value
  1117. decomposition of `V`.
  1118. If some of the singular values of `V` are so small that they are
  1119. neglected, then a `RankWarning` will be issued. This means that the
  1120. coefficient values may be poorly determined. Using a lower order fit
  1121. will usually get rid of the warning. The `rcond` parameter can also be
  1122. set to a value smaller than its default, but the resulting fit may be
  1123. spurious and have large contributions from roundoff error.
  1124. Fits using Legendre series are usually better conditioned than fits
  1125. using power series, but much can depend on the distribution of the
  1126. sample points and the smoothness of the data. If the quality of the fit
  1127. is inadequate splines may be a good alternative.
  1128. References
  1129. ----------
  1130. .. [1] Wikipedia, "Curve fitting",
  1131. https://en.wikipedia.org/wiki/Curve_fitting
  1132. Examples
  1133. --------
  1134. """
  1135. return pu._fit(legvander, x, y, deg, rcond, full, w)
  1136. def legcompanion(c):
  1137. """Return the scaled companion matrix of c.
  1138. The basis polynomials are scaled so that the companion matrix is
  1139. symmetric when `c` is an Legendre basis polynomial. This provides
  1140. better eigenvalue estimates than the unscaled case and for basis
  1141. polynomials the eigenvalues are guaranteed to be real if
  1142. `numpy.linalg.eigvalsh` is used to obtain them.
  1143. Parameters
  1144. ----------
  1145. c : array_like
  1146. 1-D array of Legendre series coefficients ordered from low to high
  1147. degree.
  1148. Returns
  1149. -------
  1150. mat : ndarray
  1151. Scaled companion matrix of dimensions (deg, deg).
  1152. Notes
  1153. -----
  1154. .. versionadded:: 1.7.0
  1155. """
  1156. # c is a trimmed copy
  1157. [c] = pu.as_series([c])
  1158. if len(c) < 2:
  1159. raise ValueError('Series must have maximum degree of at least 1.')
  1160. if len(c) == 2:
  1161. return np.array([[-c[0]/c[1]]])
  1162. n = len(c) - 1
  1163. mat = np.zeros((n, n), dtype=c.dtype)
  1164. scl = 1./np.sqrt(2*np.arange(n) + 1)
  1165. top = mat.reshape(-1)[1::n+1]
  1166. bot = mat.reshape(-1)[n::n+1]
  1167. top[...] = np.arange(1, n)*scl[:n-1]*scl[1:n]
  1168. bot[...] = top
  1169. mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*(n/(2*n - 1))
  1170. return mat
  1171. def legroots(c):
  1172. """
  1173. Compute the roots of a Legendre series.
  1174. Return the roots (a.k.a. "zeros") of the polynomial
  1175. .. math:: p(x) = \\sum_i c[i] * L_i(x).
  1176. Parameters
  1177. ----------
  1178. c : 1-D array_like
  1179. 1-D array of coefficients.
  1180. Returns
  1181. -------
  1182. out : ndarray
  1183. Array of the roots of the series. If all the roots are real,
  1184. then `out` is also real, otherwise it is complex.
  1185. See Also
  1186. --------
  1187. numpy.polynomial.polynomial.polyroots
  1188. numpy.polynomial.chebyshev.chebroots
  1189. numpy.polynomial.laguerre.lagroots
  1190. numpy.polynomial.hermite.hermroots
  1191. numpy.polynomial.hermite_e.hermeroots
  1192. Notes
  1193. -----
  1194. The root estimates are obtained as the eigenvalues of the companion
  1195. matrix, Roots far from the origin of the complex plane may have large
  1196. errors due to the numerical instability of the series for such values.
  1197. Roots with multiplicity greater than 1 will also show larger errors as
  1198. the value of the series near such points is relatively insensitive to
  1199. errors in the roots. Isolated roots near the origin can be improved by
  1200. a few iterations of Newton's method.
  1201. The Legendre series basis polynomials aren't powers of ``x`` so the
  1202. results of this function may seem unintuitive.
  1203. Examples
  1204. --------
  1205. >>> import numpy.polynomial.legendre as leg
  1206. >>> leg.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0, all real roots
  1207. array([-0.85099543, -0.11407192, 0.51506735]) # may vary
  1208. """
  1209. # c is a trimmed copy
  1210. [c] = pu.as_series([c])
  1211. if len(c) < 2:
  1212. return np.array([], dtype=c.dtype)
  1213. if len(c) == 2:
  1214. return np.array([-c[0]/c[1]])
  1215. # rotated companion matrix reduces error
  1216. m = legcompanion(c)[::-1,::-1]
  1217. r = la.eigvals(m)
  1218. r.sort()
  1219. return r
  1220. def leggauss(deg):
  1221. """
  1222. Gauss-Legendre quadrature.
  1223. Computes the sample points and weights for Gauss-Legendre quadrature.
  1224. These sample points and weights will correctly integrate polynomials of
  1225. degree :math:`2*deg - 1` or less over the interval :math:`[-1, 1]` with
  1226. the weight function :math:`f(x) = 1`.
  1227. Parameters
  1228. ----------
  1229. deg : int
  1230. Number of sample points and weights. It must be >= 1.
  1231. Returns
  1232. -------
  1233. x : ndarray
  1234. 1-D ndarray containing the sample points.
  1235. y : ndarray
  1236. 1-D ndarray containing the weights.
  1237. Notes
  1238. -----
  1239. .. versionadded:: 1.7.0
  1240. The results have only been tested up to degree 100, higher degrees may
  1241. be problematic. The weights are determined by using the fact that
  1242. .. math:: w_k = c / (L'_n(x_k) * L_{n-1}(x_k))
  1243. where :math:`c` is a constant independent of :math:`k` and :math:`x_k`
  1244. is the k'th root of :math:`L_n`, and then scaling the results to get
  1245. the right value when integrating 1.
  1246. """
  1247. ideg = pu._deprecate_as_int(deg, "deg")
  1248. if ideg <= 0:
  1249. raise ValueError("deg must be a positive integer")
  1250. # first approximation of roots. We use the fact that the companion
  1251. # matrix is symmetric in this case in order to obtain better zeros.
  1252. c = np.array([0]*deg + [1])
  1253. m = legcompanion(c)
  1254. x = la.eigvalsh(m)
  1255. # improve roots by one application of Newton
  1256. dy = legval(x, c)
  1257. df = legval(x, legder(c))
  1258. x -= dy/df
  1259. # compute the weights. We scale the factor to avoid possible numerical
  1260. # overflow.
  1261. fm = legval(x, c[1:])
  1262. fm /= np.abs(fm).max()
  1263. df /= np.abs(df).max()
  1264. w = 1/(fm * df)
  1265. # for Legendre we can also symmetrize
  1266. w = (w + w[::-1])/2
  1267. x = (x - x[::-1])/2
  1268. # scale w to get the right value
  1269. w *= 2. / w.sum()
  1270. return x, w
  1271. def legweight(x):
  1272. """
  1273. Weight function of the Legendre polynomials.
  1274. The weight function is :math:`1` and the interval of integration is
  1275. :math:`[-1, 1]`. The Legendre polynomials are orthogonal, but not
  1276. normalized, with respect to this weight function.
  1277. Parameters
  1278. ----------
  1279. x : array_like
  1280. Values at which the weight function will be computed.
  1281. Returns
  1282. -------
  1283. w : ndarray
  1284. The weight function at `x`.
  1285. Notes
  1286. -----
  1287. .. versionadded:: 1.7.0
  1288. """
  1289. w = x*0.0 + 1.0
  1290. return w
  1291. #
  1292. # Legendre series class
  1293. #
  1294. class Legendre(ABCPolyBase):
  1295. """A Legendre series class.
  1296. The Legendre class provides the standard Python numerical methods
  1297. '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
  1298. attributes and methods listed in the `ABCPolyBase` documentation.
  1299. Parameters
  1300. ----------
  1301. coef : array_like
  1302. Legendre coefficients in order of increasing degree, i.e.,
  1303. ``(1, 2, 3)`` gives ``1*P_0(x) + 2*P_1(x) + 3*P_2(x)``.
  1304. domain : (2,) array_like, optional
  1305. Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
  1306. to the interval ``[window[0], window[1]]`` by shifting and scaling.
  1307. The default value is [-1, 1].
  1308. window : (2,) array_like, optional
  1309. Window, see `domain` for its use. The default value is [-1, 1].
  1310. .. versionadded:: 1.6.0
  1311. symbol : str, optional
  1312. Symbol used to represent the independent variable in string
  1313. representations of the polynomial expression, e.g. for printing.
  1314. The symbol must be a valid Python identifier. Default value is 'x'.
  1315. .. versionadded:: 1.24
  1316. """
  1317. # Virtual Functions
  1318. _add = staticmethod(legadd)
  1319. _sub = staticmethod(legsub)
  1320. _mul = staticmethod(legmul)
  1321. _div = staticmethod(legdiv)
  1322. _pow = staticmethod(legpow)
  1323. _val = staticmethod(legval)
  1324. _int = staticmethod(legint)
  1325. _der = staticmethod(legder)
  1326. _fit = staticmethod(legfit)
  1327. _line = staticmethod(legline)
  1328. _roots = staticmethod(legroots)
  1329. _fromroots = staticmethod(legfromroots)
  1330. # Virtual properties
  1331. domain = np.array(legdomain)
  1332. window = np.array(legdomain)
  1333. basis_name = 'P'