2.2 非线性方程组求解方法

对于非线性方程组F（x）=0, x=（x1,x2,…, xn）T∈Rn，其中F=（f1,f2,…, fn）T是从Ω∈Rn到Rn的光滑映射。求解非线性方程组的解具有重要意义。我们知道，牛顿迭代法和它的变体是解非线性方程组解的经典方法，至今仍是基本而重要的方法。但是通常情况下，它的收敛半径很小，因此需要寻找到较好的初值x0。而寻求一个好的初值，本身就是一个困难的问题。在20世纪70年代，发展了一种延拓法，也称同伦延拓法，这种算法具有大范围的收敛性，初值的选取范围被显著扩大，成为求解非线性方程组最有效的方法。蔡大用、白峰杉[8]指出：“克服牛顿法的局部收敛性，寻找大范围收敛的算法可以说是几代人的梦想。”而同伦延拓法已经部分实现了这种梦想，但它用于计算非线性微分方程组的多解时，有时也可能发散。于是学者又提出了多启动的同伦延拓法[9]，相比同伦延拓法，使用多条延拓曲线，具有更好的收敛性，在某种程度上克服了这种缺点，在多解的计算研究中起到一定作用。另一种思路，不使用牛顿迭代算法，转而使用单纯形法、长方体法来搜索非线性方程组的解，也可以使用优化方法例如遗传算法、郭涛算法来求解一些非线性方程组的解。最后，我们本书后面介绍了扩展的同伦延拓法，用于求带有参数的非线性方程组在参数连续变化时，不断变化的方程组的解。

2.2.1 线性方程组高斯消元法和共轭梯度法

非线性问题通常会转化为一系列线性问题来处理，因此，求解线性方程组的方法是求解非线性方程组的基础。考虑一般的n阶线性方程组

Ax=b

其中A是n阶方阵A={aij}n×n，而x=（x1,x2,…, xn）T,b=（b1,b2,…, bn）T是n维列矢量，假定系数矩阵A非奇异，即其行列式不为零，|A|=det（A）≠0。

常用的求解方法是高斯消元法，对于对称正定的矩阵A，还可以使用共轭梯度法来求解。

2.2.1.1 高斯消元法

高斯消元法是一种经典的直接解法，实践表明，目前仍然是一种重要的有效算法，特别是在系数矩阵A非对称、也不正定的一般情形。

进行高斯消元计算时，一般应采用主元消去法，以保证计算过程的舍入误差不会被放大与扩散。许多由微分方程导出的方程组，矩阵A是对角占优的，对角线上的元素就是主元素，经过消元后仍是主元。

高斯消元法求解n阶线性方程组，总共需要做乘除法运算约n3/3次。对于规模较大的问题，需要使用计算机集群、超级计算机，甚至云计算平台。实践表明，更有效的方法是针对问题的特征，构造高效的算法来解决问题，例如，当矩阵A是对称正定矩阵时，共轭梯度法就是一种高效的算法。

2.2.1.2 共轭梯度法

M.Hestenes和E.Stiefel[10]提出用正交投影思想求解线性方程组的方法，开拓了另一种类型的快速高效算法。

当矩阵A对称正定时，对x≠0，内积（Ax,x）＞0，可定义范数‖x‖A=。若（Ax,y）=0，称x, y为A正交。仿照变分法，求解线性方程组Ax=b可以转化为极小化相应的二次函数

在几何上，F（x）=F（x0）是n-1维二次椭球曲面。梯度是F（x）上升最快的方向，二次函数F（x）在x的梯度为

gradF（x）=Ax-b=-r,r={r1,r2,…, rn}T

其中r为残量（residual）, F（x）的极小点x满足gradF（x）=0，即Ax-b=0。

先简单介绍最速下降法（steepest descent method）。

最速下降法 为了求F（x）的极小点，基本思想是：先从已有近似值xk开始，x0 是初始值，沿着下降最快的方向（即反梯度方向pk =rk =-gradF（xk））搜索下一个近似值xk+1，使得在该方向上F（x）达到极小，即取迭代

xk+1 =xk+ckpk,k=0,1,2, …

这里ck应极小化以下二次函数：

同二次函数求极值一样，我们有

于是得到最速下降法的迭代格式如下：

取任意初值x0，再

（1）求残量rk=b-Axk,k=0,1,2, …；

（2）计算修正参数ck=（rk,rk）/（Ark,rk）；

