数值分析（二续）三次样条插值二类边界完整matlab代码

原创

已于 2024-03-18 14:06:37 修改 · 3.4w 阅读

156

658

标签

#matlab #算法

于 2020-07-01 19:03:20 首次发布

本文详细介绍了三次样条插值法的两种边界条件下的MATLAB代码实现，包括自然边界和第一边界条件，通过具体实例展示了如何使用这些代码进行插值计算并绘制插值函数图像。

前言

根据上篇文章链接： 数值分析（二) 三次样条插值法matlab程序 其中只提及到了自然边界条件情况下的matlab代码，本篇文章将来填补上篇文章的其他内容给出完整的三次样条插值函数matlab代码。
注意：上篇文章所有计算原理都已讲过，本篇文章将不会重复论述上篇已有的东西，这里直接给出三种边界代码即每种边界对应一种实例题。如果还不了解原理的读者可以看上篇文章。

重新更新于（2022.5.3）新增：根据大部分同学们的询问，各段表达式显示以及各段表达式系数，本人修改代码增加了输出各段系数矩阵，以及各段表达式输出。对内容非常熟悉的可以直接看新增板块。

1. 第一边界

代码第一边界matlab代码如下：

function [D,h,A,g,M]=three1(X,Y,y0,yn)
%        自然边界条件的三次样条函数(第一种边界条件)
%        此函数为M值求值函数
%        D,h,A,g,M输出量分别为系数矩阵D，插值宽度h，差商表A，g值,M值 
         n=length(X); 
         A=zeros(n,n);A(:,1)=Y';D=zeros(n,n);g=zeros(n,1);
         for  j=2:n
            for i=j:n
                A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
            end
         end
         
         for i=1:n-1
             h(i)=X(i+1)-X(i);
         end
         for i=1:n
             D(i,i)=2; 
             D(1,2)=1;
             D(n,n-1)=1;
             if (i==1)
                 g(i,1)=6/h(i)*(A(2,2)-y0); 
             elseif (i==n) 
                     g(i,1)=6/h(i-1)*(yn-A(i,2));
             else 
                 g(i,1)=(6/(h(i-1)+h(i)))*(A(i+1,2)-A(i,2));
             end
           
         end  
         for i=1:n-2
             u(i)=h(i)/(h(i)+h(i+1));
             n(i)=1-u(i);  
             D(i+1,i+2)=n(i);
             D(i+1,i)=u(i);             %改到这里
         end
         M=D\g;
         %M=[0;M;0];         
end

function s=threesimple1(X,Y,x,y0,yn)
%        三次样条插值函数第一类型代码 
%        s函数表示三次样条插值函数插值点对应的函数值
%        根据三次样条参数函数求出的D,h,A,g,M
%        x表示求解插值点函数点，X为已知插值点        
         [D,h,A,g,M]=three1(X,Y,y0,yn)
         n=length(X); m=length(x);    
         for t=1:m
            for i=1:n-1
               if (x(t)<=X(i+1))&&(x(t)>=X(i))
                  p1=M(i,1)*(X(i+1)-x(t))^3/(6*h(i));
                  p2=M(i+1,1)*(x(t)-X(i))^3/(6*h(i));
                  p3=