非線性方程求解

1、非線性方程求解第1頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.1 化工實(shí)際問(wèn)題的提出 求解非線性方程是化工設(shè)計(jì)及模擬計(jì)算中必須解決的一個(gè)問(wèn)題。與線性方程相比，非線性方程問(wèn)題無(wú)論是從理論上還是從計(jì)算公式上，都要比線性方程復(fù)雜的多。對(duì)于一般的非線性f(x)=0，計(jì)算方程的根既無(wú)一定章程可循也無(wú)直接方法可言。例如，求解高次方程組7x6-x3+x-1.5=0的根，求解含有指數(shù)和正弦函數(shù)的超越方程ex-sin(x)=0的零點(diǎn)。解非線性方程或非線性方程組也是計(jì)算方法中的一個(gè)主題。一般地，我們用符號(hào)f(x)來(lái)表示方程左端的函數(shù)，方程的一般形式表示為f(x)=0，方程的解稱為方程的根或函數(shù)的零

2、點(diǎn)。 通常，非線性方程的根不止一個(gè)，而任何一種方法只能算出一個(gè)根。因此，在求解非線性方程時(shí)，要給定初始值或求解范圍。而對(duì)于具體的化工問(wèn)題，初值和求解范圍常?？筛鶕?jù)具體的化工知識(shí)來(lái)決定。常見(jiàn)的雷諾數(shù)和摩擦系數(shù)關(guān)系方程在雷諾數(shù)低于4000時(shí)有以下關(guān)系式： (4-1)4.14.84.74.54.34.2總目錄4.94.64.4第2頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.1 化工實(shí)際問(wèn)題的提出 這是一個(gè)典型的非線性方程。我們?cè)诠苈吩O(shè)計(jì)中經(jīng)常碰到。當(dāng)我們已知雷諾數(shù)Re，如何根據(jù)公式（4-1）求出摩擦系數(shù)，這是我們?cè)诠苈吩O(shè)計(jì)中必須首先解決的問(wèn)題。對(duì)于方程（4-1）而言，無(wú)法用解析的方法

3、求出摩擦系數(shù)，只能用數(shù)值求解的方法。如用在下面即將介紹的松弛迭代法，假設(shè)： 則利用松弛迭代公式可得： 經(jīng)11次迭代可得摩擦系數(shù)為0.07593。 同樣，在n個(gè)組分的等溫閃蒸計(jì)算中，通過(guò)物料和相平衡計(jì)算，我們可得到如下非線性方程： 4.14.84.74.54.34.2總目錄4.94.64.4第3頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.1 化工實(shí)際問(wèn)題的提出 在方程（4-3 ）中只有是未知數(shù)，ki為相平衡常數(shù)，zi為進(jìn)料組分的摩爾濃度，均為已知數(shù)。和上面的情況一樣，方程（4-3 ）也無(wú)法直接解析求解，必須利用數(shù)值的方法，借助于計(jì)算機(jī)方可精確的計(jì)算。對(duì)于這個(gè)問(wèn)題的求解，可利用我們

4、下面介紹的牛頓迭代法進(jìn)行計(jì)算，也可利用其他迭代公式進(jìn)行計(jì)算，如采用牛頓迭代公式，則可以得到如下的具體迭代公式： （4-4） 飽和蒸氣壓是我們經(jīng)常要用到的數(shù)據(jù)，雖然我們可以通過(guò)實(shí)驗(yàn)測(cè)量來(lái)獲取飽和蒸氣壓的數(shù)據(jù)，但我們通常利用前人已經(jīng)測(cè)量得到的數(shù)據(jù)或回歸的公式來(lái)獲取，這可以減輕我們大量的基礎(chǔ)實(shí)驗(yàn)工作。公式（4-5）是一種常用的飽和蒸氣壓計(jì)算公式： (4-5) 其中p為飽和蒸氣壓，單位為mmHg，T為溫度，單位為K，A、B、C、D為已知系數(shù)。要想得到某一溫度下的飽和蒸氣壓，直接利用公式（4-5）是無(wú)法得到的。因?yàn)楣剑?-5）兩邊都有未知變量，并且無(wú)法用解析的方法求解，必須用數(shù)值計(jì)算的方法求解。通過(guò)上

