[Maxima] ratinterpol always rats?

Raymond Toy toy.raymond at gmail.com
Sat Nov 8 07:27:17 CST 2008


Mario Rodriguez wrote:
>   
>>> I introduced the rat call in the lagrange function to force rational
>>> arithmetic due to the inestability of floating point calculations in
>>> high degree polynomials.
>>>       
>> I don't think using rat fixes the stability problem. :-)
>>     
>
> Hi Ray,
>
> Let me be more specific. The idea behind introducing rat was to solve
> this type of results:
>
> Let function lagrange2 be the same as lagrange, but without applying rat
> to input data.
>
>
> (%i9) 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] ]$
>
> (%i10) subst(x=140, lagrange2(m));
> (%o10)                              68192.0
>
> With some more precision I get a completely different result:
>
> (%i11) (fpprec:25, subst(x=140, lagrange2(bfloat(m))));
> (%o11)                   1.572003180533647537231445b1
>
>
> This is what we now get with the lagrange function as defined in
> interpol:
>
> (%i12) float( subst(x=140, lagrange(m)) );
> (%o12)                               15.72
>
>
>   
Neat.  But what about subst([x=140.000001], lagrange(m))?  I get 10112,
which is very far from 15.72.  Evaluating using Horner's rule doesn't
help either.  I get -5870.546.  Of course, using 140+1/1000000 gives an
answer that's much closer to 15.72.

I think this is just the nature of polynomial evaluation.  Also, isn't
it usually better (numerically) to leave the lagrange interpolating
polynomial in Lagrange form, that is, as a sum of products?
> I took this decision after a long thread in the list some time ago,
> although I am not sure this is the best solution.
>
> I'll introduce option force_rat and make documentation more clear on
> this point.
>   
I don't think force_rat is necessary.  But a note in the docs saying
floats are converted to rationals before computing the results would be
nice. :-)

Ray



More information about the Maxima mailing list