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 算法小結

這里我們介紹了求解非線性方程組的多種方法。首先，介紹了求解線性方程組的方法——高斯消元法和共軛梯度法，這是非線性方程組的求解基礎。之后，介紹了最常用的非線性方程組求解方法——牛頓法及其變體。為了解決牛頓法由于初值選擇不好而導致無法收斂的問題，介紹了具有較大范圍收斂性的同倫延拓法；為了改進同倫延拓法強烈依賴于初值的不足，介紹了多啟動的同倫延拓法，使用多條延拓曲線來繞過無解區或者奇異點區域。隨后，介紹了兩種搜索算法——長方體法和郭濤算法，前者通過不斷搜索構造的高維長方體來求解方程組的解，后者通過類似遺傳算法的優化算法來尋找方程組的解。最后，為了求解帶參數的非線性方程組的解，介紹了擴展的同倫延拓法，用于求解參數連續變化時方程組的解。

以下通過兩個例子來說明化工過程中的多穩態現象。