（3）修正xk为xk+1=xk+ckrk。

只要A正定，该方法肯定收敛。但理论分析与计算表明，最速下降法收敛速度较慢，其收敛估计为

当条件数cond（A）=λn/λ 1很大时，F（xk）的等高线是一个很扁的椭圆，其负梯度方向可能偏离真解x*很远，迭代过程的逼近效果不好。虽然如此，但是最速下降法仍然提供了一种很好的正交投影思想。

共轭梯度法（CGM）是最速下降法的改进算法，除了负梯度方向rk之外，还引进共轭关系（即A正交）选取另一个修正方向pk（称为A共轭方向），其计算步骤如下：

第1步，与最速下降法一样，对任意初值x0，求反梯度方向和修正方向

第k步，求取的修正方向pk，不是简单地取为负梯度方向rk，而是用rk-1和pk-1的线性组合来表示pk：

pk=rk-1+dkpk-1,k=2,3, …

并使pk为pk-1的共轭方向，即满足共轭关系（也称A正交）：

（pk,Apk-1）=（Apk,pk-1）=0

由此可以确定参数

然后使用线性组合xk=xk-1+ckpk，使得二次函数F（x）极小化，即可确定常数

由上述算法得到了两个向量序列{rk,pk}，具有下述正交性质：

（rk,pk）=0, k=1,2, …;（rk,rl）=0,（pk,Apl）=0, k≠l

因此{rl}是正交系，{pl}是A正交系。利用这些正交性质，两个参数ck和dk有更简单的表示方法，文献[11]中对两种正交序列的性质做了较好的论证。

综上所述，可以得到共轭梯度法更便于计算的格式：

（1）定任意初值x0，计算残量r0=b-Ax0，取p1=r0（此步与最速下降法相同）。

（2）对k=1,2,3, …采用迭代公式计算：

使用该方法，迭代m次后得到近似解xm，持续迭代计算，直到计算结果满足精度要求为止。

共轭梯度法计算比较简单，在整个计算过程中，A的形状保持不变，只要计算Apk及一些点乘积（rk,rk）及（pk,Apk）。当A是高度稀疏矩阵时，计算Apk的工作量会显著减少。

如果没有舍入误差，对n阶方程组Ax=b，理论上能保证最多n步求得精确解，因此共轭梯度法实质上是一种直接方法。但是由于计算过程中的舍入误差与机器误差，正交性在迭代计算中逐渐受到破坏，经n步迭代并未终止，因此在一段时间里共轭梯度法并没有受到特别的重视。直到后来发现上述迭代有很好的收敛性质，k次迭代有误差估计为

它的收敛速度比最速下降法好很多。许多计算表明，共轭梯度法迭代次数不是很多时，精度已经很好（而迭代次数很大时，正交性变坏，收敛速度变慢），因此，人们更愿意把共轭梯度法当作一种迭代算法使用，它是目前最常用的计算方法之一。

2.2.2 牛顿法及其变体

求解非线性方程组，通常使用迭代法，同时需要找到一个好的初值。在求解非线性方程组的历史中，牛顿迭代法是最经典的算法，至今仍是非常重要的计算方法，文献[12]中对这一方法做了详细介绍。

设方程组F（x）=0的根是x*，雅可比矩阵DF（x）在x* 的某领域N（x*）⊂Ω中非奇异。方向矢量为

V（x）=-（DF（x））-1F（x）

由于在邻域N（x*）中每点x都对应着一个方向V（x），因此它们在区域N（x*）中构成了一个向量场，这些向量曲线都向点x*汇集。设已知根x*的一个近似值x0，用多元函数的Taylor展开

F（x*）-F（x0）=DF（x0）（x*-x0）+（x*-x0）T（DF（ξ））2（x*-x0）

可以解出

x*-x0=V（x0）-（DF（x0））-1（x*-x0）T（DF（ξ））2（x*-x0）

由于x*-x0 很小，可弃去等号右边第2项，得到x*的近似值，即牛顿公式：

x1 =x0+V（x0）

与前式相减，误差表示为

x*-x1=-（DF（x0））-1（x*-x0）（DF（ξ））2（x*-x0）T

因此总误差估计为

‖x*-x1‖≤β0γ‖x*-x0‖2 xn+1 =xn+V（xn）, n=0,1,2, …

