三次样条插值算法

   在[a,b]上分为n段，共n+1个点。插入3次多项式，并使其二阶导数连续的方法称为三次样条插值算法。

思路：

1.二阶导数为线性函数。

2.插值点的函数值已知、一阶导数、二阶导数连续。

3.加上边界条件即可求解。

边界条件

1.夹持条件：已知起点和终点的速度。

2.自然边界条件：已知奇点和终点的加速度。

3.周期性条件：假设f为以b-a为周期的周期函数        

 

插值方法：

假设：第i个点的二阶导数为Mi，则在[xi-1,xi]上，其二阶导数为线性函数

则在[xi-1,xi]上进行线性插值

令 

进行二次积分得到原函数 

 

线性插值有：

 得到原函数：

对各段函数求导 

 

由连续性条件 

两边同乘 

则可写成

添加边界条件1 

 

可写成

可写成

写成矩阵形式

结合以下五个式子即可求解

 

实例matlab介绍 

%三次样条插值
x=[0 1 4 5];
y=[0 -2 -8 -4];
dy1=5/2;
dy4=19/4;
n=4;
 
%计算h
for i=2:n
    h(i)=x(i)-x(i-1);
end
 
%计算u，v，g
for i=2:n-1
    u(i)=h(i)/(h(i)+h(i+1));
    v(i)=1-u(i);
    g(i)=6/(h(i)+h(i+1))*((y(i+1)-y(i))/h(i+1)-(y(i)-y(i-1))/h(i));%明天修改
end
 
g(1)=6/h(2)*((y(2)-y(1))/h(2)-dy1);
g(4)=6/h(4)*(dy4-(y(4)-y(3))/h(4));
 
 
%系数矩阵
A=zeros(4,4);
for i=1:n
    if i==1
        A(i,1)=2;
        A(i,2)=1;
    else if i==n
            A(i,n-1)=1;
            A(i,n)=2;
        else 
            A(i,i)=2;
            A(i,i-1)=u(i);
            A(i,i+1)=v(i);
        end
    end
end
 
m=inv(A)*g';
syms t
for i=2:4
    a(i)=m(i-1)/(6*h(i));
    b(i)=m(i)/(6*h(i));
    c(i)=(y(i-1)-m(i-1)*h(i)^2/6)/h(i);
    d(i)=(y(i)-m(i)*h(i)^2/6)/h(i);
end
t1=0:0.01:1;
t2=1:0.01:4;
t3=4:0.01:5;
plot(t1,fthreesample(m,x,y,h,2,t1))
hold on
plot(t2,fthreesample(m,x,y,h,3,t2))
plot(t3,fthreesample(m,x,y,h,4,t3))

function f=fthreesample(m,x,y,h,k,t)
% 
for i=2:4
    a(i)=m(i-1)/(6*h(i));
    b(i)=m(i)/(6*h(i));
    c(i)=(y(i-1)-m(i-1)*h(i)^2/6)/h(i);
    d(i)=(y(i)-m(i)*h(i)^2/6)/h(i);
end
 
f=a(k)*(x(k)-t).^3+b(k)*(t-x(k-1)).^3+c(k)*(x(k)-t)+d(k)*(t-x(k-1))

插值结果如下：