三次样条插值及三弯矩法完整（Matlab实现）

目录

1、原理

2、案例

3、代码实现 

4、结果 

---

1、原理

     

       

       

2、案例

3、代码实现 

clear clc
x=[-1 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.5];
y=[-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31];
%%方法一
%for i=1:8
%    h(i)=x(i+1)-x(i);  %步长
%    fyjcs(i)=(y(i+1)-y(i))/h(i);  %一阶差商
%end
%for i=1:7
%     fejcs(i)=(fyjcs(i+1)-fyjcs(i))/(h(i+1)+h(i));  %二阶差商
%     % 二阶差商中“(x(i+2)-x(i))”出现8超过迭代i最大值7，故用(h(i+1)+h(i))代替
%end
%for i=2:8
%    u(8)=1;  % 边界条件代入
%    v(1)=1;  % 边界条件代入
%    dy1=5;  % 边界条件代入
%    dy9=29.16;  % 边界条件代入
%    d(1)=6/h(1)*((y(2)-y(1))/h(1)-dy1);  % 边界条件代入
%    d(9)=6/h(8)*(dy9-(y(9)-y(8))/h(8));  % 边界条件代入
%    u(i-1)=h(i-1)/(h(i-1)+h(i));  % 公式
%    v(i)=1-u(i-1);  % 公式
%    d(i)=6*fejcs(i-1);  % 公式
%end
%%三弯矩方程组：AM=d
%A=2*eye(9,9)+diag(v,1)+diag(u,-1)
%d1=d'
%AN=inv(A);
%M=AN*d1
%M2=M(2,1);
%M3=M(3,1);
%M7=M(7,1);
%M8=M(8,1);
%syms z c
%S2(z)=((((x(3)-z)^3)/(6*h(3)))*M2)+(((z-x(2))^3)/(6*h(3))*M3)+((y(2)-(M2*(((h(3))^2)/6)))*(x(3)-z)/h(3))+((y(3)-(M3*((h(3))^2/6)))*(z-x(2))/h(3));
%z=-0.02;
%fzl=vpa(S2(z))
%S7(c)=((((x(8)-c)^3)/(6*h(8)))*M7)+(((c-x(7))^3)/(6*h(8))*M8)+((y(7)-(M7*(((h(8))^2)/6)))*(x(8)-c)/h(8))+((y(8)-(M8*((h(8))^2/6)))*(c-x(7))/h(8));
%c=2.56;
%fcl=vpa(S7(c))
 
%%方法二

df1=5;df2=29.16;%边界条件
n=length(x);%计算维度
h=zeros(1,n);%初始化赋值
dy=zeros(1,n);%初始化赋值
u=zeros(1,n);%初始化赋值
lamda=zeros(1,n);%初始化赋值
d=zeros(1,n);%初始化赋值
l=zeros(1,n);%初始化赋值
r=zeros(1,n);%初始化赋值
z=zeros(1,n);%初始化赋值
M=zeros(1,n);%初始化赋值
for i=2:n
h(i)=x(i)-x(i-1);  %步长
end
for i=1:n-1
dy(i)=(y(i)-y(i+1))/(x(i)-x(i+1));  %一阶差商
end
for i=2:n-1
u(i)=h(i)/(h(i)+h(i+1));%求u
lamda(i)=1-u(i);%求lamda
d(i)=6*(dy(i)-dy(i-1))/(h(i)+h(i+1));%求d
end
d(1)=6*(dy(1)-df1)/h(2);
d(n)=6*(df2-dy(8))/h(n);
u(n)=1;  %u的最后一个元素为1
lamda(1)=1;  %lamda第1个元素为1
r(1)=2;
for i=2:n    %追赶法求解M
l(i)=u(i)/r(i-1);
r(i)=2-l(i)*lamda(i-1);
end
z(1)=d(1);
for i=2:n
z(i)=d(i)-l(i)*z(i-1);
end
M(n)=z(n)/r(n);
for i=n-1:-1:1
M(i)=(z(i)-lamda(i)*M(i+1))/r(i)
end
c=input('请输入c：');
for k=2:8 
if x(k-1)<=c & c<=x(k) 

S=M(k-1)/6/h(k)*(x(k)-c)^3+M(k)/6/h(k)*(c-x(k-1))^3+1/h(k)*(y(k)-M(k)*h(k)^2/6)*(c-x(k-1))+1/h(k)*(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-c)%第一公式

%S=((((x(k)-c)^3)/(6*h(k)))*M(k-1))+(((c-x(k-1))^3)/(6*h(k))*M(k))+((y(k)-(M(k-1)*(((h(k))^2)/6)))*(x(k)-c)/h(k))+((y(k)-(M(k)*((h(k))^2/6)))*(c-x(k-1))/h(k))%第二公式
return; 
end 
End
三次样条插值：
syms xi
x=[-1 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.5];
y=[-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31];
xi=-1:0.01:3.5;
yi=spline(x,y,xi);
plot(x,y,'o',xi,yi)%画图
title('三次样条插值');
xlabel('x');
ylabel('y');

4、结果 

注意：二阶差商中“(x(i+2)-x(i))”出现8超过迭代i最大值7，故用(h(i+1)+h(i))代替，用第一种方法与第二种方法求得M相同：

M = [ -99.8212；55.0298；-17.2683；1.7565；15.2518；-34.0445；59.0295；-90.7190；158.3059 ]

令x=-0.02，得S =-3.1058。

令x=2.56，得S =4.7834。

对比MATLAB三次样条插值和五次多项式拟合（结果见下表）：

  

    

若增加运算次数：第二个公式（程序中标了注释），怀疑在编写后程序运行时：MATLAB内部计算时步骤太多，导致放大误差，结果失真。具体见下表

令x=-0.02，得S =-2.3468。

令x=2.56，得S =7.8756。