第2章非線性方程求解

第2章非線性方程求解

2.1化工實(shí)際問題的提出2.2實(shí)根的對(duì)分法2.3直接迭代法2.4松弛迭代法2.5韋格斯坦法2.6牛頓迭代法2.7割線法2.8非線性方程組的牛頓方法2.9化工生產(chǎn)中非線性方程組求解應(yīng)用實(shí)例

目錄2.1化工實(shí)際問題的提出

求解非線性方程是化工設(shè)計(jì)及模擬計(jì)算中必須解決的一個(gè)問題。與線性方程相比，非線性方程問題無論是從理論上還是從計(jì)算公式上，都要比線性方程復(fù)雜的多。對(duì)于一般的非線性f(x)=0，計(jì)算方程的根既無一定章程可循也無直接方法可言。例如，求解高次方程組7x6-x3+x-1.5=0的根，求解含有指數(shù)和正弦函數(shù)的超越方程ex-sin(x)=0的零點(diǎn)。解非線性方程或非線性方程組也是計(jì)算方法中的一個(gè)主題。一般地，我們用符號(hào)f(x)來表示方程左端的函數(shù)，方程的一般形式表示為f(x)=0，方程的解稱為方程的根或函數(shù)的零點(diǎn)。通常，非線性方程的根不止一個(gè)，而任何一種方法只能算出一個(gè)根。因此，在求解非線性方程時(shí)，要給定初始值或求解范圍。而對(duì)于具體的化工問題，初值和求解范圍常?？筛鶕?jù)具體的化工知識(shí)來決定。常見的雷諾數(shù)和摩擦系數(shù)關(guān)系方程在雷諾數(shù)低于4000時(shí)有以下關(guān)系式：

2.12.82.72.52.32.2總目錄2.92.62.42.1化工實(shí)際問題的提出這是一個(gè)典型的非線性方程。我們?cè)诠苈吩O(shè)計(jì)中經(jīng)常碰到。當(dāng)我們已知雷諾數(shù)Re，如何根據(jù)公式（2-1）求出摩擦系數(shù)λ，這是我們?cè)诠苈吩O(shè)計(jì)中必須首先解決的問題。對(duì)于方程（2-1）而言，無法用解析的方法求出摩擦系數(shù)，只能用數(shù)值求解的方法。如用在下面即將介紹的松弛迭代法，假設(shè)：

則利用松弛迭代公式可得：經(jīng)11次迭代可得摩擦系數(shù)為0.07593。同樣，在n個(gè)組分的等溫閃蒸計(jì)算中，通過物料和相平衡計(jì)算，我們可得到如下非線性方程：

2.12.82.72.52.32.2總目錄2.92.62.42.1化工實(shí)際問題的提出在方程（2-3）中只有α是未知數(shù)，ki為相平衡常數(shù)，zi為進(jìn)料組分的摩爾濃度，均為已知數(shù)。和上面的情況一樣，方程（2-3）也無法直接解析求解，必須利用數(shù)值的方法，借助于計(jì)算機(jī)方可精確的計(jì)算。對(duì)于這個(gè)問題的求解，可利用我們下面介紹的牛頓迭代法進(jìn)行計(jì)算，也可利用其他迭代公式進(jìn)行計(jì)算，如采用牛頓迭代公式，則可以得到如下的具體迭代公式：

（2-4）

飽和蒸氣壓是我們經(jīng)常要用到的數(shù)據(jù)，雖然我們可以通過實(shí)驗(yàn)測(cè)量來獲取飽和蒸氣壓的數(shù)據(jù)，但我們通常利用前人已經(jīng)測(cè)量得到的數(shù)據(jù)或回歸的公式來獲取，這可以減輕我們大量的基礎(chǔ)實(shí)驗(yàn)工作。公式（2-5）是一種常用的飽和蒸氣壓計(jì)算公式：其中p為飽和蒸氣壓，單位為mmHg，T為溫度，單位為K，A、B、C、D為已知系數(shù)。要想得到某一溫度下的飽和蒸氣壓，直接利用公式（2-5）是無法得到的。因?yàn)楣剑?-5）兩邊都有未知變量，并且無法用解析的方法求解，必須用數(shù)值計(jì)算的方法求解。通過上面的一些例子，我們可以發(fā)現(xiàn)，如果沒有適當(dāng)?shù)氖侄魏娃k法來求解非線性方程，那么化學(xué)化工中的許多研究、設(shè)計(jì)等工作將無法展開，這勢(shì)必影響化學(xué)化工的發(fā)展，下面我們將介紹一些實(shí)用的非線性方程求解方法，并提供計(jì)算機(jī)程序。

