Natural Cubic Spline Approximation

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.

Answer 1

[Walter Roberson]

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.

[Jacob]

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??