完整word版数值分析作业三次样条插值文档格式.docx
《完整word版数值分析作业三次样条插值文档格式.docx》由会员分享,可在线阅读,更多相关《完整word版数值分析作业三次样条插值文档格式.docx(19页珍藏版)》请在冰豆网上搜索。
理论上证明三次样条插值函数的收敛性是比较困难的,通过本实验可以证明这一理论结果。
按照一定的规则分别选择等距或非等距的插值节点,并不断增加插值节点的个数。
实验要求:
(1)随着节点个数的增加,比较被逼近函数和三样条插值函数的误差变化情况,分析所得结果并与拉格朗日插值多项式比较;
(2)三次样条插值函数的思想最早产生于工业部门.作为工业应用的例子,考虑如下例子:
某汽车制造商根据三次样条插值函数设计车门曲线,其中一段数据如下:
1
4
5
6
7
8
9
10
0.79
1.53
2.19
2。
71
3。
03
3.27
89
06
3.19
3.29
0.8
算法描述:
拉格朗日插值:
其中
是拉格朗日基函数,其表达式为:
牛顿插值:
三样条插值:
所谓三次样条插值多项式Sn(x)是一种分段函数,它在节点Xi(a〈X0<
X1……<
Xn〈b)分成的每个小区间[xi—1,xi]上是三次多项式,其在此区间上的表达式如下:
式中Mi=
.
因此,只要确定了Mi的值,就确定了整个表达式,Mi的计算方法如下:
令
则Mi满足如下n—1个方程:
常用的边界条件有如下几类:
(1)给定区间两端点的斜率m0,mn,即
(2)给定区间两端点的二阶导数M0,Mn,即
(3)假设y=f(x)是以b-a为周期的周期函数,则要求三次样条插值函数S(x)也为周期函数,对S(x)加上周期条件
对于第一类边界条件有
对于第二类边界条件有
那么解就可以为
对于第三类边界条件,
由此推得
,其中
,那么解就可以为:
程序代码:
1拉格朗日插值函数
Lang。
m
functionf=lang(X,Y,xi)
%X为已知数据的横坐标
%Y为已知数据的纵坐标
%xi插值点处的横坐标
%f求得的拉格朗日插值多项式的值
n=length(X);
f=0;
fori=1:
n
l=1;
forj=1:
i-1
l=l。
*(xi—X(j))/(X(i)—X(j));
end;
forj=i+1:
l=l.*(xi—X(j))/(X(i)-X(j));
end;
%拉格朗日基函数
f=f+l*Y(i);
end
fprintf('
%d\n'
,f)
return
2牛顿插值函数
newton.m
functionf=newton(X,Y,xi)
%X为已知数据的横坐标
%Y为已知数据的纵坐标
%f求得的拉格朗日插值多项式的值
n=length(X);
newt=[X’,Y'
];
%计算差商表
forj=2:
fori=n:
-1:
ifi〉=j
Y(i)=(Y(i)-Y(i—1))/(X(i)—X(i—j+1));
elseY(i)=0;
newt=[newt,Y'
];
end
%计算牛顿插值
f=newt(1,2);
fori=2:
z=1;
fork=1:
z=(xi—X(k))*z;
f=f+newt(i—1,i)*z;
fprintf(’%d\n'
return
3三次样条插值第一类边界条件
Threch。
functionS=Threch1(X,Y,dy0,dyn,xi)
%X为已知数据的横坐标
%S求得的三次样条插值函数的值
%dy0左端点处的一阶导数
%dyn右端点处的一阶导数
n=length(X)—1;
d=zeros(n+1,1);
h=zeros(1,n-1);
f1=zeros(1,n—1);
f2=zeros(1,n—2);
fori=1:
n%求函数的一阶差商
h(i)=X(i+1)—X(i);
f1(i)=(Y(i+1)-Y(i))/h(i);
end
fori=2:
n%求函数的二阶差商
f2(i)=(f1(i)—f1(i-1))/(X(i+1)—X(i-1));
d(i)=6*f2(i);
d
(1)=6*(f1
(1)-dy0)/h
(1);
d(n+1)=6*(dyn—f1(n—1))/h(n—1);
%¸
赋初值
A=zeros(n+1,n+1);
B=zeros(1,n—1);
C=zeros(1,n-1);
fori=1:
n-1
B(i)=h(i)/(h(i)+h(i+1));
C(i)=1—B(i);
A(1,2)=1;
A(n+1,n)=1;
n+1
A(i,i)=2;
n
A(i,i-1)=B(i-1);
A(i,i+1)=C(i—1);
M=A\d;
symsx;
Sx(i)=collect(Y(i)+(f1(i)—(M(i)/3+M(i+1)/6)*h(i))*(x—X(i))..。
+M(i)/2*(x—X(i))^2+(M(i+1)—M(i))/(6*h(i))*(x-X(i))^3);
digits(4);
Sx(i)=vpa(Sx(i));
%三样条插值函数表达式
disp(’S(x)='
);
fprintf('
%s(%d,%d)\n'
,char(Sx(i)),X(i),X(i+1));
ifxi>
=X(i)&&
xi〈=X(i+1)
S=Y(i)+(f1(i)—(M(i)/3+M(i+1)/6)*h(i))*(xi—X(i))+M(i)/2*(xi—X(i))^2+(M(i+1)—M(i))/(6*h(i))*(xi-X(i))^3;
disp('
xiS’);
%d,%d\n’,xi,S);
return
4三次样条插值第二类边界条件
Threch2.m
function[Sx]=Threch2(X,Y,d2y0,d2yn,xi)
X为已知数据的横坐标
%d2y0左端点处的二阶导数
%d2yn右端点处的二阶导数
n=length(X)-1;
d=zeros(n+1,1);
h=zeros(1,n—1);
f1=zeros(1,n—1);
f2=zeros(1,n-2);
n%求一阶差商
h(i)=X(i+1)-X(i);
n%求二阶差商
f2(i)=(f1(i)—f1(i—1))/(X(i+1)-X(i-1));
d(i)=6*f2(i);
d
(1)=2*d2y0;
d(n+1)=2*d2yn;
%赋初值
A=zeros(n+1,n+1);
B=zeros(1,n-1);
n—1
B(i)=h(i)/(h(i)+h(i+1));
A(1,2)=0;
A(n+1,n)=0;
A(i,i)=2;
fori=2:
A(i,i—1)=B(i-1);
A(i,i+1)=C(i-1);
symsx;
Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x—X(i))。
.。
+M(i)/2*(x—X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);
);
%s(%d,%d)\n'
char(Sx(i)),X(i),X(i+1));
=X(i)&
&xi〈=X(i+1)
S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi—X(i))^2+(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3;
%d,%d\n'
,xi,S);
5插值节点处的插值结果
main3。
clear
clc
X=[0。
0,0。
1,0。
2,0.3,0。
4];
Y=[0.5000,0.5398,0.5793,0.6179,0.7554];
xi=0。
13;
%xi=0.36;
disp(’xi=0.13’);
%disp('
xi=0.36'
disp(’拉格朗日插值结果'
lang(X,Y,xi);
disp(’牛顿插值结果’);
newton(X,Y,xi);
disp(’三次样条第一类边界条件插值结果’);
Threch1(X,Y,0.40,0。
36,xi);
%0.4,0。
36分别为两端点处的一阶导数
disp(’三次样条第二类边界条件插值结果'
Threch2(X,Y,0,-0.136,xi);
%0,-0。
136分别为两端点处的二阶导数
6将多种插值函数即原函数图像画在同一张图上
main2.m
0,0.1,0.2,0.3,0。
4];
Y=[0。
5000,0。
5398,0。
5793,0.6179,0。
7554];
a=linspace(0,0。
4,21);
NUM=21;
L=zeros(1,NUM);
N=zeros(1,NUM);
S=zeros(1,NUM);
B=zeros(1,NUM);
NUM
xi=a(i);
L(i)=lang(X,Y,xi);
%拉格朗日插值
N(i)=newton(X,Y,xi);
%牛顿插值
B(i)=normcdf(xi,0,1);
%原函数
S(i)=Threch1(X,Y,0。
4,0。
%三次样条函数第一类边界条件
plot(a,B,’——r'
holdon;
plot(a,L,'
b’);
plot(a,N,’r’);
holdon;
plot(a,S,'
r+’);
legend(’原函数'
,'
拉格朗日插值’,’牛顿插值'
三次样条插值'
2);
holdoff
7增加插值节点观察误差变化
main4。
clear;
clc;
N=5;
%4.5第一问
Ini=zeros(1,1001);
a=linspace(—1,1,1001);
Ini=1。
/(1+25*a.^2);
3%节点数量变化次数
N=2*N;
t=linspace(-1,1,N+1);
%插值节点
ft=1./(1+25*t.^2);
%插值节点函数值
val=linspace(-1,1,101);
forj=1:
101
L(j)=lang(t,ft,val(j));
S(j)=Threch1(t,ft,0。
074,—0.074,val(j));
%三样条第一类边界条件插值
plot(a,Ini,'
k'
)%原函数图象
holdon
plot(val,L,'
r'
)%拉格朗日插值函数图像
plot(val,S,'
b’)%三次样条插值函数图像
str=sprintf('
插值节点为%d时的插值效果’,N);
title(str);
legend(’原函数’,’拉格朗日插值'
%显示图例
holdoff
figure
8车门曲线
main5.m
X=[0,1,2,3,4,5,6,7,8,9,10];
Y=[0.0,0。
79,1。
53,2。
19,2.71,3.03,3。
27,2.89,3。
06,3。
19,3.29];
dy0=0。
8;
dyn=0。
2;
n=length(X)-1;
h=zeros(1,n-1);
f2=zeros(1,n—2);
nh(i)=X(i+1)—X(i);
f1(i)=(Y(i+1)—Y(i))/h(i);
nf2(i)=(f1(i)—f1(i—1))/(X(i+1)-X(i-1));
d
(1)=6*(f1
(1)-dy0)/h
(1);
d(n+1)=6*(dyn-f1(n-1))/h(n—1);
A=zeros(n+1,n+1);
B=zeros(1,n—1);
C(i)=1-B(i);
A(1,2)=1;
A(n+1,n)=1;
A(i,i+1)=C(i-1);
M=A\d;
x=zeros(1,n);
S=zeros(1,n);
x(i)=X(i)+0.5;
S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x(i)—X(i))+M(i)/2*(x(i)-X(i))^2+(M(i+1)—M(i))/(6*h(i))*(x(i)—X(i))^3;
plot(X,Y,'
holdon;
plot(x,S,’o’);
title(’三次样条插值效果图’);
legend('
已知插值节点'
’三次样条插值’);
holdoff
实验结果:
1计算插值节点处的函数值
xi=0.13时
Xi=0.36时
2将多种插值函数即原函数图像画在同一张图上
4.5.1增加插值节点观察误差变化
从上面三张图可以看出增加插值节点并不能改善差之效果
5。
2车门曲线