1、第4章 非線性方程求解,4.1 化工實際問題的提出 4.2 實根的對分法 4.3直接迭代法 4.4松弛迭代法 4.5韋格斯坦法 4.6牛頓迭代法 4.7 割線法 4.8非線性方程組的牛頓方法 4.9 化工生產(chǎn)中非線性方程組求解應(yīng)用實例,4.1 化工實際問題的提出,求解非線性方程是化工設(shè)計及模擬計算中必須解決的一個問題。與線性方程相比，非線性方程問題無論是從理論上還是從計算公式上，都要比線性方程復(fù)雜的多。對于一般的非線性f(x)=0，計算方程的根既無一定章程可循也無直接方法可言。例如，求解高次方程組7x6-x3+x-1.5=0的根，求解含有指數(shù)和正弦函數(shù)的超越方程ex-sin(x)=0的零點(diǎn)。解

2、非線性方程或非線性方程組也是計算方法中的一個主題。一般地，我們用符號f(x)來表示方程左端的函數(shù)，方程的一般形式表示為f(x)=0，方程的解稱為方程的根或函數(shù)的零點(diǎn)。 通常，非線性方程的根不止一個，而任何一種方法只能算出一個根。因此，在求解非線性方程時，要給定初始值或求解范圍。而對于具體的化工問題，初值和求解范圍常常可根據(jù)具體的化工知識來決定。常見的雷諾數(shù)和摩擦系數(shù)關(guān)系方程在雷諾數(shù)低于4000時有以下關(guān)系式： (4-1),4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.1 化工實際問題的提出,這是一個典型的非線性方程。我們在管路設(shè)計中經(jīng)常碰到。當(dāng)我們已知雷諾

3、數(shù)Re，如何根據(jù)公式（4-1）求出摩擦系數(shù)，這是我們在管路設(shè)計中必須首先解決的問題。對于方程（4-1）而言，無法用解析的方法求出摩擦系數(shù)，只能用數(shù)值求解的方法。如用在下面即將介紹的松弛迭代法，假設(shè)： 則利用松弛迭代公式可得： 經(jīng)11次迭代可得摩擦系數(shù)為0.07593。 同樣，在n個組分的等溫閃蒸計算中，通過物料和相平衡計算，我們可得到如下非線性方程：,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.1 化工實際問題的提出,在方程（4-3 ）中只有是未知數(shù)，ki為相平衡常數(shù)，zi為進(jìn)料組分的摩爾濃度，均為已知數(shù)。和上面的情況一樣，方程（4-3 ）也無法直接解析

4、求解，必須利用數(shù)值的方法，借助于計算機(jī)方可精確的計算。對于這個問題的求解，可利用我們下面介紹的牛頓迭代法進(jìn)行計算，也可利用其他迭代公式進(jìn)行計算，如采用牛頓迭代公式，則可以得到如下的具體迭代公式： （4-4） 飽和蒸氣壓是我們經(jīng)常要用到的數(shù)據(jù)，雖然我們可以通過實驗測量來獲取飽和蒸氣壓的數(shù)據(jù)，但我們通常利用前人已經(jīng)測量得到的數(shù)據(jù)或回歸的公式來獲取，這可以減輕我們大量的基礎(chǔ)實驗工作。公式（4-5）是一種常用的飽和蒸氣壓計算公式： (4-5) 其中p為飽和蒸氣壓，單位為mmHg，T為溫度，單位為K，A、B、C、D為已知系數(shù)。要想得到某一溫度下的飽和蒸氣壓，直接利用公式（4-5）是無法得到的。因為公式（

5、4-5）兩邊都有未知變量，并且無法用解析的方法求解，必須用數(shù)值計算的方法求解。通過上面的一些例子，我們可以發(fā)現(xiàn)，如果沒有適當(dāng)?shù)氖侄魏娃k法來求解非線性方程，那么化學(xué)化工中的許多研究、設(shè)計等工作將無法展開，這勢必影響化學(xué)化工的發(fā)展，下面我們將介紹一些實用的非線性方程求解方法，并提供計算機(jī)程序。,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.2 實根的對分法,4.2.1 使用對分法的條件 4.2.2對分法求根算法 4.2.3對分法VB程序清單,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.2.1 使用對分法的條件,對分法或稱