5、面的一些例子，我們可以發(fā)現(xiàn)，如果沒(méi)有適當(dāng)?shù)氖侄魏娃k法來(lái)求解非線性方程，那么化學(xué)化工中的許多研究、設(shè)計(jì)等工作將無(wú)法展開(kāi)，這勢(shì)必影響化學(xué)化工的發(fā)展，下面我們將介紹一些實(shí)用的非線性方程求解方法，并提供計(jì)算機(jī)程序。 4.14.84.74.54.34.2總目錄4.94.64.4第4頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.2 實(shí)根的對(duì)分法 4.2.1 使用對(duì)分法的條件 對(duì)分法求根算法 對(duì)分法VB程序清單 4.14.84.74.54.34.2總目錄4.94.64.4第5頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.2.1 使用對(duì)分法的條件 對(duì)分法或稱二分法是求方程近似解的一

6、種簡(jiǎn)單直觀的方法。設(shè)函數(shù)f(x)在a,b上連續(xù)，且f(a)f(b)0，則f(x)在a,b上至少有一零點(diǎn)，這是微積分中的介值定理，也是使用對(duì)分法的前提條件。計(jì)算中通過(guò)對(duì)分區(qū)間，逐步縮小區(qū)間范圍的步驟搜索零點(diǎn)的位置。 如果我們所要求解的方程從物理意義上來(lái)講確實(shí)存在實(shí)根，但又不滿足f(a)f(b)0，這時(shí)，我們必須通過(guò)改變a和b的值來(lái)滿足二分法的應(yīng)用條件。 4.14.84.74.54.34.2總目錄4.94.64.4第6頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四對(duì)分法求根算法 計(jì)算f(x)=0的一般計(jì)算步驟如下： 1、輸入求根區(qū)間a,b和誤差控制量，定義函數(shù)f(x)。 2、判斷: 如果

7、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,xa,b，重復(fù)3（3）f(x)f(b)0：修正區(qū)間x,ba,b，重復(fù)3 4、輸出近似根x*。 右圖給出對(duì)分法的示意圖。 x3=(x0+x2)/2 x2= (x0+x1)/2 x0 x3 x1 x1 4.14.84.74.54.34.2總目錄4.94.64.4第7頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四對(duì)分法VB程序清單 Private Sub Command1_Click()Dim

8、 x1, x2, x, y1, y2, y, eer80 x1 = InputBox(x1)x2 = InputBox(x2)eer=inputbox(“eer”)y1 = f(x1)y2 = f(x2)If y1 * y2 0 Then GoTo 100Else Print please repeat input x1 and x2 GoTo 80End If100 x = (x1 + x2) / 2y = f(x)If Abs(y) = 0.001 Then Print the function root is ; x Print y=; yElse If y1 * y 0 Then x2

9、 = x y2 = y GoTo 100 Else x1 = x y1 = y GoTo 100 End IfEnd IfEnd SubPublic Function f(x)Dim yy = x 3 + x 2 - 1f = yEnd Function 4.14.84.74.54.34.2總目錄4.94.64.4第8頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四對(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,

10、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é)果見(jiàn)表4-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)的方程，我們可以通過(guò)將給定的區(qū)間a,b進(jìn)行細(xì)分，然后在細(xì)分后的區(qū)間內(nèi)用二分法分別求解，從而得到多個(gè)零點(diǎn)。例如求方程在0-30內(nèi)的所有根

