Re: lspoly fails with higher order polynomials

>
>I'm trying to run lspoly to get a good polynomial interpolation. Let's
>consider this example :
>
>------------------------
>begin
> x = (/0.,1.,2.,3.,4.,5.,6./)
> y = (/0.,1.,2.,3.,4.,5.,6./)
> c = lspoly(x,y,1,5)
> print(c)
>end
>------------------------
>
>which returns obviously the output (y = x) :
>
>------------------------
>Variable: c
>Type: float
>Total Size: 20 bytes
> 5 values
>Number of Dimensions: 1
>Dimensions and sizes: [5]
>Coordinates:
>(0) 0
>(1) 1
>(2) 0
>(3) 0
>(4) 0
>------------------------
>
>
>But, if I just change the order of the requested polynomial from 4 to 5 :
>
>------------------------
> c = lspoly(x,y,1,6)
>------------------------
>
>I get the result :
>
>------------------------
>Dimensions and sizes: [6]
>Coordinates:
>(0) 1e+20
>(1) 1e+20
>(2) 1e+20
>(3) 1e+20
>(4) 1e+20
>(5) 1e+20
>------------------------
>
>Is this a bug that could be fixed soon ? I would really need a higher
>order polynomial.
>_______________________________________________

It is not a bug per se. It is a limitation of the approach
used to estimate the polynomials. The advantage was low memory
and speed.

The original fortran subroutine documentation says:

C ACCURACY FOR LOWER ORDER POLYNOMIALS (N .LE. 5), LSPOLY
C CAN BE EXPECTED TO GIVE SATISFACTORY RESULTS.
C AS THE DEGREE OF THE FITTING POLYNOMIAL IS
C INCREASED, THE PROBLEM BECOMES NUMERICALLY
C UNSTABLE. IF A HIGHER DEGREE POLYNOMIAL FIT
C IS DESIRED, THEN ROUTINE RGRSN1 SHOULD
C BE USED.
  

===
A person is looking into this.
===

If not, perhaps you could call a NAG routine [eg: E02ADF ]
the following shows how to invoke an IMSL code.

http://www.ncl.ucar.edu/Support/talk_archives/2005/0677.html

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Jan 23 2007 - 15:35:08 MST