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!