目录

• 前言
• 1. 第一边界
• 2. 第二边界
• 3. 实例分析
• 1. 第一小问
• 2. 第二小问

4. 总结

5. 补充

前言

  根据上篇文章链接：  数值分析（二) 三次样条插值法matlab程序 其中只提及到了自然边界条件情况下的matlab代码，本篇文章将来填补上篇文章的其他内容给出完整的三次样条插值函数matlab代码。
  

注意

：上篇文章所有计算原理都已讲过，本篇文章将不会重复论述上篇已有的东西，这里直接给出三种边界代码即每种边界对应一种实例题。如果还不了解原理的读者可以看上篇文章。

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=(A(i,1)-M(i,1)/6*(h(i))^2)*(X(i+1)-x(t))/h(i);
p4=(A(i+1,1)-M(i+1,1)/6*(h(i))^2)*(x(t)-X(i))/h(i);
s(t)=p1+p2+p3+p4;
break;
else
s(t)=0;
end
end
end
end

2. 第二边界

代码

此次编写的第二边界matlab代码其中包含了自然边界，程序如下：

function [D,h,A,g,M]=three2(X,Y,y0,yn)
%        第二边界条件的三次样条函数(包含自然边界条件)
%        y0,yn表示的是S''(x0)=f''(x0)=y0，S''(xn)=f''(xn)=yn,自然边界即条件值为0
%        此函数为M值求值函数
%        D,h,A,g,M输出量分别为系数矩阵D，插值宽度h，差商表A，g值,M值
n=length(X);
A=zeros(n,n);A(:,1)=Y';D=zeros(n-2,n-2);g=zeros(n-2,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-2
D(i,i)=2;
if (i==1)
g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2))-h(i)/(h(i)+h(i+1))*y0;
elseif (i==(n-2))
g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2))-(1-h(i)/(h(i)+h(i+1)))*yn;
else
g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2));
end
end
for i=2:n-2
u(i)=h(i)/(h(i)+h(i+1));
n(i-1)=h(i)/(h(i-1)+h(i));
D(i-1,i)=n(i-1);
D(i,i-1)=u(i);
end
M=D\g;
M=[y0;M;yn];
end

function s=threesimple2(X,Y,x,y0,yn)
%        第二边界条件函数
%        s函数表示三次样条插值函数插值点对应的函数值
%        根据三次样条参数函数求出的D,h,A,g,M
%        x表示求解插值点函数点，X为已知插值点
[D,h,A,g,M]=three2(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=(A(i,1)-M(i,1)/6*(h(i))^2)*(X(i+1)-x(t))/h(i);
p4=(A(i+1,1)-M(i+1,1)/6*(h(i))^2)*(x(t)-X(i))/h(i);
s(t)=p1+p2+p3+p4;
break;
else
s(t)=0;
end
end
end
end

3. 实例分析

  题目选自：《数值分析第五版李庆扬》书P49页第20题 题目如下

1. 第一小问

  属于第一边界值问题，将数值写入，调用上文先写好的matlab代码：

x=[0.25 0.3 0.39 0.45 0.53];
y=[0.5 0.5477 0.6245 0.6708 0.7280];
y0=1.000 ;          %  S'(x0)=f'(x0)=y0
yn=0.6868;          %  S'(xn)=f'(xn)=yn
x0=0.25:0.01:0.53;
s=threesimple1(x,y,x0,y0,yn)
plot(x0,s)        %绘制第一边界条件插值函数图像
hold on
grid on
plot(x,y,'o')
axis([0.2 0.55 0.4 0.75])
xlabel('自变量 X'), ylabel('因变量 Y')
title('插值点与三次样条函数')
legend('三次样条插值点坐标','插值点')

运行后的结果如下：

2. 第二小问

  属于第二边界值（自然边界）问题，将数值写入，调用上文先写好的matlab代码：

x=[0.25 0.3 0.39 0.45 0.53];
y=[0.5 0.5477 0.6245 0.6708 0.7280];
y0=0;          %  S''(x0)=f''(x0)=y0
yn=0;          %  S''(xn)=f''(xn)=yn
x0=0.25:0.01:0.53;
s=threesimple2(x,y,x0,y0,yn)
plot(x0,s)        %绘制第二边界条件插值函数图像
hold on
grid on
plot(x,y,'o')
axis([0.2 0.55 0.4 0.75])
xlabel('自变量 X'), ylabel('因变量 Y')
title('插值点与三次样条函数')
legend('三次样条插值点坐标','插值点')

运行后的结果如下：

4. 总结

  看上去代码大同小异，实则里面还是改变了几个逻辑点，需要读者慢慢去理解，希望这些代码配上实例可以帮助读者更好的理解公式。其实一些想要写C/C++的读者也可以认真的看看代码，主要是逻辑都是一样的，只是某些地方写法不同。感谢读者的耐心阅读，如果觉得写的好的话能否给编者点个赞非常感谢。

5. 补充

  本人又对代码经行改进优化及封装，完整的资料点击链接: 优化及完整代码