【Matlab】三次样条插值实现

Myth丶恋晨 2022-11-20 03:56 442阅读 0赞

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

--------------------

### 文章目录 ###

*   *  题目：三次样条插值
     *  程序1：Gauss列主元消去法
     *  程序2：三次样条插值
     *  程序3：主函数及输入数据
     *  运行结果

## 题目：三次样条插值 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70_pic_center]

## 程序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：三次样条插值 ##

思路：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70_pic_center 1]

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

## 运行结果 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70_pic_left]  
![在这里插入图片描述][20201113194958634.png_pic_left]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70_pic_center]: /images/20221120/09afe0861b1e47e28d805b2f49efd027.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221120/b07887d80a9e4f7690be0942e61cc627.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70_pic_left]: /images/20221120/e148e9c12e1b447dba9122e050666360.png
[20201113194958634.png_pic_left]: https://img-blog.csdnimg.cn/20201113194958634.png#pic_left