西北农林科技大学数值分析报告数值法实验报告材料.docx
《西北农林科技大学数值分析报告数值法实验报告材料.docx》由会员分享，可在线阅读，更多相关《西北农林科技大学数值分析报告数值法实验报告材料.docx（21页珍藏版）》请在冰点文库上搜索。
西北农林科技大学数值分析报告数值法实验报告材料
数值法实验报告
专业班级：
信息与计算科学121：
金辉学号：
2012014280
1）实验目的
本次实验的目的是熟练《数值分析》第二章“插值法”的相关容，掌握三种插值方法：
牛顿多项式插值，三次样条插值，拉格朗日插值，并比较三种插值方法的优劣。
本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并在MATLAB软件中去实现。
2）实验题目
实验一：
已知函数在下列各点的值为
xi
0.2
0.4
0.6
.0.8
1.0
f（xi）
0.98
0.92
0.81
0.64
0.38
试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值。
用图给出{（xi，yi），xi=0.2+0.08i，i=0，1,11,10}，P4（x）及S（x）。
实验二：
在区间[-1,1]上分别取
用两组等距节点对龙格函数
作多项式插值及三次样条插值，对每个
值，分别画出插值函数即
的图形。
实验三：
下列数据点的插值
x
0
1
4
9
16
25
36
49
64
y
0
1
2
3
4
5
6
7
8
可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9各点作8次多项式插值L8（x）.
（2）用三次样条（自然边界条件）程序求S（x）。
从结果看在[0,64]上，那个插值更精确；在区间[0,1]上，两种哪个更精确？
3）实验原理与理论基础
《数值分析》第二章“插值法”的相关容，包括：
牛顿多项式插值，三次样条插值，拉格朗日
4）实验容
实验一：
已知函数在下列各点的值为
xi
0.2
0.4
0.6
.0.8
1.0
f（xi）
0.98
0.92
0.81
0.64
0.38
试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值。
用图给出{（xi，yi），xi=0.2+0.08i，i=0，1,11,10}，P4（x）及S（x）。
（1）首先我们先求牛顿插值多项式，此处要用4次牛顿插值多项式处理数据。
已知n次牛顿插值多项式如下：
Pn=f（x0）+f[x0,x1]（x-x0）+f[x0,x1,x2]（x-x0）（x-x1）+···+f[x0,x1，···xn]（x-x0）···（x-xn-1）
我们要知道牛顿插值多项式的系数，即均差表中得部分均差。
在MATLAB的Editor中输入程序代码，计算牛顿插值中多项式系数的程序如下：
functionvarargout=newtonliu（varargin）
clear,clc
x=[0.20.40.60.81.0];
fx=[0.980.920.810.640.38];
newtonchzh（x,fx）;
functionnewtonchzh（x,fx）
%由此函数可得差分表
n=length（x）;
fprintf（'*****************差分表*****************************\n'）;
FF=ones（n,n）;
FF（:
1）=fx';
fori=2:
n
forj=i:
n
FF（j,i）=（FF（j,i-1）-FF（j-1,i-1））/（x（j）-x（j-i+1））;
end
end
fori=1:
n
fprintf（'%4.2f',x（i））;
forj=1:
i
fprintf（'%10.5f',FF（i,j））;
end
fprintf（'\n'）;
end
由MATLAB计算得：
xi
f（xi）
一阶差商
二阶差商
三阶差商
四阶差商
0.20
0.980
0.40
0.920
-0.30000
0.60
0.810
-0.55000
-0.62500
0.80
0.640
-0.85000
-0.75000
-0.20833
1.00
0.380
-1.30000
-1.12500
-0.62500
-0.52083
所以有四次插值牛顿多项式为：
P4（x）=0.98-0.3（x-0.2）-0.62500（x-0.2）（x-0.4）-0.20833（x-0.2）（x-0.4）（x-0.6）-0.52083（x-0.2）（x-0.4）（x-0.6）（x-0.8）
（2）接下来我们求三次样条插值函数。
用三次样条插值函数由上题分析知,要求各点的M值：
三次样条插值函数计算的程序如下：
functiontgsanci（n,s,t）%n代表元素数，s,t代表端点的一阶导。
x=[0.20.40.60.81.0];
y=[0.980.920.810.640.38];
n=5
forj=1:
1:
n-1
h（j）=x（j+1）-x（j）;
end
forj=2:
1:
n-1
r（j）=h（j）/（h（j）+h（j-1））;
end
forj=1:
1:
n-1
u（j）=1-r（j）;
end
forj=1:
1:
n-1
f（j）=（y（j+1）-y（j））/h（j）;
end
forj=2:
1:
n-1
d（j）=6*（f（j）-f（j-1））/（h（j-1）+h（j））;
end
d
（1）=0
d（n）=0
a=zeros（n,n）;
forj=1:
1:
n
a（j,j）=2;
end
r
（1）=0;
u（n）=0;
forj=1:
1:
n-1
a（j+1,j）=u（j+1）;
a（j,j+1）=r（j）;
end
b=inv（a）
m=b*d'
p=zeros（n-1,4）;%p矩阵为S（x）函数的系数矩阵
forj=1:
1:
n-1
p（j,1）=m（j）/（6*h（j））;
p（j,2）=m（j+1）/（6*h（j））;
p（j,3）=（y（j）-m（j）*（h（j）^2/6））/h（j）;
p（j,4）=（y（j+1）-m（j+1）*（h（j）^2/6））/h（j）;
end
p
得到m=（0-1.6071-1.0714-3.10710）T
即M0=0;M1=-1.6071;M2=-1.0714;M3=-3.1071;M4=0
则根据三次样条函数定义，可得：
S（x）=
接着，在CommandWindow里输入画图的程序代码，
下面是画牛顿插值以及三次样条插值图形的程序：
x=[0.20.40.60.81.0];
y=[0.980.920.810.640.38];
plot（x,y）
holdon
fori=1:
1:
5
y（i）=0.98-0.3*（x（i）-0.2）-0.62500*（x（i）-0.2）*（x（i）-0.4）-0.20833*（x（i）-0.2）*（x（i）-0.4）*（x（i）-0.6）-0.52083*（x（i）-0.2）*（x（i）-0.4）*（x（i）-0.6）*（x（i）-0.8）
end
k=[011011]
x0=0.2+0.08*k
fori=1:
1:
4
y0（i）=0.98-0.3*（x（i）-0.2）-0.62500*（x（i）-0.2）*（x（i）-0.4）-0.20833*（x（i）-0.2）*（x（i）-0.4）*（x（i）-0.6）-0.52083*（x（i）-0.2）*（x（i）-0.4）*（x（i）-0.6）*（x（i）-0.8）
end
plot（x0,y0,'o',x0,y0）
holdon
y1=spline（x,y,x0）
plot（x0,y1,'o'）
holdon
s=csape（x,y,'variational'）
fnplt（s,'r'）
holdon
gtext（'三次样条自然边界'）
gtext（'原图像'）
gtext（'4次牛顿插值'）
运行上述程序可知：
给出的{（xi，yi），xi=0.2+0.08i，i=0，1,11,10}点，S（x）及P4（x）图形如下所示：
实验二：
在区间[-1,1]上分别取
用两组等距节点对龙格函数
作多项式插值及三次样条插值，对每个
值，分别画出插值函数即
的图形。
我们先求多项式插值：
在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令（如牛顿插值公式）：
function[C,D]=newpoly（X,Y）
n=length（X）;
D=zeros（n,n）
D（:
1）=Y'
forj=2:
n
fork=j:
n
D（k,j）=（D（k,j-1）-D（k-1,j-1））/（X（k）-X（k-j+1））;
end
end
C=D（n,n）;
fork=（n-1）:
-1:
1
C=conv（C,poly（X（k）））
m=length（C）;
C（m）=C（m）+D（k,k）;
end
当n=10时，我们在CommandWindow中输入以下的命令：
clear,clf,holdon;
X=-1:
0.2:
1;
Y=1./（1+25*X.^2）;
[C,D]=newpoly（X,Y）;
x=-1:
0.01:
1;
y=polyval（C,x）;
plot（x,y,X,Y,'.'）;
gridon;
xp=-1:
0.2:
1;
z=1./（1+25*xp.^2）;
plot（xp,z,'r'）
得到插值函数和f（x）图形：
当n=20时，我们在CommandWindow中输入以下的命令：
clear,clf,holdon;
X=-1:
0.1:
1;
Y=1./（1+25*X.^2）;
[C,D]=newpoly（X,Y）;
x=-1:
0.01:
1;
y=polyval（C,x）;
plot（x,y,X,Y,'.'）;
gridon;
xp=-1:
0.1:
1;
z=1./（1+25*xp.^2）;
plot（xp,z,'r'）
得到插值函数和f（x）图形：
下面再求三次样条插值函数，在MATLAB的Editor中建立一个多项式的M-file,
输入下列程序代码：
functionS=csfit（X,Y,dx0,dxn）
N=length（X）-1;
H=diff（X）;
D=diff（Y）./H;
A=H（2:
N-1）;
B=2*（H（1:
N-1）+H（2:
N））;
C=H（2:
N）;
U=6*diff（D）;
B
（1）=B
（1）-H
（1）/2;
U
（1）=U
（1）-3*（D
（1））;
B（N-1）=B（N-1）-H（N）/2;
U（N-1）=U（N-1）-3*（-D（N））;
fork=2:
N-1
temp=A（k-1）/B（k-1）;
B（k）=B（k）-temp*C（k-1）;
U（k）=U（k）-temp*U（k-1）;
end
M（N）=U（N-1）/B（N-1）;
fork=N-2:
-1:
1
M（k+1）=（U（k）-C（k）*M（k+2））/B（k）;
end
M
（1）=3*（D
（1）-dx0）/H
（1）-M
（2）/2;
M（N+1）=3*（dxn-D（N））/H（N）-M（N）/2;
fork=0:
N-1
S（k+1,1）=（M（k+2）-M（k+1））/（6*H（k+1））;
S（k+1,2）=M（k+1）/2;
S（k+1,3）=D（k+1）-H（k+1）*（2*M（k+1）+M（k+2））/6;
S（k+1,4）=Y（k+1）;
end
当n=10时，我们在CommandWindow中输入以下的命令：
clear,clc
X=-1:
0.2:
1;
Y=1./（25*X.^2+1）;
dx0=0.0739644970414201;dxn=-0.0739644970414201;
S=csfit（X,Y,dx0,dxn）
x1=-1:
0.01:
-0.5;y1=polyval（S（1,:
）,x1-X
（1））;
x2=-0.5:
0.01:
0;y2=polyval（S（2,:
）,x2-X
（2））;
x3=0:
0.01:
0.5;y3=polyval（S（3,:
）,x3-X（3））;
x4=0.5:
0.01:
1;y4=polyval（S（4,:
）,x4-X（4））;
plot（x1,y1,x2,y2,x3,y3,x4,y4,X,Y,'.'）
结果如图：
当n=20时，我们在CommandWindow中输入以下的命令：
clear,clc
X=-1:
0.1:
1;
Y=1./（25*X.^2+1）;
dx0=0.0739644970414201;dxn=-0.0739644970414201;
S=csfit（X,Y,dx0,dxn）
x1=-1:
0.01:
-0.5;y1=polyval（S（1,:
）,x1-X
（1））;
x2=-0.5:
0.01:
0;y2=polyval（S（2,:
）,x2-X
（2））;
x3=0:
0.01:
0.5;y3=polyval（S（3,:
）,x3-X（3））;
x4=0.5:
0.01:
1;y4=polyval（S（4,:
）,x4-X（4））;
plot（x1,y1,x2,y2,x3,y3,x4,y4,X,Y,'.'）
结果如图：
实验三：
下列数据点的插值
x
0
1
4
9
16
25
36
49
64
y
0
1
2
3
4
5
6
7
8
可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9各点作8次多项式插值L8（x）.
（2）用三次样条（自然边界条件）程序求S（x）。
从结果看在[0,64]上，那个插值更精确；在区间[0,1]上，两种哪个更精确？
L8（x）可由公式Ln（x）=
得出。
三次样条可以利用自然边界条件。
写成矩阵：
其中
j=
，
i=
，dj=6f[xj-1,xj,xj+1]，
n=
0=0d0=dn=0
l0（x）=
l1（x）=
l2（x）=
l3（x）=
l4（x）=
l5（x）=
l6（x）=
l7（x）=
l8（x）=
L8（x）=l1（x）+2l2（x）+3l3（x）+4l4（x）+5l5（x）+6l6（x）+7l7（x）+8l8（x）
求三次样条插值函数由MATLAB计算：
可得矩阵形式的线性方程组为：
在MATLAB中的Editor中输入程序代码，
以下是三次样条函数的程序代码：
functiontgsanci（n,s,t）%n代表元素数，s,t代表端点的一阶导。
y=[012345678];
x=[01491625364964];
n=9
forj=1:
1:
n-1
h（j）=x（j+1）-x（j）;
end
forj=2:
1:
n-1
r（j）=h（j）/（h（j）+h（j-1））;
end
forj=1:
1:
n-1
u（j）=1-r（j）;
end
forj=1:
1:
n-1
f（j）=（y（j+1）-y（j））/h（j）;
end
forj=2:
1:
n-1
d（j）=6*（f（j）-f（j-1））/（h（j-1）+h（j））;
end
d
（1）=0
d（n）=0
a=zeros（n,n）;
forj=1:
1:
n
a（j,j）=2;
end
r
（1）=0;
u（n）=0;
forj=1:
1:
n-1
a（j+1,j）=u（j+1）;
a（j,j+1）=r（j）;
end
b=inv（a）
m=b*d'
t=a
p=zeros（n-1,4）;%p矩阵为S（x）函数的系数矩阵
forj=1:
1:
n-1
p（j,1）=m（j）/（6*h（j））;
p（j,2）=m（j+1）/（6*h（j））;
p（j,3）=（y（j）-m（j）*（h（j）^2/6））/h（j）;
p（j,4）=（y（j+1）-m（j+1）*（h（j）^2/6））/h（j）;
end
p
解得:
M0=0;M1=-0.5209;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=--0.0009;
M8=0，则三次样条函数：
S（x）=
下面进行画图，在CommandWindow中输入画图的程序代码：
%画图形比较那个插值更精确的函数：
x0=[01491625364964];
y0=[012345678];
x=0:
64;y=lagr1（x0,y0,x）;
plot（x0,y0,'o'）
holdon
plot（x,y,'r'）;
holdon;
pp=csape（x0,y0,'variational'）
fnplt（pp,'g'）;
holdon;
plot（x0,y0,':
b'）;holdon
%axis（[0201]）;%看[01]区间的图形时加上这条指令
gtext（'三次样条插值'）
gtext（'原图像'）
gtext（'拉格朗日插值'）
%下面是求拉格朗日插值的函数
functiony=lagr1（x0,y0,x）
n=length（x0）;m=length（x）;
fori=1:
m
z=x（i）;
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=k
p=p*（z-x0（j））/（x0（k）-x0（j））;
end
end
s=p*y0（k）+s;
end
y（i）=s;
end
拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三次样条插值曲线，蓝色虚线条为x=y2的曲线，另外一条红色线条为拉格朗日插值曲线。
图3-1为[01]的曲线，图3-2为[064]区间上的曲线。
图3—1
图3—2
由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线，在[0，1]朗格朗日插值更精确。
而在区间[0，64]上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的，而红色线条在[30，70]之间有很大的振荡，所在在区间[0，64]三次样条插值更精确写。
5）实验结果
单个多项式高次插值效果并不理想，有龙格现象，偏差大，没有使用价值。
而分段低次插值则精确度较高，拟合效果较好，而三次样条插值具有良好的收敛性与稳定性，与分段低次插值相比较光滑度更高，而且提供的信息也相对少一些。
我们可以看到，在以上的三道实验题里，我们可以从图形中看出，三次样条的拟合程度是三种插值函数里最好的。
6）实验结果分析与小结
通过此次实验，我对牛顿多项式插值，三次样条插值，拉格朗日插值有了更进一步的了解，知道了三次样条的拟合程度在高次的情况下更高，在理论上和应用上都有重要意义，在利用计算机编程软件进行高次插值的时候，我们可以多考虑利用三次样条进行插值。
升级会员 