Natural Cubic Spline Approximation

1 次查看（过去 30 天）

Jacob 2012-4-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/35098-natural-cubic-spline-approximation

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));