How to get an equation of spline interpolation for a set of data X and Y?
8 views (last 30 days)
Show older comments
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);
2 Comments
Answers (2)
John D'Errico
on 18 Jun 2020
Edited: John D'Errico
on 18 Jun 2020
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.
0 Comments
Fayyaz Ahmad
on 17 Oct 2021
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
2 Comments
John D'Errico
on 20 Oct 2021
Edited: John D'Errico
on 20 Oct 2021
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)
bk
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
on 20 Oct 2021
@ 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

See Also
Categories
Find more on Spline Postprocessing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!