数值分析上机实习报告.doc

数值分析 上机实习报告学 号: 姓 名: 专 业: 联系电话:任课老师: 二零一二年十二月数值分析上机实习报告 第III页序 言数值分析在现代科学发展中有着重要的作用，而随着科学的发展进步，越来越多的数值分析问题不能够光靠人力计算，这就要借助计算机进行计算。而在利用计算机解决实际问题时，要根据具体情况作出可靠的理论分析，才能够写出比较可靠的程序。现在面向数值分析问题的计算机软件有：C、C+、MATLAB、Python、Fortran等。C+是笔者在本科学过的唯一一门编程语言，但是由于学习时间较短，而且在学习时不精，再加上时间已久远，对这门编程语言课程已经几乎没有多少印象了。Python是一种面向对象的解释性的计算机程序设计语言，也是一种功能强大而完善的通用型语言，已经具有十多年的发展历史，成熟且稳定。Python 具有脚本语言中最丰富和强大的类库，足以支持绝大多数日常应用。Fortran为“公式翻译器”，它是世界上最早出现的计算机高级程序设计语言，广泛应用于科学和工程计算领域。Fortran语言以其特有的功能在数值、科学和工程计算领域发挥着重要作用。MATLAB（矩阵实验室）是一个功能强大的软件，是一种数值计算环境和编程语言。在当今世界流行的30多个数学类软件中，MATLAB语言处于数值计算型软件的主导地位，适用范围涵盖了工程数学的各个方面。它的有点主要有：1、matlab是以矩阵为基础的工具，若是编一些对速度没有要求的，进行数值计算或者信号处理的小程序，可以用matlab，且简单。2、matlab除具备卓越的数值计算能力外，它还提供有专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。3、matlab的基本数据单位是矩阵，它的指令表达式与数学，工程中常用的形式十分相似，所以用matlab来解算问题要比用C、FORTRAN等语言完成相同的事情简捷得多。在新版本中也加入了对C、FORTRAN、c+、JAVA的支持，使用时可以直接调用，也可将编写的实用程序导入到matlab函数库中方便以后使用时调用。本次编程所用的软件为MATLAB，希望通过这次作业，能够对它有了初步的认识，为以后的学习和工作奠定一定的基础。目 录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算法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，如果取则迭代格式如下：所以具体迭代格式过程为：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) 其输入方法同上，不再赘述。这里只给出运行结果：Gauss消元法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迭代法求解第三小题时，出现不明确的数值结果，说明在使用这一方法时应该注意方程组的形式； 数值分析上机实习报告 第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步之后的计算不再有舍入误差，只是由k引起的扰动m，都有m0, 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(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=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*(-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,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),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b); while toleps k=k+1; h=h/2; Q=0; for i=1:M x=a+h*(2*i-1); Q=Q+subs(sym(f),findsym(sym(f),x); end T(k+1,1)=T(k,1)/2+h*Q; M=2*M; for j=1:k T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j)/(4j-1); end tol=abs(T(k+1,j+1)-T(k,j); end I=T(k+1,k+1); step=k;end选做题二：newton算法function y = newton(x)%UNTITLED5 Summary of this function goes here% Detailed explanation goes hereformat long;x1=x-func2_1(x)/func2_1_1(x);if(abs(x1)1e-6|abs(func2_1(x)1e-6) x=x1; x1=x-func2_1(x)/func2_1_1(x); if(abs(x1)1.5) delt=abs(x1-x); else delt=abs(x1-x)/x1); end if func2_1(x1)=0 break end endy=x;en