6、二分法是求方程近似解的一種簡單直觀的方法。設(shè)函數(shù)f(x)在a,b上連續(xù)，且f(a)f(b)0，則f(x)在a,b上至少有一零點(diǎn)，這是微積分中的介值定理，也是使用對分法的前提條件。計算中通過對分區(qū)間，逐步縮小區(qū)間范圍的步驟搜索零點(diǎn)的位置。 如果我們所要求解的方程從物理意義上來講確實存在實根，但又不滿足f(a)f(b)0，這時，我們必須通過改變a和b的值來滿足二分法的應(yīng)用條件。,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.2.2對分法求根算法,計算f(x)=0的一般計算步驟如下： 1、輸入求根區(qū)間a,b和誤差控制量，定義函數(shù)f(x)。 2、判斷: 如果f(a

7、)f(b)0則轉(zhuǎn)下，否則，重新輸入a和b的值。 3、計算中點(diǎn) x=(a+b)/2以及f(x)的值 分情況處理 （1）|f(x)|：停止計算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*。 右圖給出對分法的示意圖。,x3=(x0+x2)/2,x2= (x0+x1)/2,x0 x3 x1 x1,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.2.3對分法VB程序清單,Private Sub Command1_Click() Dim x1, x2, x, y1

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

9、 x y2 = y GoTo 100 Else x1 = x y1 = y GoTo 100 End If End If End Sub Public Function f(x) Dim y y = x 3 + x 2 - 1 f = y End Function,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.2.3對分法求解實例,用對分法求 在區(qū)間1,2之間的根。 解： (1) f(1)= -2.8，f(2)=0.3，由介值定理可得有根區(qū)間a,b=1,2。 (2) 計算x2=(1+2)/2=1.5，f(1.5)= -0.45，有根區(qū)間a,b=1.5,2。

10、 (3) 計算x3=(1.5+2)/2=1.75，f(1.75)=0.078125，有根區(qū)間a,b=1.5,1.75。 一直做到|f(xn)|(計算前給定的精度)或|a-b|時停止。詳細(xì)計算結(jié)果見表4-1。 對分法的算法簡單，然而，若f(x)在a,b上有幾個零點(diǎn)時，如不作特殊處理只能算出其中一個零點(diǎn)；另一方面，即使f(x)在a,b上有零點(diǎn)，也未必有f(a)f(b)0。這就限制了對分法的使用范圍。對分法只能計算方程f(x)=0的實根。 對于多個零點(diǎn)的方程，我們可以通過將給定的區(qū)間a,b進(jìn)行細(xì)分，然后在細(xì)分后的區(qū)間內(nèi)用二分法分別求解，從而得到多個零點(diǎn)。例如求方程在0-30內(nèi)的所有根。需要對二分法進(jìn)

11、行以下處理：即先給定一個a，本例中為0，然后不斷增加，直到找到一個b，使f(a)f(b)0，調(diào)用二分法，計算在a，b范圍內(nèi)的根，然后將b作為a，重復(fù)上面的工作，直到計算范圍超出30為止。,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,VB調(diào)用,4.3直接迭代法,對給定的方程f(x)=0，將它轉(zhuǎn)換成等價形式： 。給定初值x0，由此來構(gòu)造迭代序列 ，k=1,2,，如果迭代收斂，即 有 ,則就是方程f(x)=0的根。在計算中當(dāng) 小于給定的精度控制量時，取 為方程的根。 例如，代數(shù)方程x3-2x-10=0的三種等價形式及其迭代格式如下： 對于方程 構(gòu)造的多種迭代格式

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

