非线性方程组数值解法——Newton法遇奇异点的处理

杨家岭

(中国矿业大学 理学院,江苏 徐州 221008)

设映像F：D⊂Rn→Rn,非线性方程组F(x)=0的数值解法——Newton法，xk+1=xk-[F'(xk)]-1F(xk),k=0,1,…中，当F'(xk)奇异或病态时迭代就不能进行，以下给出三种避免F'(xk)奇异的方法.

1 经典的引入阻尼因子λk法与KaHTOPOBHY定理

1.1 引入阻尼因子λk

引入参数λk使F'(xk)+λkI非病态，此时得到迭代序列

xk+1=xk-[F'(xk)+λkI]-1F(xk),k=0,1,…

(1)

λk称阻尼因子，只要选取λk足够大使得F' (xk)+λkI成为对角占优就可以消除奇异性，且关于(1)的局部收敛定理如下.

定理1 设x*是非线性方程组F(x)=0的解：F:D⊂Rn→Rn在x*的邻域S0⊂D上的连续可微且F' (x*)非奇异，又设μ1,…，μn为矩阵F' (x*)的特征值，令

若不存在Reμi>0或Reμi<0的特征值可分取β=+∞或η=+∞，则对任意λ∈(-β,η)，由(1)迭代产生的序列

{xk}⊂S0,

且收敛于x*.

证明 若对固定的λ∈(-β,η)，在S0上定义

A(x)=F'(x)+λI,

由于F'(x)在S0上连续，故A(x)在x*处连续，下面证明A(x*)非奇异，假定A(x*)奇异，则必存在μj使得μj+λ=0,λ≠0,若λ>0，因λ∈(-β,η)故有

这是不可能的，若λ<0，则μj>0，故有

这也不可能，这矛盾说明A(x*)非奇异,于是映像

G(x)=x-[F'(x)+λI]-1F(x)

在球S=S(x*,δ)⊂S0中适定且在x*处可导，并有

G'(x*)=I-[F'(x*)+λI]-1F'(x*),

即

这表明λ∈(-β,η)时，ρ(G'(x*))<1，从而只要(1)中λk∈(-β,η),则它产生的序列{xk}是收敛于x*的[1].

此定理的不足之处就是它需要事先知道方程组的解x*才能求出λ的范围,但在实际问题中x*大都是未知的.

1.2 依据KaHTOPOBHY定理

定理2(KaHTOPOBHY) 假设F:D⊂Rn→Rn在凸集D0⊂D上F可导，并且对任何x,y∈D0，有常数r>0使

‖F'(x)-F'(y)‖≤r‖x-y‖,

2 Newton法与简化Newton法相结合的方法

Newton迭代程序

xk+1=xk-[F'(xk)]-1F(xk)，k=0,1,2,…

(2)

若出现F'(xk)为奇异矩阵，则用F'(xk-1)代替F'(xk)，即在此时用一步简化Newton法：xk+1=xk-[F'(xk-1)]-1F(xk).替代(2)式,还若F'(xk+1)再奇异继续用F'(xk-1)代替F'(xk+1)，即用了两步简化Newton法，往后以此类推，若F'(xk+1)为非奇异矩阵，则再回到Newton迭代程序(2).

2.1 算法

步1 给出初始近似值x0及计算精度ε1和ε2；

步2 假定已经进行了k次迭代，已求出xk及F(xk)，计算

F'(xk)=Ak,并记bk=F(xk);

步3 若|F'(xk)|=0，则F'(xk)→F'(xk-1)，转4，否则直接转4；

步4 解线性方程组AkΔxk=-bk得Δxk;

步5 求xk+1=xk+Δxk及F(xk+1);

步6 若‖ΔxK‖≤ε1‖xk‖或‖F(xk+1)‖≤ε2，则置x*=xk+1打印x*，‖F(x*)‖及‖Δxk‖，转步7，否则k+1→k,xk+1→xk,F(xk+1)→F(xk)转步2;

步7 结束.

2.2 收敛性证明

整个程序过程是：Newton法+若干步简化Newton法(奇异处) +Newton法+若干步简化Newton法(奇异处)+…，而Newton法是超线性收敛到方程组的解x*的[3],简化Newton法是线性收敛到x*的[2]，所以整体得到的迭代序列是收敛到x*的.

其中Newton法的收敛域问题可详见[4].这种避免奇异性的方法也可以用在其他Newton型迭代中，简便易行.

参考文献：

[1]李庆杨,莫孜聪,补力群.非线性方程组的数值解法[M].北京:科学出版社,1992.

[2]冯果忱.非线性方程组迭代解法[M].上海:上海科学技术出版社,1989.

[3]黄象鼎,曾钟钢,马亚南.非线性数值分析的理论与方法[M].武汉:武汉大学出版社,2004.

[4]王兴华，关于牛顿法的收敛域[J].科学通报数理化专辑，1980(1).