clc;clear
syms x
%输入目标函数f、区间x、节点横坐标xi（均可任意输入）
fx=1/(1+x^2);
x=linspace(-5,5,1000);
xi=-5:1:5;
%绘制目标函数图像
figure('numbertitle','off','name','三次样条插值法')
plot(x,eval(fx),'k','linewidth',1);
xlabel('x','FontSize',20)
ylabel('y','FontSize',20)
title('三次样条插值函数与目标函数图像的比较','FontSize',15)
hold on
grid on
%求解矩阵h,u,lambda（r）,d,进而求解线性方程组，得到多项式系数矩阵M
n=length(xi);
dfx=diff(fx);
for i=1:n
x=xi(i);
y(i)=eval(fx);
dy(i)=eval(dfx);
end
for i=1:n-1
h(i)=xi(i+1)-xi(i);
end
for i=1:n-2
r(i+1)=h(i+1)/(h(i)+h(i+1));
u(i)=1-r(i+1);
d(i+1)=6*deviation(i,i+2,xi,y);
end
r(1)=1;
u(n-1)=1;
d(1)=6*(deviation(1,2,xi,y)-dy(1))/h(1);
d(n)=6*(dy(n)-deviation(n-1,n,xi,y))/h(n-1);
C=2*eye(n)+diag(r,1)+diag(u,-1);
M=C\d';
%用求得的多项式系数矩阵M通过循环画出每个区间的插值函数图像
for i=1:n-1
t=xi(i):0.01:xi(i+1);
plot(t,M(i)*(xi(i+1)-t).^3/(6*h(i))+M(i+1)*(t-xi(i)).^3/(6*h(i))+...
(y(i)-M(i)*h(i).^2/6)*(xi(i+1)-t)/h(i)+(y(i+1)-M(i+1)*h(i).^2/6)*(t-xi(i))/h(i));
hold on
end
%调用求解f[x(i-1),x(i),x(i+1)]的函数
function f=deviation(i,j,x,y)
if i<j
f=(deviation(i+1,j,x,y)-deviation(i,j-1,x,y))/(x(j)-x(i));
elseif i==j
f=y(i);
end
end