13、價形式 。,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,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）將迭代格式寫為 迭代格式不能保證收斂,但并不一定不收斂。 VB程序界面：,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,【例4-3】（V312）試編程計算0.01mol/L的 HAc(Ka1.75E-5）溶液的pH值。 在新建窗體上添加一個Command1命令框（在屬性窗口將Command

14、1的屬性Caption改為計算），對命令框?qū)ο?Command1的事件 Click、Private Function f(X）和 Private Sub ITR（XO，E，x）通用過程輸入如下代碼。在窗體“通用聲明”段用 Dim定義 KA，C能被本模塊的 Private Function f（x）過程存取。,Dim KA，C Private Sub Command1_Click） VB調(diào)用 C=0.01 KA0.0000175 XO=Sqr（K*C） E1.0E-5 Call ITR（X0,E,x） Print pH=”；Int（Log（x）Log（10）*10005）100 End Sub

15、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,4.4松弛迭代法,有些非線性方程或方程組當(dāng)用上一節(jié)中的直接迭代法求解時,迭代過程是發(fā)散的。這時可引入松弛因子,利用松弛迭代法。通過選擇合適的松弛因子,就可以使迭代過程收斂。松弛法的迭代公式如下: （4-7） 由上式可知，當(dāng)松弛因子等于1時，松弛迭代變?yōu)橹苯拥．?dāng)松弛因子大于1時松弛法使迭代步長加大,可加速迭代，但有可能使原來收斂的迭代變成發(fā)散。當(dāng)01時

16、, 松弛法使迭代步長減小,這適合于迭代發(fā)散或振蕩收斂的情況,可使振蕩收斂過程加速。當(dāng)0時,將使迭代反方向進(jìn)行,可使一些迭代發(fā)散過程收斂。 松弛法是否有效的關(guān)鍵因子是松弛因子的值能否正確選定。如果值選用適當(dāng),能使迭代過程加速,或使原來不收斂的過程變成收斂。但如果值選用不合適，則效果相反,有時甚至?xí)乖瓉硎諗康倪^程變得不收斂。松弛因子的數(shù)值往往要根據(jù)經(jīng)驗選定,但選用較小的松弛因子,一般可以保證迭代過程的收斂。,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.4松弛迭代法,例:用松弛迭代法求解下面非線性方程組,并分析松弛因子對迭代次數(shù)及收斂過程的影響。已知迭代初值

17、x和y均為0,收斂精度=0.001 解: 取以下迭代表達(dá)式: 若取松弛因子為1.1則其迭代過程如左表。 若改變松弛因子,迭代過程及迭代所需的次數(shù)亦將發(fā)生變化，詳見右表。 由右表數(shù)據(jù)可知,當(dāng)松弛因子小于1時,增大松弛因子,可加速迭代過程,減少 迭代次數(shù),但當(dāng)松弛因子大于1時,迭代次數(shù)反而增加,當(dāng)松弛因子達(dá)到1.56 時,迭代過程分散。,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,此法是一種迭代加速方法。在圖4-3中,從x0開始,曲線y=(x)和直線y=x之間的階梯形折線為單變量的直接迭代過程。其迭代步長為每個階梯的水平距離。 若利用圖4-3中曲線(x)上的兩個

18、點(diǎn)作一直線, 通過求該直線和y=x的交點(diǎn)C的橫坐標(biāo)xW作為新的迭代點(diǎn),這樣的迭代計算過程就是韋格斯坦法。即韋格斯坦法需要在已知兩點(diǎn)的前提下迭代計算第三點(diǎn)。一般第一點(diǎn)為人為設(shè)定,第二點(diǎn)利用直接迭代計算而得,則第三點(diǎn)就可以用韋格斯坦法迭代求解。已知A點(diǎn)的坐標(biāo)為(x1,y1),B點(diǎn)的坐標(biāo)為(x2,y2)。則過A、B兩點(diǎn)的直線方程為： 求該直線和y=x的交點(diǎn)可得: 韋格斯坦法的一般計算通式為,4.5韋格斯坦法,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.5韋格斯坦法,由上述公式可知, 韋格斯坦法也是一種松弛法,其松弛因子為 一般情況下,當(dāng)1k0時,迭代過程為單調(diào)

