第5章非线性方程(组) 非线性方程(组)是科学技术与工程计算中的重要问题之一，求解方法远比线性方程(组)复杂，通常都要使用数值迭代技术。本章主要介绍二分法、牛顿法、拟牛顿法和不动点法，并讨论它们的收敛性和计算复杂性。这些迭代技术在实践中非常重要，应用十分广泛。阅读本章的数学基础是初等微分学。 5.1非线性方程 5.1.1二分法 定义5.1.1如果f(ξ)=0，则称ξ是方程f(x)=0的一个根(解)，或者说是函数y=f(x)的一个零点。若f(x)是非线性函数，则称f(x)=0为非线性方程。 二分法是求解非线性方程最简单的实用方法，其理论基础源于微分学中的介值定理和实数连续性的闭区间套定理。 介值定理： 如果f(x)在闭区间［a，b］上连续，则对f(a)和f(b)之间的任何一个数y，都存在ξ∈［a，b］使得f(ξ)=y，即f(x)能取到f(a)和f(b)之间的所有数。 闭区间套定理： 设［an，bn］［an+1，bn+1］，n=1，2，…，是一个闭区间套，若limn→∞(bn-an)=0，则存在唯一ξ∈［an，bn］，n=1，2，…，即长度趋于零的闭区间套的交是单点集，{ξ}=∩∞n=1［an，bn］。 若f(x)是闭区间［a，b］上的连续函数，且在端点异号f(a)·f(b)<0，则0必在f(a)和f(b)之间，根据介值定理存在ξ∈［a，b］使得f(ξ)=0，即方程f(x)=0在［a，b］上有根。 在［a，b］上找根： 首先取中点c=(a+b)/2，如果f(c)=0则找到一个根c； 否则比较f(c)和f(a)的符号： 若f(a)·f(c)<0，令a1=a，b1=c，否则a1=c，b1=b。于是，得到一个新区间［a1，b1］［a，b］，且 b1-a1=b-a2，f(a1)·f(b1)<0 在［a1，b1］上重复前面的操作，并将这个过程继续下去，要么得到一个根c； 要么得到一个区间套 ［a，b］［a1，b1］…［an，bn］… bn-an=b-a2n→0(n→∞) f(an)·f(bn)<0 在第二种情况下，由闭区间套定理存在ξ∈［an，bn］，且f(ξ)在f(an)和f(bn)之间，因此必有f(ξ)=0。不难看出 an+bn2-ξ≤b-a2n+1(5.1) 因此，根是闭区间套中点的极限 ξ= limn→∞an+bn2 这就是二分法的寻根过程，首先找有根区间，然后二分区间将寻根范围缩小到原区间的二分之一，继续这个过程直至达到期望精度，得到近似根ξ≈(an+bn)/2。在寻根过程中，不需要观察函数曲线的性态，只须计算函数值yn=f((an+bn)/2)，算法极为简单。以下是二分法的概述。 二分法 给定初始区间［a，b］，f(a)·f(b)<0； 控制误差δ (1) (b-a)/2>δ，置c=(a+b)/2。若f(c)=0，输出根c； 否则 若f(a)·f(c)<0，令b=c，否则a=c； (2) (b-a)/2≤δ，输出近似根c=(a+b)/2，否则返回(1)。 例5.1.1利用二分法求方程f(x)=x3-x-1=0在［1，2］上的根。 解f(1)·f(2)=-5<0，因此在［1，2］上方程有根，表5.1给出了二分法前9次迭代的计算结果，图5.1是函数f(x)=x3-x-1曲线和最初4次迭代区间端点的图示。 表5.1二分法寻方程x3-x-1=0的根： 前9次迭代 nanbncn=(an+bn)/2f(cn) 01.02.01.50.125 11.0 1.51.25-0.609375 21.251.51.375-0.291016 31.3751.51.4375-0.0959473 41.43751.51.468750.0112 51.43751.468751.45313-0.0431938 61.453131.468751.46094-0.0162034 71.460941.468751.46484-0.00255352 81.464841.468751.46680.00431024 91.464841.46681.465820.00087512 图5.1曲线f(x)=x3-x-1和二分法前4次迭代的区间端点 二分法第n次迭代的近似根cn=(an+bn)/2，与真值之间的误差为 |cn-ξ|≤(b-a)/2n+1(5.2) 迭代过程中只需计算(n+1)次函数值，因而计算量极小。每次迭代都以常数因子1/2减小误差，因此二分法是线性收敛的，或者说二分法线性收敛。 如果取控制误差δ=0.5×10-p，即期望根精确到小数点后p位，则从(b-a)/2n+1<0.5×10-p估算出所需要的迭代次数 n>plog210+log2(b-a)≈3.3219p+log2(b-a)(5.3) 实际迭代次数要小于这个估计数。例如： 在例5.1.1中，如果期望根精确到小数点后3位，则迭代次数n>3.3219p+log21=9.99657，即迭代10次后可以使根精确到小数点后3位。事实上，这个方程的根是1.46557(精确到小 图5.2二分法失效 数点后5位)，从表5.1可以看出第9次迭代根就已经精确到小数点后3位。 二分法的困难部分也许是确定端点异号的初始区间，一旦这种区间被找到二分法就保证迭代一定收敛到根。值得指出的是，二分法不适合重根的情况，对图像如y=x2m的一类问题失效，如图5.2所示。下面介绍的牛顿法可解决这类问题，收敛速度也比二分法快。 5.1.2牛顿法 牛顿法，也称牛顿拉夫逊(NewtonRaphson)方法，它的收敛速度比二分法快，求根的迭代过程有明显的几何意义，其原理是用一系列线性方程的根逼近非线性方程的根。 为求方程f(x)=0的根，先给定一个初始猜测值(简称初始值)x0，然后过点(x0，f(x0))作曲线y=f(x)的切线 y=f(x0)+f′(x0)(x-x0) 它与x轴的交点x1，如图5.3所示，是线性方程f(x0)+f′(x0)(x-x0)=0的解 x1=x0-f(x0)f′(x0) 图5.3牛顿迭代几何 再以x1代替x0重复上面的操作，得到 x2=x1-f(x1)f′(x1) 将这个过程继续下去，就得到牛顿迭代点列 xk+1=xk-f(xk)f′(xk)，k=0，1，2，… 可见，它是一系列线性方程f(xk)+f′(xk)(x-xk)=0，k=0，1，2，…的解。 牛顿法 解方程f(x)=0 给定初始值： x0 迭代： xk+1=xk- f(xk)f′(xk),k=0，1，2，… 由Taylor公式也可导出牛顿法。由Taylor定理，在x和x0之间存在ηx使得 f(x)=f(x0)+f′(x0)(x-x0)+f″(ηx)2(x-x0)2 当x在x0附近时，二次项很小可忽略不计，即 f(x)≈f(x0)+f′(x0)(x-x0) 当x0在方程f(x)=0的根附近时，自然希望下次迭代x1更靠近方程的根，因此希望x1满足 f(x0)+f′(x0)(x1-x0)=0(5.4) 于是 x1=x0-f(x0)f′(x0) 满足式(5.4)的x1只是比x0更接近于根ξ，但未必有f(x1)=0。为了找到更精确的解，需要继续这个过程，因此就产生了逐渐趋近于根ξ的牛顿迭代，xk+1=xk-f(xk)/f′(xk)，k=0，1，2，…。 和二分法一样，牛顿法也是通过迭代逼近方程的根，不同的是牛顿法将导数用在迭代过程中，从而充分利用了函数的局部变化信息(导数是函数的瞬时变化率)，因此收敛速度比二分法更快。很可惜，牛顿迭代不总是收敛的。 例5.1.2从给定的初始值开始对下述方程作8次牛顿迭代，并观察收敛性： (i) x2-2=0，x0=2.5 (ii) sinx=0，x0=0.8 (iii) x3-x-3=0，x0=0.0 (iv) x1/3=0，x0=1.0 解这四个方程的牛顿迭代分别为 (i) xk+1=xk-2x2k-12xk，x0=2.5 (ii) xk+1=xk-sinxkcosxk，x0=0.8 (iii) xk+1=xk-x3k-xk-33x2k-1，x0=0.0 (iv) xk+1=xk-x1/3kx-2/3k/3，x0=1.0 表5.2给出了8次迭代结果，可以看出： 对于x2-2=0和sinx=0牛顿法收敛，且第4次迭代就分别给出精确到小数点后8位和24位的根。对于方程x3-x-3=0，从初始值x0=0.0开始，迭代3次后进入循环状态，因此对给定的初始值牛顿迭代不收敛。对于x1/3=0，牛顿迭代对给定的初始值是发散的。 表5.2牛顿迭代 x2-2=0sinx=0x3-x-3=0x1/3=0 2.50.80.000001.00000 1.65-0.2296385570503640-3.00000-2.00000 1.4310606060606060.00412357916974798-1.961534.00000 1.414312727593564-0.0000000233724753-0.0065716.0000 1.4142135658496033.308722450212×10-24-3.00038-32.0000 1.4142135623730950.-1.9618164.0000 1.4142135623730950.-1.14743-128.000 1.4142135623730950.-0.00725256.000 1.4142135623730950.-3.00047-512.000 方程x3-x-3=0精确到小数点后10位的根是ξ=1.6716998816，选择充分靠近根ξ的点作为初始值，比如x0=1.0，你会发现牛顿迭代收敛到根ξ，即从初始值x0=1.0开始牛顿法收敛。这说明牛顿法的收敛性依赖于初始值的选择，有时在充分靠近根的情况下才收敛，也就是说牛顿法局部收敛。最特别的是方程x1/3=0，不管如何选取初始值，无论多么靠近方程的根(只要不等于根ξ)，经过若干次迭代后将逐渐偏离方程的根，也就是说牛顿法对这类方程失效。 1. 二阶收敛性 牛顿法何时收敛，收敛速度如何？这是下面要讨论的问题，先给出收敛速度的定义。 定义5.1.2假定ξ是方程f(x)=0的一个根，{xk: k=1，2，…}是收敛到ξ的迭代点列。如果存在常数c使得 limk→∞|xk+1-ξ||xk-ξ|p=c(5.5) 则称迭代最终收敛速度是p阶的，或者说迭代p阶收敛。特别地，p=1，称为线性收敛； p=2，称为平方收敛； 1cξ时必存在δ使得， x，η∈(ξ-δ，ξ+δ)，f″(η)2f′(x)0，收敛到2； 取x0<0，收敛到-2。 牛顿法不像二分法那样，可以通过估计绝对误差来终止迭代，在实际问题中需要设置终止条件。给定容差ε，常用的终止条件有 绝对终止条件： |xk+1-xk|≤ε。 相对终止条件： |xk+1-xk|xk+1≤ε (根不在0附近)。 绝对/相对终止条件： |xk+1-xk|max{|xk+1|，θ}≤ε，其中θ>0，常用于根在0附近。 2. 线性收敛性 在定理5.1.1中，f′(ξ)≠0是牛顿法收敛的充分条件，但不是必要条件。例如方程xp=0，零是它的p重根，即f(0)=f′(0)=…=f(p-1)(0)=0，f(p)(0)=p≠0，牛顿迭代 xk+1=xk-xpkpxp-1k=p-1pxk 以(p-1)/p的因子减小，因此收敛到方程的根。从(xk+1-0)/(xk-0)=(p-1)/p知，此时牛顿法的收敛速度是线性的，不再有二阶收敛性。下述定理概括了这类问题。 定理5.1.2假定f(x)有连续的p+1阶导数，ξ是f(x)=0的p+1重根，即 f(ξ)=f′(ξ)=…=f(p)(ξ)=0，f(p+1)(ξ)≠0 若x0充分接近ξ，则牛顿迭代收敛于ξ，且 limk→∞xk+1-ξxk-ξ=pp+1(5.8) 即牛顿迭代的收敛速度是线性的。 证明在牛顿迭代xk+1=xk-f(xk)f′(x1)的两边同减去ξ，得到 xk+1-ξ=xk-ξ-f(xk)f′(xk)(5.9) 由重根条件，对f(xk)和f′(xk)在点ξ作Taylor展开 f(xk)=f(p+1)(ηk)(p+1)!(xk-ξ)(p+1)，f′(xk)=f(p+1)(η′k)p!(xk-ξ)p 其中ηk，η′k位于ξ，xk之间。将上式代入式(5.9)，得到 xk+1-ξ=(xk-ξ)1-1p+1·f(p+1)(ηk)f(p+1)(η′k)(5.10) 由于f(x)是连续p+1阶可导的，且f(p+1)(ξ)≠0，所以取0<ε<1/(p+1)，一定存在δ>0使得当η，η′∈(ξ-δ，ξ+δ)时， 1-1p+1·f(p+1)(η)f(p+1)(η′)<1-ε(5.11) 其中0<1-ε<1。因x0充分接近ξ，不妨假定x0∈(ξ-δ，ξ+δ)，根据式(5.10)和式(5.11)得到 |x1-ξ|=(x0-ξ)1-1p+1·f(p+1)(ηk)f(p+1)(η′k)<(1-ε)|x0-ξ|<δ 从而x1∈(ξ-δ，ξ+δ)，同理可推知所有xk∈(ξ-δ，ξ+δ)，因此 |xk+1-ξ|<(1-ε)|xk-ξ|<…<(1-ε)k+1|x0-ξ| 因(1-ε)k+1→0(k→∞)，必有xk+1-ξ→0。再由式(5.10)， limk→∞xk+1-ξxk-ξ= limk→∞1-1p+1·f(p+1)(ηk)f(p+1)(η′k)=pp+1 证毕。 注释对于p+1重根牛顿法线性收敛，如果利用下面的迭代 xk+1=xk-(p+1)f(xk)f′(xk)，k=1，2，… 则恢复到牛顿法的二阶收敛性。因事先难以确定根的重数，这种迭代方式仅有理论意义而无实用价值，对此不作详细讨论。 5.1.3拟牛顿法 牛顿法对单根二阶收敛，对重根线性收敛，但在迭代过程中必须计算导数。实际问题的函数可能非常复杂，或者根本没有解析表达式，此时牛顿法需要花费很大代价来计算导数，自然希望有不需计算导数的迭代方法来解决此类问题。拟牛顿法是避免计算导数的一类方法，它的一般迭代形式为 xk+1=xk-f(xk)Tk 这里Tk是f′(xk)的近似Tk≈f′(xk)。选择不同的Tk导致不同的拟牛顿法，其中最著名的是下面将要介绍的割线法，其收敛速度介于线性与二阶之间。 1. 割线法 众所周知，导数是切线的斜率，而切线在几何上是割线的极限，因此导数可用割线的斜率来近似。通过(xk-1，f(xk-1))和(xk，f(xk))两点的割线斜率 图5.4割线法迭代几何 Tk=f(xk)-f(xk-1)xk-xk-1=f［xk-1，xk］ 替代牛顿法中f′(xk)，得到割线法的迭代 xk+1=xk-f(xk)f［xk-1，xk］(5.12) 事实上，xk+1是割线与x轴的交点，如图5.4所示。与牛顿法不同，实施割线法需要两个初始值。 割线法 解方程f(x)=0 给定初始值： x0，x1 迭代： xk+1=xk-f(xk)f［xk-1,xk］,k=1,2,… 在讨论割线法的收敛性和收敛速度之前，先看一个例子。继续考察方程x2-2=0，此次使用割线法，初始值取x0=2.5，x1=1.5。按迭代公式(5.12)，最初6次迭代的近似值及其误差如下： x2=1.4375e2=0.0232864 x3=1.41489361e3=0.0006800 x4=1.41429114e4=5.551×10-6 x5=1.41421356e5=1.334×10-9 x6=1.41421356e6=2.664×10-15 观察误差可以发现，割线法的收敛速度快于线性收敛速度，但比二阶收敛速度慢。根据下面的定理，准确收敛阶p=(1+5)/2 ≈1.618，位于线性与二阶中间偏上一些。 定理5.1.3假定f(x)有连续的二阶导数，ξ是方程f(x)=0的根，f′(ξ)≠0。若x0，x1充分接近ξ，则割线法收敛到ξ，且 limk→∞ek+1ekek-1=cξcξΔf″(ξ)2f′(ξ)(5.13) 证明由式(5.12)，得到 xk+1-ξ=xk-ξ-f(xk)(xk-ξ)-(xk-1-ξ)f(xk)-f(xk-1) 因此 ek+1=ek-f(xk)ek-f(xk)ek-1f(xk)-f(xk-1)=ekek-11ek-1-f(xk)ek-1-f(xk)ekf(xk)-f(xk-1) =ekek-1·f(xk)/(xk-ξ)-f(xk-1)/(xk-1-ξ)xk-xk-1A(xk-1，xk)·xk-xk-1f(xk)-f(xk-1)B(xk-1，xk)(5.14) 当xk，xk-1→ξ时 B(xk，xk-1)→1f′(ξ)(5.15) 考虑A(xk，xk-1)在xk，xk-1→ξ时的极限： 定义h(x)=f(x)x-ξ，则 A(xk，xk-1)=h(xk)-h(xk-1)xk-xk-1 当xk，xk-1→ξ时，A(xk，xk-1)→h′(ξ)。下面计算h′(ξ)： h′(x)=f′(x)(x-ξ)-f(x)(x-ξ)2 根据求极限的洛必达法则 h′(ξ)= limx→ξf′(x)(x-ξ)-f(x)(x-ξ)2= limx→ξf″(x)(x-ξ)2(x-ξ)=f″(ξ)2 因此 A(xk，xk-1)→f″(ξ)2(5.16) 由式(5.15)和式(5.16)， A(xk，xk-1)·B(xk，xk-1)→f″(ξ)2f′(ξ)=cξ(5.17) 于是，给定ε>0必存在δ>0使得当xk，xk-1∈(ξ-δ，ξ+δ)时， A(xk，xk-1)·B(xk，xk-1)1 2(Δ)=［0.5848，0.7937］Δ，maxx∈Δ|′2(x)|=0.9746<1 3(Δ)=［0.6097，0.8］Δ，maxx∈Δ|′3(x)|=0.6495<1 4(Δ)=［0.6823，0.7142］Δ，maxx∈Δ|′3(x)|=0.1756<1 可见，收敛的迭代函数都有共同的质性(Δ)Δ，max|′(x)|<1，发散的迭代函数没有这个特性； max|′(x)|越小收敛越快。这一现象可由后面的定理得到进一步证实。 1. 整体收敛性 定义5.1.4如果存在常数ρ<1使得 |(x)-(y)|≤ρ|x-y|，x，y∈［a，b］ 则称(x)是［a，b］上压缩映射。 如果MΔmaxa≤x≤b|′(x)|<1，则(x)必是［a，b］上压缩映射。压缩映射缩小任意两点间的距离，因此在［a，b］上一致连续。 定理5.1.4(整体收敛的充分条件)如果(x)是［a，b］上的压缩映射，且(［a，b］)［a，b］，则 (i) 在［a，b］上(x)存在唯一的不动点ξ (ii) 对任意x0∈［a，b］为初始值的迭代xk+1=(xk)∈［a，b］，且 |xk+1-ξ|≤ρk+11-ρ|x1-x0|(5.21) (iii) xk→ξ(k→∞) 证明令Φ(x)=x-(x)。因(［a，b］)［a，b］，即x∈［a，b］，a≤(x)≤b，所以Φ(a)≤0，Φ(b)≥0。由压缩映射的连续性和介值定理，在［a，b］上存在ξ使Φ(ξ)=0，即ξ是(x)的不动点ξ。如果还存在一个不动点η≠ξ，则 |ξ-η|=|(ξ)-(η)|≤ρ|ξ-η|ρ≥1 与ρ<1矛盾，故不动点是唯一的。 |xk+1-ξ|=|(xk)-(ξ)|≤ρ|xk-ξ|≤ρ(|xk+1-xk|+|xk+1-ξ|) 由此可推知， |xk+1-ξ|≤ρ1-ρ|xk+1-xk|(5.22) 此外 |xk+1-xk|=|(xk)-(xk-1)|≤ρ|xk-xk-1||xk+1-xk|≤ρk|x1-x0| 代入式(5.22)， |xk+1-ξ|≤ρk+11-ρ|x1-x0| 最后，因ρ<1，得到limk→∞|xk+1-ξ|=0。证毕。 称式(5.21)和式(5.22)为迭代误差的事先估计和事后估计。事先估计是迭代开始时对第k+1迭代误差的一个粗略估计，以便为终止迭代提供信息； 事后估计是当前迭代的误差估计，比事先估计准确。给定精度ε，由事先估计 ρk+11-ρ|x1-x0|=εk+1=ρ-1ln1-ρ|x1-x0| 得到需要迭代的次数。当然，这必须以已知ρ为前提，在ρ未知情况下仍使用5.1.2节介绍的终止条件。 定理5.1.5(整体发散的充分条件)假定存在常数ρ>1使得|(x)-(y)|≥ρ|x-y|(x，y∈［a，b］)，或者mina≤x≤b|′(x)|>1。若初始点x0不是不动点，则迭代发散。 证明略，留给读者。 例5.1.5下述三个方程在［2，3］上都与方程x3-2x-5=0同解，试分析不动点迭代的收敛性： (i) x=x3-x-5 (ii) x=(2x+5)1/3 (iii) x=(2+5/x)1/2 解 (i) (x)=x3-x-5，′(x)=3x2-1，min2≤x≤3|′(x)|=|′(2)|=11>1，由定理5.1.5，迭代发散。 (ii) (x)=(2x+5)1/3，′(x)=2(2x+5)-2/3/3 min2≤x≤3(x)=(2)=91/3>2，max2≤x≤3(x)=(3)=111/3<3  (［2，3］)［2，3］ max2≤x≤3|′(x)|=′(2)=2·9-2/3/3≈0.15408<1 根据定理5.1.4，在区间［2，3］上迭代整体收敛。 (iii) (x)=(2+5/x)1/2，′(x)=-5(2+5/x)-1/2/2x min2≤x≤3(x)=(3)=(11/2)1/2<2，max2≤x≤3(x)=(2)=(9/2)1/2<3  (［2，3］)［2，3］ 因此，不满足定理5.1.4条件。如果将区间缩小到［2，2.5］，就有 (［2，2.5］)［2，2.5］ max2≤x≤2.5|′(x)|=′(2)=5·2-1/2/12<1 由定理5.1.4，迭代在区间［2，2.5］上整体收敛。事实上，在区间［2，3］上迭代也是整体收敛的，这与定理5.1.4不矛盾，因它给出的仅是收敛性的充分条件。 2. 局部收敛性 如果初始值充分靠近不动点，则迭代的收敛性仅依赖于迭代函数的局部性质，对此有如下定理： 定理5.1.6(局部收敛的充分条件)若(x)在不动点 ξ的某个邻域有连续的导数，且|′(ξ)|<1，则迭代局部收敛于 ξ，也就是说当x0充分靠近不动点时，迭代xk+1=(xk)必收敛到 ξ。 证明因′(x)在包含 ξ的邻域内连续且′(ξ)<1，所以存在δ>0使得 ′(x)<1，x∈［ξ-δ，ξ+δ］ 记 ρ=max{|′(x)|:x∈［ξ-δ，ξ+δ］}<1 由微分中值定理，x，y∈［ξ-δ，ξ+δ］， (x)-(y)=′(η)(x-y)|(x)-(y)|≤ρ|x-y| 因此 |(x)-ξ|=|(x)-(ξ)|≤ρ|x-ξ|<δ(［ξ-δ，ξ+δ］)［ξ-δ，ξ+δ］ 根据定理5.1.4，迭代xk+1=(xk)在［ξ-δ，ξ+δ］上收敛到 ξ。证毕。 定理5.1.7(局部发散的充分条件)假定(x)在不动点 ξ的某个邻域有连续的导数且|′(ξ)|>1，则若x1≠x0(即x0不是不动点)，则迭代局部发散。 例5.1.6不难验证： 迭代函数1(x)=x2-1，2(x)=2x-(1+5)/2，3(x)=1+x有同一个不动点 ξ=(1+5)/2≈1.61803 由于′1(ξ)=1+5>1，′2(ξ)=2>1，|′3(ξ)|=1/2<1，根据定理5.1.7，对于1(x)，2(x)迭代都是局部发散的； 根据定理5.1.6对于3(x)迭代收敛。两个发散迭代的性态不相同，对于1(x)迭代最终进入周期状态； 而对于2(x)随着迭代次数的增加离不动点越来越远。表5.5给出了最初15次的迭代结果。 表5.5不动点局部迭代(初始值： 1.60) xk+1=1(xk)xk+1=2(xk)xk+1=3(xk) 1.5600000000000001.5819660112501051.612451549659710 1.4336000000000011.5458980337503151.616307999627456 1.0552089600000041.4737620787507361.617500540843018 0.1134659492642901.3294901687515771.617869135883066 -0.9871254783575531.0409463487532601.617983045610511 -0.0255832899773710.4638587087566261.618018246377497 -0.999345495273933-0.6903165712366411.618029124081979 -0.001308581075696-2.9986671312231771.618032485484138 -0.999998287615568-7.6153682511962501.618033524215162 -3.4247659311814×10-6-16.848770491142391.618033845200761 -0.999999999988270-35.315574971034681.618033944390772 -2.3458124331909×10-11-72.249183930819271.618033975042172 -1.-146.11640185038841.618033984513975 0. -293.85083768952671.618033987440923 -1.-589.31970936780341.618033988345400 最后，以不动点迭代的收敛阶定理结束本节。 定理5.1.8如果(x)在不动点 ξ的某个邻域有连续的p阶导数，且 ′(ξ)=″(ξ)=…=(p-1)(ξ)=0，p(ξ)≠0 则迭代p阶局部收敛于 ξ。 证明由′(ξ)=0，迭代必局部收敛于 ξ。对(x)在点 ξ应用Taylor定理： xk=(xk-1)=(ξ)+′(ξ)2(xk-1-ξ)+…+(p-1)(ξ)(p-1)!(xk-1-ξ)p-1+ (p)(ηk-1)p!(xk-1-ξ)p 于是 xk=ξ+(p)(ηk-1)p!(xk-1-ξ)p 因此 |xk-ξ||xk-1-ξ|p=(p)(ηk-1)p!→(p)(ξ)p!(k→∞) 故迭代p阶收敛。证毕。 5.2非线性方程组 非线性方程组是指由n个n元非线性方程构成的方程组，一般形式为 f1(x1，x2，…，xn)=0 f2(x1，x2，…，xn)=0  fn(x1，x2，…，xn)=0 (5.23) 记 f(x)=（f1(x)，f2(x)，…，fn(x)）T，x=(x1，x2，…，xn)T 它是Rn到Rn的映射，称为向量值函数。利用这个记号，将非线性方程组写成简洁的向量形式 f(x)=0(5.24) 求解非线性方程组远比求解一元非线性方程困难，本节的主要目的是将上节中的牛顿法、拟牛顿法和不动点迭代法拓广到非线性方程组。 5.2.1多元牛顿法 1. 向量函数的微分 为了将牛顿法推广到非线性方程组，需要类似一元函数导数的向量值函数的导数概念，尤其是向量值函数的一阶Taylor展开。 如果对每个1≤i≤n，函数fi(x)在点x0可微，则称向量值函数f(x)在点x0可微。矩阵 Df(x0)=f1x1(x0)f1x2(x0)…f1xn(x0) f2x1(x0)f2x2(x0)…f2xn(x0)  fnx1(x0)fnx2(x0)…fnxn(x0) 称为f(x)在点x0的Jacobi矩阵，在向量值函数分析中它相当于一元函数的导数。由Jacobi矩阵定义的线性映射 L(d)=Df(x0)d(5.25) 称为f(x)在点x0的微分，它是f(x0+d)-f(x0)的局部线性近似。也就是说，在x0的邻域内 f(x0+d)-f(x0)≈L(d)=Df(x0)d 即 f(x0+d)≈f(x0)+Df(x0)d(5.26) 例如： 令f(x)=x21+4x22-4 4x21+x22-4，(f1(x)=x21+4x22-4，f2(x)=4x21+x22-4)，它在点x的Jacobi矩阵是 Df(x)=f1x1(x)f1x2(x) f2x1(x)f2x2(x)=2x18x2 8x12x2 在点x0=(1，1)T的微分 L(d)=Df(1，1)d=28 82d1 d2 由式(5.26)，得到f(1.1，1.1)的近似值 f(1.1，1.1)≈f(1，1)+28 820.1 0.1=2 2 而精确值f(1.1，1.1)=(2.05，2.05)T。可见，f(x)在点x0邻域内可由该点微分得到很好的线性近似。 回忆n元函数在点x0的一阶Taylor展开 fi(x)=fi(x0)+Tfi(x0)(x-x0)+ri(x0，x) 其中 fi(x0)=fix1(x0)，fix2(x0)，…，fixn(x0)T 是fi(x)在点x0的梯度； ri(x0，x)=o(‖x0-x‖2) 是‖x0-x‖2的无穷小量，即 lim‖x-x0‖2→0ri(x0，x)‖x-x0‖2=0 因此 f(x)=f(x0)+ Tf1(x0) Tf2(x0)  Tfn(x0)(x-x0)+r1(x0，x) r2(x0，x)  rn(x0，x) =f(x0)+Df(x0)(x-x0)+r(x0，x)(5.27) 其中r(x0，x)=o(‖x0-x‖2)是‖x0-x‖2的无穷小向量，即 lim‖x-x0‖2→0‖r(x0，x)‖2‖x-x0‖2=0 公式(5.27)称为向量值函数的一阶Taylor展开式。 2. 多元牛顿法 求解非线性方程组的牛顿法是一元方程牛顿思想的自然推广。用Jacobi矩阵Df(xk)替代一元牛顿法中的导数f′(xk)，得到多元牛顿迭代 xk+1=xk-D-1f(xk)f(xk)，k=0，1，2，…(5.28) 多元牛顿法 解方程f(x)=0 给定初始值： x0 迭代： Df(xk)dk=-f(xk) xk+1=xk+dk，k=0，1，2,… 已知第k次迭代xk，自然希望下次迭代xk+1更靠近方程组f(x)=0的解。由一阶Taylor展开式 0≈f(xk+1)=f(xk)+Df(xk)(xk+1-xk)+r(xk+1，xk) ≈f(xk)+Df(xk)(xk+1-xk) 于是，我们希望 f(xk)+Df(xk)(xk+1-xk)=0 故，第k+1次迭代 xk+1=xk-D-1f(xk)f(xk) 在求解大型方程组中，为避免计算Jacobi矩阵的逆，可先解线性方程组(称为牛顿方程) Df(xk)dk=-f(xk)（dk=xk+1-xk） 得到dk，再实现牛顿迭代xk+1=xk+dk。牛顿法在每步迭代的计算量是： 计算n个分量函数的值，n2个偏导数(Jacobi矩阵)和解一次牛顿方程，因此有O(n3)次运算。 例5.2.1用牛顿法解方程组 f(x)=x21+4x22-4 4x21+x22-4=0 0(5.29) 解计算f(x)在点x的Jacobi矩阵 Df(x)=2x18x2 8x12x2 初始值取x0=6.5 5.5，Df(x0)=1344 5211，f(x0)=159.25 195.25。从牛顿方程Df(x0)d0=-f(x0)得到 d0=-829/260 -101/110 因此 x1=x0+d0=6.5 5.5+-829/260 -101/110=3.311538461538461 2.822727272727272 表5.6给出了最初10次迭代结果。这个方程组有四个解x1=±2/5，x2=±2/5，分别以初始值x0=(-6.5，5.5)、(-6.5，-5.5)和(6.5，-5.5)进行迭代可得到另外三个解，但这样做没有必要，因为根据方程组的平方项特性，从一个解就可确定其他三个解。 表5.6多元牛顿迭代 kxk 0{6.5，5.5} 1{3.3115384615384613， 2.8227272727272723} 2{1.7765590100955955， 1.5530705606792563} 3{1.1134338611524544， 1.0340895796272630} 4{0.9159659295118101， 0.9038584739280965} 5{0.8946804303627222， 0.8944763962130312} 6{0.8944272268396262， 0.8944271923533070} 7{0.8944271909999166， 0.8944271909999160} 8{0.8944271909999159， 0.8944271909999159} 9{0.8944271909999159， 0.8944271909999159} 10{0.8944271909999159， 0.8944271909999159} 从表5.6可以看出，第6次迭代得到精确到小数点后6位的解。此后，收敛速加快，第8次迭代精度准确到小数点后16位，呈现出局部二阶收敛性特征。 事实上，如果f(x)可微且在解处的Jacobi矩阵Df(x*)可逆，则牛顿法超线性收敛； 如果Df(x)在解的局部邻域内还满足李普希兹条件‖Df(x)-Df(y)‖2