[Scilab-users] Cubic spline

[Claus Futtrup]

Hi Rafael Guerra, et al.

Thank you for clarifying this to me - that interp1 uses not-a-knot as 
default. IMHO Scilab documentation should be clear about this.

 >What do you mean by your result is “OK, but it's not perfect, compared 
to a Fortran script”?

I can try to make a snapshot ... what kind do you need? Like a 
difference in the resulting vector, a piece of source code ... or maybe 
a graph plot illustrating the difference?

Since the default acc. to your script is not-a-knot ... returning to my 
original question, is it true I can alter this to 'natural' cubic spline 
by writing:

a(:,1) = _interp1_(log(f3),a1,log(f),'spline','natural');

?

... The documentation for interp1 spline isn't clear whether 'natural' 
is an option for the extrapolation. The documentation in Scilab ( 
https://help.scilab.org/docs/6.0.0/en_US/interp1.html ) only says:

The|extrapolation|parameter sets the evaluation rule for extrapolation, 
i.e for|xp(i)|not in [x1,xn] interval

"extrap"

    the extrapolation is performed by the defined method.
    yp=interp1(x,y,xp,method,"extrap")


... And the documentation provides no examples of adding this "extrap" 
parameter to the spline fitting (or any other fitting).


???

Anyway, the Fortran code uses 'natural' cubic spline and (as mentioned) 
DGTTRF and DGTTRS - are there any Scilab equivalent functions? ... I 
might wish to try to code (exactly) the same spline functionality, just 
to eliminate a potential error here.

DGTTRF computes an LU factorization of a real tridiagonal matrix A
* using elimination with partial pivoting and row interchanges.

DGTTRS solves one of the systems of equations
* A*X = B or A'*X = B,
* with a tridiagonal matrix A using the LU factorization computed
* by DGTTRF.

( Source: 
https://www.gfd-dennou.org/arch/ruby/products/ruby-lapack/doc/dgt.html )

P.S. Sorry for the not-so-fast response time, 1) I'm thinking a lot 
about it, 2) I shall try not to make too many false statements ... but 
it's a bit difficult because I don't understand why I get the 
differences that I observe - and I feel a bit like I'm searching for the 
reason (the solution) while being blind-folded.

Best regards,
Claus

On 11-09-2016 18:48, Rafael Guerra wrote:
>
> Hi,
>
> The code line using *interp1* and the spline method seems to perform 
> "not_a_knot" cubic spline interpolation, as demonstrated by the simple 
> test here below:
>
> x0= [0 1 2 3 4];
> y0= [0 -1 0 2 1];
> x= 0:0.05:4;
> y1= _interp1_(x0,y0,x,'spline');
> d= splin(x0, y0,"not_a_knot");
> e= splin(x0, y0,"natural");
> y2= interp(x, x0, y0, d);
> y3= interp(x, x0, y0, e);
> _clf_();
> _plot_(x,y1,'blue',x,y2,'--green',x,y3,'red',x0,y0,'Xblack')
>
> What do you mean by your result is “OK, but it's not perfect, compared 
> to a Fortran script”?
>
> Could you provide a snapshot?
>
> Regards,
>
> Rafae;
>
> *From:*Claus Futtrup
> *Sent:* Friday, September 09, 2016 8:34 PM
> *To:* International users mailing list for Scilab. 
> <users at lists.scilab.org>
> *Subject:* [Scilab-users] Cubic spline
>
> Hi there
>
> In Scilab I've used interp1 to calculate a cubic spline interpolation, 
> like this:
>
> a(:,1) = _interp1_(log(f3),a1,log(f),'spline');
>
> a = amplitude (magnitude). f3 is a frequency (27000 linear spaced 
> data), a1 is the original data, f is the resampled 1200 frequencies 
> (log-spaced), where I need the spline to interpolate some data for me.
>
> Above should work OK, but it's not perfect, compared to a Fortran 
> script. The fortran script calculates with its own cubic spline 
> routine, utilizing LAPACK (DGTTRF and DGTTRS) to solve for polynomial 
> coefficients.
>
> The question is - above code line with the interp1 spline, which kind 
> of spline is it?
>
> Digging into interp, I see multiple options. Digging into splin, I 
> also see multiple options.
>
> I looks like the interp1 is using "natural" spline - is this correct?
>
> It's strange because the splin help documentation doesn't recommend 
> this. It says: Don't use the natural type unless the underlying 
> function have zero second end points derivatives.
>
> This might be my problem.
>
> The Scilab help for interp1 doesn't give any examples, but does 
> mention I can add an "extrap" method. Could this be any of the 
> suggestions in the splin documentation. For example, could I write:
>
> a(:,1) = _interp1_(log(f3),a1,log(f),'spline','not-a-knot');
>
> ?
>
> P.S. Since not-a-knot is mentioned as the default for the splin 
> function, I think it should also be made the default for interp1 ... 
> just my two cents.
>
> /Claus
>
>
>
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20160920/d54527b7/attachment.htm>