《西南交大数值分析上机作业.doc》由会员分享，可在线阅读，更多相关《西南交大数值分析上机作业.doc（23页珍藏版）》请在三一文库上搜索。
1、数值分析上机作业姓 名：XXX 学 号：XXX 专 业：XXX 联系电话：XXX 数值分析上机作业 II序言长期以来，在数值分析教学中，数值分析作为一门纯粹的理论课程进行授课，但结果便忽略数值分析上机的试验环节，忽略了将编程语言与教材中给出的算法结合起来，以解决实际问题。而对于C语言和Fortran语言，只有基本熟练掌握后才具有一定得编程能力，而且此两种语言编程效率较低，不适于数值分析编程。本次数值分析上机实习所用的语言是MATLAB语言。MATLAB语言是著名的矩阵计算数学软件，除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。MATLAB语
2、言主要特点：1）高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；2）具有完备的图形处理功能，实现计算结果和编程的可视化；3）友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；4）功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。编程环境:MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。简单易用：MATLAB是一个高级的矩阵/阵列语言，
3、它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。而且这种语言可移植性好、可拓展性极强。强处理能力：MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。函数中所使用的算法都是科研和工程计算中的最新研究成果，而且经过了各种优化和容错处理。图形处理：MATLAB自产生之日起就具有方便的数据可视化功能，以将向量和矩阵用图形表现出来，并且可以对图形进行标注和打印。高层次的作图包括二维和三维的可视化、图象处理、动画
4、和表达式作图。可用于科学计算和工程绘图。程序接口：MATLAB可以利用MATLAB编译器和C/C+数学库和图形库，将自己的MATLAB程序自动转换为独立于MATLAB运行的C和C+代码。工具箱是MATLAB函数的子程序库，每一个工具箱都是为某一类学科专业和应用而定制的，主要包括信号处理、控制系统、神经网络、模糊逻辑、小波分析和系统仿真等方面的应用。此数值分析上机实习报告包含两道大题，其中第一题包含两个小题，为必做题，第二题为选作题，从五个题中选做两个。每道题分为四个部分，分别为题目分析、编写程序、计算结果和结果分析，对每个题的不同情况进行了一定的分析，做到理论与实际结合，以解决实际问题。 数值
5、分析上机作业 III目录序言I1.必做题一11.1 Gauss消元法11.1.1 第一小题11.1.2 第二小题21.1.3 第三小题21.2高斯塞德尔21.2.1 第一小题21.2.2 第二小题31.2.3 第三小题31.3 计算结果对比分析42. 必做题二52.1 步长为0.0552.2 步长为0.162.3 步长为0.272.4 计算结果对比分析73. 选做题题三93.1 选取初值为293.2 选取初值为393.3 选取初值为493.4 结果分析104. 选做题题四114.1 复化梯形公式114.2 复化Simpson公式114.3 Romberg算法114.4 结果分析11总结12附录
6、13 数值分析上机作业 191.必做题一写出对一般的线性方程组通用的Gauss消元, Gauss-Seidel迭代程序。并以下面的线性方程组为例进行计算，讨论所得到的计算结果是否与理论一致。（1）（2）（3）1.1 Gauss消元法1.1.1 第一小题1.1.1.1右端项为b1A=6,2,-1;1,4,-2;-3,1,4;b=-3,2,4;x=Gauss(A,b)计算结果： x = -0.7273 0.8081 0.25251.1.1.2右端项为b2A=6,2,-1;1,4,-2;-3,1,4;b=100,-200,345;x=Gauss(A,b)计算结果： x = 36.3636 -2.07
7、07 114.04041.1.2 第二小题1.1.2.1右端项为b1A=1,0.8,0.8;0.8,1,0.8;0.8,0.8,1;b=3,2,1;x=Gauss(A,b)计算结果：x = 5.7692 0.7692 -4.23081.1.2.2右端项为b2A=1,0.8,0.8;0.8,1,0.8;0.8,0.8,1;b=5,0,-10;x=Gauss(A,b)计算结果： x = 32.6923 7.6923 -42.30771.1.3第三小题A=1,3;-7,1;b=4,6;x=Gauss(A,b)计算结果： x = -0.6364 1.54551.2高斯塞德尔1.2.1第一小题1.3.1
8、.1右端项为b1A=6,2,-1;1,4,-2;-3,1,4;b=-3,2,4;x=Gauss_seidel(A,b)计算结果：x = -0.7273 0.8081 0.25251.2.1.2右端项为b2A=6,2,-1;1,4,-2;-3,1,4;b=100,-200,345;x=Gauss_seidel(A,b)计算结果： x = 36.3636 -2.0707 114.04041.2.2第二小题1.3.2.1右端项为b1A=1,0.8,0.8;0.8,1,0.8;0.8,0.8,1;b=3,2,1;x=Gauss_seidel(A,b)计算结果 x = 5.7692 0.7692 -4.
9、23071.2.2.2右端项为b2A=1,0.8,0.8;0.8,1,0.8;0.8,0.8,1;b=5,0,-10;x=Gauss_seidel(A,b)计算结果： x = 32.6923 7.6923 -42.30771.2.3第三小题A=1,3;-7,1;b=4,6;x=Gauss_seidel(A,b)计算结果： x = 0.175310131 1.22681.3 计算结果对比分析表1.1 Gauss消元法与Gauss-Seidel迭代法对比计算结果方法（1）（2）（3）b1b2b1b2bGauss消元x1-0.727336.36365.769232.6923-0.6364x20.80
10、81-2.07070.76927.69231.5455x30.2525114.0404-4.2308-42.3077Gauss-Seidelx1-0.727336.36365.769232.69230.175310131x20.8081-2.07070.76927.69231.2268x30.2525114.0404-4.2307-42.3077表1.2 高斯赛德尔迭代法敛散性分析题目序列号特征值谱半径敛散性（1）、（2）1= 02= 0.0469 + 0.3504i3= 0.0469 - 0.3504i（B）=0.3536（B）1收敛（3）、（4）1=02= 0.7040 + 0.1280i
11、3= 0.7040 - 0.1280i（B）= 0.7155（B）1发散小结：从表1.2可知，对矩阵方程组，利用高斯赛德尔迭代法计算，当矩阵的谱半径（B）1时，结果发散。而且（B）值越小，收敛速度越快。计算结果收敛与否只取决于迭代矩阵的谱半径，与初始量及方程组的右端项无关，只与迭代矩阵有关。前两个小题高斯-塞德尔迭代方法在b不同时其计算结果与理论值很接近，但第三小题结果却与理论值有很大差异，原因是此种情况下高斯-塞德尔迭代法不收敛。通过计算可知： =Max=211所以此时高斯-塞德尔迭代算法计算结果与理论值差异很大。2.必做题二给定初值问题 （精确解为），用Runge-Kutta 4阶算法按步
12、长求解，分析其中遇到的现象及问题。2.1 步长为0.05fun=(x,y)-20*y;x,y=runge_kutta(fun,0,2,1,0.05)计算结果：N =40x = Columns 1 through 6 0 0.0500 0.1000 0.1500 0.2000 0.2500 Columns 7 through 12 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 Columns 13 through 18 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 Columns 19 through 24 0.9000
13、 0.9500 1.0000 1.0500 1.1000 1.1500 Columns 25 through 30 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 Columns 31 through 36 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 Columns 37 through 41 1.8000 1.8500 1.9000 1.9500 2.0000y = Columns 1 through 6 1.0000 0.3750 0.1406 0.0527 0.0198 0.0074 Columns 7 thro
14、ugh 12 0.0028 0.0010 0.0004 0.0001 0.0001 0.0000 Columns 13 through 18 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 19 through 24 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 25 through 30 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 31 through 36 0.0000 0.0000 0.0000 0.0000 0.0000 0.
15、0000 Columns 37 through 41 0.0000 0.0000 0.0000 0.0000 0.00002.2 步长为0.1fun=(x,y)-20*y;x,y=runge_kutta(fun,0,2,1,0.1)计算结果：N =20x = Columns 1 through 6 0 0.1000 0.2000 0.3000 0.4000 0.5000 Columns 7 through 12 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 Columns 13 through 18 1.2000 1.3000 1.4000 1.5000
16、1.6000 1.7000 Columns 19 through 21 1.8000 1.9000 2.0000y = Columns 1 through 6 1.0000 0.3333 0.1111 0.0370 0.0123 0.0041 Columns 7 through 12 0.0014 0.0005 0.0002 0.0001 0.0000 0.0000 Columns 13 through 18 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 19 through 21 0.0000 0.0000 0.00002.3 步长为0.
17、2fun=(x,y)-20*y;x,y=runge_kutta(fun,0,2,1,0.2)计算结果：N =10x = Columns 1 through 6 0 0.2000 0.4000 0.6000 0.8000 1.0000 Columns 7 through 11 1.2000 1.4000 1.6000 1.8000 2.0000y = Columns 1 through 5 1 5 25 125 625 Columns 6 through 10 3125 15625 78125 390625 1953125 Column 11 97656252.4计算结果对比分析表2.1 步长为
18、0.05时的Runge-Kutta四阶算法K001.00001010.05000.37500.3678790.0071220.10000.14060.1353350.0052630.15000.05270.0497870.0029140.20000.01980.0183160.0014850.25000.00740.0067380.0006660.30000.00280.0024790.0003270.35000.00100.0009128.810-580.40000.00040.0003356.510-590.45000.00010.0001232.3410-5100.50000.00014
19、.5410-55.510-5110.55000.00001.6710-51.6710-5表2.2 步长为0.1时的Runge-Kutta四阶算法K00 1.00001010.10000.33330.1353350.19796520.20000.11110.0183160.09278430.30000.03700.0024790.03452140.40000.01230.0003350.01196550.50000.00414.5410-50.00405560.60000.00146.1410-60.00139470.70000.00058.3210-70.00049980.80000.0002
20、1.1310-70.000290.90000.00011.5210-8110-4101.00000.00002.0610-92.0610-9表2.3 步长为0.2时的Runge-Kutta四阶算法K0011010.200050.0183164.98168420.4000250.00033524.9996630.60001256.1410-612540.80006251.1310-762551.000031252.0610-9312561.2000156253.7810-111562571.4000781256.9110-137812581.60003906251.2710-1439062591
21、.800019531252.3210-161953125102.000097656254.2510-189765625小结：（1）用标准的四阶Runge-Kutta算法计算常微分方程时，不同的步长对算法的稳定性有影响。不够稳定的步长下面的计算，误差会越来越大，结果失真严重。（2）一般情况下，步长越小，标准的四阶Runge-Kutta算法的稳定性越高，精度也越高。3. 选做题题三利用牛顿法求方程于区间的根，考虑不同初值下牛顿法的收敛情况。3.1 选取初值为2计算结果： 3.146193220620583 8 1.918995658978622 1.837015121183504 3.486711
22、770340445 3.097944090589694 3.145068876896332 3.146197295211114 3.146193220281240 3.146193220620583Elapsed time is 3.721679 seconds.3.2 选取初值为3计算结果：3.146193220620600 6 2.982091433074349 2.962152133230459 3.148673175461638 3.146157564897870 3.146193214077455 3.146193220620600Elapsed time is 3.250174 s
23、econds.3.3 选取初值为4计算结果：3.146193220620586 7 4.155299956602238 4.348876621999810 3.199581233044784 3.149700130114401 3.146206862616837 3.146193224159589 3.146193220620586Elapsed time is 3.548942 seconds.3.4 结果分析表3 不同初值下牛顿法的收敛情况234计算结果3.1461932206205833.1461932206206003.146193220620586迭代次数867所用时间（秒）3.72
24、16793.2501743.548942迭代过程n11.9189956589786222.9820914330743494.15529995660223821.8370151211835042.9621521332304594.34887662199981033.4867117703404453.1486731754616383.19958123304478443.0979440905896943.1461575648978703.14970013011440153.1450688768963323.1461932140774553.14620686261683763.146197295211
25、1143.1461932206206003.14619322415958973.1461932202812403.14619322062058683.146193220620583由上表可知，当初值离根越近，牛顿法收敛的速度越快，迭代次数越少，计算机的计算时间越少；相反，初值离根越远，收敛速度越慢，迭代次数越多，计算机的计算时间越多，占用的内存也就越高。4. 选做题题四已知，利用复化梯形公式、复化Simpson公式和Romberg算法求的近似值；并观察实际计算结果，比较它们的收敛速度。4.1 复化梯形公式计算结果：I =3.1399259889071594.2 复化Simpson公式计算结果：
26、I =3.1415926139392154.3 Romberg算法计算结果：I =3.1415926652777174.4 结果分析的精确值为3.1415926535897932384根据以三种的计算结果可知，收敛速度最快的是Romberg算法，复化Simpson算法次之，复化梯形算法最慢，原因如下表所示。表4 Romberg算法、复化Simpson算法、复化梯形算法比较求积公式代数精度误差量级Tn(f)1O(h2)Sn(f)3O(h4)Rn(f)7O(h8)其中上表中h为步长，O(hm)表示同阶无穷小。因此它们的收敛速度依次为Romberg算法复化Simpson算法复化梯形算法.总结通过本次
27、的编程作业，我收获很多。逐渐学会了MATLAB编程语言的初步使用，在编程过程中加强了对高斯-塞德尔方法、初值问题、数值积分收敛速度等方面的理解，更加培养了对于数值分析学习的乐趣，对我今后的学习工作都将会产生积极的影响。编程要重视最基本的原理，即这些方法的最基本的步骤，写出合理的算法，能够有效地提高编程速度。同一个题目，可能有程序方法实现，各种方法都有各自的有点，要合理选择合适的方法。在大型计算中，计算速度非常重要，因此，我们要想办法优化程序，提高运速度。对于实际的工程而言，比如土木、机械行业，在建立合适的力学模型后，绝大多数情况下，解析解难以求出，运用编程的方法求出数值解，对于解决工程实际问题
28、，意义重大。因此，我们要重视课本上的基本知识，基本迭代原理，加强编程训练，以利于以后的工作学习中实际问题的解决。附录程序一1.Gauss法调用M文件程序function x=Gauss(A,b)%线性方程组的高斯消去法%A为系数矩阵，b为方程组的右端项n,m=size(A);nb=length(b);if n=merror(系数矩阵必须是方的);endif m=nberror(b的维数与方程组的行数不匹配);endfor k=1:n-1a_max=0;for i=k:nif abs(A(i,k)a_maxa_max=abs(A(i,k);r=i;endendif a_maxkfor j=k:n
29、;z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%消元过程for i=k+1:nm=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endend%回带过程if abs(A(n,n)1e-15error(系数矩阵奇异，无法求解方程组);endx=zeros(size(b);for k=n:-1:1for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end2. Gauss-seid
30、elfunction x=Gauss_seidel(A,b)%线性方程组的Gauss-seidel迭代法%A为系数矩阵，b为方程组的右端项n,m=size(A);nb=length(b);if n=merror(系数矩阵必须是方的);endif m=nberror(b的维数与方程组的行数不匹配);endx=zeros(nb,1);k=1;B=-tril(A)triu(A,1);f=tril(A)b;while k100y=B*x+f;if norm(y-x,inf)1e-5break;endx=y;k=k+1;end程序二Runge-Kutta调用M文件程序function x,y=runge
31、_kutta(fun,a,b,y0,h)%四阶runge_kutta算法%fun一阶微分方程的函数%a,b为区间端点；y0为初始值%h为步长，N为区间的等分数N=(b-a)/hx=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=y0;for n=1:Nx(n+1)=x(n)+h;k1=h*feval(fun,x(n),y(n);k2=h*feval(fun,x(n)+1/2*h,y(n)+1/2*k1);k3=h*feval(fun,x(n)+1/2*h,y(n)+1/2*k2);k4=h*feval(fun,x(n)+h,y(n)+k3);y(n+1)=y(
32、n)+1/6*(k1+2*k2+2*k3+k4);end程序三牛顿法调用M文件程序function x,iter,X=newton(fun,x0,eps,maxiter)%x非线性方程的近似根%iter迭代次数%X每步迭代的结果%fun输入函数%x0初始迭代点%eps计算精度，默认值为1e-6%maxiter最大迭代次数，默认值为1e4tic %保存当前时间if nargin2,error(输入参数至少2个),endif nargin3|isempty(eps),eps=1e-6;endif narginepsk=k+1;fx0,dfx0=feval(fun,x0);if dfx0=0erro
33、r(fx在x0处的导数为0，终止运算)endx1=x0-fx0/dfx0;err=x1-x0;x0=x1;X(k)=x1;endif k=maxitererror(迭代次数超限，迭代失败)endx=x1;iter=k;X=X;format long; %双精度disp(x);disp(iter);disp(X);toc %记录程序完成时间%f.m函数描述文件，主要求解函数的表达式值及导数值function y,dy=f(x);y=x-2-log(x);ff=sym(x-2-log(x);dy=diff(ff);dy=subs(dy,x);程序四1. 复化梯形公式调用M文件程序function
34、I=T(x,y)%复化梯形求积%x为自变量，y为节点处的函数值n=length(x);m=length(y);if n=merror(x与y的维数不相同);endh=(x(n)-x(1)/(n-1);a=1 2*ones(1,n-2) 1;I=h/2*sum(a.*y);2. 复化Simpson公式调用M文件程序function I=simpson(x,y)%复化simpson求积公式%x为自变量，y为节点处的函数值n=length(x);m=length(y);if n=merror(x与y的维数不相同);endif rem(n-1,2)=0I=T(x,y);return;endN=(n-1
35、)/2;h=(x(n)-x(1)/N;a=zeros(1,n);for k=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4;a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);2. Romberg算法调用M文件程序function I=Romberg(fun,a,b)%Romberg求积公式%fun被积函数，a，b积分区间端点m=1;h=b-a;I=h/2*(feval(fun,a)+feval(fun,b);T(1,1)=I;while 1N=2(m-1);h=h/2;I=I/2;for i=1:NI=I+h*feval(fun,a+(2*i-1)*h);endT(m+1,1)=I;M=2*N;k=1;while M1;T(m+1,k+1)=(4k*T(m+1,k)-T(m,k)/(4k-1);M=M/2;k=k+1;endif abs(T(k,k)-T(k-1,k-1)1e-5break;endm=m+1;endI=T(k,k);