三次样条插值---三弯矩法（Matlab实现）

目录

1、概述

2、代码

3、结果

---

1、概述

      在前面一章我们已经总结完知识点了，我们这篇文章直接用Matlab进行编程，下面直接上代码和运行结果。

2、代码

function Q=san(ssss,p) 
Q=zeros(2,1); 
x=[1;2;3;4;5;6;7;8;9;10]; 
y=[0;0.69314718;1.0986123;1.3862944;1.6094378;1.7917595;
1.9459101;2.079445;2.1972246;2.3025851]; 
h=zeros(10,1); 
d=zeros(10,1); 
u=zeros(10,1); 
v=zeros(10,1); 
r=zeros(10,1); 
l=zeros(10,1); 
z=zeros(10,1); 
m=zeros(10,1); 
for t=1:1:9; 
    h(t)=x(t+1)-x(t); 
end 
d(1)=6/h(1)*((y(2)-y(1))/h(1)-1); 
d(10)=6/h(9)*(0.1-(y(10)-y(9))/h(9)); 
for t=1:1:8 
	u(t+1)=h(t)/(h(t)+h(t+1)); 
v(t+1)=1-u(t+1); 
d(t+1)=6/(h(t)+h(t+1))*((y(t+2)-y(t+1))/(x(t+2)-x(t+1))-(y(t+1)-y (t))/(x(t+1)-x(t))); 
end 
u(10)=1;v(1)=1;r(1)=d(1); 
for t=2:1:10 
	l(t)=u(t)/r(t-1); 
	r(t)=d(t)-l(t)*v(t-1); 
end 
z(1)=d(1); 
for t=2:1:10 
	z(t)=d(t)-l(t)*z(t-1); 
end 
m(10)=z(10)/r(10); 
for t=9:-1:1 
	m(t)=(z(t)-v(t)*m(t+1))/r(t); 
end 
for t=1:1:10 
	if p>=t&&p<(t+1) 
		Q(:,1)=feval(ssss,p,t,x,m,h,y);break 
	end
end 
function Q=ssss(p,t,x,m,h,y) 
Q=zeros(2,1); 
Q(1,1)=((power((x(t+1)-p),3)*m(t)+power((p-x(t)),3)*m(t+1))/6+(y( t)-m(t)*h(t)*h(t)/6)*(x(t+1)-p)+(y(t+1)-m(t+1)*h(t)*h(t)/6)*(p-x( t)))/h(t);
Q(2,1)=(-(power((x(t+1)-p),2)*m(t)+power((p-x(t)),2)*m(t+1))/2+(y (t)-m(t)*h(t)*h(t)/6)+(y(t+1)-m(t+1)*h(t)*h(t)/6))/h(t); 
End

3、结果