2.12.82.72.52.32.2總目錄2.92.62.42.2實(shí)根的對(duì)分法

2.2.1使用對(duì)分法的條件

2.2.2對(duì)分法求根算法

2.2.3對(duì)分法VB程序清單

2.12.82.72.52.32.2總目錄2.92.62.42.2.1使用對(duì)分法的條件對(duì)分法或稱二分法是求方程近似解的一種簡(jiǎn)單直觀的方法。設(shè)函數(shù)f(x)在[a,b]上連續(xù)，且f(a)f(b)<0，則f(x)在[a,b]上至少有一零點(diǎn)，這是微積分中的介值定理，也是使用對(duì)分法的前提條件。計(jì)算中通過對(duì)分區(qū)間，逐步縮小區(qū)間范圍的步驟搜索零點(diǎn)的位置。如果我們所要求解的方程從物理意義上來講確實(shí)存在實(shí)根，但又不滿足f(a)f(b)<0，這時(shí)，我們必須通過改變a和b的值來滿足二分法的應(yīng)用條件。

2.12.82.72.52.32.2總目錄2.92.62.42.2.2對(duì)分法求根算法

計(jì)算f(x)=0的一般計(jì)算步驟如下：

1、輸入求根區(qū)間[a,b]和誤差控制量ε，定義函數(shù)f(x)。

2、判斷:如果f(a)f(b)<0則轉(zhuǎn)下，否則，重新輸入a和b的值。

3、計(jì)算中點(diǎn)x=(a+b)/2以及f(x)的值分情況處理（1）|f(x)|<ε：停止計(jì)算x*=x，轉(zhuǎn)向步驟4（2）f(a)f(x)<0：修正區(qū)間[a,x]→[a,b]，重復(fù)3（3）f(x)f(b)<0：修正區(qū)間[x,b]→[a,b]，重復(fù)34、輸出近似根x*。

右圖給出對(duì)分法的示意圖。

x3=(x0+x2)/2

x2=

(x0+x1)/2

x0

x3

x1

x1

2.12.82.72.52.32.2總目錄2.92.62.42.2.3對(duì)分法VB程序清單

PrivateSubCommand1_Click()Dimx1,x2,x,y1,y2,y,eer80x1=InputBox("x1")x2=InputBox("x2")y1=f(x1)y2=f(x2)Ify1*y2<0Then

GoTo100ElsePrint"pleaserepeatinputx1andx2"

GoTo80EndIf100x=(x1+x2)/2y=f(x)IfAbs(y)<=0.001ThenPrint"thefunctionrootis";xPrint"y=";yElseIfy1*y<0Thenx2=xy2=yGoTo100Elsex1=xy1=yGoTo100EndIfEndIfEndSubPublicFunctionf(x)Dimyy=x^3+x^2-1f=yEndFunction2.12.82.72.52.32.2總目錄2.92.62.42.2.3對(duì)分法求解實(shí)例用對(duì)分法求在區(qū)間[1,2]之間的根。解：(1)f(1)=-2.8，f(2)=0.3，由介值定理可得有根區(qū)間[a,b]=[1,2]。

(2)計(jì)算x2=(1+2)/2=1.5，f(1.5)=-0.45，有根區(qū)間[a,b]=[1.5,2]。

(3)計(jì)算x3=(1.5+2)/2=1.75，f(1.75)=0.078125，有根區(qū)間[a,b]=[1.5,1.75]。

一直做到|f(xn)|<ε(計(jì)算前給定的精度)或|a-b|<ε時(shí)停止。詳細(xì)計(jì)算結(jié)果見表2-1。對(duì)分法的算法簡(jiǎn)單，然而，若f(x)在[a,b]上有幾個(gè)零點(diǎn)時(shí)，如不作特殊處理只能算出其中一個(gè)零點(diǎn)；另一方面，即使f(x)在[a,b]上有零點(diǎn)，也未必有f(a)f(b)<0。這就限制了對(duì)分法的使用范圍。對(duì)分法只能計(jì)算方程f(x)=0的實(shí)根。

