How to get an equation of spline interpolation for a set of data X and Y?

Hi , i have a set of data which is x and y and lets say for example

x=[1 2 3 4 5 6 7 8 9 10 ]

y=[100 200 300 400 500 600 700 800 900 1000]

How to get the equation of spline interpolation for this data, as i used the command

OT = spline(x,y);

[KSSV]

Check OT, it has piece wise polynomials. To get the values use ppval.

[zaid jamil]

ok , OT is a struct that have many variables wich contains coeffecants and breaks and many other things... how i can use it to make an equation so that by using it in anther program it can work perfectly

Answer 1

[John D'Errico]

Congratulations! You win the door prize as the person to ask this question the millionth time!

Seriously, this is a question I get asked so many times I have lost count. And there really is no easy formula you can write down, because the answer is a bit of a mess. I do offer a code to provide what you want. But reember that a spline is a segmented thing, composed of segments of cubic polynomials, defined in a piecewise sense.

x = linspace(0,2*pi,9)
x =
            0       0.7854       1.5708       2.3562       3.1416        3.927       4.7124       5.4978       6.2832
y = sin(x);
spl = spline(x,y);

As I said, spl is a cubic spline. In my SLM toolbox, I provide a code that will extract the polynomial segments.

funs = slmpar(spl,'symabs')
funs =
  2×8 cell array
  Columns 1 through 7
    {1×2 double}    {1×2 double}    {1×2 double}    {1×2 double}    {1×2 double}    {1×2 double}    {1×2 double}
    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }
  Column 8
    {1×2 double}
    {1×1 sym   }

So in the half open interval [0,pi/4), the cubic polynomial is given as

>> funs{1,1}
ans =
            0       0.7854
>> vpa(funs{2,1},10)
ans =
- 0.08493845396*x^3 - 0.1356173501*x^2 + 1.059224243*x

In the half open interval [pi, pi + p/4) there is a different cubic.

>> funs{1,5}
ans =
       3.1416        3.927
>> vpa(funs{2,5},10)
ans =
0.1568856928*x^3 - 1.47861282*x^2 + 3.648107873*x - 1.731986498

With 9 data points, we will have 8 cubic segments. So the final segment has this form:

>> funs{1,8}
ans =
       5.4978       6.2832
>> vpa(funs{2,8},10)
ans =
- 0.08493845396*x^3 + 1.736669488*x^2 - 10.70470091*x + 19.76765782

I could have extracted the polynomial segments in a variety of forms.

Answer 2

[Fayyaz Ahmad]

clc;

clear all;

close all;

format short g

x = linspace(0,2*pi,9);

y = sin(x);

spl = spline(x,y);

bk = spl.breaks;

cf = spl.coefs;

[m,n] = size(cf);

syms z;

