数值分析上机实习报告

1、.数值分析 上机实习报告学 号: 姓 名: 专 业: 联系电话:任课老师: 二零一二年十二月；数值分析上机实习报告 第III页序 言数值分析在现代科学发展中有着重要的作用，而随着科学的发展进步，越来越多的数值分析问题不能够光靠人力计算，这就要借助计算机进行计算。而在利用计算机解决实际问题时，要根据具体情况作出可靠的理论分析，才能够写出比较可靠的程序。现在面向数值分析问题的计算机软件有：C、C+、MATLAB、Python、Fortran等。C+是笔者在本科学过的唯一一门编程语言，但是由于学习时间较短，而且在学习时不精，再加上时间已久远，对这门编程语言课程已经几乎没有多少印象了。Python是一

2、种面向对象的解释性的计算机程序设计语言，也是一种功能强大而完善的通用型语言，已经具有十多年的发展历史，成熟且稳定。Python 具有脚本语言中最丰富和强大的类库，足以支持绝大多数日常应用。Fortran为“公式翻译器”，它是世界上最早出现的计算机高级程序设计语言，广泛应用于科学和工程计算领域。Fortran语言以其特有的功能在数值、科学和工程计算领域发挥着重要作用。MATLAB（矩阵实验室）是一个功能强大的软件，是一种数值计算环境和编程语言。在当今世界流行的30多个数学类软件中，MATLAB语言处于数值计算型软件的主导地位，适用范围涵盖了工程数学的各个方面。它的有点主要有：1、matlab是以

3、矩阵为基础的工具，若是编一些对速度没有要求的，进行数值计算或者信号处理的小程序，可以用matlab，且简单。2、matlab除具备卓越的数值计算能力外，它还提供有专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。3、matlab的基本数据单位是矩阵，它的指令表达式与数学，工程中常用的形式十分相似，所以用matlab来解算问题要比用C、FORTRAN等语言完成相同的事情简捷得多。在新版本中也加入了对C、FORTRAN、c+、JAVA的支持，使用时可以直接调用，也可将编写的实用程序导入到matlab函数库中方便以后使用时调用。本次编程所用的软件为MATLAB，希望通过这次作业，能够对它

4、有了初步的认识，为以后的学习和工作奠定一定的基础。目 录1.第一题11.1.题目11.2.Gauss消元法及Guass-Seidel迭代法11.3问题的求解11.4 方法总结32.第二题42.1题目42.2 Runge-Kutta法的基本思想42.3 问题的求解42.4方法总结53.选做题一63.1.题目63.2.基本理论63.3.问题求解83.4.方法总结84.选做题二104.1.题目104.2.基本理论104.3.问题求解104.4.方法总结12总结13附件14第一题：gauss.m文件14gauseidel.m文件14第二题：四阶Runge-Kutta 算法15选做题一：Romberg算

5、法16选做题二：newton算法16西南交通大学数值分析上机实习报告 第3页1. 第一题1.1. 题目写出对一般的线性方程组通用的Gauss消元, Gauss-Seidel迭代程序。并以下面的线性方程组为例进行计算，讨论所得到的计算结果是否与理论一致。(1) (2)(3)1.2. Gauss消元法及Guass-Seidel迭代法1.2.1 Gauss消元法高斯消元法的原理是：若用初等行变换将增广矩阵化为行阶梯阵，则AX = B与CX = D是同解方程组。所以我们可以用初等行变换把增广矩阵转换为行阶梯阵，然后回代求出方程的解。1.2.2 Guass-Seidel迭代法A分解为A=D-L-U，如果

6、取则迭代格式如下：所以具体迭代格式过程为：1.3问题的求解1) 在命令行中分别调用函数文件gauss.m及gauseidel.m文件。命令如下A=6 2 -1;1 4 -2;-3 1 4; b=-3 2 4'X=gauss(A,b)x,n=gauseidel(A,b,0 0 0',1e-6)运行结果如下：X = -0.7273 0.8081 0.2525x = -0.7273 0.8081 0.2525n = 15其理论值计算，在matlab命令行中输入Ab运行结果ans = -0.7273 0.8081 0.25252) 其输入方法同上，不再赘述。这里只给出运行结果：Gaus

7、s消元法5.76920.7692-4.2308Gauss-Seidel法5.76920.7692-4.2308理论值5.76920.7692-4.23083) 其输入方法同上，不再赘述。这里只给出运行结果:Gauss消元法-0.63641.5455Gauss-Seidel法NaNNaN理论值-0.63641.54551.4 方法总结在用Gauss消元法及Guass-Seidel迭代法时，可以得出以下结论：1、Gauss消元法及Guass-Seidel迭代法在求解线性方程组时通常具有较高的精度，其与理论计算结果一致；2、在使用Guass-Seidel迭代法求解第三小题时，出现不明确的数值结果，说