對(duì)于多個(gè)零點(diǎn)的方程，我們可以通過將給定的區(qū)間[a,b]進(jìn)行細(xì)分，然后在細(xì)分后的區(qū)間內(nèi)用二分法分別求解，從而得到多個(gè)零點(diǎn)。例如求方程在0-30內(nèi)的所有根。需要對(duì)二分法進(jìn)行以下處理：即先給定一個(gè)a，本例中為0，然后不斷增加，直到找到一個(gè)b，使f(a)f(b)<0，調(diào)用二分法，計(jì)算在[a，b]范圍內(nèi)的根，然后將b作為a，重復(fù)上面的工作，直到計(jì)算范圍超出30為止。

2.12.82.72.52.32.2總目錄2.92.62.4VB調(diào)用2.3直接迭代法

對(duì)給定的方程f(x)=0，將它轉(zhuǎn)換成等價(jià)形式：。給定初值x0，由此來構(gòu)造迭代序列，k=1,2,…，如果迭代收斂，即有,則就是方程f(x)=0的根。在計(jì)算中當(dāng)小于給定的精度控制量時(shí)，取為方程的根。例如，代數(shù)方程x3-2x-10=0的三種等價(jià)形式及其迭代格式如下：對(duì)于方程構(gòu)造的多種迭代格式，怎樣判斷構(gòu)造的迭代格式是否收斂？收斂是否與迭代的初值有關(guān)？根據(jù)數(shù)學(xué)知識(shí)，我們可以直接利用以下收斂條件：

1、

當(dāng)有

2、

在[a,b]上可導(dǎo)，并且存在正數(shù)L<1，使任意的，有則在[a,b]上有唯一的點(diǎn)滿足，稱為的不動(dòng)點(diǎn)。而且迭代格式對(duì)任意初值均收斂于的不動(dòng)點(diǎn)，并有下面誤差估計(jì)式：

(2-6)

要構(gòu)造滿足收斂條件的等價(jià)形式一般比較困難。事實(shí)上，如果為f(x)的零點(diǎn)，若能構(gòu)造等價(jià)形式，而，由的邊疆性，一定存在的鄰域，其上有，這時(shí)若初值迭代也就收斂了。由此構(gòu)造收斂迭代格式有兩個(gè)要素，其一，等價(jià)形式應(yīng)滿足；其二，初值必須取自的充分小鄰域，其大小決定于函數(shù)f(x)，及做出的等價(jià)形式。2.12.82.72.52.32.2總目錄2.92.62.42.3直接迭代法例：求代數(shù)方程x3-2x-5=0，在x0=2附近的零點(diǎn)。

解：1）x3=2x+5構(gòu)造的迭代序列收斂。取x0=2，則

準(zhǔn)確的解是x=2.09455148150。2）將迭代格式寫為

迭代格式不能保證收斂,但并不一定不收斂。VB程序界面：

2.12.82.72.52.32.2總目錄2.92.62.42.4松弛迭代法

有些非線性方程或方程組當(dāng)用上一節(jié)中的直接迭代法求解時(shí),迭代過程是發(fā)散的。這時(shí)可引入松弛因子,利用松弛迭代法。通過選擇合適的松弛因子,就可以使迭代過程收斂。松弛法的迭代公式如下:

（2-7）由上式可知，當(dāng)松弛因子ω等于1時(shí)，松弛迭代變?yōu)橹苯拥?。?dāng)松弛因子ω大于1時(shí)松弛法使迭代步長(zhǎng)加大,可加速迭代，但有可能使原來收斂的迭代變成發(fā)散。當(dāng)0<ω<1時(shí),松弛法使迭代步長(zhǎng)減小,這適合于迭代發(fā)散或振蕩收斂的情況,可使振蕩收斂過程加速。當(dāng)ω<0時(shí),將使迭代反方向進(jìn)行,可使一些迭代發(fā)散過程收斂。松弛法是否有效的關(guān)鍵因子是松弛因子ω的值能否正確選定。如果ω值選用適當(dāng),能使迭代過程加速,或使原來不收斂的過程變成收斂。但如果ω值選用不合適，則效果相反,有時(shí)甚至?xí)乖瓉硎諗康倪^程變得不收斂。松弛因子的數(shù)值往往要根據(jù)經(jīng)驗(yàn)選定,但選用較小的松弛因子,一般可以保證迭代過程的收斂。2.12.82.72.52.32.2總目錄2.92.62.42.4松弛迭代法例:用松弛迭代法求解下面非線性方程組,并分析松弛因子對(duì)迭代次數(shù)及收斂過程的影響。已知迭代初值x和y均為0,收斂精度ε=0.001

