西北农林科技大学数值分析数值法实验报告Word文档下载推荐.docx

1、1491625364964y235678可以得到平方根函数的近似，在区间0,64上作图。（1）用这9各点作8次多项式插值L8（x）.（2）用三次样条（自然边界条件）程序求S（x）。从结果看在0,64上，那个插值更精确；在区间0,1上，两种哪个更精确？3）实验原理与理论基础数值分析第二章“插值法”的相关内容，包括：牛顿多项式插值，三次样条插值，拉格朗日4）实验内容（1）首先我们先求牛顿插值多项式，此处要用4次牛顿插值多项式处理数据。已知n次牛顿插值多项式如下：Pn=f（x0）+fx0,x1（x-x0）+ fx0,x1,x2（x-x0） （x-x1）+ fx0,x1，xn（x-x0） （x-xn-

2、1）我们要知道牛顿插值多项式的系数，即均差表中得部分均差。在MATLAB的Editor中输入程序代码，计算牛顿插值中多项式系数的程序如下：function varargout=newtonliu（varargin）clear,clcx=0.2 0.4 0.6 0.8 1.0;fx=0.98 0.92 0.81 0.64 0.38;newtonchzh（x,fx）;function newtonchzh（x,fx）%由此函数可得差分表n=length（x）;fprintf（*差分表*n）;FF=ones（n,n）;FF（:,1）=fx;for i=2:n for j=i: FF（j,i）=（FF

3、（j,i-1）-FF（j-1,i-1）/（x（j）-x（j-i+1）; endendfor i=1: fprintf（%4.2f,x（i）; for j=1:i%10.5f,FF（i,j）;n由MATLAB计算得：xi f（xi） 一阶差商二阶差商三阶差商四阶差商0.200.9800.400.920-0.300000.600.810-0.55000-0.625000.800.640-0.85000-0.75000-0.208331.000.380-1.30000-1.12500-0.52083所以有四次插值牛顿多项式为：P4（x）=0.98-0.3（x-0.2）-0.62500 （x-0.2）

4、（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值：三次样条插值函数计算的程序如下：function tgsanci（n,s,t） %n代表元素数，s,t代表端点的一阶导。y=0.98 0.92 0.81 0.64 0.38; n=5for j=1:1:n-1 h（j）=x（j+1）-x（j）;for j=2: r（j）=h（j）/（h（j）+h（j-1）; u（j）=1-r（j）; f（j）=（y（j+1）-y（

5、j）/h（j）; d（j）=6*（f（j）-f（j-1）/（h（j-1）+h（j）; d（1）=0 d（n）=0 a=zeros（n,n）; a（j,j）=2; r（1）=0; u（n）=0; a（j+1,j）=u（j+1）; a（j,j+1）=r（j）;b=inv（a）m=b*dp=zeros（n-1,4）; %p矩阵为S（x）函数的系数矩阵 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）; p得到m=（0 -1

6、.6071 -1.0714 -3.1071 0）T 即M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0则根据三次样条函数定义，可得：S（x）= 接着，在Command Window里输入画图的程序代码，下面是画牛顿插值以及三次样条插值图形的程序：plot（x,y）hold on 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

7、）*（x（i）-0.8）k=0 1 10 11x0=0.2+0.08*ky0（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）plot（ x0,y0,o,x0,y0 ）hold ony1=spline（x,y,x0）plot（x0,y1,）s=csape（x,y,variational fnplt（s,rgtext（三次样条自然边界原图像4次牛顿插值运行

8、上述程序可知：给出的（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10点，S（x）及P4（x）图形如下所示：我们先求多项式插值：在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令（如牛顿插值公式）：function C,D=newpoly（X,Y）n=length（X）;D=zeros（n,n） D（:,1）=Y for k=j: D（k,j）=（D（k,j-1）- D（k-1,j-1）/（X（k）-X（k-j+1）;C=D（n,n）;for k=（n-1）:-1: C=conv（C,poly（X（k） m=length（C）; C（m）= C（m）+