其中β0=‖（DF（x0））-1‖, γ=‖（DF（y））2‖。一般地，可将已求得xn作为初值，使用牛顿迭代公式

计算出下一步的新值xn+1。其误差估计为

‖x*-xn‖≤βn-1γ‖x*-xn-1‖2 ≤β0β1…βn-1γn‖x*-x0‖2n

其中βj=‖（DF（xj））-1‖。由此可以看到，在根x*附近，牛顿迭代法具有二次收敛性，在局部的范围下，它是沿着向量场

V（x）=-（DF（x））-1F（x）

的方向前进的，这是最佳的选择方向，因此牛顿法至今都是求解非线性问题最基本最重要的方法。

然而，牛顿法有一个严格的限制，初值x0必须在根x*附近，才有上述收敛效果。也就是说，有一个以x*为中心，以r为半径的球K（x*,r），取其中的任意一点为初值，使用牛顿迭代法都会快速收敛到x*。r为它的收敛半径，或称此收敛球为解x的吸引域。20世纪50年代，苏联的Meisovskix和Kantorovicz研究了这个问题[13]，估计出它的收敛半径，一般都很小，具体如下。

Kantorovicz定理 设F是从G⊂Rd 到Rd 的连续函数，在凸集G0⊂G上可导，对任意x,y∈G0，存在常数γ＞0，使得

‖DF（x）-DF（y）‖≤γ‖x-y‖

假设存在x0∈G0，使得‖（DF（x0））-1‖≤β, ‖（DF（x0））-1F（x0）‖≤η，且满足h=2βγη＜1，那么使用牛顿迭代法产生的序列{xn}均在球K（x0,r1）内，并收敛到唯一解x*∈K（x0,r1）∩G0，其中，半径r1=2η/b,b=1+。误差估计为

由此看到，要求h=2βγη≤1是一个非常严格的限制，因为对于多维问题，β、γ可能会很大，此时要求η很小，即要求‖F（x0）‖很小。虽然后续对牛顿迭代的收敛域还有许多研究，但都不能对此缺点有本质的改变。

要选择初值x0 很接近零点x*一般是很困难的，因为我们事先并不知道x*在哪里。一旦x0 初值选得不好，第1次迭代结果x1 将会落到收敛球之外，以后的迭代将不会收敛，即第1次迭代不好，可能导致整个迭代的发散。实际计算时，也经常遇到这种迭代不收敛的现象。

为了克服牛顿迭代可能落到收敛域之外的不足，提出了一种单调牛顿迭代格式：

xn+1 =xn+hV（xn）, 0＜h＜1，

只要选择适当的0＜h＜1，可保证单调下降性‖F（xn+1）‖＜‖F（xn）‖。这种思想很重要，后文介绍的（多启动）延拓法也用到了这个思想。

牛顿法的另一个缺点是计算工作量很大。因为每次迭代，都要计算d×d雅可比矩阵DF（xn）并求逆。简化的方法之一是采用修正技巧：每计算一次矩阵Bn=DF（xn）-1，保持不变，作m次（如3～5次）内迭代，而这些内迭代只改变F（x）的值。即先定义xn,0=xn，逐个计算

xn,j =xn,j-1-BnF（xn,j-1）, j=1,2, …, m

并定义内迭代终值为xn+1,0=xn,m。这种方法称为拟牛顿法，或改进的牛顿法。

这里虽然内迭代m次的工作量增加了，但一个循环达到了m+1次收敛性：

‖xn,m-x*‖≤Cm‖xn-x*‖m+1

适当选择内迭代次数m，可将整个计算工作量显著减小。

2.2.3 同伦延拓法

为了解决由于初值选择不好而导致的牛顿迭代发散的问题，同伦延拓法（Homotopy continuation）[3～5],[8],[14]提供了一种大范围收敛的迭代算法。同伦延拓法起初是研究算子方程的一种理论工具，后来应用于求解非线性方程与方程组。经过众多学者的不断研究，同伦延拓法目前已经成为大范围求解非线性问题的最有效方法。对于任意给定的初值，使用同伦延拓法通常可以求得问题的一个解。

考虑F:D⊂Rn→Rn，求F（x）=0的解x*∈D。同伦延拓算法的思想是引入一个参数t，同时引入一个函数G（x），满足G（x0）=0。构造一个延拓函数H（t,x），例如

H（t,x）=tF（x）+（1-t）G（x）=0, 0≤t≤1