8、明在使用这一方法时应该注意方程组的形式； 数值分析上机实习报告 第16页2. 第二题2.1题目给定初值问题 （精确解为），用Runge-Kutta 4阶算法按步长求解，分析其中遇到的现象及问题。2.2 Runge-Kutta法的基本思想Runge-Kutta法的基本思想是通过f (x, y)某些点函数值的适当线性组合替换Euler 法中的f(xk,yk)，可能使得方法的精确度更高。由于有了这种思想，因此标准的四阶Runge-Kutta的算法如下(h为求解步长）：K1=hf(xk,yk)如果yk某种方法第k步的近似值，y(xk)是其准确值，其绝对误差为k，即有k=y(xk)-yk。假定第k步之后

9、的计算不再有舍入误差，只是由k引起的扰动m，都有m<k，则称此方法是绝对稳定的。数值解法的稳定性一般与步长h以及f(x,y)有关。2.3 问题的求解在MATLAB 中按照标准四阶Runge-Kutta 的算法编写程序代码，并且运行程序。解：首先探讨不同步长下算法稳定的问题，取不同步长得到求解的微分方程的解同精确解的误差分别如表6、表7、表8。表6 h=0.05下的Runge-Kutta四阶算法0.10000000000000 0.13533954843050 0.13533528323661 0.00000426519389 0.20000000000000 0.018316793369

10、37 0.01831563888873 0.00000115448064 0.30000000000000 0.00247898654331 0.00247875217667 0.00000023436664 0.40000000000000 0.00033550491934 0.00033546262790 0.00000004229143 0.50000000000000 0.00004540708428 0.00004539992976 0.00000000715452 0.60000000000000 0.00000614537428 0.00000614421235 0.000000

11、00116193 0.70000000000000 0.00000083171218 0.00000083152872 0.00000000018346 0.80000000000000 0.00000011256355 0.00000011253517 0.00000000002838 0.90000000000000 0.00000001523430 0.00000001522998 0.00000000000432 1.00000000000000 0.00000000206180 0.00000000206115 0.00000000000065 表7 h=0.1下的Runge-Kut

12、ta四阶算法xkyky(xk)k0.10000000000000 0.333333333333320.135335283236610.197998050096710.20000000000000 0.111111111111100.018315638888730.092795472222370.30000000000000 0.037037037037030.002478752176670.034558284860370.40000000000000 0.012345679012340.000335462627900.012010216384440.50000000000000 0.00411

13、5226337450.000045399929760.004069826407690.60000000000000 0.001371742112480.000006144212350.001365597900130.70000000000000 0.000457247370830.000000831528720.000456415842110.80000000000000 0.000152415790280.000000112535170.000152303255100.90000000000000 0.000050805263430.000000015229980.0000507900334

14、51.00000000000000 0.000016935087810.000000002061150.00001693302665表8 h=0.2下的Runge-Kutta四阶算法xk（*10-3）yk（*10-3）y(xk) （*10-3）k（*10-3）0.000200000000000.005000000000000.000018315638890.004981684361110.000400000000000.025000000000000.000000335462630.024999664537370.000600000000000.125000000000010.00000000

15、6144210.124999993855790.000800000000000.625000000000040.000000000112540.624999999887500.001000000000003.125000000000250.000000000002063.124999999998192.4方法总结（1）用标准的四阶 Runge-Kutta 算法计算常微分方程时，不同的步长算法的稳定性有影响。不够稳定的步长下面的计算，误差会越来越大，结果失真严重。（2）一般情况下，步长越小，标准的四阶 Runge-Kutta 算法的稳定性越高，精度也越高。3. 选做题一3.1. 题目已知，利用复

16、化梯形公式、复化Simpson公式和Romberg算法求的近似值；并观察实际计算结果，比较它们的收敛速度。3.2. 基本理论将积分区间a,b划分n等分，步长，求积节点，在每个小区间上应用梯形公式：然后将它们累加求和，作为所求积分I的近似值.记     式为复化梯形求积公式，下标n表示将区间n等分。将积分区间a,b划分n等分，记子区间的中点为在每个小区间上应用辛普森公式，则有其中记                          

17、60;式为复化辛普森求积公式。准备初值，计算且（为等份次数） 2按梯形公式的递推关系，计算 3按龙贝格公式计算加速值 4精度控制。对给定的精度，若则终止计算，并取作为所求结果；否则，重复24步，直到满足精度为止。3.3. 问题求解1) 复化梯形公式程序如下：x=0:0.1:1;Y=4./(1+x.2);T=trapz(x,Y);运行结果：T = 3.13992) 复化Simpson公式程序如下： fun=inline('4./(1+x.2)');Q=quad(fun,0,1);运行结果：Q = 3.14163) Romberg算法编写了Rg.m文件，在命令行进行调用，调用命令如

18、下：q,s=Rg('4/(1+x2)',0,1)其运行结果：q = 3.1416s = 43.4. 方法总结通过比较其运算结果，复化梯形公式求得的结果精度稍差一些、复化Simpson公式和Romberg算法求得的结果较为一致；但是从其收敛速度来看，复化梯形公式最慢，复化Simpson公式次之，Romberg算法的收敛速度最快。4. 选做题二4.1. 题目利用牛顿法求方程于区间的根，考虑不同初值下牛顿法的收敛情况。4.2. 基本理论首先，选择一个接近函数零点的，计算相应的和切线斜率（这里表示函数的导数）。然后我们计算穿过点并且斜率为的直线和轴的交点的坐标，也就是求如下方程的解：我

