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