可以看出，对于延拓函数H（t,x）, x0=x（0）是H（0, x）=0的根。当t从0连续变化，逐渐增加到1时，存在一条相应的延拓曲线x（t），其终点x（1）=x*是H（1, x）=F（x）=0的解。延拓过程如图2-1所示。

从数值上求解上式，需要对参数t进行划分：t0=0＜t1＜t2＜…＜tm=1，其中，要求步长hj=tj-tj-1较小，但是不一定相等。在实际计算中，步长将根据迭代计算的收敛情况来逐渐调整，也就是说，如果连续多步的迭代计算都很快收敛，那么可以适当增大步长，反之，如果迭代不容易收敛，那么就相应减小步长。

通过划分参数t，我们需要逐个求解以下子问题（j=1,2, …, m）：

Pj:H（tj,x）=tjF（x）+（1-tj）G（x）=0，初值x（0）=x（tj-1）

我们选取前一个子问题Pj-1的解x（tj-1），作为该子问题Pj 的初值，特别地，当j=1时，x（0）=x0。因为从tj-1到tj 的变化很小，通常情况下，初值x（tj-1）是解x（t）在（tj-1,tj）中一个较好的近似值，可以使用牛顿法求解这个子问题的解。因为步长很小，而且只需要近似地求解中间过渡的每个子问题，所以，对于每个子问题，只需迭代3至5次即可得到满意的近似解。而当t=1时，需要计算准确的解，此时增加迭代次数，直到计算结果满足精度要求为止。

延拓函数的构造方法有很多，常见的延拓函数有D型、F型[9]，构造方法如下：

D型：H（t,x）=tF（x）+（1-t）DF（x0）（x-x0）

F型：H（t,x）=F（x）-（1-t）F（x0）

虽然延拓法在大多数情况下对解决一些困难的非线性方程组很有效，但是有些时候仍然存在一些问题：

（1）如果选取的初值x0 远离方程组的解x*，那么DF（x0）与DF（x*）差异较大，此时同伦曲线x（t）存在很大的弯曲，不能保证x（t）会逼近目标解x*。因此初值的选择依然有一定的限制，为了确保延拓算法收敛到x*，需要确保初值x0 属于x*的某个吸引域D（x*）。吸引域的大小与延拓函数的构造有关。

（2）当t接近1时，x（t）将逼近目标解x*。虽然F（x）与1-t的数值都变得很小，但G（x*）的数值可能会很大，这使得求解H（t,x）=0变得困难，即牛顿迭代可能很难收敛。为了克服这一困难，我们需要构造G（x*）数值尽可能小的函数。

综上所述，虽然同伦延拓法是求解非线性方程组最有效的方法之一，极大地扩展了初值的收敛范围，但是对于一些特殊情形也会失效。

2.2.4 多启动延拓法

如前所述，牛顿迭代xj+1=xj+V（xj）（j=1,2, …）在接近零点时具有2次收敛性，其中V（x）=-（DF（x））-1F（x）是理想的方向场，但是当选取的初值不好时，它的第1次迭代就可能导致发散。另一方面，同伦延拓法具有较好的大范围收敛性，但它强烈地依赖于选择的初值，有些情况下，会使延拓曲线x（t）偏离理想的向量场V（x），严重地弯曲，以致绕到了无解区，或遇到奇点，或收敛到别的解而使迭代计算失败。为克服这些缺点，同时充分利用这两种方法的优点，有学者[9]提出了多启动延拓法（MSEM），它总是沿着矢量场V（x）方向迭代，并使整个连通的非奇异域Dj 都成为吸引域。多启动延拓法具体如下：

构造多启动延拓函数

H（t,x）=tF（x）+（1-t）G（x,x′）=0, 0≤t≤1

其中，G（x,x′）=DF（x′）（x-x′）或G（x,x′）=F（x）-F（x′）。

第1步 取0到1的k-1个划分，记为tj0=0＜tj1＜…＜1（j=1,2, …, k-1），但对每个第j次延拓只解决第1个子问题，即

Pj1:H（t,x）=tF（x）+（1-t）G（x,x′）=0，

t=tj1,x=x（tj1）, 1≤j≤k-1

当j=1时，取初值x′=x0，当1＜j＜k时，取上次的结果为初值x′=x（tj-1,1）。