19、们将新求得的点的坐标命名为，通常会比更接近方程的解。因此我们现在可以利用开始下一轮迭代。迭代公式可化简为如下所示：已经证明，如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那么牛顿法必定收敛。 并且，如果不为0, 那么牛顿法将具有平方收敛的性能. 粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。4.3. 问题求解编制了newton.m函数文件，在命令行中进行调用：y=newton(2)运行结果：y =3.146194257027137考虑不同初值情况下分别计算了其解，并代会原式，计算了其误差。见表4-1所示。表4-1 不同初值下计算

20、结果X0Y代入原式真实值1InfInf223.1461942570271372.000000706990555233.1461934407979102.000000150195192243.1461932212421842.000000000424029253.1461932843192922.000000043452430263.1461937809908402.000000382259707273.1461932206209122.000000000000225283.1461932206226392.000000000001403293.1461932206285502.00000000

21、00054352103.1461932206434282.0000000000155842113.1461932206739442.0000000000364012其收敛情况见图4-2所示。图4-2 收敛情况比较表4.4. 方法总结当初始值为1时，计算结果不收敛。初始值为36时收敛不佳。当初始值为711时，收敛结果比较良好。 总结通过此次数值分析上机实习，认识到了计算机编程在数值分析课程中的重要性。也让我意识到了学习编程语言的重要性。要想要将数学应用于实际工程中，特别是对于工科的学生来讲，这是一门非常重要的实践课程。在此次报告中，首次接触了Matlab这门软件，之所以选这门，是因为很多工程中对

22、数据的处理需要使用，这对我本身的专业是非常重要的。虽然之前学过C+，但是对于这门新的计算语言还是很头痛的，例如对函数的定义、调用以及很多形式方面跟之前的比较，还是有了很大的差异，这不得不使我在编程上面再次下功夫了。本学期学习的数值分析方法，让我在对数学的理解上又有了新的认识。在一连串看似杂乱无章的数据中，用数值分析进行处理，有了很多的收获。在学习的过程中，意识到了这门课对我专业的重要性，让我对自己的专业有了更好的发挥。对我写论文以及解决实际工程问题有很大的帮助。最后，感谢老师给我这样的机会去接触这门语言，虽然只了解了皮毛，可是仍然收获颇多。由于初次接触这门软件，在报告中仍难免会有不完善甚至错误

23、的地方，望谅解！附件第一题：gauss.m文件function X=gauss(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, return end if RA=RB if RA=n X=zeros(n,1);C=zeros(1,n+1); for p= 1:n-1 for 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); end end b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A

24、(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); end else end endendgauseidel.m文件function y,n=gauseidel(A,b,x0,eps)%UNTITLED Summary of this function goes here% Detailed explanation goes hereif nargin=3 eps=1e-6;elseif nargin<3 error returnendD=diag(diag(A);L=-tril(A,-1);U=-triu(A

25、,1);G=(D-L)U;f=(D-L)b;y=G*x0+f;n=1;while norm(y-x0)>=eps x0=y; y=G*x0+f; n=n+1;endend第二题：四阶Runge-Kutta 算法h=0.1,0.2k=1for xk=0.1:0.1:1 %取不同的x进行计算for i=1:1:4 %取不同的步长计算每个xhj=h(i)yk=1for xkp=hj:hj:xk %四阶标准Runge-Kutta迭代算法k1=hj*(-20*yk);format longk1;k2=hj*(-20*(yk+k1*0.5);k3=hj*(-20*(yk+k2*0.5);k4=hj*

26、(-20*(yk+k3);ykm=yk+(k1+2*k2+2*k3+k4)*(0.16666666666667);yk=ykm;endyk=ykm;yeal=exp(-20*xkp) ; %求出精确值ck=abs(yk-yeal); %精确值和计算值之间的绝对误差ykeal(k)=yeal;ckk(k)=ck;yxk(k)=yk;xkk(k)=xk;k=k+1;endendD=yxk' ykeal' ckk' %将迭代得到的结果及误差存到矩阵D中for i=1:1:10 %将同一个步长下不同X得到的结果集中到一个矩阵for j=1:1:3Ah1(i,j)=D(4*i-3

27、,j);Ah2(i,j)=D(4*i-2,j);Ah3(i,j)=D(4*i-1,j);Ah4(i,j)=D(4*i,j);endAxk(i)=xkk(4*i-3);endCh1=Axk' Ah1Ch2=Axk' Ah2Ch3=Axk' Ah3Ch4=Axk' Ah4Cheng=Ch1(1,2 3 4);Ch2(1,2 3 4); Ch3(1,2 3 4);Ch4(1,2 3 4)选做题一：Romberg算法function I,step = Rg(f,a,b,eps)%UNTITLED4 Summary of this function goes here% Detailed explanation goes here if nargin=3 eps=1e-4; end M=1; tol=10; k=0; T=zeros(1,1); h=b-a; T(1,1)=(h/2)*(subs(sym(f),fin