西北农林科技大学数值分析数值法实验报告.docx

1、西北农林科技大学数值分析数值法实验报告数值法实验报告专业班级：信息与计算科学121 姓名：金辉 学号：20120142801）实验目的本次实验的目的是熟练数值分析第二章“插值法”的相关容，掌握三种插值法：牛顿多项式插值，三次样条插值，拉格朗日插值，并比较三种插值法的优劣。本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并在MATLAB软件中去实现。2）实验题目实验一：已知函数在下列各点的值为xi0.20.40.6.0.81.0f（xi）0.980.920.810.640.38试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值。用图给出（x

2、i，yi），xi=0.2+0.08i，i=0，1, 11, 10，P4（x）及S（x）。实验二：在区间-1,1上分别取用两组等距节点对龙格函数作多项式插值及三次样条插值，对每个值，分别画出插值函数即的图形。实验三：下列数据点的插值x01491625364964y012345678可以得到平根函数的近似，在区间0,64上作图。（1）用这9各点作8次多项式插值L8(x).（2）用三次样条（自然边界条件）程序求S（x）。从结果看在0,64上，那个插值更精确；在区间0,1上，两种哪个更精确？3）实验原理与理论基础数值分析第二章“插值法”的相关容，包括：牛顿多项式插值，三次样条插值，拉格朗日4）实验容实

3、验一：已知函数在下列各点的值为xi0.20.40.6.0.81.0f（xi）0.980.920.810.640.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)+fx0,x1(x-x0)+ fx0,x1,x2(x-x0) (x-x1)+ fx0,x1，xn(x-x0) (x-xn-1)我们要知道牛顿插值多项式的系数，即均差表中得部分均差。在M

4、ATLAB的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:n FF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-

5、i+1); endendfor i=1:n fprintf(%4.2f,x(i); for j=1:i fprintf(%10.5f,FF(i,j); end fprintf(n);end由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.62500-0.52083所以有四次插值牛顿多项式为：P4（x）=0.98-0.3(x-0.2)-0.62500 (

6、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值：三次样条插值函数计算的程序如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。x=0.2 0.4 0.6 0.8 1.0;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);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+

7、h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfor j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=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);for j=1:1:n a(j,j)=2;end r(1)=0; u(n)=0;for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,j+1)=r(j);endb=inv(a)m=b*dp=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:

8、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.1071 0)T 即M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0则根据三次样条函数定义，可得：S(x)= 接着，在Command Window里输入画图的程序代码，下面是画牛顿插值以及三次样条插值图形的程序：x=0.2 0.4 0.6 0.8 1.

9、0;y=0.98 0.92 0.81 0.64 0.38;plot(x,y)hold on for i=1:1:5y(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)endk=0 1 10 11x0=0.2+0.08*kfor i=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4) -

10、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)endplot( x0,y0,o,x0,y0 )hold ony1=spline(x,y,x0)plot(x0,y1,o)hold ons=csape(x,y,variational) fnplt(s,r)hold ongtext(三次样条自然边界)gtext(原图像)gtext(4次牛顿插值)运行上述程序可知：给出的（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10点，S（x）及P4（x）图形

11、如下所示：实验二：在区间-1,1上分别取用两组等距节点对龙格函数作多项式插值及三次样条插值，对每个值，分别画出插值函数即的图形。我们先求多项式插值：在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令（如牛顿插值公式）：function C,D=newpoly(X,Y)n=length(X);D=zeros(n,n) D(:,1)=Y for j=2:n for k=j:n D(k,j)=(D(k,j-1)- D(k-1,j-1)/(X(k)-X(k-j+1); endendC=D(n,n);for k=(n-1):-1:1 C=conv(C,poly(X(k) m=le

12、ngth(C); C(m)= C(m)+D(k,k);end当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:1; y=polyval(C,x); plot(x,y,X,Y,.); grid on; xp=-1:0.2:1; z=1./(1+25*xp.2); plot(xp,z,r)得到插值函数和f（x）图形：当n=20时，我们在Command Window中输入以下的命令：clear,clf,hold on; X=-1:0.

13、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,.); grid on; xp=-1:0.1:1; z=1./(1+25*xp.2); plot(xp,z,r)得到插值函数和f（x）图形：下面再求三次样条插值函数，在MATLAB的Editor中建立一个多项式的M-file,输入下列程序代码：function 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

14、: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=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/