第2步 最后对第k次延拓，取一个完全的0到1的划分Zk:tk0=0＜tk1＜tk2＜…＜tkm=1，并取k-1次延拓的结果为初值x′=x（tk-1,1），逐一求解此延拓的m个子问题：

Pki:H（t,x）=tF（x）+（1-t）G（x,x′）=0，

t=tki,x=x（tki）, i=1,2, …, m

于是x*=x（tkm）就是我们所希望的解。

可以看到，新算法的第1步能够提供一条从起点x0通往目标解x*的较好路径，使得到的解曲线x（t）基本上沿着从x0发出的理想矢量场V（x）=-（DF（x））-1F（x）连续地前进，且初值的影响变得越来越小。当接近解x*时，在第2步采用经典的延拓法或牛顿法来求解。对于计算量而言，第1步花了许多计算，但由于每次增加的步长很小，用牛顿迭代计算3～5次就可以得到满意的解，总工作量不一定增加。

为解决第1步，有两种简单的方式来选择它们的第一步长，具体如下。

（1）增加型。对1≤j≤k-1及适当大的L，取第1个划分点为

tj1=1/（k+L-j）, 1≤j≤k-1

若认为E=（0,1）是整个考虑的区间，用步长h1=t11，求解第1个子问题P11后，剩下的长度为r1=1-t11=（k+L-1）。对j=2，它的实际步长仅仅是h2=r1× t21=1/（k+L-1），而剩下r2=r1-h2=（k+L-3）/（k+L-1）。于是最后一个子问题Pk-1,1有相同的实际步长hk-1=1/（k+L-1），余下区间的长度是rk-1=（L-1）/（k+L-1）。例如取k=41, L=10，我们将有相同的实际步长hj=1/50，最后剩下的区间长为rk-1=10/50=1/5。只在这个小区间上最后完成第2步的延拓计算。

（2）减小型。对1≤j≤k-1取第1个划分点为

tj1=c/k=h0,c＞0

解第1个P11后，剩下r1=1-c/k。一般地，对第j个子问题Pj1，我们将有实际步长hj=rj-1×h0，并剩下rj=（1-h0）j。最后对j=k-1，剩下长度为rk-1=（1-c/k）k-1≈e-c=γ（c）。例如，对c=3、4、5、6，剩下长度为γ（c）≈0.05、0.018、0.0067、0.0025，已很小很小，最后在第2步容易用延拓法或牛顿法解决。

上述两种选择，在解各种问题时有不同的效果。若初值接近某个奇点，前若干步应是小步长，不得不采用（增加型）划分1。在两种情形，当改变参数L、k（或c）、局部步长，及在第j个问题Pj中增加延拓步数，都能适当改变延拓曲线的路径。

多启动延拓法的第1步，实质上是求解一个特殊的Davidenko方程：

应当指出，多启动延拓法用于计算非线性问题的多解是有效的。在有些情形，若经典的D型与F型延拓法失败了，用多启动延拓法有时能收敛到所希望的解。原因是多启动延拓法可能改变了解曲线x（t）的路径，致使它避开了某些奇点，或者绕开了同伦延拓法无解的地方。

2.2.5 单纯形算法和长方体算法

前面介绍的所有方法都假设已知一个较好的初值，然后迭代计算非线性方程组的解，并没有涉及如何来选取较好的初值。在实际应用中，如何搜索与计算非线性方程组的实根（特别是求所有实根），是一个很困难的问题。本节首先介绍如何求解方程组一个根的单纯形法，然后介绍求所有根的长方体法。

2.2.5.1 单纯形算法

为了搜索非线性方程组的一个较好的初值，Scarf[15]在1967年提出了单纯形算法，用此方法一般可以找到一个解，文献[12],[16]对此做了较详细的介绍与评述。

由n+1个节点构成的n维多面体称为n维单纯形，如一维的线段，二维的三角形，三维的四面体等，它是n维空间中最简单的几何图形。单纯形法的基本思想是：在一个n维有界域Ω 内，如多维长方体区域[a,b]n中，寻求方程组F（x）=0的根x*，先将Ω划分为有限个小的单纯形e∈Jh，如果找到了一个e，在它的n+1个顶点上，每个分量值Fj（z）（j=1,2, …, n）都至少改变了1次符号，则在此e上可能存在一个点x*，使得所有的Fj（x*）=0。于是在e中此根可用某种迭代法求得，因此关键的问题是如何找到这种单纯形。为此，在每个节点z上，定义唯一的标号l（z）如下：

