三次插值函数的Matlab实现

科学计算课的上机作业，留下来供日后参考。

mSpline文件

clear,clc

​
%定义全局变量，方便在函数ms中调用
​
global x;
​
global y;
​
global M;
​
global n;
​
global df0;
​
global dfn;
​
%输入数据
​
x=[0:10];
​
y=[2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80];
​
df0=0.8;
​
dfn=0.2;
​
%点的个数
​
nt=size(x);
​
n=nt(2)-1;
​
%计算h,f
​
h=zeros(1);
​
f=zeros(1);
​
for k=1:n
​
   h(k)=x(k+1)-x(k);
   f(k)=(y(k+1)-y(k))./h(k);
​
end
​
%计算lambda,mu,e
​
lambda=zeros(1);
​
mu=zeros(1);
​
me=zeros(1);
​
for k=2:n
​
   lambda(k)=h(k)./(h(k-1)+h(k));
   mu(k)=h(k-1)./(h(k-1)+h(k));
   me(k)=3*(lambda(k)*f(k-1)+mu(k)*f(k));
​
end
​
%构造矩阵AM=B
​
A=zeros(n-1,n-1);
​
A(1,1)=2;
​
A(1,2)=mu(2);
​
for i=2:n-2
​
   A(i,i-1)=lambda(i+1);
   A(i,i)=2;
   A(i,i+1)=mu(i+1);
​
end
​
A(n-1,n-2)=lambda(n);
​
A(n-1,n-1)=2;
​
M=zeros(n-1,1);
​
B=zeros(n-1,1);
​
B(1)=me(2)-lambda(2)*df0;
​
for i=2:n-2
​
   B(i)=me(i+1);
​
end
​
B(n-1)= me(n) - mu(n)*dfn;
​
%求解M
​
M=A\B;
​
mvx=0:0.01:10;
​
%调用ms生成插值函数
​
   mvy=0;
​
for i=1:1001
​
   mvy(i)=ms(mvx(i));
​
end
​
plot(x,y,'o',mvx,mvy,'black');
​
legend('原始数据','三次样条插值')

ms.m文件

function z = ms(mx)

%MS 计算s
​
%声明全局变量
global M;
global x;
global y;
global n;
global df0;
global dfn;
​
%这是添加了初始条件后的M矩阵
mM=[df0;M;dfn];
​
%这里的k表示书中的k+1
z=0;
for k=1:n+1
   z=z+(y(k)*malpha(k,mx)+mM(k)*mbeta(k,mx));
end
​
end
​
%实现的alpha函数
function my=malpha(k,mx)
global x;
global y;
global n;
​
if k==1
   if x(k)<= mx & mx <= x(k+1)
       my=(1+2*(mx-x(k))/(x(k+1)-x(k)))*((mx-x(k+1))/(x(k)-x(k+1)))^2;
   else
       my=0;
   end
   return
end
​
if 2<= k & k<= n
   if x(k-1) <= mx & mx <= x(k)
       my=(1+2*(mx-x(k))/(x(k-1)-x(k)))*((mx-x(k-1))/(x(k)-x(k-1)))^2;
   else
       if x(k) < mx & mx <= x(k+1)
         my=(1+2*(mx-x(k))/(x(k+1)-x(k)))*((mx-x(k+1))/(x(k)-x(k+1)))^2; 
       else
           my=0;
           return
       end
   end
end
​
if k==n+1
   if x(1)<= mx & mx < x(k-1)
       my=0;
   else
       if x(k-1) <= mx & mx <= x(k)
           my=(1+2*(mx-x(k))/(x(k-1)-x(k)))*((mx-x(k-1))/(x(k)-x(k-1)))^2; 
       else
           my=0;
           return
       end
   end
end
end
​
function my=mbeta(k,mx)
global x;
global y;
global n;
​
if k==1
   if x(k)<= mx & mx <= x(k+1)
       my=(mx-x(k))*((mx-x(k+1))/(x(k)-x(k+1)))^2;
   else
       my=0;
   end
   
   return
end
​
if 2<= k & k<= n
   if x(k-1) <= mx & mx <= x(k)
       my=(mx-x(k))*((mx-x(k-1))/(x(k)-x(k-1)))^2;
   else
       if x(k) < mx & mx <= x(k+1)
         my=(mx-x(k))*((mx-x(k+1))/(x(k)-x(k+1)))^2; 
       else
           my=0;
           return
       end
   end
end
​
if k==n+1
   if x(1)<= mx & mx < x(k-1)
       my=0;
   else
       if x(k-1) <= mx & mx <= x(k)
           my=(mx-x(k))*((mx-x(k-1))/(x(k)-x(k-1)))^2; 
       else
           my=0;
           return
       end
   end
end
end