Natural Cubic Spline Approximation
1 次查看(过去 30 天)
显示 更早的评论
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.
0 个评论
采纳的回答
Walter Roberson
2012-4-8
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Splines 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!