当i=1,2, …, n时，如果对i有Fi（z）＞0，且对所有j＜i还有Fj（z）≤0（对i=1无此要求），那么规定l（z）=i；如果对所有j=1,2, …, n都有Fj（z）≤0，那么规定l（z）=0。

若e的n+1个顶点上有完全的标号{0,1,2, …, n}，则称为全标号单纯形。如果具有标号为{0,1,2, …, n-1}，则称为几乎全标号单纯形。全标号单纯形是我们要找的单纯形。因为在全标号单纯形上，每个分量Fi（x）都变号1次，其中可能有F（x）的零点。

单纯形法中包含着一个搜索算法，具体如下：搜索从区域的边界（它是n-1维域）开始，首先在边界上找到一个n-1维的全标号单纯形，它的侧面连接着许多其他的单纯形（在Ω内），从中寻找n维全标号单纯形或几乎全标号单纯形，依此类推，直到找到一个n维全标号单纯形为止。

单纯形法很优美，是拓扑学与计算数学的结合，并且首次对不动点定理给出了一个构造性的证明。但是在具体的实施中还有一定困难：

（1）一般只给出了一个根，无法求解所有的根。

（2）计算工作量很大，当n较大时，即使在n-1维边界上搜索n-1维全标号单纯形的工作量也很大。

（3）划分一个n维立方体为许多单纯形不容易，也不直观。单纯形看似简单，但将一个区域划分为许多单纯形却很复杂（例如将1个三维立方体划分为6个四面体就比较困难），其单元与节点的编号就非常复杂，而要找到这种全标号单纯形的程序就更加复杂。由于这些原因，单纯形法并未普遍应用。文献[16]第120页写道：“由于该算法程序实现比较困难，目前还处于实际应用人员的视野之外。”

2.2.5.2 长方体算法

单纯形法虽然很难应用，但是它却为如何判断在一个小单元中是否存在一个根提供了启发。有学者[9]汲取其中有益的思想，提出了长方体算法（或立方体，多方体算法）。长方体算法简单实用，当未知数不太多时，能搜索所有的解。

设Ω是一个n维长方体，每个方向用m+1个节点，将它划分为mn个n维长方体单元e，共有（m+1）n个节点。整个计算通过逐层单元循环计算来完成。

长方体算法由以下3部分组成：

（1）判断单元中是否有根。对每个节点z及每个分量Fi（z），规定一个指标d（z,i）如下：

若Fi（z）＞0，定义d（z,i）=1；

若Fi（z）＜0，定义d（z,i）=0；

若Fi（z）=0，定义d（z,i）=0.5。

长方体单元e的2n个顶点集合记为T（e）。对每个分量Fi（z），记所有顶点z∈T（e）上的指标d（z,i）之和为g（i）=d（z,i）。

若满足0＜g（i）＜2n，则记dd（i）=1，即Fi（x）在这些顶点上至少变号1次，否则记dd（i）=0，分量Fi（z）没有变号。若在单元e上，所有节点的和dd（i）=2n，即每个分量都至少变号1次，则在此单元e中可能有1个根x*。通过这个过程，找到可能有根的单元。

（2）在单元e中求根。一旦找到了含根的单元e，即转入在此单元内求根的过程。我们取此单元的中心y作为初始点（当然，一个单元中有时可能有很多根，为了防止漏掉根，也可取更多的点作初始点，由于这些单元的数量相对较少，计算工作量并不大），用牛顿法迭代（最好用延拓法或多启动延拓法迭代，因为它们的数值逼近比较稳定），得到近似根yj（j=0,1,2, …），直到满足精度要求。

（3）逻辑判断。为了计算和识别这个根，还要完成某些逻辑判断。例如，由于在e中计算的根很可能落到其他单元，一个根也可能被计算多次，因此要比较y*与前面已计算的根，如果y*与已有的所有根不同，那么记录y*是一个新的根。重复以上过程，直到最后得到所有不同的实根为止。

长方体搜索法在理论上没有单纯形法完美，在单纯形法中n个方程用n+1个点上的符号就可判断是否有根，而长方体法需要计算所有2n个顶点才能确定。长方体法的优点是：划分采用n维长方体，结构简单直观，编程也简单，而且一次计算能得到所有的不同实根。缺点是计算量很大，需要对所有节点进行全面搜索。

2.2.6 郭涛算法