9、D（k,k）;当n=10时，我们在Command Window中输入以下的命令：clear,clf,hold on; X=-1:0.2:1; Y=1./（1+25*X.2）; C,D=newpoly（X,Y）; x=-1:0.01: y=polyval（C,x）; plot（x,y,X,Y,. grid on; xp=-1: z=1./（1+25*xp.2）; plot（xp,z,得到插值函数和f（x）图形：当n=20时，我们在Command Window中输入以下的命令：0.1:下面再求三次样条插值函数，在MATLAB的Editor中建立一个多项式的M-file,输入下列程序代码：funct

10、ion S=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）;for k=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）;for k=

11、N-2: M（k+1）=（U（k）-C（k）*M（k+2）/B（k）;M（1）=3*（D（1）-dx0）/H（1）-M（2）/2;M（N+1）=3*（dxn-D（N）/H（N）-M（N）/2;for k=0: 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）;X=-1:Y=1./（25*X.2+1）;dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit（X,Y,

12、dx0,dxn） x1=-1:-0.5;y1=polyval（S（1,:）,x1-X（1）;x2=-0.5:0;y2=polyval（S（2,:）,x2-X（2）;x3=0:0.5; y3=polyval（S（3,:）,x3-X（3）;x4=0.5:y4=polyval（S（4,:）,x4-X（4）;plot（x1,y1,x2,y2,x3,y3,x4,y4, X,Y,结果如图：L8（x）可由公式Ln（x）=得出。三次样条可以利用自然边界条件。写成矩阵：其中j=，i=，dj=6fxj-1,xj,xj+1，n=0=0 d0=dn=0l0（x）=l1（x）= l2（x）= l3（x）= l4（x）=

13、 l5（x）= l6（x）= l7（x）= l8（x）= L8（x）= l1（x）+2 l2（x）+3 l3（x）+4 l4（x）+5 l5（x）+6 l6（x）+7 l7（x）+8 l8（x）求三次样条插值函数由MATLAB计算：可得矩阵形式的线性方程组为：在MATLAB中的Editor中输入程序代码，以下是三次样条函数的程序代码：y=0 1 2 3 4 5 6 7 8;x=0 1 4 9 16 25 36 49 64; n=9t=a解得:M0=0;M1=-0.5209;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=-0.00

14、09;M8=0，则三次样条函数：下面进行画图，在Command Window中输入画图的程序代码：%画图形比较那个插值更精确的函数：x0=0 1 4 9 16 25 36 49 64;y0=0 1 2 3 4 5 6 7 8;x=0:64;y=lagr1（x0,y0,x）;plot（x0,y0,plot（x,y,hold on;pp=csape（x0,y0,fnplt（pp,g:b%axis（0 2 0 1）; %看0 1区间的图形时加上这条指令三次样条插值拉格朗日插值%下面是求拉格朗日插值的函数function y=lagr1（x0,y0,x） n=length（x0）; m=length（

15、x）; for i=1:m z=x（i）; s=0.0; for k=1: p=1.0; if j=k p=p*（z-x0（j）/（x0（k）-x0（j）; s=p*y0（k）+s; y（i）=s;拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三次样条插值曲线，蓝色虚线条为x=y2的曲线，另外一条红色线条为拉格朗日插值曲线。图3-1为0 1的曲线，图3-2为0 64区间上的曲线。图31图32由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线，在0，1朗格朗日插值更精确。而在区间0，64上从图3-2中可以看出蓝色虚线条和绿色线条

16、是几乎重合的，而红色线条在30，70之间有很大的振荡，所在在区间0，64三次样条插值更精确写。5）实验结果 单个多项式高次插值效果并不理想，有龙格现象，偏差大，没有使用价值。而分段低次插值则精确度较高，拟合效果较好，而三次样条插值具有良好的收敛性与稳定性，与分段低次插值相比较光滑度更高，而且提供的信息也相对少一些。我们可以看到，在以上的三道实验题里，我们可以从图形中看出，三次样条的拟合程度是三种插值函数里最好的。6）实验结果分析与小结通过此次实验，我对牛顿多项式插值，三次样条插值，拉格朗日插值有了更进一步的了解，知道了三次样条的拟合程度在高次的情况下更高，在理论上和应用上都有重要意义，在利用计算机编程软件进行高次插值的时候，我们可以多考虑利用三次样条进行插值。