解:取以下迭代表達(dá)式:

若取松弛因子為1.1則其迭代過程如左表。若改變松弛因子,迭代過程及迭代所需的次數(shù)亦將發(fā)生變化，詳見右表。由右表數(shù)據(jù)可知,當(dāng)松弛因子小于1時(shí),增大松弛因子,可加速迭代過程,減少迭代次數(shù),但當(dāng)松弛因子大于1時(shí),迭代次數(shù)反而增加,當(dāng)松弛因子達(dá)到1.56

時(shí),迭代過程分散。

2.12.82.72.52.32.2總目錄2.92.62.4此法是一種迭代加速方法。在圖2-3中,從x0開始,曲線y=φ(x)和直線y=x之間的階梯形折線為單變量的直接迭代過程。其迭代步長(zhǎng)為每個(gè)階梯的水平距離。若利用圖2-3中曲線φ(x)上的兩個(gè)點(diǎn)作一直線,通過求該直線和y=x的交點(diǎn)C的橫坐標(biāo)xW作為新的迭代點(diǎn),這樣的迭代計(jì)算過程就是韋格斯坦法。即韋格斯坦法需要在已知兩點(diǎn)的前提下迭代計(jì)算第三點(diǎn)。一般第一點(diǎn)為人為設(shè)定,第二點(diǎn)利用直接迭代計(jì)算而得,則第三點(diǎn)就可以用韋格斯坦法迭代求解。已知A點(diǎn)的坐標(biāo)為(x1,y1),B點(diǎn)的坐標(biāo)為(x2,y2)。則過A、B兩點(diǎn)的直線方程為：求該直線和y=x的交點(diǎn)可得:

韋格斯坦法的一般計(jì)算通式為

2.5韋格斯坦法

A

y1

y0

φ(x)Bx0x1x2x*

xw圖2.-3韋格斯坦法迭代示意圖2.12.82.72.52.32.2總目錄2.92.62.42.5韋格斯坦法由上述公式可知,韋格斯坦法也是一種松弛法,其松弛因子為

一般情況下,當(dāng)1>k>0時(shí),迭代過程為單調(diào)收斂過程。當(dāng)-1<k<0時(shí),迭代過程為振蕩收斂過程,但當(dāng)k=1時(shí),收斂將發(fā)散,故在編程計(jì)算時(shí)應(yīng)注意當(dāng)k=1時(shí)則取k=0進(jìn)行計(jì)算。韋格斯坦法求解方程x3-2x2+5x-4=0根的QB程序見課本2.12.82.72.52.32.2總目錄2.92.62.42.6牛頓迭代法

2.6.1牛頓法的理論推導(dǎo)2.6.2牛頓法的幾何意義2.12.82.72.52.32.2總目錄2.92.62.4對(duì)方程f(x)=0可構(gòu)造多種迭代格式，牛頓迭代法是借助于對(duì)函數(shù)f(x)=0的泰勒展開而得到的一種迭代格式。將f(x)=0在初始值x0做泰勒展開得：取展開式的線性部分作為的近似值，則有：設(shè)則令類似地，再將f(x)=0在x1作泰勒展開并取其線性部分得到：一直做下去得到牛頓法的迭代格式：牛頓迭代格式對(duì)應(yīng)于f(x)=0的等價(jià)方程為：若b是f(x)的單根時(shí)，，則有，只要初值x0充分接近b，牛頓迭代都收斂。牛頓迭代是二階迭代方法?？梢宰C明，b為f(x)的a重根時(shí)，迭代也收斂，但這是一階迭代，收斂因子為，若這時(shí)取下面迭代格式，它仍是二階方法：

2.6.1牛頓法的理論推導(dǎo)