11、。需要對(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為止。 4.14.84.74.54.34.2總目錄4.94.64.4VB調(diào)用第9頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.3直接迭代法 對(duì)給定的方程f(x)=0，將它轉(zhuǎn)換成等價(jià)形式： 。給定初值x0，由此來(lái)構(gòu)造迭代序列 ，k=1,2,，如果迭代收斂，即 有 ,則就是方程f(x)=0的根。在計(jì)算中當(dāng) 小于給定的精度控制量時(shí)，取 為方程的根。 例如，代數(shù)方程x3-2x-10=0的三種等價(jià)

12、形式及其迭代格式如下： 對(duì)于方程 構(gòu)造的多種迭代格式 ，怎樣判斷構(gòu)造的迭代格式是否收斂？收斂是否與迭代的初值有關(guān)？根據(jù)數(shù)學(xué)知識(shí)，我們可以直接利用以下收斂條件： 1、 當(dāng) 有 2、 在a,b上可導(dǎo)，并且存在正數(shù)L1，使任意的 ，有 則在a,b上有唯一的點(diǎn) 滿足 ， 稱 為 的不動(dòng)點(diǎn)。而且迭代格式對(duì)任意初值均收斂于的不動(dòng)點(diǎn)，并有下面誤差估計(jì)式： (4-6) 要構(gòu)造滿足收斂條件的等價(jià)形式一般比較困難。事實(shí)上，如果 為f(x)的零點(diǎn)，若能構(gòu)造等價(jià)形式 ，而 ，由 的邊疆性，一定存在的鄰域 ，其上有 ，這時(shí)若初值 迭代也就收斂了。由此構(gòu)造收斂迭代格式有兩個(gè)要素，其一，等價(jià)形式 應(yīng)滿足； 其二，初值必須取

13、自 的充分小鄰域，其大小決定于函數(shù)f(x)，及做出的等價(jià)形式 。4.14.84.74.54.34.2總目錄4.94.64.4第10頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四2.3直接迭代法例：求代數(shù)方程x3-2x-5=0，在x0=2附近的零點(diǎn)。 解：1）x3=2x+5 構(gòu)造的迭代序列收斂。取x0=2，則 準(zhǔn)確的解是x=2.09455148150。2）將迭代格式寫(xiě)為 迭代格式不能保證收斂,但并不一定不收斂。VB程序界面： 4.14.84.74.54.34.2總目錄4.94.64.4第11頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四【例4-3】（V312）試編程計(jì)算0.0

14、1mol/L的 HAc(Ka1.75E-5）溶液的pH值。 在新建窗體上添加一個(gè)Command1命令框（在屬性窗口將Command1的屬性Caption改為計(jì)算），對(duì)命令框?qū)ο?Command1的事件 Click、Private Function f(X）和 Private Sub ITR（XO，E，x）通用過(guò)程輸入如下代碼。在窗體“通用聲明”段用 Dim定義 KA，C能被本模塊的 Private Function f（x）過(guò)程存取。 Dim KA，C Private Sub Command1_Click） VB調(diào)用 C=0.01 KA0.0000175 XO=Sqr（K*C） E1.0E-5

15、 Call ITR（X0,E,x） Print pH=”；Int（Log（x）Log（10）*10005）100 End Sub Private Sub ITR（X0,E,X） xXO Do XOX X=f(XO） LOOp While Abs（XXO）E End Sub Private Function f(x） f= Sqr（KA*（CX） End Function第12頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.4松弛迭代法 有些非線性方程或方程組當(dāng)用上一節(jié)中的直接迭代法求解時(shí),迭代過(guò)程是發(fā)散的。這時(shí)可引入松弛因子,利用松弛迭代法。通過(guò)選擇合適的松弛因子,就可以使迭代過(guò)程

16、收斂。松弛法的迭代公式如下: （4-7） 由上式可知，當(dāng)松弛因子等于1時(shí)，松弛迭代變?yōu)橹苯拥?。?dāng)松弛因子大于1時(shí)松弛法使迭代步長(zhǎng)加大,可加速迭代，但有可能使原來(lái)收斂的迭代變成發(fā)散。當(dāng)01時(shí), 松弛法使迭代步長(zhǎng)減小,這適合于迭代發(fā)散或振蕩收斂的情況,可使振蕩收斂過(guò)程加速。當(dāng)k0時(shí),迭代過(guò)程為單調(diào)收斂過(guò)程。當(dāng)-1ke End sub第22頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.7 割線法 在牛頓迭代格式中： ，用差商代導(dǎo)數(shù) ，并給定初始值x0和x1 ，那么迭代格式可寫(xiě)成如下形式： 上式稱為割線法。 用割線法迭代求根，每次只需計(jì)算一次函數(shù)值，而用牛頓迭代法每次要計(jì)算一次函數(shù)值和

17、一次導(dǎo)數(shù)值。但割線收斂速度稍慢于牛頓迭代法，割線法為1.618階迭代方法。 做過(guò)兩點(diǎn)(x0,f(x0)和(x1,f(x1)的一條直線（弦），該直線與x軸交點(diǎn)就是生成的迭代點(diǎn)x2，再做過(guò)(x1,f(x1)和(x2,f(x2)的一條直線，x3是該直線與x軸的交點(diǎn)，繼續(xù)做下去得到方程的根f(a)=0，如圖4.4所示。例4.5：用割線法求方程 的根，取x0=1.5，x1=4.0。解： 計(jì)算結(jié)果列于表 VB程序4.14.84.74.54.34.2總目錄4.94.64.4第23頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四【例 4-6】(V316）計(jì)算 500mL、1.0mol/L Na2S2O

18、3溶液可以溶解多少克 AgBr。 在新建窗體上添加兩個(gè)標(biāo)簽框Label1、Label2（在屬性窗口將Label1的屬性Caption改為割線法，屬性Autosize設(shè)為T(mén)rue；在屬性窗口將Label2的屬性Caption置為空）和一個(gè)Command1命令框（在屬性窗口將Command1的屬性Caption改為計(jì)算），對(duì)命令框?qū)ο驝ommand1的事件Click、Private Sub GX（XO，H，E，XZ）和 Private Function f(x)通用過(guò)程輸入如下代碼，在窗體“通用聲明”段用Dim定義c，k變量為單精度型能被本模塊的PrivateFunction f（x）過(guò)程存取。

19、VB程序第24頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四Dim c,k as singlePrivate sub command1_click()H=0.01C=1#K1=29000000000000#Ksp=0.00000000000077K=k1*kspX0=val(inputbox(“x0=“)Call gx(x0,h,e,x)X=int(x*100+0.5)/100S=500*x/1000*188W$=chr(13)+chr(10)Label2.caption=“x=“&x&w$&”s=“&s&”g”P(pán)rintEnd subPrivate function f(x)F=(

20、4*k-1)*x*x-4*k*c*x+k*c*cEnd function第25頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四Private Sub gx(x0, h, e, x2)bx0 = f(x0)If bx0 0 ThenX1 = x0 - hElseX1 = x0 + hEnd IfDox2 = X1 - f(X1) / (f(X1) - f(x0) * (X1 - x0)If Abs(x2 - X1) / x2 e ThenExit Dox0 = X1X1 = x2LoopEnd Subnext第26頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.8非線性方程組

21、的牛頓方法 為了敘述的簡(jiǎn)單，我們以解二階非線性方程組為例演示解題的方法和步驟，類似地可以得到解更高階非線性方程組的方法和步驟。 設(shè)二階方程組 其中x,y為自變量。為了方便起見(jiàn)，將方程組寫(xiě)成向量形式： 將 在（x0,y0）附近進(jìn)行二元泰勒展開(kāi)，并取其線性部分，得到下面方程組： 令 則有 4.14.84.74.54.34.2總目錄4.94.64.4第27頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.8非線性方程組的牛頓方法如果再將原方程組在u1處進(jìn)行二元泰勒展開(kāi)，并取其線性部分，得到下面方程組： 解出 得出 繼續(xù)做下去，每一次迭代都是一個(gè)方程組 4.14.84.74.54.34.2總

22、目錄4.94.64.4即 為止。 第28頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.8非線性方程組的牛頓方法例4.7：求解下面非線性方程組 取初始值 解： 解方程得 繼續(xù)做下去，直到 時(shí)停止。 4.14.84.74.54.34.2總目錄4.94.64.4第29頁(yè)，共32頁(yè)，2022年，5月20日，3點(diǎn)40分，星期四4.9 化工生產(chǎn)中非線性方程組求解應(yīng)用實(shí)例 在化工生產(chǎn)中,為了求解反應(yīng)前后各物料的濃度,常常要聯(lián)立求解一些非線性方程組,這些方程組難以用常規(guī)的解析方法求解,一般只能利用數(shù)值求解的方法加以求解。下面是在合成氨生產(chǎn)中利用非線性方程組求解方法求解烴類蒸氣轉(zhuǎn)化反應(yīng)前后各物料濃度的實(shí)例。 例 4.8 在合成氨生產(chǎn)中,烴類蒸氣發(fā)生以下轉(zhuǎn)化反應(yīng): 已知進(jìn)料甲