Problem with expand in lagrange interpolation

I think that your problem is that you (or this routine) are using floating
point numbers to compute values that should be computed more accurately if
you want the numbers to be correct.

The difference between 150^10   and 150.0^10  is about 91136

You might find a divided-difference formula to be better.
RJF


> -----Original Message-----
> From: maxima-bounces at math.utexas.edu 
> [mailto:maxima-bounces at math.utexas.edu] On Behalf Of Mario Rodriguez
> Sent: Wednesday, June 13, 2007 8:54 AM
> To: maxima at math.utexas.edu
> Subject: Problem with expand in lagrange interpolation
> 
> I have found a problem with lagrange interpolation.
> 
> Please, look at this output:
> 
> (%i11) m:[[140, 15.72], [141, 15.53], [142, 15.19],
>           [143, 16.56], [144, 16.21], [145, 17.39],
>           [146, 17.36], [147, 17.42], [148, 17.60],
>           [149, 17.75], [150, 18.95] ]$
> (%i12) load(interpol)$
> (%i13) poly:lagrange(m)$
> (%i14) subst(x=140,poly);
> (%o14)                              68192.0
> 
> 
> The result should be 15.72. The origin of this incorrect output is the
> last sentence in the lagrange routine:
> 
> 
> .........
>    for i:1 thru n do(
>       prod: 1,
>       for k:1 thru n do
>          if k#i then prod: prod * (defaults[1]-tab[k][1]) /
> (tab[i][1]-tab[k][1]),
>       sum: sum + prod * tab[i][2] ),
>    expand(sum))$
> 
> If I don't call expand and just simply throw sum, the output is now
> correct:
> 
> (%i15) load(interpol)$
> (%i16) poly:lagrange(m)$
> (%i17) subst(x=140,poly);
> (%o14)                               15.72
> 
> The error above happens sometimes, mostly when first 
> coordinates are big
> numbers.
> 
> When lagrange returns sum without postprocessing, the output 
> is not very
> nice (that's why I introduced expand), but accurate:
> 
> (%i15) poly;
> (%o15) 5.222111992945325E-6 (x - 149) (x - 148) (x - 147) (x 
> - 146) (x -
> 145)
>  (x - 144) (x - 143) (x - 142) (x - 141) (x - 140)
>  - 4.891424162257496E-5 (x - 150) (x - 148) (x - 147) (x - 146) (x -
> 145)
>  (x - 144) (x - 143) (x - 142) (x - 141) (x - 140)
>  + 2.1825396825396828E-4 (x - 150) (x - 149) (x - 147) (x - 146) (x -
> 145)
>  (x - 144) (x - 143) (x - 142) (x - 141) (x - 140)
>  - 5.760582010582012E-4 (x - 150) (x - 149) (x - 148) (x - 146) (x -
> 145)
>  (x - 144) (x - 143) (x - 142) (x - 141) (x - 140)
>  + 0.00100462962962963 (x - 150) (x - 149) (x - 148) (x - 
> 147) (x - 145)
>  (x - 144) (x - 143) (x - 142) (x - 141) (x - 140)
>  - .001207638888888889 (x - 150) (x - 149) (x - 148) (x - 
> 147) (x - 146)
>  (x - 144) (x - 143) (x - 142) (x - 141) (x - 140)
>  + 9.380787037037037E-4 (x - 150) (x - 149) (x - 148) (x - 147) (x -
> 146)
>  (x - 145) (x - 143) (x - 142) (x - 141) (x - 140)
>  - 5.476190476190477E-4 (x - 150) (x - 149) (x - 148) (x - 147) (x -
> 146)
>  (x - 145) (x - 144) (x - 142) (x - 141) (x - 140)
>  + 1.8836805555555555E-4 (x - 150) (x - 149) (x - 148) (x - 147) (x -
> 146)
>  (x - 145) (x - 144) (x - 143) (x - 141) (x - 140)
>  - 4.2796516754850084E-5 (x - 150) (x - 149) (x - 148) (x - 147) (x -
> 146)
>  (x - 145) (x - 144) (x - 143) (x - 142) (x - 140)
>  - 4.3320105820105814E-6 (141 - x) (x - 150) (x - 149) (x - 148) (x -
> 147)
>  (x - 146) (x - 145) (x - 144) (x - 143) (x - 142)
> 
> 
> 
> I have also tested ratsimp(sum), but the output doesn't look 
> nicer. As a
> conclusion, I think lagrange should return simply sum; or do 
> you have a
> better idea?
> 
> -- 
> Mario Rodriguez Riotorto
> www.biomates.net
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>