MATLAB: Natural Cubic Spline Approximation

natural cubic spline

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.

Best Answer

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.
Related Question