【Matlab】三次样条插值实现

敲了一下午终于把代码敲完了，留念。

---

题目：三次样条插值

程序1：Gauss列主元消去法

function [X] = Gauss_elimination(A, Y)
    %% 参数初始化
    matrix = [A, Y];
    [n,m] = size(matrix);         % 矩阵大小
    l = zeros(n,n);          % 比例系数矩阵
    X = zeros(1,n);          % 方程组解
  
    
    %% 列主元消元
    for k=1:n-1     % 1~~n-1列,1~~n-1行
        [~,m] = max(matrix(:,k));
        change = matrix(k,:);
        matrix(k,:)=matrix(m,:);
        matrix(m,:)=change;
        if matrix(k,k)==0 
            flag=0;break;
        end 
        for i=k+1:n
           l(i,k)=matrix(i,k)/matrix(k,k); 
           for j=k:n+1
              matrix(i,j)=matrix(i,j)-l(i,k)*matrix(k,j);
           end
        end
    end
    
    % 回代
    X(1,n)=matrix(n,n+1)/matrix(n,n);
    for k=n-1:-1:1
        sum=matrix(k,end);
        for j=k+1:n            
            sum=sum-matrix(k,j)*X(1,j);
        end
        X(1,k)=sum/matrix(k,k);
    end     
    matrix;
end

程序2：三次样条插值

思路：

function yy = Interpolation_Spline(x, y, dy0, dyn)
    %% 计算三阶差商表（不含导数）
    [~,length] = size(y);
    max_order = length - 2;
    diffQuo = zeros(length, 4);
    diffQuo(:,1) = x';
    diffQuo(:,2) = y';
    if(max_order < 3)
        printf("没有足够的数据计算三阶差商")
    else
        for i=3:4 % 第三列到第四列
            for j=1:length -i+2  % 行
                diffQuo(j,i) = (diffQuo(j,i-1) - diffQuo(j+1, i-1)) / (diffQuo(j,1) - diffQuo(j+i-2, 1));
            end
        end
    end
    
    % 改写(将导数加入差商表）
    dd0 = (diffQuo(1,3) - dy0) / (diffQuo(2,1) - diffQuo(1,1));
    ddn = (diffQuo(length-1,3) - dyn) / (diffQuo(length-1,1) - diffQuo(length,1));
    diffQuo1 = [x(1),y(1),dy0,dd0];
    diffQuon = [x(length),y(length),0,0];
    diffQuo = [diffQuo1;diffQuo;diffQuon];
    diffQuo(length+1,3) = dyn;
    diffQuo(length,4) = ddn;

    %% 求三弯矩方程组
    H = zeros(1,length-1);
    mu = zeros(1,length-2);
    lamda = zeros(1,length-2);
    matrix = zeros(length,length);
    matrix(1,2) = 1;
    matrix(length,length-1) =  1;
    for i=1:length
        matrix(i,i) = 2;
    end
    for i=1:length-1
        H(1,i) = x(i+1) - x(i);
    end
    for i=1:length-2
        mu(1,i) = H(1,i) / (H(1,i) + H(1,i+1));
        lamda(1,i) = 1 -mu(1,i);
        matrix(1+i,i) = mu(1,i);
        matrix(1+i,i+2) = lamda(1,i);
    end
    
    D = diffQuo(1:length,4)';
    for i=2:length-1
        D(1,i) = diffQuo(i-1,4);
    end
    
    %% 列主元消去法求解
    M = Gauss_elimination(matrix, D');
    
    %% 计算结果
    output = zeros(1,10);
    xx = x;
    yy = y;
    digits(5)  % 控制显示及计算精度
    for i=1:length-1
        syms x
        l1(x) = (diffQuo(i+1,3) - (1/3*M(i) + 1/6 * M(i+1)) * H(i)) * (x - xx(i));
        l2(x) = 1 / 2 * M(i) * (x - xx(i)) ^ 2;
        l3(x) = 1 / 6/ H(i) * (M(i+1) - M(i)) * (x - xx(i)) ^ 3;
        fprintf("------------------------------------------------------------------------------------------\n")
        Sx(x) = yy(i) + l1 + l2  +l3;
        Sx(x) = vpa(Sx);    % 控制显示及计算精度,与digits配合使用
        fprintf('S(x)=%s \t x∈(%d,%d)\n',char(Sx),i-1,i);
    end
    fprintf("\n\n\n")
    
    % 计算0.5~9.5的数值
    for i=1:length-1
        syms x
        l1(x) = (diffQuo(i+1,3) - (1/3*M(i) + 1/6 * M(i+1)) * H(i)) * (x - xx(i));
        l2(x) = 1 / 2 * M(i) * (x - xx(i)) ^ 2;
        l3(x) = 1 / 6/ H(i) * (M(i+1) - M(i)) * (x - xx(i)) ^ 3;
        Sx(x) = yy(i) + l1 + l2  +l3;

        Sx(x) = vpa(Sx);    % 控制显示及计算精度,与digits配合使用
        x_i = i - 0.5;
        y_i = Sx(x_i);
        output(i) = y_i;
        fprintf('%d\t %f\t %f\n',i,x_i,y_i); 
    end
    output;
end

程序3：主函数及输入数据

clc;clear all;
%% 给定条件
x = 0:1:10;
y = [2.51, 3.30, 4.04, 4.70, 5.22, 5.54, 5.78, 5.40, 5.57, 5.70, 5.80];
dy0 = 0.8;
dyn = 0.2;

%% 求三次样条插值
Interpolation_Spline(x, y, dy0, dyn);

运行结果