数值分析实验指导书

1、数值分析实 验 指 导 书潍坊学院数学与信息科学学院2012年04月I / 28II / 28目 录目录I实验一 插值与曲线拟合的最小二乘法1实验二 数值积分4实验三 解线性方程组的直接法9实验四 解线性方程组的迭代法11实验五 非线性方程的数值解法13实验六 常微分方程数值解法17I / 28实验一 插值与曲线拟合的最小二乘法一、实验目的：1了解拉格朗日插值法、牛顿插值法、曲线拟合最小二乘法的基本原理和方法；2掌握拉格朗日插值多项式牛顿插值多项式的用法；3掌握最小二乘原理，会求拟合函数及超定方程组的最小二乘解。二、实验内容：1用拉格朗日插值公式和牛顿插值公式确定函数值；2对函数f (x)进行

2、拉格朗日插值和牛顿插值；3利用Polyfit拟合幂函数，利用Polyfit拟合多项式。三、实验过程：1给定函数四个点的数据如下:，试用插值公式确定函数在处的函数值。MATLAB程序如下：X=1.1,2.3,3.9,5.1; Y =3.877,4.726,4.651 ,2.117;p1=poly(X(1); p2=poly(X(2);p3=poly(X(3); p4=poly(X(4); l01= conv ( conv (p2, p3), p4)/( X(1)- X(2)* ( X(1)- X(3) * ( X(1)- X(4), l11= conv ( conv (p1, p3), p4)/

3、( X(2)- X(1)* ( X(2)- X(3) * ( X(2)- X(4),l21= conv ( conv (p1, p2), p4)/( X(3)- X(1)* ( X(3)- X(2) * ( X(3)- X(4),l31= conv ( conv (p1, p2), p3)/( X(4)- X(1)* ( X(4)- X(2) * ( X(4)- X(3),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31),P = l01* Y(1)+ l11* Y(2) + l21* Y(3) +

4、l31* Y(4),运行后输出的基函数l0，l1，l2和l3为1 / 28l0 =-1/24*x3+1/8*x2-1/12*x，l1 =1/4*x3-1/4*x2-x+1l2 =-1/3*x3+4/3*x，l3 =1/8*x3+1/8*x2-1/4*x输入程序>> L=poly2sym (P),x=2.101; Y = polyval(P,x)运行后输出插值多项式和插值为L=-629/5376*x3+31433/53760*x2-63029765850741/281474976710656*x+2010616283501353/562949953421312Y =4.5969输入程

5、序>> L=poly2sym (P),x=4.234; Y = polyval(P,x)运行后输出插值多项式和插值为L=-629/5376*x3+31433/53760*x2-63029765850741/281474976710656*x+2010616283501353/562949953421312Y =4.2244L=145616387951645/9007199254740992*x3-2517512191700115/4503599627370496*x2+14477/6000*x+2007/1000Y = 3.4290输入程序>> syms M; x=2.

6、101; R3=M*abs(x-X(1)*(x-X(2) *(x-X(3) *(x-X(4)/24运行后输出误差限为R3 =435065974692861/36028797018963968*M2在区间上取结点数，等距间隔的节点为插值点，对于函数进行拉格朗日插值。MATLAB程序如下t=-5: 1:5;ft=(1+t.*t).5;t1=-5:1:5;ft1=(1+t1.*t1).5;2 / 28y1=Lagran(t1,ft1,t);plot(t,ft,'b:',t,y1,'g+');xlabel('x');ylabel('y')

7、;对一组数据做拉格朗日的M文件如下Lagranmfunction fi =Lagran(x,f,xi)fi=zeros(size(xi);npl=length(f);for i=1:nplz=ones(size(xi);for j=1: npl if i=j z=z.*(xi-x(j)/(x(i)-x(j); endendfi=fi+z*f(i);end3给出节点数据，作三阶牛顿插值多项式，计算，并估计其误差。MATLAB程序如下syms M,X=-4,0,1,2; Y =27,1,2,17; x=-2.345; y,R,A,C,P=newdscg(X,Y,x,M)运行后输出y = 22.32

8、11R =1323077530165133/562949953421312*M（即R =2.3503*M）A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167C =0.9167 4.2500 -4.1667 1.0000P =11/12*x3+17/4*x2-25/6*x+13 / 28其中求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序如下：function y,R,A,C,L=newdscg(X,Y,x,M)n=length(X); m=length(x)

9、;for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y's=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k);d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q

10、1/c1;L(k,:)=poly2sym(C);4给定数据如下 试用幂函数拟合以上数据。Matlab程序如下x=0.15,0.4,0.6,1.01,1.5,2.2,2.4,2.7,2.9,3.5,3.8,4.4,4.6,5.1,6.6,7.6y=4.4964,5.1284,5.6931,6.2884,7.0989,7.5507,7.5106,8.0756,7.8708,8.2403,8.5303,8.7394,8.9981,9.1450,9.5070,9.9115c=polyfit(x, y,1)5用以下数据求二次多项式的系数，并做出拟合曲线：4 / 28 Matlab程序如下x=0.1,0.

11、4,0.5,0.7,0.7,0.9;y=0.61,0.92,0.99,1.52,1.47,2.03;cc=polyfit(x,y,2)xx=x(1):0.1:x(length(x);yy=polyval(cc,xx);plot(xx,yy);hold onplot(x,y,'x')axis(0,1,0,3)xlabel('x')ylabel('y')四、实验要求：1学会Lagrange插值方法和牛顿插值方法并应用上述方法于实际问题；2将拟合的结果与拉格朗日插值的结果比较。3归纳总结数值实验结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中

12、选择方法应注意的问题。5 / 28实验二 数值积分一、实验目的：1理解数值积分的概念，掌握各种数值积分方法，包括梯形公式、抛物线公式等，通过实际计算体会各种数值积分方法的精确度；2掌握复化梯形公式、复化Simpson公式及其截断误差的分析；3了解Romberg公式及Romberg积分法。二、实验内容：1利用低阶求积公式求积分的近似值；利用复化梯形公式和复化Simpson公式计算定积分；2利用复化梯形公式计算无穷限广义积分；3利用Romberg积分法计算定积分。三、实验过程：1用求截断误差公式的MATLAB主程序，求计算定积分d的近似值的阶牛顿科茨公式的截断误差公式。Matlab程序如下n=1,

13、 RNC1=ncE(n), n=2, RNC2=ncE(n), n=3, RNC3=ncE(n)n=4, RNC4=ncE(n), n=8, RNC8=ncE(n)运行后屏幕显示结果如下n= 1 RNC1 =1/12*(b-a)3*fxn1n = 2 RNC2 =1/90*(1/2*b-1/2*a)5*fxn2n = 3 RNC3=3/80*(1/3*b-1/3*a)5*fxn1n = 4 RNC4=8/945*(1/4*b-1/4*a)7*fxn2n = 8 RNC8=35065906189543/6926923254988800*(1/8*b-1/8*a)11*fxn2其中计算n阶牛顿-科

14、茨的公式的截断误差公式的MATLAB主程序function RNC=ncE(n)6 / 28suk=1; p=n/2-fix(n/2);if p=0for k=1:n+2suk=suk*k;endsuk; syms t a b fxn2,su=t2; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNC= (b-a)/n)(n+3)*fxn2*abs(y)/ suk;elsefor k=1:n+1suk=suk*k;endsuk; syms t a b fxn1,su=t; for u=1:nsu=su*(t-u)

15、;endsu; intf=int(su,t,0,n); y=double(intf);RNC= (b-a)/n)(n+2)*fxn1*abs(y)/ suk;end2分别利用复化梯形公式和复化Simpson公式计算定积分 取，精确值为。复化梯形公式MATLAB程序如下clear;Iexact=4.006994;a=0;b=2;fprintf('n n I Errorn');n=1; for k=1:4n=2*n;7 / 28h=(b-a)/n;i=1:n+1;x=a+(i-1)*h;f=sqrt(1+exp(x);I=travez_v(f,h);fprintf('%3.

16、0f %10.5f %10.5fn',n,I,Iexact-I);end其中复化梯形求积法的M文件如下trapez_v.mfunction I=trapez_v(f,h)I=h*sum(f)-(f(1)+f(length(f)/2;复化Simpson求积公式MATLAB程序如下：clear;Iexact=4.006994;a=0;b=2;fprintf('n n I Errorn');n=1;for k=1:4,n=2*nh=(b-a)/n;i=1:n+1;x=a+(i-1)*h;f=sqrt(1+exp(x); I=(h/3)*(f(1)+4*sum(f(2:2:n)

17、+f(n+1); if n>2 I=I+(h/3)*2*sum(f(3:2:n); end fprintf('%3.0f %10.5f %10.5fn',n,I,Iexact-I);end3计算无穷限广义积分，其精确值。用-10,10替换积分限，则有8 / 28利用复化梯形公式，调用trapez_v.m文件，取n=10,20的MATLAB程序如下：I=1;a=-10;b=10;fprintf('n h n I1n');n=0;for k=1:2;n=n+20;h=(b-a)/n;i=1:n+1;x=a+(i-1)*h;f=exp(-x.2);I=trape

18、z_v(f,h);I1=1/sqrt(pi)*I;fprintf('%3.1f %10.0f %10.8fn',h,n,I1);end4取精度为，分别用和作为计算停止的条件，用Romberg程序计算d，取精度为,并与精确值比较。然后取精度为，用作为计算停止的条件，观察用龙贝格求积公式计算的结果与精确值的绝对误差是否满足。MATLAB程序如下F=inline('1./(1+x)'); RT,R,wugu,h=romberg(F,0,1.5,1.e-8,13)syms x fi=int(1/(1+x),x,0,1.5); Fs=double(fi), wR=doub

19、le(abs(fi-R), wR1= wR - wugu其中龙贝格积分的MATLAB程序function RT,R,wugu,h=romberg(fun,a,b, wucha,m)n=1;h=b-a; wugu=1; x=a;k=0; RT=zeros(4,4); RT(1,1)=h*(feval(fun,a)+feval(fun,b)/2;while(wugu>wucha)&(k<m)|(k<4) k=k+1; h=h/2; s=0; for j=1:n x=a+h*(2*j-1); s=s+feval(fun,x);end9 / 28RT(k+1,1)= RT(k

20、,1)/2+h*s; n=2*n;for i=1:kRT(k+1,i+1)=(4i)*RT(k+1,i)-RT(k,i)/(4i-1);endwugu=abs(RT(k+1,k)-RT(k+1,k+1);endR=RT(k+1,k+1);四、实验要求：1写出梯形公式、Simpson公式、柯特斯公式求数值积分的算法；2写出复化梯形公式和复化Simpson公式Romberg方法求数值积分的算法；3进一步加深对数值积分的理解。10 / 28实验三 解线性方程组的直接法一、实验目的：1了解求线性方程组的直接法的有关理论和方法；2会编制Gauss顺序消去法、列主元消去法的程序；3熟悉Gauss顺序消去法

21、、列主元消去法求解线性方程组的过程.。二、实验内容： 用Gauss顺序消去法、列主元消去法求线性方程组的解。三、实验过程：1用顺序Gauss消去法求解线性方程组：Gauss消去法的MATLAB程序如下：a=-0.04 0.04 0.12 3; 0.56 -1.56 0.32 1; -0.24 1.24 -0.28 0;x=0,0,0;temp=a(2,:); a(2,:) = a(1,:); a(1,:)=temp;aa(2,:)=a(2,:)- a(1,:)*a(2,1)/a(1,1);a(3,:)=a(3,:)- a(1,:)*a(3,1)/a(1,1);atemp=a(3,:); a(3

22、,:) = a(2,:); a(2,:)=temp;aa(3,:)=a(3,:)- a(2,:)*a(3,2)/a(2,2);ax(3)=a(3,4)/a(3,3);x(2)=(a(2,4)-a(2,3)*x(3)/a(2,2);x(1)=(a(1,4)-a(1,2)*x(2)-a(1,3)*x(3)/a(1,1);x2用列主元消元法解线性方程组的MATLAB程序解方程组列主元消元法的MATLAB程序如下：A=0 -1 -1 1;1 -1 1 -3;2 -2 -4 6;1 -2 -4 1; 11 / 28b=0;1;-1;-1; RA,RB,n,X=liezhu(A,b)运行后输出结果RA =

23、 4，RB = 4，n = 4，X =0 -0.5 0.5 0其中列主元消元法解线性方程组的MATLAB主程序如下function RA,RB,n,X=liezhu(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if zhica>0,dispreturnendif RA=RB if RA=ndisp X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for

24、k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse dispendend四、实验要求：1写出高斯顺序消去法、列主元消去法解线性方程组的算法，编写程序上机调试出结果；13 / 282进一步加深对高斯消去法的理解。13 / 28实验四 解线性方程组的迭代法一、实验目的：1熟悉解线性方程组的迭

25、代法的有关理论和方法；2掌握用迭代法求解线性方程组的基本思想和计算步骤；3熟悉逐次超松弛迭代法求解线性方程组的过程；4注意所用方法的收敛性及其收敛速度问题。二、实验内容：用高斯-塞德尔迭代法、逐次超松弛迭代法求线性方程组的解。三、实验过程：1用高斯-塞德尔迭代法解下列线性方程组,取初始值，要求当时，迭代终止。（a） （b）用高斯-塞德尔迭代定义解线性方程组Ax=b的MATLAB主程序如下：function X=gsdddy(A,b,X0,P,wucha,max1)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); dD=det(D);if dD=0dispel

26、sedispiD=inv(D-L); B2=iD*U;f2=iD*b;jX=Ab; X=X0; n m=size(A);for k=1:max1X1= B2*X+f2; djwcX=norm(X1-X,P);xdwcX=djwcX/(norm(X,P)+eps);if (djwcX<wucha)|(xdwcX<wucha) return else k,X1,k=k+1;X=X1;end14 / 28endif (djwcX<wucha)|(xdwcX<wucha) dispelsedispX=X;jX=jX'endendX=X;D,U,L,jX=jX (a)的M

27、ATLAB程序如下：A=10 -1 -2;-1 10 -2;-1 -1 0.5; b=7.2;8.3;4.2; X0=0;0;0;X=gsdddy(A,b,X0,inf, 0.001,100)(b)的MATLAB程序如下：A=3 4 -5 7;2 -8 3 -2;4 51 -13 16;7 -2 21 3;b=5;2;-1;21;X0=0;0;0;0;X=gsdddy(A,b,X0,inf,0.001,100)2用逐次超松弛迭代法解由一电路得到的方程组解的精度为0.001。逐次超松弛法的MATLAB程序如下：a(1,1)=1/2+1/4+1/3;a(1,2)=-1/4;a(1,3)=-1/3;

28、a(2,1)=a(1,2);a(2,2)=1/4+1/3+1/5;a(2,3)=-1/5;a(3,1)=a(1,3);a(3,2)=a(2,3);a(3,3)=1/3+1/3+1/5;y(1)=20/2;y(2)=0;y(3)=5/3;x=zeros(1,3);w=1.2;for it=1:50 error=0; for i=1:315 / 28 s=0;xb=x(i); for j=1:3 if i=j,s=s+a(i,j)*x(j);end end x(i)=w*(y(i)-s)/a(i,i)+(1-w)*x(i); error=error+abs(x(i)-xb); end fprint

29、f('it.no.=%3.0f,error=%7.2en',it,error) if error/3<0.001,break;endendx四、实验要求：1掌握常用的几种迭代方法；2掌握迭代法收敛性及误差估计；3学会用高斯-塞德尔迭代法、逐次超松弛迭代法求解线性方程组。16 / 28实验五 非线性方程的数值解法一、实验目的：1掌握非线性方程的各种解法，包括迭代法、牛顿法，并通过编程练习与上机运算，体会迭代法、牛顿法的不同特点；2掌握解非线性方程的弦截法，并与牛顿法作比较；3了解各种方法的收敛性。二、实验内容：用二分法、牛顿法、弦截法求非线性方程的根。三、实验过程：1用迭代

30、法求方程的一个正根。迭代法的MATLAB程序function k,piancha,xdpiancha,xk=diedai1(x0,k)x(1)=x0; for i=1:k x(i+1)=fun1(x(i); piancha= abs(x(i+1)-x(i); xdpiancha=piancha/( abs(x(i+1)+eps); i=i+1;xk=x(i);(i-1) piancha xdpiancha xkendif (piancha >1)&(xdpiancha>0.5)&(k>3) disp return; end if (piancha < 0

31、.001)&(xdpiancha< 0.0000005)&(k>3) disp return; endp=(i-1) piancha xdpiancha xk;用不同迭代法求解的MATLAB程序如下：function m=fun1 (x)m=(10-x*x)/2;17 / 28k,piancha,xdpiancha,xk= diedai1(2,5)运行后输出用迭代公式的结果k,piancha,xdpiancha,xk=1.00000000000000 1.00000000000000 0.33333333333333 3.000000000000002.000000

32、00000000 2.50000000000000 5.00000000000000 0.50000000000000 3.00000000000000 4.37500000000000 0.89743589743590 4.87500000000000 4.00000000000000 11.75781250000000 1.70828603859251 -6.882812500000005.00000000000000 11.80374145507813 0.63167031671297 -18.68655395507813请用户注意：此迭代序列发散k=5，piancha = 11.803

33、74145507813，xdpiancha = 0.63167031671297，xk = -18.68655395507813由以上运行后输出的迭代序列与根=2.316 624 790 355 40相差越来越大，即迭代序列发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于0.5。请重新输入新的迭代公式function m=fun1 (x)m=10/(x+2);k,piancha,xdpiancha,xk= diedai1(2,5) 用迭代公式运行后输出的结果k,piancha,xdpiancha,xk=1.00000000000000 0.5000

34、0000000000 0.20000000000000 2.50000000000000 2.00000000000000 0.27777777777778 0.12500000000000 2.222222222222223.00000000000000 0.14619883040936 0.06172839506173 2.36842105263158 4.00000000000000 0.07926442612555 0.03462603878116 2.28915662650602 5.00000000000000 0.04230404765128 0.01814486863115 2

35、.33146067415730k=5,piancha =0.04230404765128,xdpiancha = 0.01814486863115,xk = 2.33146067415730可见，偏差piancha和偏差的相对误差xdpiancha的值逐渐变小，且第5次的迭代值xk =2.331 460 674 157 30与根=2.316 624 790 355 40接近，则迭代序列19 / 28收敛，但收敛速度较慢，此迭代法较为成功。请重新输入新的迭代公式function m=fun1 (x)m=x-(x*x+2*x-10)/(2*x+2);k,piancha,xdpiancha,xk=

36、diedai1(2,5)用迭代公式运行后输出的结果k,piancha,xdpiancha,xk=1.00000000000000 0.33333333333333 0.14285714285714 2.333333333333332.00000000000000 0.01666666666667 0.00719424460432 2.31666666666667 3.00000000000000 0.00004187604690 0.00001807631822 2.31662479061977 4.00000000000000 0.00000000026437 0.0000000001141

37、2 2.316624790355405.00000000000000 0 0 2.31662479035540k = 5; piancha = 0; xdpiancha = 0; y = 2.31662479035540.可见，偏差piancha和偏差的相对误差xdpiancha的值越来越小，且第5次的迭代值xk=2.316 624 790 355 40与根=2.316 624 790 355 40相差无几，则迭代序列收敛，且收敛速度很快，此迭代法成功。2取初值，用牛顿法解非线性方程，要求精度为0.00001.牛顿法的MATLAB程序如下function x=Newt_n(f_name,x0)

38、x0=10;x=x0,xb=x-999;n=0;del_x=0.01;while abs(x-xb)>0.00001 n=n+1;xb=x; if n>300,break;end y=feval('f_name',x);y_driv=(feval('f_name',x+del_x)-y)/del_x;19 / 28 x=xb-y/y_driv; fprintf('n=%3.0f,x=%12.5e,y=%12.5e',n,x,y) fprintf('yd=%12.5en',y_driv)endfprintf('n

39、 final answer=%12.5en',x);其中调用的函数f_name(x)定义如下function y=f_name(x)y=x2-155;3用弦截法解由一物理问题抽象出的物理方程 弦截法的MATLAB程序如下function y=zlj(x0,x1)x2=x1-fc(x1)*(x1-x0)/(fc(x1)-fc(x0);n=1;while(abs(x1-x0)>1.0e-4)&(n<=10000000) x0=x1;x1=x2; x2=x1-fc(x1)*(x1-x0)/(fc(x1)-fc(x0); n=n+1; fprintf('n=%3.0

40、fn',n-1) fprintf('fv=%12.5en',x0)endnx2四、实验要求：掌握解非线性问题的二分法、牛顿法、弦截法，并会用牛顿法、弦截法解非线性方程。20 / 28实验六 常微分方程数值解法一、实验目的：1熟悉求解常微分方程初值问题的有关方法和理论；2熟悉显式欧拉公式、二龙格-库塔公式和四阶龙格-库塔公式的使用；3会编制上述方法的计算程序；4通过对各种求解方法的计算实习，体会各种解法的功能、优缺点及适用场合，会选取适当的求解方法。二、实验内容：1列出常微分方程并利用显式欧拉公式求其数值解；2利用二龙格-库塔公式和四阶龙格-库塔公式求微分方程的数值解。三、实验过程：1某跳伞者在t=0时刻从飞机上跳出，假设初始时刻的垂直速度