[数值分析拟合]Matlab三次样条插值拟合数据

           三次样条插值是一种运用极为广泛的工程插值算法，本文章编写的函数默认使用端点处的导数值代替给定的两端点的导数值使用三转角构造法进行插值（该函数也可传入端点导数数值进行分析），对数据进行方便而迅速的拟合（但是目前没有三弯矩构造法）

        一、三次样条插值的基本原理：

        首先，三次样条插值将曲线的每一段通过三次函数进行拟合。曲线会通过每一个数据点，并且在每个数据点处都有二阶连续条件（节点处的一阶导数值和二阶导数值均相同），当已知n+1个数据点时，可以设前n个节点的三次函数分别为  （i=1,2,...n）其中有4n个未知数，需要建立4n个方程来解决（下面是自己整理的总体思路，因为也是第一次写比较多的公式，如有错误和不妥之处请在评论区指正）

1) 首先有  (i=1,2,...n),(通过值条件)可建立n个方程，即

2) 由 (i=1,2,...n) (连续性条件)建立n个方程

3) 使用内节点处导数连续条件，即（i=1,2,...n-1）,由求导可得，

记,

则有

故可利用上式得出

,

这样就建立了n-1个节点处导数连续的方程

3)同样地，解出二阶导数，并使得到

其中设为该点处的二阶导数值，可建立n-1个相应的方程组

于是，这样对一条插值曲线的4n个未知数就产生了4n-2个约束，另外两个约束由边界条件确定

通过边界条件等4n个方程将Ai,Bi,Ci,Di用mi表示，建立方程组解出m，从而确定共4n个未知数的值，构造出每一段的三次曲线方程

        二、简介三转角构造法和三弯矩构造法

1)   三转角构造法：采用两端的导数值作为边界条件建立2个补充方程

2）三弯矩构造法：给出两端点的二阶导数值，作为补充的2个边界条件

其余的还有自然边界法等方法，这个方法会在下方的文章链接里面有详细解释

        三、详细的推导过程可以参考下面这篇文章：

****************************************

三次样条插值的原理和C语言实现

****************************************

另外这篇文章的一部分内容也是有一些比较难以推导的部分和误区，在这里补充如下：

在上述文章的【端点条件】中（端点条件即为使用三转角构造法），有推导过程

他的第二式b0=A是这样推导的：

1）利用前面的式子b.       (1)

可获得b0为上式值

2）所给出的边界条件中，这一张里面下标应当为bn=B（是最后一个节点的），通过这个推导下方的式子hn-1.....这里要用到导数连续条件     (2)

（有,）

 将（1）式中的i换成n-1代入（2）式（i=n）中,即可得到图中的第二个式子

        四、代码部分

下面附上完整的程序：

1）首先要在该函数的文件夹下创建一个函数len.m以便于后期取长度时调用

%返回一维数组长度(行数或列数其中大的一个，必须是一位数组)
function [length]=len(A)
 
    s=size(A);
    if s(2)==1
        length=s(1); %取行
    end
    if s(1)==1
        length=s(2);
    end
    if s(1)~=1 && s(2)~=1
        disp('必须是单行向量或单列向量')
        return
    end
            
end

2)Spline_3.m   (三次样条插值主函数)

