非线性方程求根、解线性方程组的迭代法.doc

福建农林大学计算机与信息学院课程实习报告课程名称：数值分析课程实习实习题目：非线性方程求根、解线性方程组的迭代法姓 名：系：数学与应用数学专 业：数学与应用数学年 级：2008级学 号：指导教师：职 称：讲师2011年6月 3日 福建农林大学计算机与信息学院数学类课程实习报告结果评定评语：成绩：指导教师签字：评定日期： 目 录 一、非线性方程求根11.实习的目的和任务12.实习要求13.实习地点14.主要仪器设备（实验用的软硬件环境）15.实习内容15.1二分法：15.2牛顿迭代法45.3弦截法66.问题讨论与分析97.结束语11参考文献12二、解线性方程组的迭代法131. 实习的目的和任务132. 实习要求133.实习地点134.主要仪器设备（实验用的软硬件环境）135.实习内容135.1 雅可比迭代法135.2 高斯-塞德尔迭代法166.问题讨论与分析187.结束语19参考文献20一、非线性方程求根1. 实习的目的和任务a) 掌握非线性方程（组）的各种解法，包括二分法、牛顿迭代法等，并通过编程与上机运算，体会二分法与牛顿迭代法的不同特点；b) 掌握解非线性方程的弦截法，并与牛顿迭代法作比较；c) 了解各种方法的收敛性。2. 实习要求a) 能够熟练应用所学的数值计算方法；b) 能够熟练使用MATLAB软件；c) 对数值分析及其计算方法有进一步了解，通过实例，能够对所学非线性方程求根方法的优缺点有所体会。3. 实习地点 田家炳实验楼5134. 主要仪器设备（实验用的软硬件环境） 方正商祺PC，MATLAB7.05. 实习内容5.1二分法：用二分法求方程在附近的根。5.1.1二分法求解原理二分法又称为对分法，是求解非线性方程根最简单是方法，用它求解的基本思想是由介值定理的大的。介值定理：若函数在区间上连续，且，则存在一点，使得。介值定理相当于方程根的存在性定理，称满足定理条件的区间为有根区间。考察有根区间，取中点将它分为两半，然后进行根的搜索，即检查与是否同号，如果确系同号，说明所求的根在的右侧，这时令；否则在的左侧，这时令。不管出现哪一种情况，新的有根区间的长度仅为的一半。对压缩了的有根区间又可施行同样的手续，即用中点将区间再分为两半，然后通过根的搜索判定所求的根在的哪一侧，从而又确定一个新的有根区间，其长度是的一半。如此反复二分下去，即可得出一系列有根区间，其中每个区间都是前一个区间的一半，因此的长度当时趋于零。就是说，如果二分过程无限地继续下去，这些区间最终必将收缩于一点，该点显然就是所求的根。每次二分后，设取有根区间的中点作为根的近似值，则在二分过程中可以获得一个近似根的序列，该序列必以根为极限。不过在实际计算时，不可能完成这个无限过程，其实也没有这种必要，因为数值分析的结果允许带有一定的误差。由于，只要二分足够多次（即充分大），便有，这里为预定的精度。5.1.2二分法求解的计算步骤步骤一：准备 计算在有根区间端点处的值；步骤二：二分 计算在区间中点处的值；步骤三：判断 若，则即是根，计算过程结束。否则作如下检验： 若与异号，则根位于区间内，这时以代替； 若与同号，则根位于区间内，这时以代替；反复执行步骤二和步骤三，直到区间长度缩小到允许误差范围之内，此时区间中点即可作为所求的根。5.1.3二分法的方程求解1）利用MATLAB软件，建立m文件dichotomy.mfunction k,x,e,y=dichotomy(fun,a,b,precision)%二分法%输入a,b为区间左右端点，precision是预选设定的精度%输出k为二分法的次数，x是方程在(a,b)内的实根x*的近似值，if nargin4 | narginprecision ya=subs(fun,a);yb=subs(fun,b); if ya*yb0 disp(警告：区间端点的函数值同号，请重新输入区间端点！),return end x=(a+b)/2;y=subs(fun,x); if y=0,return elseif y*yb0,b=x; else a=x; end e=abs(a-b); k=k+1; k,xendx=(a+b)/2;e=abs(a-b);y=subs(fun,x);2）在MATALB窗口输入如下指令： k,x,e,y=dichotomy(x3+x2-3*x-3,1,2,1e-6)3）运行结果：ans =1.000000000000000 1.5000000000000002.000000000000000 1.7500000000000003.000000000000000 1.6250000000000004.000000000000000 1.6875000000000005.000000000000000 1.7187500000000006.000000000000000 1.7343750000000007.000000000000000 1.7265625000000008.000000000000000 1.7304687500000009.000000000000000 1.73242187500000010.000000000000000 1.73144531250000011.000000000000000 1.73193359375000012.000000000000000 1.73217773437500013.000000000000000 1.73205566406250014.000000000000000 1.73199462890625015.000000000000000 1.73202514648437516.000000000000000 1.73204040527343817.000000000000000 1.73204803466796918.000000000000000 1.73205184936523419.000000000000000 1.73204994201660220.000000000000000 1.732050895690918k = 20x = 1.732050418853760e = 9.536743164062500e-007y = -3.678838435661191e-0064）结果分析：该求解过程进行了20次迭代，得到最后结果x =1.7320504188537605.2牛顿迭代法用牛顿迭代法求方程在附近的根，并由迭代次数分析结果。5.2.1 Newton法求解原理对于方程，如果是线性函数，则对它求根是容易的。法实质上是一种线性化方法，其基本思想是将非线性方程逐步归结为某种线性方程来求解。把函数在某一初值附近展开成Taylor级数，有取其线性部分，近似地代替函数可得方程的近似式设，解该近似方程可得可作为方程式的近似解。重复以上过程，得迭代公式法有明显的几何解释。如图（1）所示方程的根可解释为曲线与轴的交点的横坐标。设是根的某个近似值，过曲线上横坐标为的点引切线，并将该切线与轴的横坐标作为的新的近似值。注意到切线方程为图（1）著名的公式为5.2.2Newton法求解计算步骤步骤一：准备 选型初始近似值，计算，；步骤二：迭代 按公式迭代一次，得新的近似值，计算；步骤三：控制 如果满足或，则终止迭代，以作为所求的根；否则转步4。此处是允许误差，而其中是取绝对误差或相对误差的控制常数，一般可取;步骤四：修改 如果迭代次数达到预先指定的次数或者，则方法失败，否则以代替转步骤二继续迭代。5.2.3 Newton法的方程求解1）利用MATLAB软件，建立m文件Newton.mfunction x_fin,index,it=Newton(fun,x,ep,it_max)if nargin 4 it_max=100;endif nargin 3 ep=1e-5;endindex=0;k=1;while k = it_max x1=x;f=feval(fun,x); if abs(f(2) ep break ;end x=x-f(1)/f(2); if abs(x-x1) fun=inline(x3+x2-3*x-3,3*x2+2*x-3);x_fin,index,it=Newton(fun,1.5)3）运行结果：ans =2.000000000000000 1.7777777777777783.000000000000000 1.7333606669400034.000000000000000 1.732051929409472x_fin =1.732050807569701index =1it = 44）结果分析：该求解过程进行了4次迭代，得到最后结果x_fin =1.7320508075697015.3弦截法用弦截法求方程在附近的根，并由迭代次数分析结果。5.3.1弦截法求解原理设是的近似根，利用构造一次差值多项式，并用的根作为的新的近似根，由于（5.5.1.1）因此有（5.5.1.2）这样导出的迭代公式（5.5.1.2）可以看做Newton公式中的导数用差商取代的结果。其迭代的几何意义，如图（2）所示曲线上的横坐标的点分别记作，则弦线的斜率等于差商值，其方程是图（2）因此，按照式（5.5.1.2）求得的实际上是弦线与x轴的交点的横坐标。这种算法因此称为弦截法。其主要思想为：任取两个数，求得对应的函数值。如果两函数值同号，则重新取数，直到这两个函数值异号为止。连接与这两点形成的直线与轴相交于一点，求得对应的，判断其与f中的哪个值同号。如与同号，则为新的。将新的与连接，如此循环。体现的是极限的思想。5.3.2弦截法求解计算步骤步骤一：准备 选取初始近似值，计算相应的函数值步骤二：迭代 按公式迭代得到新的近似值，计算。步骤三：控制 如果满足或者，则认为过程收敛，终止迭代过程而以作为所求的根；否则执行步骤四，因此此处的，是允许的误差，而其中C是预期指定的控制数。步骤四：修改 如果迭代次数达到预先指定的次数N，则认为过程不收敛，计算失败；否则以，分别代替，而后转步骤二继续迭代。5.3.3 弦截法的方程求解1）利用MATLAB软件，建立m文件secant.mfunction k,x,e,y=secant(fun,x0,x1,precision,N)%切线法,弦截法%输入初始值x0,x1，要求的近似根xk的函数值f(xk)的误差限precision，迭代次数的最大值N及其函数fun%输出迭代序列xk、迭代次数、迭代k次得到的迭代值xk、函数值f(xk)、相邻两次迭代的偏差e=|y-y2|（迭代次数超过最大值N时，输出信息警告：迭代次数超过给定的最大值，过程不收敛。）。if nargin5 | nargin3 error(错误的输入参数个数！)endswitch nargin case 3,precision=1e-6;N=100; case 4,N=100;endk=1;format long;while k=N f=0;y0=subs(fun,x0);y1=subs(fun,x1); x=x1-y1*(x1-x0)/(y1-y0);y=subs(fun,x);e=abs(y-y1); if e k,x,e,y=secant(x3+x2-3*x-3,1,2,1e-5)3）运行结果：ans =2.000000000000000 1.5714285714285713.000000000000000 1.7054108216432874.000000000000000 1.7351357706607395.000000000000000 1.7319963707826996.000000000000000 1.732050697785584k =6x =1.732050807572790e =1.039037202943405e-006y =3.702993467413762e-0114）结果分析：该求解过程进行了6次迭代，得到最后结果x =1.732050807572790；结果显示用牛顿迭代法进行了4次迭代，用弦截法用了6次迭代结果x_fin =1.732050807569701，由此可知用牛顿迭代法求解方程收敛性更好。6. 问题讨论与分析1、何谓二分法？二分法的优点是什么？如何估计误差？二分法又称为对分法，是考察有根区间，取中点将它分为两半，然后进行根的搜索，即检查与是否同号，如果确系同号，说明所求的根在的右侧，这时令；否则在的左侧，这时令。不管出现哪一种情况，新的有根区间的长度仅为的一半。对压缩了的有根区间又可施行同样的手续，即用中点将区间再分为两半，然后通过根的搜索判定所求的根在的哪一侧，从而又确定一个新的有根区间，其长度是的一半。如此反复二分下去，即可得出一系列有根区间，其中每个区间都是前一个区间的一半，因此的长度当时趋于零。就是说，如果二分过程无限地继续下去，这些区间最终必将收缩于一点，该点显然就是所求的根。每次二分后，设取有根区间的中点作为根的近似值，则在二分过程中可以获得一个近似根的序列，该序列必以根为极限。二分法的优点是算法简单，而且收敛性总能得到保证。不过在实际计算时，不可能完成这个无限过程，其实也没有这种必要，因为数值分析的结果允许带有一定的误差。由于，只要二分足够多次（即充分大），便有，这里为预定的精度。2、何谓迭代法？它的收敛条件、误差估计式是什么？迭代法也称辗转法，是一种不断用变量的旧值递推新值的过程,跟迭代法相对应的是直接法（或者称为一次解法），即一次性解决问题。迭代法又分为精确迭代和近似迭代。收敛条件：设为方程的根，在的邻近连续且，则迭代过程在邻近具有局部收敛性。误差估计式：3、 比较迭代法收敛的快慢？迭代法收敛快慢的可以用迭代法收敛速度来比较。收敛速度：设序列收敛于，令若存在某实数及正常数, 使则称序列为阶收敛, 如果序列 是由x n生的, 且 阶收敛, 则称这种迭代过程是 阶收敛的。当, 且 时, 称为线性收敛；当=2, 称为平方收敛(或二次收敛)；当 时, 称为超线性收敛。4、 法与牛顿法的优劣牛顿法要求初值充分接近根以保证局部收敛性；牛顿迭代法的主要优点是收敛较快，是平方收敛的缺点是公式中需要求 的导数；若比较复杂，则使用牛顿公式就大为不便。牛顿法用单点迭代法，在计算时只用到前一步的值；弦截法多点迭代法，在求时要用到前两步的值和。与牛顿法相比，弦截法的收敛速度也是比较快的。可以证明，弦截法具有超线性收敛速度，收敛阶为7. 结束语这次的解非线性方程中，了解了二分法、牛顿法、弦截法的不同和优劣。了解了如何比较不同迭代法收敛快慢的比较。在计算上，弦截法只计算了一个公式，并不需要计算求导公式，比牛顿法减少了很多计算量。而计算最简单的是二分法，且其估计误差控制迭代次数也比其他两种方法简单。21参考文献1 张德丰.Matlab数值分析与应用（第2版）M.北京：北京国防工业出版社，2009.162-1772 李庆扬等.数值分析（第4版）M.武汉：华中科技大学出版社，2006.3 胡洁. 牛顿法和弦截法在方程中的应用比较研究J. 武汉：武汉理工大学理学院，20064 张菁, 张丽梅. 迭代法收敛速度的比较J. 大连：大连水产学院,，2007二、解线性方程组的迭代法1. 实习的目的和任务（a）熟悉迭代法的有关理论和方法；（b）会编制雅可比迭代法、高斯-塞德尔迭代法的程序；（c）注意所用方法的收敛性及其收敛速度问题。2. 实习要求a) 能够熟练应用所学的数值计算方法；b) 能够熟练使用MATLAB软件；c) 对数值分析及其计算方法有进一步了解，通过实例，能够对所学解线性方程组的迭代法优缺点有所体会。3. 实习地点 田家炳实验楼5134. 主要仪器设备（实验用的软硬件环境） 方正商祺PC，MATLAB7.05. 实习内容5.1 雅可比迭代法用雅可比迭代法解方程组.注意：若用高斯-塞德尔迭代法则发散。这里取。5.1.1雅可比迭代法原理设有方程组，即式中，非奇异，且；先将所给的方程组化成如下形式的同解方程组：然后给出初始向量 并按迭代公式：因此可以得到：其中相应的迭代公式为：称上式的迭代公式为雅可比迭代。雅可比迭代法的收敛性：求解出雅可比迭代法的迭代矩阵，记，则，方程组改写为：相应的矩阵形式的迭代公式为则该迭代法的迭代矩阵为：求出迭代矩阵的特征值，判断给最大特征值的绝对值是否小于1，如果小于1，说明该迭代法收敛，否则发散。5.1.2雅可比迭代法算法（1） 取初始点，置，精度要求和最大迭代次数N。（2） 用上述迭代公式计算。（3） 若，则停止计算（作为线性方程组的解）。若，则停止计算（输出某些信息）；否则置，转（2）。5.1.3雅可比迭代法解方程组1）编写Jacobi迭代法的M文件Jacobi.mfunction x, k, index=Jacobi(A, b, ep, it_max)% 求解线性方程组的Jacobi迭代法，其中% A - 方程组的系数矩阵% b - 方程组的右端项% ep - 精度要求。省缺为1e-5% it_max - 最大迭代次数，省缺为100% x - 方程组的解% k - 迭代次数% index - index=1表示迭代收敛到指定要求；% index=0表示迭代失败if nargin 4 it_max=100; endif nargin 3 ep=1e-5; endn=length(A); k=0;x=zeros(n,1);y=zeros(n,1);index=1;while 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-A(i,j)*x(j); end end if abs(A(i,i)1e-10 | k=it_max index=0; return; end y(i)=y(i)/A(i,i); end if norm(y-x,inf) A=1 2 -2;1 1 1;2 2 1;b=7;2;5;ep=1e-5;x,k,index=Jacobi(A,b,ep)3）运行结果：x =1 2 -1k =3index =14）结果分析：求解得到结果经过3次迭代，方程组的解为5.2 高斯-塞德尔迭代法用高斯-塞德尔迭代法解方程组.注意：若用雅可比迭代法则发散。这里取。 5.2.1高斯-塞德尔迭代法原理Gauss-Seidel迭代法是Jacobi迭代法的一种改进形式，仔细研究Jacobi迭代公式就会发现，当计算到时，分量，已都计算出，但在计算时却将它们放置一边，而仍然使用旧分量，。从直观上看，新计算出的分量要比旧的准确些，因此可以设想，一旦新分量已求出就马上用新分量进行计算，这样得到的计算结果要优于原来的计算结果。即在Jacobi迭代法中用，替代，就是Gauss-Seidel迭代法的基本思想。Gauss-Seidel迭代法的迭代公式如下： （1）下面写出迭代公式的矩阵形式，分量形式的迭代公式（1）等价于因此得到 （2）写成 （3）式中，。称迭代公式（1）或（2）为Gauss-Seidel迭代，称式（3）中的为Gauss-Seidel迭代矩阵。5.2.2高斯-塞德尔迭代法算法：（1） 取初始点，置，精度要求和最大迭代次数（2） 用式（1）或式（2）计算（3） 若，则停止计算（作为线性方程组的解）若,则停止计算（输出某些信息）；否则置，转（2）。5.2.3高斯-塞德尔迭代法解方程组1）编写Gauss_seidel迭代法的M文件Gauss_Seidel.mfunction x,k,flag=Gauss_Seidel(A,b,delta,max1)% 求解线性方程组的迭代法，其中% A为方程组的系数矩阵% b为方程组的右端项% delta为精度要求，缺省值为1e-5% max1为最大迭代次数，缺省值为100% x为方程组的解% k为迭代次数% flag为指标变量 flag=OK!表示迭代收敛到指标要求% flag=fail!表示迭代失败 if nargin4 max1=100;endif nargin3 delta=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);flag=OK!;while 1 y=x; for i=1:n z=b(i); for j=1:n if j=i z=z-A(i,j)*x(j); end end if abs(A(i,i)1e-10|k=ma