Natural Cubic Spline Approximation
Show older comments
Hi All, I am doing a natural cubic spline problem and at some point I have to divide the second derivative function with the original vector but still I can not go through. Obviously, the second derivative of the function is 2 values less than the original function that's why I get the error that,"Index exceeds matrix dimensions", please I need help on this, anyone with ideas:
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=2;
d2x = [d2x 0 0];
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
Here is what the program tells me:
??? 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 54
c1=(d2x(i)/6)/(x(i+1)-x(i));
Answers (0)
Categories
Find more on Common Operations 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!