数值分析实验报告1

1、实验报告实验项目名称 插值法 实验室数学实验室 所属课程名称 数值逼近 实 验 类 型 算法设计 实 验 日 期 班 级 学 号 姓 名 成 绩 实验概述：【实验目的及要求】本次实验的目的是熟练数值分析第二章“插值法”的相关内容，掌握三种插值方法：牛顿多项式插值，三次样条插值，拉格朗日插值，并比较三种插值方法的优劣。本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并在MATLAB软件中去实现。【实验原理】数值分析第二章“插值法”的相关内容，包括：牛顿多项式插值，三次样条插值，拉格朗日插值的相应算法和相关性质。【实验环境】（使用的软硬件）软件：MATLAB 2012a硬件：

2、电脑型号：联想 Lenovo 昭阳E46A笔记本电脑操作系统：Windows 8 专业版 处理器：Intel（R）Core（TM）i3 CPU M 350 2.27GHz 2.27GHz实验内容：【实验方案设计】第一步，将书上关于三种插值方法的内容转化成程序语言，用MATLAB实现；第二步，分别用牛顿多项式插值，三次样条插值，拉格朗日插值求解不同的问题。【实验过程】（实验步骤、记录、数据、分析）实验的主要步骤是：首先分析问题，根据分析设计MATLAB程序，利用程序算出问题答案，分析所得答案结果，再得出最后结论。实验一：已知函数在下列各点的值为xi0.20.40.6.0.81.0f（xi）0.9

3、80.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)我们要知道牛顿插值多项式的系数

4、，即均差表中得部分均差。在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:n FF(j,i)=(FF(j

5、,i-1)-FF(j-1,i-1)/(x(j)-x(j-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.

6、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值：三次样条插值函数计算的程序如下：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=

7、1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+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);

8、endb=inv(a)m=b*d'p=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=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.1071 0)T 即M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0则根据三次样条函数定义，可得：S(x)= 接着，在Comm

9、and Window里输入画图的程序代码，下面是画牛顿插值以及三次样条插值图形的程序：x=0.2 0.4 0.6 0.8 1.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

10、=1:1:4y0(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)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')ho

11、ld ongtext('三次样条自然边界')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

12、,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=length(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

13、(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.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(x

14、p,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: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 tem

15、p=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)/2;for k=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)

16、/6; S(k+1,4)=Y(k+1);end当n=10时，我们在Command Window中输入以下的命令：clear,clcX=-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:

17、1;y4=polyval(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,'.')结果如图：当n=20时，我们在Command Window中输入以下的命令：clear,clcX=-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(

18、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,'.')结果如图：实验三：下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间0,64上作图。（1）用这9各点作8次多项式插值L8(x).（2）用三次样条（自然边界条件）程序求S（x）。从结果看在0,64上，那个插值更精确；在区间0,1上，两种哪个更精确？L8(x)可由公式Ln(x)=得出。 三次

19、样条可以利用自然边界条件。写成矩阵：其中j=，i=，dj=6fxj-1,xj,xj+1，n=0=0 d0=dn=0l0(x)=l1(x)= l2(x)= l3(x)= l4(x)= 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中输入程序代码，以下是三次样条函数的程序代码：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶

20、导。y=0 1 2 3 4 5 6 7 8;x=0 1 4 9 16 25 36 49 64; n=9for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+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

21、)=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*d't=ap=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=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=

22、-0.0029;M6=-0.0008;M7=-0.0009;M8=0，则三次样条函数：S(x)= 下面进行画图，在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,'o')hold onplot(x,y,'r');hold on;pp=csape(x0,y0,'variational')fnplt(pp,'g');hold on;plot(x0,y0,':b');hold on%axis(0 2 0 1); %看0 1区间的图形时加上这条指令gtext('三次样条插值')gtext('原图像')gtext('拉格朗日插值'