2.12.82.72.52.32.2總目錄2.92.62.4以為斜率作過(x0，f(x0))點(diǎn)的直線，即作f(x)在x0的切線方程：令y=0，則在x1處的切線與x軸的交點(diǎn)x1，即：再作f(x)在x1處的切線，得交點(diǎn)x2，逐步逼近方程的根b。如圖2-4所示。在區(qū)域[x0,x0+h]的局部“以直代曲”是處理非線性問題的常用手法。在泰勒展開中，截取函數(shù)展開的線性部分替代f(x)。

yxx2x1x02.6.2牛頓法的幾何意義

圖2-4牛頓切線法示意圖

2.12.82.72.52.32.2總目錄2.92.62.42.6.2牛頓法的幾何意義例：用牛頓迭代法求方程f(x)=x3-7.7x2+19.2x-15.3，在x0=1附近的零點(diǎn)。解：

計(jì)算結(jié)果列于表2-4中。

比較表2-1和表2-4的數(shù)值，可以看到牛頓迭代法的收斂速度明顯快于對(duì)分法。牛頓迭代法也有局限性。在牛頓迭代法中，選取適當(dāng)?shù)跏贾祒0是求解的前提，當(dāng)?shù)某跏贾祒0在某根的附近時(shí)迭代才能收斂到這個(gè)根，有時(shí)會(huì)發(fā)生從一個(gè)根附近跳向另一個(gè)根附近的情況，尤其在導(dǎo)數(shù)數(shù)值很小時(shí)，如圖2-5所示。x2yx1xx0x2x0x3x1如果f(x)=0沒有實(shí)根，初始值x0是實(shí)數(shù)，則迭代序列不收斂。圖2-6給出迭代函數(shù)f(x)=2+x2，初始值x0=2的發(fā)散的迭代序列。圖2-5圖2-62.12.82.72.52.32.2總目錄2.92.62.42.7割線法

在牛頓迭代格式中：

，用差商代導(dǎo)數(shù)，并給定初始值x0和x1

，那么迭代格式可寫成如下形式：上式稱為割線法。用割線法迭代求根，每次只需計(jì)算一次函數(shù)值，而用牛頓迭代法每次要計(jì)算一次函數(shù)值和一次導(dǎo)數(shù)值。但割線收斂速度稍慢于牛頓迭代法，割線法為1.618階迭代方法。做過兩點(diǎn)(x0,f(x0))和(x1,f(x1))的一條直線（弦），該直線與x軸交點(diǎn)就是生成的迭代點(diǎn)x2，再做過(x1,f(x1))和(x2,f(x2))的一條直線，x3是該直線與x軸的交點(diǎn)，繼續(xù)做下去得到方程的根f(a)=0，如圖2.4所示。例2.5：用割線法求方程的根，取x0=1.5，x1=4.0。解：

計(jì)算結(jié)果列于表

VB程序2.12.82.72.52.32.2總目錄2.92.62.42.8非線性方程組的牛頓方法

為了敘述的簡(jiǎn)單，我們以解二階非線性方程組為例演示解題的方法和步驟，類似地可以得到解更高階非線性方程組的方法和步驟。設(shè)二階方程組

其中x,y為自變量。為了方便起見，將方程組寫成向量形式：將在（x0,y0）附近進(jìn)行二元泰勒展開，并取其線性部分，得到下面方程組：令則有

2.12.82.72.52.32.2總目錄2.92.62.42.8非線性方程組的牛頓方法如果再將原方程組在u1處進(jìn)行二元泰勒展開，并取其線性部分，得到下面方程組：解出得出

繼續(xù)做下去，每一次迭代都是一個(gè)方程組

2.12.82.72.52.32.2總目錄2.92.62.4即為止。2.8非線性方程組的牛頓方法例2.6：求解下面非線性方程組取初始值解：解方程得繼續(xù)做下去，直到時(shí)停止。

2.12.82.72.52.32.2總目錄2.92.62.42.9化工生產(chǎn)中非線性方程組求解應(yīng)用實(shí)例

在化工生產(chǎn)中,為了求解反應(yīng)前后各物料的濃度,常常要聯(lián)立求解一些非線性方程組,這些方程組難以用常規(guī)的解析方法求解,一般只能利用數(shù)值求解的方法加以求解。下面是在合成氨生產(chǎn)中利用非線性方程組求解方法求解烴類蒸氣轉(zhuǎn)化反應(yīng)前后各物料濃度的實(shí)例。例2.7

在合成氨生產(chǎn)中,