19、收斂過程。當(dāng)-1k0時,迭代過程為振蕩收斂過程,但當(dāng)k=1時,收斂將發(fā)散,故在編程計算時應(yīng)注意當(dāng)k=1時則取k=0進(jìn)行計算。,韋格斯坦法求解方程x3-2x2+5x-4=0根的QB程序見課本,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.6牛頓迭代法,4.6.1牛頓法的理論推導(dǎo) 4.6.2牛頓法的幾何意義,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,對方程f(x)=0可構(gòu)造多種迭代格式 ，牛頓迭代法是借助于對函數(shù)f(x)=0的泰勒展開而得到的一種迭代格式。 將f(x)=0在初始值x0做泰勒展開得： 取展開式的線性部分作為的

20、近似值，則有： 設(shè) 則 令 類似地，再將f(x)=0在x1作泰勒展開并取其線性部分得到： 一直做下去得到牛頓法的迭代格式： 牛頓迭代格式對應(yīng)于f(x)=0的等價方程為： 若b是f(x)的單根時， ，則有 ，只要初值x0充分接近b，牛頓迭代都收斂。牛頓迭代是二階迭代方法?？梢宰C明，b為f(x)的a重根時，迭代也收斂，但這是一階迭代，收斂因子為 ，若這時取下面迭代格式，它仍是二階方法：,4.6.1牛頓法的理論推導(dǎo),4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,以 為斜率作過(x0，f(x0)點(diǎn)的直線，即作f(x)在x0的切線方程： 令y=0，則在x1處的切線與x軸

21、的交點(diǎn)x1，即： 再作f(x)在x1處的切線，得交點(diǎn)x2，逐步逼近方程的根b。如圖4-4所示。 在區(qū)域x0,x0+h的局部“以直代曲”是處理非線性問題的常用手法。在 泰勒展開中，截取函數(shù)展開的線性部分替代 f(x)。,4.6.2牛頓法的幾何意義,圖4-4牛頓切線法示意圖,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.6.2牛頓法的幾何意義,例：用牛頓迭代法求方程 f(x)=x3-7.7x2+19.2x-15.3，在x0=1附近的零點(diǎn)。 解： 計算結(jié)果列于表4-4中。,比較表4-1和表4-4的數(shù)值，可以看到牛頓迭代法的收斂速度明顯快于對分法。牛頓迭代法也有局

22、限性。在牛頓迭代法中，選取適當(dāng)?shù)跏贾祒0是求解的前提，當(dāng)?shù)某跏贾祒0在某根的附近時迭代才能收斂到這個根，有時會發(fā)生從一個根附近跳向另一個根附近的情況，尤其在導(dǎo)數(shù)數(shù)值很小時，如圖4-5所示。,如果f(x)=0沒有實根，初始值x0是實數(shù)，則迭代序列不收斂。圖4-6給出迭代函數(shù)f(x)=2+x2，初始值x0=2的發(fā)散的迭代序列。,圖4-5,圖4-6,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,【例 4-4】(V214）1mL 0.2 mol/L HCl溶液中含Cu2+5 mg，若在室溫及1.013 X 105Pa下通入 H2S氣體達(dá)飽和（0.lmol/L)

23、則析出 CuS沉淀，達(dá)平衡時，求殘留在溶液中的Cu2+？ 在新建窗體上添加一個Commandl命令框（在屬性窗口將Commandl的屬性Caption改為計算），對命令框?qū)ο?Commandl的事件 Click、Pivate Sub NT（X0，E，x）、Private Function f(x)和Private Function g(x)通用過程輸入如下代碼。在窗體“通用聲明”段用Dim定義C1，C2,a,K變量為雙精度型能被本模塊的 Private Function f(x)和 Private Function g(x)過程存取。 VB程序,Dim C1，C2，a，K As Double

