Yaroslav Bulatov wrote:

I'm trying to approximate a function using an expansion over few dozen of Legendre polynomials, but I'm running into a large numerical error.

For instance the following special.legendre(30).integ()(1)-special.legendre(30).integ()(-1) evaluates to 2e-4, instead of 0

So my question is -- is such numerical error normal? Is it realistic to try to approximate a function using 30 Legendre polynomials?

The problem, I think, is that special.legendre returns a poly1d object. The coefficients are reasonably well computed, as this comparison (for P_30) with Mathematica shows: Python: In [31]: p30coefs Out[31]: [-0.14446444809436668, 1.5679191456676283e-14, 67.175968363880514, -2.180105233807014e-12, -5172.5495640188101, 7.7713427093987573e-11, 156900.67010857066, 4.5831526974104543e-09, -2487996.3402931187, 3.887890884521149e-07, 23718898.444125757, 3.1658444632926091e-06, -147344672.15291381, -1.590865831210945e-05, 626619649.70533192, -0.00028178808880156549, -1879858949.11516, -0.00073104158643097776, 4042311073.5892482, 3.7885858888136625e-05, -6254944503.3432274, 0.0010140695072848572, 6904808867.3253412, 0.00041825387498683678, -5303693767.6564827, 6.9383401445695963e-05, 2692644528.1947165, 8.5971433363986932e-06, -812067397.39206791, 2.6902196481091132e-07, 110142474.588809] Mathematica: N[LegendreP[30, x], 20] \!\(\(-0.14446444809436798095703125`20. \) + 67.17596836388111114501953125`20. \ x\^2 - 5172.54956401884555816650390625`20. \ x\^4 + 156900.67010857164859771728515625`20. \ x\^6 - 2.48799634029306471347808837890625`20.*^6\ x\^8 + 2.371889844412721693515777587890625`20.*^7\ x\^10 - 1.4734467215291149914264678955078125`20.*^8\ x\^12 + 6.2661964970523901283740997314453125`20.*^8\ x\^14 - 1.87985894911571703851222991943359375`20.*^9\ x\^16 + 4.04231107358869872987270355224609375`20.*^9\ x\^18 - 6.25494450334251277148723602294921875`20.*^9\ x\^20 + 6.90480886732615046203136444091796875`20.*^9\ x\^22 - 5.30369376765631847083568572998046875`20.*^9\ x\^24 + 2.69264452819474630057811737060546875`20.*^9\ x\^26 - 8.1206739739206634461879730224609375`20.*^8\ x\^28 + 1.1014247458880899846553802490234375`20.*^8\ x\^30\) But the problem is that this is an unstable polynomial to compute directly, because of the alternation in signs. Since poly1d doesn't know anything about this fact, it happily computes an unstable quantity. You'll also notice that the odd coefficients, which are analytically zero (for even order Legendre polys), are not quite zero in this construction, another source of error. If you need accurate results with Legendre polynomials as the order increases, your best bet is to use the standard 3-term recursion, which is stable. A braindead implmenentation of the algorithm, as taken from any textbook, gives this: In [47]: p30_end = utils.plegnv(array([-1.0,1.0]),30) In [48]: for n,ends in enumerate(p30_end): ....: print ends[0]-((-1)**n)*ends[-1] ....: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 This is checking all 30 by subtracting the endpoints, accounting for even/oddness. Here's the same calculation using scipy: In [33]: for n in range(30): ....: pn = scipy.special.legendre(n) ....: print pn(-1.0)-((-1)**n)*pn(1.0) ....: 0.0 0.0 0.0 0.0 -1.11022302463e-15 -4.55191440096e-15 -1.40998324127e-14 -4.4408920985e-15 -4.41868763801e-14 -3.41948691585e-14 -2.05613304161e-13 -1.99495975295e-12 2.8537172625e-12 -3.97726296342e-12 -6.78088696304e-11 1.88647986121e-10 3.10480530175e-10 -9.75348690702e-11 -5.69281288776e-09 -6.65902566421e-09 -7.10117353808e-08 1.19151672973e-07 7.9948702647e-07 1.12660495866e-06 -3.17713138309e-06 1.19201396691e-05 8.83580242022e-05 -4.29196349421e-05 0.000923915025271 0.000167352152224 As you can see, the 3-term recursion fares much, much better. HTH, f