毕业设计（论文）-基于Newton法对几类非线性方程组求解的研究.pdf

基于 Newton 法对几类非线性方程组 求解的研究 学生姓名学生姓名：钟良艳钟良艳 学学号：号：200664090202200664090202 班班级：级：数学数学 06020602 班班 所在院所在院( (系系) ): :数学与计算科学学院数学与计算科学学院 指导教师指导教师：张宏伟张宏伟 完成日期完成日期: :20102010 年年 6 6 月月 CHANGSHACHANGSHA UNIVERSITYUNIVERSITY OFOF SCIENCESCIENCE simplified Newton iteration，modified Newton iteration, approximate Newton iteration and a higher order Taylor expansion about the iterative formula. In addition, through some examples we compare the iterative formula of the Newton method with those improved iterative formula,and we find those improved methods show a lots of advantages KeyKeyKeyKey words:words:words:words:the Newton method ; second order convergence ; the nonlinear equation 基于 Newton 法对几类非线性方程组求解的研究 目 录 1引言1 2Newton 法1 2.1 Newton 法的基本思想2 2.2 Newton 法的几何意义2 2.3 Newton 法的局部收敛性3 2.4 Newton 法的算法4 2.5 实例计算4 3求解非线性方程组的 Newton 迭代法5 3.1 解非线性方程组的 Newton 迭代法5 3.2 基于 Newton 法对几类非线性方程组求解的研究7 3.2.1 关于 Newton 法对含n阶多项式方程组的研究7 3.2.2 关于 Newton 法对含xsin,xcos 方程组的研究9 3.2.3 关于 Newton 法对含 x e方程组的研究10 4Newton 法的改进11 4.1 简化 Newton 迭代法11 4.2 修正 Newton 迭代法12 4.3 拟 Newton 迭代法16 4.4 近似程度更高的二阶泰勒展开式的迭代公式17 5Newton 法的优缺点19 5.1 Newton 法的优点19 5.2 Newton 法的缺点19 6结束语20 参考文献21 致谢22 附录23 基于 Newton 法对几类非线性方程组求解的研究 第 1页共 33 页 基于 Newton 法对几类非线性方程组求解的研究 1引言 在科学研究和工程计算中常常碰到非线性方程（组）求解问题。非线性方 程0)(=xf的解一般不能解析求出。所以数值解法显得非常重要，而数值解法在 实际中的实现则更为重要。 在实际应用中，只要能获得给定误差限度内的近似值 就可以了，这就需要方程近似求解。方程近似求解的方法很多，比如二分法、 简 单迭代法、Newton 迭代法等等。而 Newton 迭代法是求解非线性方程（组）零点 的一种非常重要的方法，其应用范围较广，可解代数方程和超越方程，也可解非 线性方程组，既可求方程实根，也可求复根。事实上，虽然求方程的根的方法很 多，但是其中牛顿法计算简单，收敛速度也好，为一般科技人员所采用。 在本论文中， 主要采用理论与实例相联系的分析方法，借助现代互联网这一 信息平台，获得了大量的有关这方面的具体文献，从而使问题的研究有据可依。 为了全面地分析探究这个课题，本文是先将用于求解线性问题的 Newton 法推广 至求解非线性问题，然后基于 Newton 法对含n阶多项式、xxcos,sin及 x e等非线 性方程组解进行研究。为了更深一步对牛顿法的掌握，还对 Newton 法迭代公式 进行了改进和推广：简化 Newton 迭代法、修正 Newton 迭代法、拟 Newton 法、 近似程度更高的二阶泰勒展开式的迭代公式，这样能使计算简便，迭代公式的收 敛速度加快，从而具有广大的应用范围，更能够提高创造的能力。并且对改进的 迭代公式都通过例子与 Newton 迭代公式进行了比较，显示其优越性。 2Newton 法 解非线性方程0)(=xf的牛顿(Newton) 法，就是将非线性方程线性化的一 种方法。它是解代数方程和超越方程的有效方法之一。 基于 Newton 法对几类非线性方程组求解的研究 第 2页共 33 页 2.1 Newton 法的基本思想 设函数)(xf在有限区间ba,上二阶连续可微， 0 x是根的某个近似值， 将)(xf在 0 x处作Taylor 展开，则有 2 0 0 00 0 )( 2 )( )()()(xx f xxxfxfxf+= （2-1） 其中 0 介于x和 0 x之间，如果用其线性主部近似)(xf，即 )()()( 00 0 xxxfxfxf+（2-2） 则将非线性方程0)(=xf近似为线性方程 0)()( 00 0 =+xxxfxf(2-3） 若0)( 0 xf，将其解记为 )( )( 0 0 01 xf xf xx=（2-4） 则得到根的新的近似值 1 x，一般地，在 k x附近线性化方程为 0)()( =+ kkk xxxfxf（2-5） 若0)( k xf，将其解记为 )( )( 1 k k kk xf xf xx= + （2-6） 迭代公式（2-6）称为 Newton 法迭代公式。 2.2 Newton 法几何意义 Newton法有明显的几何意义 （图2.1） 。 函数0)(=xf的根是曲线)(xfy= 与x轴的交点的横坐标。过点)(,( kkk xfxM作曲线的切线 k T，则切线方程为 )()( kkk xxxfxfy+=，切线 k T与x轴的交点的横坐标为 )( )( 1 k k kk xf xf xx= + 。 因 此 Newton 法又叫做切线法。 基于 Newton 法对几类非线性方程组求解的研究 第 3页共 33 页 图 2.1 牛顿法的几何意义图解 2.3 Newton 法的局部收敛性 对于一种迭代过程，为了保证它是有效的，需要肯定它的收敛性，同时考 察它的收敛速度。 所谓迭代过程的收敛速度，是指在接近收敛的过程中迭代误差 的下降速度。 定义定义 设迭代过程)( 1kk xx= + 收敛于方程)(xx=的根，如果迭代误差 = kk xe，当k时成立下列渐进关系式： C e e p k k = +1 lim)0(为常数C（2-7） 则称该迭代过程是p阶收敛的。特别地，1=p时称为线性收敛，1p时称为超 线性收敛，2=p时称为平方收敛。 定理定理 对于迭代过程)( 1kk xx= + ，如果)( )( x p 在所求根的邻近连续，并且 = , 0)( , 0)()()( )( )1( p p （2-8） 则该迭代过程在点邻近是p阶收敛的。 由上述定理可知，迭代过程的收敛速度依赖于迭代函数)(x的选取。对于 Newton 公式（2-6） ，其迭代函数为 )( )( )( xf xf xx=， 2 )( )()( )( xf xfxf x=。假定 是)(xf的一个单根，即0)(=f，0)( f，则由上式知0)( =，于是依据 定理可以断定，Newton 法在根的邻近是平方收敛的。 基于 Newton 法对几类非线性方程组求解的研究 第 4页共 33 页 2.4 Newton 法的算法 解非线性方程0)(=xf的 Newton 法的算法如下： 输入：迭代初值 0 x；精度要求；最大迭代次数 N； 输出：近似解a，或方法失败的信息； Step1置 0 , 0xxk=； Step2)(),( xffxff=； Step3若0 =f，转 Step8； Step4置 f f x=； Step5若，转 Step8; Step7置1,+=kkx，转 Step2; Step8输出“方法失败”停机（此时可另选初值或增大 N 后再重新计算） ； 2.5 实例计算 例 1用 Newton 迭代法求方程01= x xe在 0.5 附近的根，要求精度为 9 10=。 解：Newton 迭代格式为 k x k k x k x x k kk x ex x exe ex xx k kk k + = + = + 1 1 1 ，2 , 1 , 0=k 取5 . 0 0 =x，计算结果见表 2-1 表 2-1 Newton 迭代法计算结果 k k x 1 kk xx 00.5 10.5710204390.071020439 基于 Newton 法对几类非线性方程组求解的研究 第 5页共 33 页 20.5671555680.003864871 30.5671432900.000012277 40.5671432900.000000000 经过 4 次迭代就得到结果，可取567143290 . 0 4 =x。 3求解非线性方程组的 Newton 迭代法 3.1解非线性方程组的 Newton 迭代法 讨论非线性方程组 = = = 0),( 0),( 0),( 21 212 211 nn n n xxxf xxxf xxxf （3-1） 的求解问题。为了便于讨论，引进向量记号，令 = n x x x x 2 1 ， = )( )( )( )( 2 1 xf xf xf xf n ， = 0 0 0 0 则上述方程组可写成向量值函数 0)(=xf（3-2） 类似单个方程的 Newton 迭代法一样，采用逐次线性化的方法，构造解方程 组的的 Newton 迭代法。在某个近似解 k x处，将向量函数)(xf做 Taylor 展开， 则有 )()()( kkk xxxfxfxf+（3-3） 从而得到方程组的近似方程组 0)()( =+ kkk xxxfxf（3-4） 其中 基于 Newton 法对几类非线性方程组求解的研究 第 6页共 33 页 k xx n n n nn n n k x f x f x f x f x f x f x f x f x f xf = = 1 2 2 2 1 2 1 2 1 1 1 )(（3-5） 称为)(xf的 Jacobi 矩阵。若)(xf的 Jacobi 矩阵非奇异，则方程组(3-4） 有唯一解，记为 , 2 , 1 , 0),()( 1 1 = + kxfxfxx kkkk (3-6) 这就是解非线性方程组（3-2）的 Newton 迭代法。 在应用Newton 法解非线性方程组(3-1）的实际计算过程中, 每一步计算 1+k x时, 一般不直接采用(3-6)计算)( k xf的逆矩阵 1 )( k xf，而是对(3-6)做 一些处理。于是( 3-6) 也可写成： 2 , 1 , 0, 0)()( 1 =+ + kxfxxxf kkkk (3-7) 因此， 1+k x便是线性方程组: 0)()( 1 =+ +kkkk xfxxxf(3-8) 的解。 由（3-8） ，令 kkk xxx= +1 ，将迭代公式改写为： = += + )()( 1 kkk kkk xfxxf xxx , 2 , 1 , 0=k 每一步迭代均需解方程组： )()( kkk xfxxf= 这是一个线性方程组，可用高斯消去法求解。 通常可用fStartValue; coutfAccuracy; coutiMax; coutfStartXValuefStartYValue; coutfAccuracy; coutiMax; cout=fAccuracy if(kiMax) return false; cout= 1 ierr=1; sol=xold; disp(increase in residual) disp(outstat) return; end % end while end sol=x; if debug=1 disp(outstat) 基于 Newton 法对几类非线性方程组求解的研究 第 30页共 33 页 end % % on failure, set the error flag % if fnrm stop_tol ierr = 1; end disp(iterration number: ,num2str(itc); function l, u =diffjac(x, f, f0) % compute a forward difference Jacobian f(x), return lu factors % % uses dirder.m to compute the columns % % C.T.Kelley, November 25, 1993 % % This code comes with no guarantee or warranty of any kind. % % % inputs: %x, f = point and function %f0= f(x), preevaluated % n=length(x); for j=1:n zz=zeros(n,1); zz(j)=1; jac(:,j)=dirder(x,zz,f,f0); end l, u = lu(jac); function z = dirder(x,w,f,f0) % Finite difference directional derivative % Approximate f(x) w % % C.T.Kelley, November 25, 1993 % % This code comes with no guarantee or warranty of any kind. % % function z = dirder(x,w,f,f0) % % inputs: %x, w = point and direction %f = function %f0 = f(x), in nonlinear iterations 基于 Newton 法对几类非线性方程组求解的研究 第 31页共 33 页 %f(x) has usually been computed %before the call to dirder % Hardwired difference increment. epsnew=1.d-7; %n=length(x); % scale the step %if norm(w) = 0 z=zeros(n,1); return end epsnew = epsnew/norm(w); if norm(x) 0 epsnew=epsnew*norm(x); end % % del and f1 could share the same space if storage % is more important than clarity % del=x+epsnew*w; f1=feval(f,del); z = (f1 - f0)/epsnew; / NewtonComp.cpp : Defines the entry point for the console application. / #include #include #include using namespace std; bool NewtonIteration(double fStartValue,double fAccuracy,int iMax); void main() double fStartValue;/初始值 x0 double fAccuracy;/精度 intiMax;/max 迭代次数 coutfAccuracy; coutiMax; couta，我们用下面式子代替 （12.1） (12.4) ( ) = = ZjxvVU nZnjUaU jjj n jx n jt , , 0, 0 其中前面的 xt 和表示差分系数，所以差分方程可理解为 h UU a k UU n j n j n j n j = + + 1 1 这时我们引入网格比，且当0, 0kh时，网络比)( 是常数CC h k =。 于是我们可以通过在 1+ = nn tt处定义 （12.5）0,)( 1 = + nZnjUaUEU n jj n k n j ， 从而将（12.4）当作是一个有限差分格式。 如果我们定义 Rx ，则在网络点 j xx=处可以写成 (12.6)RxxUahxUaxUExU nnn k n += + ),()1 ()()()( 1 通过迭代法，我们在 n tt=可以找到近似解。 RxxvExU n x n =),)()( 对于热传导方程来说，情况是类似的，即在最大范数里如果 1a，那么 k E是稳定的。因为 k E的各项系数都大于零，且相加等于 1。所以有： c c k vvE 因此有： c c n k c n vvEU= 很容易看到对于稳定性来说，1a这个条件是必须的。如前面所 讲的，稳定性就意味着是收敛的。 定理定理 12.112.112.112.1如果 n U、u由式（12.6）、（12.1）确定且假设 10a，那么则有 0 2 n c n nn tvhCtcuU，当 证明： 首先我们引入截断误差： （12.17）)()()(xuaxux n x n t n = 通过在差分方程精确解u处的泰勒展开式，我们可以得到在 ),( 1+ = nnn ttI 处有： （12.18）),()(),()()( nx n xnt n t n txauxuatxuxux+ 2 )(.,)(.,)( max C xxtt It vChtutukhC n + （其中：运用了hk， xxtt auu=， 2 )(., cc xx vtu） （12.7）也可以改写成为： )()()( 1 xkxuExu nn k n += + Rx 不妨设 nnn uUz=，因此有： nn k n kzEz= +1 或者，通过不断应用得： = = 1 0 10 n j jjn k n k n EkzEz 又因为0 000 =vuuUz，于是由稳定性条件和（12.8）推出 2 1 0 c n j c j c n vCnkhkz = ，当0 n t 到此，定理证明完毕。 另一方面，我们考虑如果0a，那么用有限差分（12.9）代替 （12.4），即： n jx n jt UaU=)0( n 或者代替（12.5），即： )0,()1 ()( 1 1 时当+= + nZjUaUaUEU n j n jj n k n j 那么，这时稳定性条件是10a。因为在（12.4）和（12.9）中都用 到流动方向，所以我们将此差分格式称之为迎风格式迎风格式。 下面，考虑更一般的显示有限差分格式，即： （12.10）)0,()( 1 时，当= + nznjUaUEU p n pjpj n k n j （其中：)( pp aa=、常数= h k 、Rx） (12.11)0,)()()( 1 = + nRxphxUaxUExU n p p n k n )()( 0 RxxvxU= 如果当0h时，有)()( 11 hOuEuk n k nn = + ，那么我们称这种方法为 Accurate of order r。u是（12.1）的精确解且常数= h k 。 如在 9.1 节时，我们做vEk的傅立叶变换，得到有： )( )( )( )()( = p ip pk eaEvhEvE 其中： 确切说，与抛物型方程一样，在 2 L范数里差分格式稳定的充分必 要条件为：)( 1)( RE 这个条件称为 von Neumann 条件。 由（12.5）有： aeaE i +=1)( 当变化时，in the complex plane with center at)( Ebelongs to a circle 对于（12.12）来说，1+=aE（以外除在m=） 于是，我们得到对于任意的，这种方法都是不稳定的。 在（12.14）中我们用平均值)()( 2 1 hxUhxU nn +代替其左边部 分，则有 h hxUhxU a k hxUhxUxU nn nnn 2 )()( )()( 2 1 )( 1 + = + + 或是 )()1 ( 2 1 )()1 ( 2 1 )()( 1 hxUahxUaxUExU nnn k n += + 这是 Friedrichs 格式的特殊形式，在下面我们将进一步探讨。 若sincos)( iaE+=，则)( 1sincos)( 2222 2 时当RaE+=当且 仅当1a时才成立。 Friedrichs 格式也可以用以下式子表示，即： )()(2)( 2 1 2 )()()()( 1 hxUxUhxU kh hxUhxU a k xUxU nnn nnnn + + = + 或 （12.15） n x x n x n t U h UaU += 2 1 这个方程可近似看作为传播因子为 h 2 1 的抛物型方程。在初始公式 （12.14）中引入人工传播后，可使得差分格式的稳定性得到保证。 （在流体动力学这也被称为人工粘性。） 注意到对于 Friedrichs 格式，)(xU n 以初始数据形式表示为： = = n nj nj n jhxvaxU)()( 所以)(xv在间隔区间 +=+ nn t x t xnhxnhx,里取值。当v的值为atx n + 时， n tt=处的精确解可由（12.2）式确定。显然，如果差分格式的依赖 区域比如 + nn t x t x,不包含精确解的依赖区域即点atx n +，那么这样的 差分方法是不可能成功的。条件如果简化为11a，那就是通常所 说的旧的稳定性条件。 对于一般形如（12.10）的式子，我们可以将其稳定性条件称为 Courant-Friedrichs-Lewy 条件（又称为 CFL 条件），即：为了使差分格 式稳定，那么在),(tx处有限差分格式的依赖区域必须包含微分方程的依 赖区域。 在（12.4）例子中，我们可以看到有限差分格式依赖区域为 n txx+,，此时 CFL 条件要求10a。另外，旧的稳定性条件是 1a，而且当0a时（12.4）不成立。类似地，在（12.4）中若0a就 不能用向后差分系数代替向前差分系数。但是，在0a的情况下，如 果1a，那么用 x 代替 x 后的（12.4）是稳定的。 有限差分格式（12.14）与 Friedrichs 格式有着相同的依赖区域，但 是由（12.14）得出的 CFL 条件对于稳定性条件来说，是不充分的，因 为它不能使得对于所有的都是稳定的。 如有限差分格式（12.14）一样，Friedrichs 格式也是有着一阶精确 解。如果（12.1）的精确解u是非常规范的，那么通过（12.15）的表现 则有： )( 2 1 2 1 2 1 2 hOu h aukuuu h uau n xx n x n tt n t n xx n x n t += )()( 2 1 22 hOuu h n xx n tt += )() 1( 2 1 222 hOua h n xx+ = 时当0h 由此可知除了特殊取值 a 1 =以外，误差估计都是一阶的，所以 在（12.13）中近似解等同于精确解。 下面，开始讨论如下形式的二阶精确解，即： )()()()()( 101 1 hxUaxUahxUaxUExU nnnn k n += + 由（12.13）可表示其稳定性条件为： )( 3 101 Oeeaaea iaii +=+ (当0时) 或者通过泰勒展开式有： )( 2 1 1 )( 2 1 )()( 3222 1111101 Oaia aaaaiaaa += + (当0时) 于是推出：1 101 =+ aaa aaa= 11 22 11 aaa=+ 即有：)( 2 1 22 1 aaa+= 22 0 1aa= )( 2 1 22 1 aaa= 所以：)()( 2 1 )()1 ()()( 2 1 )( 222222 hxUaaxUahxUaaxUE nnnn k += 如果有： sincos)( 2222 aiaaIE+= 那么通过简单的计算，可得到： （12.16） 22222 2 )cos1)(1 (1)( =aaE 于是按 2 L范数如果1 22 a，那么此差分格式是稳定的（参考习题 12.2），对于稳定性条件来说，这也与 CFL 必要条件相适应。 后一种差分格式称之为 Lax-Wendroff 格式。我们注意到这是一个 按最大范数不稳定而按 L2 范数却是稳定的例子。事实上，我们可以这 样表示： 如果 10 22 tRx 的相关问题，其中 T N uuu),( 1 =表示N维向量，而A是一个n阶矩阵。 例如：Friedrichs 格式采用这种形式表示： )()( 2 1 )()( 2 1 )()( 1 hxUAIhxUAIxUExU nnn k n += + Lax-Wendroff 格式 )()( 2 1 )()()()( 2 1 )()( 2222221 hxUAAxUAIhxUAAxUExU nnnn k n += + 于是这些算子变成了矩阵周期函数，分别对应为： sincos)( AiIE+= sincos)( 2222 AiAAIE+= 通过对矩阵A作相同的正交变换，很容易看出 Friedrichs 格式和 Lax- Wendroff 格式稳定的条件是1A（其中：在 N R里，A是指矩阵范数 除以 Euclidean 范数即： j Nj vv Av A , 1 0 maxsup = =，其中 j 是指A的特征 值）。 考虑以下这样一个方程组的初值问题 (12.19) 2 2 2 2 2 x w a t w = + RRtx),( 0 )0(.,ww= 1 )0(.,w t w = Rx 通过令 x w au = 1 ， t w u = 2 ，则二阶双曲型方程可化简。 这些函数满足 x u a t u = 21 ， x u a t u = 12 ， 即：令 T uuu),( 21 =, 则有： (12.20) x u a a t u = 0 0 + RRtx),( = 1 0 )0(., w aw uRx 相反地，（12.19）的解由（12.20）的解确定。 现任选差分格式（12.17）和（12.18）运用到该方程组中。由于 （12.20）的矩阵A的特征值为a，所以可确定稳定性条件为1a。 下面推广到多维空间，讨论对称双曲型方程组（或者是Friedrichs 方程组），即： （12.21） = = d j j j x u A t u 1 + RRtx),( vu=)0(., d Rx （其中：u是指N维向量函数， j A表示n阶对称矩阵） 回顾 11.4 节，其相关的初值问题都是按 2 L范数讨论的. 下面我们探讨相关的有限差分算子 （12.22）)()()( 1 hxUaxUExU nn k n = + （其中：),( 1d =表示整数向量， a指N维常数矩阵） 当 d d R=),( 1 时，符号矩阵为： = i eaE)( （ dd += 11 ） 经过傅立叶变换后，有： )()( )()( 1 nn UhEU= + 0n 所以得出： )( )( )()(vhEU n = 于是很容易得到，差分格式按 2 L范数稳定性的充分必要条件是符号矩 阵范数满足： （12.23）CE n )( (当 d Rn, 0时) 与一阶标量方程对比，矩阵范数CAn（0n）表示1A是不正 确的。比如： （12.24） = n nn n a naa a a 00 1 1 ，但是1a。 当然，如果是（12.23）的形式，则对于)( E的每一个特征值)(j来 说，有C n j )(（0,nRd），从而 1)(j（ d R） 我们将这称之为 von Neumann 的稳定性条件，这也是有限差分格式按 2 L范数稳定性的必要条件。但是不是充分条件，比如（12.24）中当 1=a时所阐述的一样。 有限差分格式按 2 L范数稳定性的充分条件是： 1)( E（ d R） 为了能够建立一个稳定的有限差分格式，需要以下结论。 引理引理 12.112.112.112.1假设 a是一个半正定对称矩阵，且1= a 则有： 1)( = i eaE（ d R） 证明：令 = = N j jjv uvu 1 ,（ N Cvu,） 由于 a是一个半正定实对称矩阵，则双线性形式vua, 满足： uvavua, =（0,uua N Cvu,） 运用这些特性，很容易证明一般情况下的Cauchy-Schwarz 不等式： 2 1 2 1 ,vvauuavua （标准的 Cauchy-Schwarz 不等式的证明：这个不等式只是一个概括， 因为 2 1 ,uua只是一个半范数。）因此，应用不等式 22 2baab+得到： vvauuavua, 2 1 , 2 1 , + 从而：wveawvE i ,)(