function Matrix4=D2Yspline3(X,Y,d2Y,m,IsPlot)
%% 说明
%X 插值点
%Y 插值点对应的函数值
%dY 边界条件
%m 精度
%% 使用示例
% clear;clc;
% X=[1 2 4 5];
% Y=[1 3 4 2];
% d2Y=[0 0];
% m=9;
% mm=D2Yspline3(X,Y,d2Y,m,true)
%% IsPlot
if nargin==4
IsPlot=false;
end
%% 初值
N=size(X,2);%插值点个数
%导数边界条件
s0=d2Y(1);
sn=d2Y(2);
%画图插值精度
interval=0.01*(X(2)-X(1));
%计算间隔h
h=X(2:N)-X(1:N-1);
%% 计算方程系数
%计算矩阵次对角线系数
mu=zeros(1,N-1);
md=zeros(1,N-1);
md(1:N-2)=h(2:N-1)./(h(2:N-1)+h(1:N-2));
md(N-1)=1;
mu(1)=1;
mu(2:N-1)=1-md(1:N-2);
%计算方程右端系数d
d=zeros(1,N);
d(2:N-1)=3*(mu(2:N-1).*(Y(3:N)-Y(2:N-1))./h(2:N-1)+md(1:N-2).*(Y(2:N-1)-Y(1:N-2))./h(1:N-2));
d(1)=3*(Y(2)-Y(1))/h(1)-h(1)/2*s0;
d(N)=3*(Y(N)-Y(N-1))/h(N-1)-h(N-1)/2*sn;
%% 以下采用追赶法计算
%主对角线化为1,q为上次对角线的系数
p=zeros(1,N);
q=zeros(1,N);
p(1)=2;
q(1)=mu(1)/2;
for i=2:N-1
p(i)=2-md(i-1)*q(i-1);
q(i)=mu(i)/p(i);
end
p(N)=2-md(N-1)*q(N-1);
%y为调整系数后,右端项的系数
y=zeros(1,N);
y(1)=d(1)/2;
for i=2:N
y(i)=(d(i)-md(i-1)*y(i-1))/p(i);
end
%x为相应方程的解
x=zeros(1,N);
x(N)=y(N);
for i=N-1:-1:1
x(i)=y(i)-q(i)*x(i+1);
end
%以下用M替换x
M=x;
%m的作用,设定精度,小数点后几位
digits(m);
%% 以下为插值函数的方程,未知数为t
syms t;
pp(1:N-1)=(1-2*(X(1:N-1)-t)./h(1:N-1)).*((X(2:N)-t)./h(1:N-1)).*((X(2:N)-t)./h(1:N-1)).*Y(1:N-1)+(1-2*(t-X(2:N))./h(1:N-1)).*((t-X(1:N-1))./h(1:N-1)).*((t-X(1:N-1))./h(1:N-1)).*Y(2:N)+(t-X(1:N-1)).*((X(2:N)-t)./h(1:N-1)).*((X(2:N)-t)./h(1:N-1)).*M(1:N-1)+(t-X(2:N)).*((t-X(1:N-1))./h(1:N-1)).*((t-X(1:N-1))./h(1:N-1)).*M(2:N);
pp(1:N-1)=simplify(pp(1:N-1));%化简
%以下计算每段多项式的系数
Matrix4=zeros(N-1,4);
for i=1:N-1
%提取多项式系数
coeff=sym2poly(pp(i));
%如果多项式的最高次项不是3,则为了通用,前面需要补0
len=length(coeff);
if len~=4
tt=coeff(1:len);
coeff(1:4)=0;
coeff(4-len+1:end)=tt;
end
Matrix4(i,:)=coeff;
%% 输出画图 可以考虑注释掉
if(IsPlot==true)
if(i==1)
fprintf('S(x)的表达式为:\n');
end
val=X(i):interval:X(i+1);
k=length(val);
fval(1:k)=coeff(1).*val(1:k).^3+coeff(2).*val(1:k).^2+coeff(3).*val(1:k)+coeff(4);
if mod(i,2)==1
plot(val,fval,'r-');
else
plot(val,fval,'b-');
end
hold on;grid on;
clear val fval;
temp=sym(coeff,'d');
temp=poly2sym(temp,'t');
fprintf('在区间[%.4f,%.4f]内,',X(i),X(i+1));
fprintf('三次样条函数S(%d)=',i);
disp(temp);
end
end
end