习题一：三弯矩法（Matlab实现）

目录

1、算法简介

 2、题目

 3、MATLAB代码

---

1、算法简介

 2、题目

 3、MATLAB代码

function yi=cubic_spline1(x,y,ydot,xi) 
% Cubic spline interpolation formula 
% X is the vector, all interpolation nodes； 
% Y is the vector and the function value at the interpolation node; 
% Ydot is the derivative value at the endpoint of the vector; 
% Xi is scalar, independent variable x;
% Yi is the function estimate of Xi;
n=length(x);ny=length(y); 
if n~=ny 
error('X and y are the same length'); 
return; 
end 
if isempty(ydot)==1 
ydot=[(y(2)-y(1))/(x(2)-x(1)) (y(n)-y(n-1))/(x(n)-x(n-1))];

end 
h=zeros(1:n);lambda=ones(1,n);mu=ones(1,n); 
M=zeros(n,1);d=zeros(n,1); 
for k=2:n 
h(k)=x(k)-x(k-1); 
if abs(h(k))<eps 
error('the DATA is error'); 
return; 
end 
end 
for k=2:n-1 
lambda(k)=h(k+1)/(h(k)+h(k+1));mu(k)=1-lambda(k); 
d(k)=6/(h(k)+h(k+1))... 
*((y(k+1)-y(k))/h(k+1)-(y(k)-y(k-1))/h(k)); 
end 
d(1)=6/h(2)*((y(2)-y(1))/h(2)-ydot(1)); 
d(n)=6/h(n)*(ydot(2)-(y(n)-y(n-1))/h(n)); 
A=diag(2*ones(1,n)); 
for i=1:n-1 
A(i,i+1)=lambda(i);A(i+1,i)=mu(i+1); 
end 
M=A\d; 
for k=2:n 
if x(k-1)<=xi & xi<=x(k) 
yi=M(k-1)/6/h(k)*(x(k)-xi)^3 
+M(k)/6/h(k)*(xi-x(k-1))^3 
+1/h(k)*(y(k)-M(k)*h(k)^2/6)*(xi-x(k-1)) 
+1/h(k)*(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-xi); 
return; 
end 
end