matlab差值 c代码,三次样条插值用C/C++如何得到如图所示的曲线函数-CSDN博客

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