%三次样条插值的matlab实现
%本例为三转角构造法，需要给出边界处的导数值
%本函数默认使用边界两点的斜率拟合最后一点的导数值
%一般默认调用方法为：[x2,y2] = Spline_3(x,y) ,  若直接使用 Spline_3(x,y)则绘图
%指定端点导数值调用方法： [x2,y2] = Spline_3(x,y,'dx',[1,-1])
function [x2,y2]=Spline_3(x,y,String,dx)  
%传入 x y 矩阵和初始参数,(默认使用三转角构造，若传入三个参数，则指定两段的导数值，
%四个参数指定两端的二阶导数值（此时第三个参数无效）（目前没有实现三弯矩构造）
n=len(x);
if n<=3
    disp('数组长度过短!')
    return
end
if n~=len(y)
    disp('两个数组长度必须相同！')
    return
end
 
if nargin==2
    dx=zeros(1,2);
    dx(1)=(y(2)-y(1))/(x(2)-x(1));
    dx(2)=(y(n)-y(n-1))/(x(n)-x(n-1));
    %指定末端点的导数值
else
    if nargin ==4
        if String == 'dx'
            %什么也不做
        end
    end
end
%========预分配内存空间===============================
h=zeros(n-1,1);  %步长
delta_y=zeros(n-1,1);   %增量
%解对应的系数
A=zeros(n-1,1);
B=zeros(n-1,1);
C=zeros(n-1,1);
D=zeros(n-1,1);
%G和f为需要解的n*n矩阵及其数值
G=zeros(n,n);  %是 n*n 矩阵,右端有n个二阶导数值
f=zeros(n,1);  % f是右边的矩阵
%=========================
 
for i=1:1:n-1
    h(i)=x(i+1)-x(i);
    delta_y(i)=y(i+1)-y(i); %每次y的增量
end
%   A是y(i)而B是导数
 
%if nargin==2     %这个留在后续可能要加三弯矩时使用
%采用三转角构造方式
    G(1,1)=2*h(1);
    G(1,2)=h(1);
    %n个数据共有n个导数值，这里取第n个来列方程（b(n)是倒数第一个）
    %参考资料的b(n-1)有误，应当是b(n)
    G(n,n-1)=h(n-1); 
    G(n,n)=2*h(n-1);
    a=dx(1);
    b=dx(2);
%=============构造矩阵=====================
    %填写G矩阵中其余的参数()
    for i=2:1:n-1  %从2行到n-1行的数据
        G(i,i-1)=h(i-1);
        G(i,i)=2*(h(i)+h(i-1));
        G(i,i+1)=h(i);
    end
    
    f(1)=6*( delta_y(1)/h(1) - a  ) ;    
    %在默认情况下，采用最前方的斜率代替导数值时，f(1)必然为0
    f(n)=6*(  b-delta_y(n-1)/h(n-1)  );
    %同理默认情况下  f（n-1）=0
    for i = 2:1:n-1
        f(i)=6*(delta_y(i)/h(i)-delta_y(i-1)/h(i));
    end
    %m为二阶导数值
    m=G\f;   
    %是G\f（实际上是inv(G)*f  ）解出m 的值（这个看错检查我好长时间QAQ）
    m=m';
    
%=====================================================
    %计算每段样条插值的系数：
    for i=1:1:n-1
        A(i)=y(i);
        B(i)=  (y(i+1)-y(i))/h(i) -h(i)*m(i)/2 - h(i)/6*(m(i+1)-m(i));
        C(i)=m(i)/2;
        D(i)=(m(i+1)-m(i))/(6*h(i));
    end
    
%*******注意:从这里设置插值后绘制函数的精度*****可以有不同的要求,
%这里设为0.001
    x_2=min(x):0.001:max(x);  %如果用x2会有问题
    x_2=x_2';   %行向量
    %绘制插值曲线
    y2=zeros(len(x_2),1);
    j=1;
    hold on
        
    for i=1:1:len(x_2)
        while x_2(i)>x(j+1)  %变换插值函数
            j=j+1;
        end
        y2(i)=A(j)   +B(j)*(x_2(i)-x(j))  +C(j)*(x_2(i)-x(j))^2  +  D(j)*(x_2(i)-x(j))^3;
    end
 
%end    %三转角构造方式
x2=x_2;
 
if nargout==0
    plot(x2,y2)
end
 
end

用法示例1:三次样条插值cos

%这个是一般的最简便调用格式
 
t=0:1:10;
y=cos(t);
Spline_3(t,y)  %绘制三次样条插值曲线
 
%下面是可以返回需要的t2和y2值
 
t=0:1:10;
y=cos(t);
[t2,y2]=Spline_3(t,y)
plot(t2,y2)

结果：

 用法2：传入端点处的一阶导数值

t=0:1:10;
y=cos(t);
Spline_3(t,y,'dx',[10,-10]);  %这样可以指定传入首尾端点的导数值，按照首尾的导数去插值
 
%最后一种就是带返回值，懂的都懂AWA

结果（修改了cosx插值的端点导数）:

 (这当然还是满足1,2,等过相应的值)

 五、实际应用

主要用于处理和模拟实验数据是极有效的模拟方式，我会最近发布一篇文章《三次样条插值整理实验数据》来对此做一步应用，有兴趣者可以参考

[matlab实践应用]matlab实现读取xls表格并三次样条插值拟合压杆稳定实验数据