不同于使用长方体法的全面搜索，我们也可以使用优化算法来搜索寻找方程组的根。本节介绍郭涛算法，文献报道，该方法可用于寻找到方程组所有的根。

郭涛算法与遗传算法相关，首先对遗传算法进行简介。遗传算法（genetic algorithms, GA）是由美国Michigan大学的John Holland教授提出的，是建立在自然选择和群体遗传学机理基础上的随机迭代和进化、具有广泛适用性的搜索方法，具有很强的全局优化搜索能力。它模拟了自然选择和自然遗传过程中发生的繁殖、交配和变异现象，根据适者生存、优胜劣汰的自然法则，利用遗传算子（选择、交叉和变异）逐代产生优选个体（即候选解），最终搜索到较优的个体。

郭涛算法是基于遗传算法思想提出的一类基于子空间搜索和群体爬山法相结合的群体随机搜索算法，也称为多父体杂交算法。对于优化问题，从理论上说，可以一次求得多个最优解。利用该特性，构造与非线性方程组等价的优化问题，也可以一次求解多个方程组的解。郭涛算法简介如下[17]。

考虑求解如下带约束的函数优化问题：

使其满足不等式gi（X）≤0, i=1,2, …，其中X=（x1,x2,…, xn）T∈Rn,D={X|Ii≤xi≤Ui,i=1,2, …, n}, f（X）为目标函数。记D中的m个点为X′j=（xj1,xj2,…, xjn）T,j=1,2, …, m。记它们所张成的子空间为

其中，ai满足条件ai=1, -0.5≤ai≤1.5。

记hi（X）=和H（X）=hi（X）。

定义逻辑函数

来表示X1优于X2。

郭涛算法如下：

    初始化P={X1,X2,…, XN}; Xi∈D
    t=0;
    Xbest=arg better（Xi,X）, ∀X∈P;
    Xworst=arg better（X,Xi）, ∀X∈P;
    while（better（Xworst,Xbest）==FALSE）{
      从P中随机选取m个点X′1,X′2,…, X′m;
      从V中随机选取点X′;
      if（better（X′, Xworst）==TRUE）{Xworst=X′};
      t=t+1;
      Xbest=arg better（Xi,X）, ∀X∈P;
      Xworst=arg better（X,Xi）, ∀X∈P;
    }
    输出t,P;
    end

其中N为群体P的大小，m为子空间的维数（如果m个向量线性无关）, t为迭代次数，Xbest=arg better（Xi,X）, ∀X∈P，表示把Xi（i=1,2, …, N）中最好的那些变量（个体）之一记作Xbest。

对于郭涛算法，以下几点需要进一步说明：

（1）初始化是随机从解空间D中选取的N个点（个体）形成的初始群体P。

（2）N的选取可根据问题的维数n与f（X）场景的复杂性而定，当n较大且场景复杂时，N的取值较大，反之，则可以适当减小。一般取20≤N≤150。

（3）m的选取，根据经验取m=7、8、9或10比较合适。

对郭涛算法进一步的分析如下：

（1）算法采用了演化计算中的群体搜索策略，保证了搜索空间的全局性。

（2）算法采用随机子空间中的随机搜索（多父体重组）策略，特别是子空间中随机搜索的非凸性，即

使算法搜索的子空间可覆盖多父体的凸组合空间，保证了随机搜索的遍历性，即解空间中不存在算法搜索不到的死角。

（3）算法采用了优胜劣汰策略，每次只把群体中适应性最差（目标函数值最大）的个体淘汰出局，淘汰压力最小，既保证了群体的多样性，也保证了适应性最好（目标函数值最小）的个体可以长期存活。这种群体爬山策略，保证了整个群体最后集体达到最深的谷底。当最优解不唯一时，该算法可能一次同时找到多个最优解。

根据算法的上述3条基本特性，如果把迭代次数t当作群体P 的生存代，则构成一个马尔可夫链。使用有关演化算法的收敛性分析方法[18]进行研究，可以得知郭涛算法能够确保计算过程的收敛性。

综上所述，从理论上来说，对于由非线性方程组转化而来的等价的优化问题，使用郭涛算法，可以同时找到多个方程解。

2.2.7 扩展的同伦延拓法

在实践中，有时需要求解一系列变化参数下的非线性方程组的解。针对这个问题，我们提出了一种基于同伦延拓法的非线性问题求解方法，称为扩展的同伦延拓法。用该方法可以方便地求出参数连续变化过程中方程组的解的变化情况。