v = z.^((n-1:-1:0)');

eq = cf*v

------------------------------------------------------

output

----------------------------------------------------------

eq =

(4770321904107807*z)/4503599627370496 - (610766247492719*z^2)/4503599627370496 - (6120460633987311*z^3)/72057594037927936

(6206087117424961*z)/9007199254740992 + 2^(1/2)/2 - (6048313895780885*z^2)/18014398509481984 - (6120460633987189*z^3)/72057594037927936

(635448956098857*z^3)/9007199254740992 - (2413390700397705*z^2)/4503599627370496 + (159898229681981*z)/36028797018963968 + 1

(1413100694996111*z^3)/9007199254740992 - (3329540071636805*z^2)/9007199254740992 - (6365985347106939*z)/9007199254740992 + 2^(1/2)/2

(2826201389992229*z^3)/18014398509481984 - (9*z^2)/9007199254740992 - (4490500002164345*z)/4503599627370496 + 4967757600021511/40564819207303340847894502572032

(1270897912197719*z^3)/18014398509481984 + (6659080143273605*z^2)/18014398509481984 - (6365985347106939*z)/9007199254740992 - 2^(1/2)/2

- (765057579248403*z^3)/9007199254740992 + (2413390700397707*z^2)/4503599627370496 + (159898229681977*z)/36028797018963968 - 1

- (6120460633987189*z^3)/72057594037927936 + (6048313895780881*z^2)/18014398509481984 + (3103043558712477*z)/4503599627370496 - 2^(1/2)/2

[John D'Errico]

@Fayyaz Ahmad

Be careful, as those polynomials cannot be used to evaluate the spline directly. (TRY IT!) You have made a mistake, if you think they can.

clc;
clear all;
close all;
format short g 
x = linspace(0,2*pi,9);
y = sin(x);
spl = spline(x,y);
bk = spl.breaks;
cf = spl.coefs;
[m,n] = size(cf);
syms z;
v = z.^((n-1:-1:0)');
eq = cf*v;

For example...

eq(2)

ans = 

bk

bk = 1×9

0 0.7854 1.5708 2.3562 3.1416 3.9270 4.7124 5.4978 6.2832

So, on the interval pi/4 to pi/2, you would imply that this polynomial should approximate sin(x), correct? Try it!

P2 = matlabFunction(eq(2));
fplot(@sin,[pi/4,pi/2])
hold on
fplot(P2,[pi/4,pi/2])
legend('Sin(x)','Polynomial 2')
xlabel 'X'
grid on

Gosh. It does not look correct to me. Not even close.

Your mistake lies in a misunderstanding on your part of what spline produces, and probably why there is a problem at all. But if you will post this as an answer, you need to understand why it fails to produce something reasonable.

[Fayyaz Ahmad]

@ John D'Errico thanks for comments, you are right. Actually the coefficients of polynomial generated by spline command are on the interval [0,bk(i+1)-bk(i)]. We have corrected the code.

-----------------------

clc ;

clear all ;

close all ;

format short g

x = linspace(pi,2*pi,9);

y = sin(x);

spl = spline(x,y);

bk = spl.breaks ;

cf = spl.coefs ;

[m,n] = size(cf );

syms z ;

v = z.^((n-1:-1:0 )');

eq = cf*v;

for i=1:m

eq(i) = simplify(subs(eq(i),z,z-bk(i)));

end

figure(1)

hold on

q = length(bk);

for i=1:q-1

s1 = linspace(bk(i),bk(i+1),10);

p = matlabFunction(eq(i));

fplot(@sin,[bk(i),bk(i+1)])

plot(s1,p(s1),'.-r')

pause

end

eq

--------------------------------------

output is

--------------------------------------

eq =

(2260863087144877*pi)/2251799813685248 - (2260863087144877*z)/2251799813685248 + (661015214105797*(z - pi)^2)/36028797018963968 + (651968472240021*(z - pi)^3)/4503599627370496 + 4967757600021511/40564819207303340847894502572032

(9349213407344919*pi)/9007199254740992 - (1038801489704991*z)/1125899906842624 + (1701418325597523*(z - (9*pi)/8)^2)/9007199254740992 + (2607873888959953*(z - (9*pi)/8)^3)/18014398509481984 - 6893811853601121/18014398509481984

(15927176730973595*pi)/18014398509481984 - (3185435346194719*z)/4503599627370496 - 2^(1/2)/2 + (6475165695337053*(z - (5*pi)/4)^2)/18014398509481984 + (6610892248307199*(z - (5*pi)/4)^3)/72057594037927936

(9475875933393265*pi)/18014398509481984 - (861443266672115*z)/2251799813685248 + (4211117090838325*(z - (11*pi)/8)^2)/9007199254740992 + (4785409777295905*(z - (11*pi)/8)^3)/144115188075855872 - 2080391759176529/2251799813685248

(27*pi)/144115188075855872 - (9*z)/72057594037927936 + (570433996317983*(z - (3*pi)/2)^2)/1125899906842624 - (2392704888647973*(z - (3*pi)/2)^3)/72057594037927936 - 1

(3445773066688457*z)/9007199254740992 - (44795049866949941*pi)/72057594037927936 + (8422234181676639*(z - (13*pi)/8)^2)/18014398509481984 - (6610892248307199*(z - (13*pi)/8)^3)/72057594037927936 - 8321567036706117/9007199254740992

(6370870692389435*z)/9007199254740992 - (44596094846726045*pi)/36028797018963968 - 2^(1/2)/2 + (809395711917129*(z - (7*pi)/4)^2)/2251799813685248 - (5215747777920001*(z - (7*pi)/4)^3)/36028797018963968

(4155205958819961*z)/4503599627370496 - (62328089382299415*pi)/36028797018963968 + (6805673302390011*(z - (15*pi)/8)^2)/36028797018963968 - (5215747777919999*(z - (15*pi)/8)^3)/36028797018963968 - 3446905926800567/9007199254740992