Natural Cubic Spline Approximation

3 Ansichten (letzte 30 Tage)
Jacob
Jacob am 8 Apr. 2012
Hi All, am writing a code to approximate the function using natural Cubic Spline. I am having a problem with the matrix dimensions. I am wondering if one may show me where I am doing wrong. Any help appreciated. Below is the code I have written and the response I get after running it:
clear all
clc
format short e
syms x r
fx = sin(exp(x)-2);
d2x = diff(diff(fx));
x = linspace(0,2,11);
y = subs(fx,x);
xu = 2.0;
yu = subs(fx,xu);
n=10;
%%Call Tridiagonal function
f(1)=2*(x(3)-x(1));
g(1)=(x(3)-x(2));
r(1)=6/(x(3)-x(2))*(y(3)-y(2));
r(1)=r(1)+6/(x(2)-x(1))*(y(1)-y(2));
for i=2:n-2
e(i)=(x(i)-x(i-1));
f(i)=2*(x(i+1)-x(i-1));
g(i)=(x(i+1)-x(i));
r(i)=6/(x(i+1)-x(i))*(y(i+1)-y(i));
r(i)=r(i)+6/(x(i)-x(i-1))*(y(i-1)-y(i));
end
e(n-1)=(x(n-1)-x(n-2));
f(n-1)=2*(x(n)-x(n-2));
r(n-1)=6/(x(n)-x(n-1))*(y(n)-y(n-1));
r(n-1)=r(n-1)+6/(x(n-1)-x(n-2))*(y(n-2)-y(n-1));
%%Call Decomposition function
for k=2:n-1
e(k)=e(k)/f(k-1);
f(k)=f(k)-e(k)*g(k-1);
end
%%Call Forward Substitution
for k=2:n-1
r(k)=r(k)-e(k)*r(k-1);
end
%%Call Back Substitution
x(n-1)=r(n-1)/f(n-1);
for k=n-1:1:-1
x(k)=(r(k)-g(k)*x(k+1))/f(k);
end
%%Call Interpolation function
flag = 0;
i=1;
while(1)
if xu>=x(i)&&xu<=x(i+1)
c1=(d2x(i)/6)/(x(i+1)-x(i));
c2=d2x(i+1)/6/(x(i+1)-x(i));
c3=y(i)/(x(i+1)-x(i))-d2x(i)*(x(i+1)-x(i))/6;
c4=y(i)/(x(i+1)-x(i))-d2x(i+1)*(x(i+1)-x(i))/6;
t1=c1*(x(i+1)-xu)^3;
t2=c2*(xu-x(i))^3;
t3=c3*(x(i+1)-xu);
t4=c4*(xu-x(i));
yu=t1+t2+t3+t4;
t1=-3*c1*(x(i+1)-xu)^2;
t2=3*c2*(xu-x(i))^2;
t3=-c3;
t4=c4;
dy=t1+t2+t3+t4;
t1=6*c1*(x(i+1)-xu);
t2=6*c2*(xu-x(i));
d2y=t1+t2;
flag=1;
else
i=i+1;
end
if i==n+1||flag==1,break
end
end
if flag==0
disp('Outside Range')
end
And here is the response which I get:
??? Error using ==> mupadmex
Error in MuPAD command: Index exceeds matrix dimensions.
Error in ==> sym.sym>sym.subsref at 1366
B = mupadmex('mllib::subsref',A.s,inds{:});
Error in ==> qn2 at 53
c1=(d2x(i)./6)./(x(i+1)-x(i));
Thanks for time and help.

Melden Sie sich an, um diese Frage zu beantworten.

Akzeptierte Antwort

Walter Roberson
Walter Roberson am 8 Apr. 2012
diff(x) is one element shorter than x. diff(diff(x)) would be two elements shorter than x. So d2x is two elements shorter than x, and you would have a problem if your "i" ever reaches length(x)-1 as d2x(i) would then be trying to access the element just past the end of d2x whereas x(i+1) and x(i) would both be fine at that "i" value.
  1 Kommentar
Jacob
Jacob am 8 Apr. 2012
Okay, I see your point. I know this is kind of strange but is there a way of extending the size of d2x like adding zeros so I may be able to do the division by making the d2x have the same size like x??

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Splines finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by