[SciPy-user] Numerical stability and Legendre polynomials

[Fernando Perez]

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