问题可以描述为：对于带参数的非线性方程组F（x,λ）=0，当已知λ=λ0时方程组的解x0，如何求取不同λ数值下的解x，其中是λ具有物理意义的参数，而且连续变化，定义域为[λl,λu]。

参照同伦延拓法的思想，我们可以将λ划分为多个小段，设xi是当λ=λi时方程组F（x,λi）=0的解，其中λl≤λi≤λu。那么，问题进一步转化为，已知xi、λi的情况下，如何确定λi+1并且求解F（x,λi+1）=0的解，其中λl≤λi+1≤λu。即已知默认步长Δλ，最优迭代次数Nopt（最优迭代次数是由结果的精度和迭代算法确定的，它是一个统计结果）、参数下限λl、参数上限λu。λi是第i步的参数值，ηi是第i步的步长伸缩因子，xi是当λ=λi时方程组F（x,λi）=0的解。算法中有三个关键步骤：步长选取，预估初值，迭代求解精确值。

由于非线性问题的特殊性，在求解过程中可能存在奇异点，此时方程组F（x,λi）=0的雅可比矩阵Fx奇异，预估及校正算法需要特别处理。为此，在计算过程中，需要不断判断矩阵Fx是否奇异。

计算过程描述如下：

计算F（x,λi）=0的雅可比矩阵Fx在点（xi,λi）的矩阵是否奇异。如果矩阵非奇异，即det（Fx）≠0，那么按照下列过程计算。

（1）使用λi、ηi、Δλ确定λi+1，即λi+1=λi+ηi ·Δλ。

（2）使用xi、xi-1、λi、λi+1预估F（x,λi+1）=0的初值，即=xi+（xi-xi-1）。

（3）使用作为初值迭代求解F（x,λi+1）=0的解xi+1。记录迭代收敛的次数Np，计算步长伸缩因子，即。

重复以上步骤直到计算出λ=λu时F（x,λ）=0的解xu。计算过程的示意图如图2-2所示。

如果矩阵奇异，即det（Fx）=0（数值计算中，当det（Fx）≤ε时，认为矩阵奇异），那么按照下列过程计算。具体过程描述如下：

（1）确定λi+1，这里取λi+1=λi-1。

（2）使用xi、xi-1预估F（x,λi+1）=0的初值，即=2xi-xi-1。

（3）使用作为初值迭代求解F（x,λi+1）=0的解xi+1。记录迭代收敛的次数Np，计算步长伸缩因子，即。

（4）将默认步长改为Δλ的相反数。转入正常点的计算过程。

奇异点附近的处理过程如图2-3所示。

由图2-3可知，预估过程是求出（xi-1,λi-1）关于奇异点（xi,λi）在该点切线的法线的对称点（, λi-1）的近似值作为初值，在校正过程中迭代收敛到方程的解（xi+1,λi-1）。

在扩展的同伦延拓算法中，每步的计算用到了前一步的计算结果以及迭代次数，算法提高了计算的有效性。该方法与同伦延拓法相比，计算中的每一步的结果都是受关注的，而同伦延拓法只有t=1时的数值是受关注的。此外，扩展的同伦延拓法也可用于求非线性方程组的解，即将待求解的非线性方程组F（x）=0，看作带参数非线性方程组F（x,λ）=0中当λ=λ0时的特殊情况。

2.2.8 算法小结

这里我们介绍了求解非线性方程组的多种方法。首先，介绍了求解线性方程组的方法——高斯消元法和共轭梯度法，这是非线性方程组的求解基础。之后，介绍了最常用的非线性方程组求解方法——牛顿法及其变体。为了解决牛顿法由于初值选择不好而导致无法收敛的问题，介绍了具有较大范围收敛性的同伦延拓法；为了改进同伦延拓法强烈依赖于初值的不足，介绍了多启动的同伦延拓法，使用多条延拓曲线来绕过无解区或者奇异点区域。随后，介绍了两种搜索算法——长方体法和郭涛算法，前者通过不断搜索构造的高维长方体来求解方程组的解，后者通过类似遗传算法的优化算法来寻找方程组的解。最后，为了求解带参数的非线性方程组的解，介绍了扩展的同伦延拓法，用于求解参数连续变化时方程组的解。

以下通过两个例子来说明化工过程中的多稳态现象。