24、Private Sub Command1_Click（） C10.07874 C20.1 a0.2 E0.00001 K11330000000000 Call NT（X0,E,x） Print x Print Format(x,”0.000e-00”） end sub Private Function f(x) f= 4 * x * x - (4 * (a + 2 * C1) + K * C2) * x + (4 * C1 * C1 + 4 * a * C1 + a * a) End Function Private Function g(x) g=8*X-(4*(a+2*C1)+K*C2)

25、End Function Private Sub NT（XO，E，x） XXO DO X0=x x=x0-f(x0)/g(x0) loop while Abs(x-x0)e End sub,4.7 割線法,在牛頓迭代格式中： ，用差商 代導(dǎo)數(shù) ，并給定初始值x0和x1 ，那么迭代格式可寫成如下形式： 上式稱為割線法。 用割線法迭代求根，每次只需計算一次函數(shù)值，而用牛頓迭代法每次 要計算一次函數(shù)值和一次導(dǎo)數(shù)值。但割線收斂速度稍慢于牛頓迭代法，割 線法為1.618階迭代方法。 做過兩點(diǎn)(x0,f(x0)和(x1,f(x1)的一條直線（弦），該直線與x軸交點(diǎn)就是 生成的迭代點(diǎn)x2，再做過(x1,f(

26、x1)和(x2,f(x2)的一條直線，x3是該直線與x軸 的交點(diǎn)，繼續(xù)做下去得到方程的根f(a)=0，如圖4.4所示。 例4.5：用割線法求方程 的根， 取x0=1.5，x1=4.0。 解：,計算結(jié)果列于表 VB程序,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,【例 4-6】(V316）計算 500mL、1.0mol/L Na2S2O3溶液可以溶解多少克 AgBr。 在新建窗體上添加兩個標(biāo)簽框Label1、Label2（在屬性窗口將Label1的屬性Caption改為割線法，屬性Autosize設(shè)為True；在屬性窗口將Label2的屬性Caption置為空

27、）和一個Command1命令框（在屬性窗口將Command1的屬性Caption改為計算），對命令框?qū)ο驝ommand1的事件Click、Private Sub GX（XO，H，E，XZ）和 Private Function f(x)通用過程輸入如下代碼，在窗體“通用聲明”段用Dim定義c，k變量為單精度型能被本模塊的Private Function f（x）過程存取。 VB程序,Dim c,k as single Private sub command1_click() H=0.01 C=1# K1=29000000000000# Ksp=0.00000000000077 K=k1*ksp

28、X0=val(inputbox(“x0=“) Call gx(x0,h,e,x) X=int(x*100+0.5)/100 S=500*x/1000*188 W$=chr(13)+chr(10) Label2.caption=“x=“&x&w$&”s=“&s&”g” Print End sub Private function f(x) F=(4*k-1)*x*x-4*k*c*x+k*c*c End function,Private Sub gx(x0, h, e, x2) bx0 = f(x0) If bx0 0 Then X1 = x0 - h Else X1 = x0 + h End I

29、f Do x2 = X1 - f(X1) / (f(X1) - f(x0) * (X1 - x0) If Abs(x2 - X1) / x2 e Then Exit Do x0 = X1 X1 = x2 Loop End Sub next,4.8非線性方程組的牛頓方法,為了敘述的簡單，我們以解二階非線性方程組為例演示解題的方法和步驟，類似地可以得到解更高階非線性方程組的方法和步驟。 設(shè)二階方程組 其中x,y為自變量。為了方便起見，將方程組寫成向量形式： 將 在（x0,y0）附近進(jìn)行二元泰勒展開，并取其線性部分，得到下面方程組： 令 則有,4.1,4.8,4.7,4.5,4.3,4.2,總目錄,4.9,4.6,4.4,4.8非線性方程組的牛頓方法,如果 再將原方程組在u1處進(jìn)行二元泰勒展開，并取其線性部分，得到下面方程 組： 解出 得出 繼續(xù)做下去，每一次迭代都是一個方程組,4.1,4.8,4.7,4.