例:已知一组数据点,编写一程序求解三次样条插值函数满足
并针对下面一组具体实验数据
求解,其中边界条件为.
1)三次样条插值自然边界条件源程序:
function s=spline3(x,y,dy1,dyn)
%x为节点,y为节点函数值,dy1,dyn分别为x=0.25,0.53处的二阶导
m=length(x);n=length(y);
if m~=n
error('x or y输入有误')
return
end
h=zeros(1,n-1);
h(n-1)=x(n)-x(n-1);
for k=1:n-2
h(k)=x(k+1)-x(k);
v(k)=h(k+1)/(h(k+1)+h(k));
u(k)=1-v(k);
end
g(1)=3*(y(2)-y(1))/h(1)-h(1)/2*dy1;
g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)/2*dyn;
for i=2:n-1
g(i)=3*(u(i-1)*(y(i+1)-y(i))/h(i)+v(i-1)*(y(i)-y(i-1))/h(i-1)); end
for i=2:n-1;
A(i,i-1)=v(i-1);
A(i,i+1)=u(i-1);
end
A(n,n-1)=1;
A(1,2)=1;
A=A+2*eye(n);
M=zhuigf(A,g); %调用函数,追赶法求M
fprintf('三次样条(三对角)插值的函数表达式\n');
syms X;
for k=1:n-1
fprintf('S%d--%d:\n',k,k+1);
s(k)=(h(k)+2*(X-x(k)))./h(k).^3.*(X-x(k+1)).^2.*y(k)...
+(h(k)-2*(X-x(k+1)))./h(k).^3.*(X-x(k)).^2.*y(k+1)... +(X-x(k)).*(X-x(k+1)).^2./h(k).^2*M(k)+(X-x(k+1)).*... (X-x(k)).^2./h(k).^2*M(k+1);
end
s=s.';
s=vpa(s,4);
%画三次样条插值函数图像
for i=1:n-1
X=x(i):0.01:x(i+1);
st=(h(i)+2*(X-x(i)))./(h(i)^3).*(X-x(i+1)).^2.*y(i)...
+(h(i)-2.*(X-x(i+1)))./(h(i)^3).*(X-x(i)).^2.*y(i+1)...
+(X-x(i)).*(X-x(i+1)).^2./h(i)^2*M(i)+(X-x(i+1)).*...
(X-x(i)).^2./h(i)^2*M(i+1);
plot(x,y,'o',X,st);
hold on
End
plot(x,y);
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%调用的函数:
%追赶法
function M=zhuigf(A,g)
n=length(A);
L=eye(n);
U=zeros(n);
for i=1:n-1
U(i,i+1)=A(i,i+1);
end
U(1,1)=A(1,1);
for i=2:n
L(i,i-1)=A(i,i-1)/U(i-1,i-1);
U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i);
end
Y(1)=g(1);
for i=2:n
Y(i)=g(i)-L(i,i-1)*Y(i-1);
end
M(n)=Y(n)/U(n,n);
for i=n-1:-1:1
M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i);
end
2)在命令窗口输入x,y,dy1,dyn,得到三次样条函数:
>>x=[0.25,0.3,0.39,0.45,0.53];
>>y=[0.5,0.5477,0.6245,0.6708,0.7280];
>>dy1=0;
>>dyn=0;
>>s=spline3(x,y,dy1,dyn)
运行结果:
三次样条(三对角)插值的函数表达式
S1--2:
S2--3:
S3--4:
S4--5:
.5000*(-3600.+.1600e5*X)*(X-.3000)^2+.5477*(5200.-.1600e5*X)*(X-.2500)^2+394.8*(X-.2500)*(X-.3000)^2+355.2*(X-.3000)*(X-.2500)^2
.5477*(-699.6+2743.*X)*(X-.3900)^2+.6245*(1193.-2743.*X)*(X-.3000)^2+109.6*(X-.300 0)*(X-.3900)^2+96.77*(X-.3900)*(X-.3000)^2
.6245*(-3333.+9259.*X)*(X-.4500)^2+.6708*(4444.-9259.*X)*(X-.3900)^2+217.7*(X-.39 00)*(X-.4500)^2+207.6*(X-.4500)*(X-.3900)^2
.6708*(-1602.+3906.*X)*(X-.5300)^2+.7280*(2227.-3906.*X)*(X-.4500)^2+116.8*(X-.45 00)*(X-.5300)^2+109.2*(X-.5300)*(X-.4500)^2
如将三次样条函数加以整理,可用如下程序:
s=collect(s);
则输出结果为:
s=
.4595000000000-13.200*X^3+9.9000000*X^2-1.48800000000*X
.37459306800000-4.2924*X^3+3.66332200*X^2-.137976090000*X
.5180681700000-3.3917*X^3+3.90576600*X^2-.733130370000*X
-.408614400000e-1+2.5768*X^3-4.03188800*X^2+2.868770320000*X。