化工應用數學-2第二章-數學模型基礎

1、第二章 化工數學模型基礎,任課老師：程道建 副教授E-mail: chengdj@mail.buct.edu.cn,本章內容,2.1 擬合建模方法2.2 迭代法求解非線性方程2.3 求解線性方程組2.3 求微分方程數值解,2.1 擬合建模方法,示例 實驗測得二甲醇（DME）的飽和蒸氣壓和溫度的關系（見表），其具體關系是什么樣的呢？,,化工實際問題的提出,,1、線性擬合（一元一次）,給定一組數據（xi, yi）,i=1,

2、 2 , …, m ，作擬合直線 p (x)=a + bx , 均方誤差為 ：,按二元函數求極值的理論，Q (a , b)的極小值需滿足方程組：,,解此聯立方程：,2.1 擬合建模方法,,2、非線性擬合（一元二次）,給定數據（xi ,yi), i=1, 2 , …, m ,用二次多項式函數擬合這組數據。 設 ，作出擬合函數與數據序列的均方誤差表達式,由數學知識可知，Q(

3、a0 ,a1 ,a2 )的極小值滿足 ：,2.1 擬合建模方法,,2、非線性擬合（一元二次）,整理得二次多項式函數擬合的滿足條件方程：,解此方程得到在均方誤差最小意義下的擬合函數p ( x )。上式稱為多項式擬合的法方程。法方程的系數矩陣是對稱的。當擬合多項式n > 5時，法方程的系數矩陣是病態(tài)的，使用直接迭代法求解線性方程時會發(fā)散，要采用一些特殊算法（ Newton迭代法） 。,2.1 擬合建模方法,示例解決,,1）采用線性擬合

4、（一元一次） ：由表的數據觀測可得，DME的飽和蒸氣壓和溫度有正相關關系，如果以函數p=a+bt來擬合,則擬合函數是一條直線。通過計算均方誤差Q ( a , b )最小值而確定直線方程。,擬合得到直線方程為： 相關系數R為0.97296， 平均絕對偏差SD為0.0707。,2.1 擬合建模方法,示例解決,,2）非線性擬合（一元二次）

5、：,通過計算下述均方誤差 擬合得二次方程為 相關系數R為0.99972，平均絕對偏差SD為0.00815,具體擬合曲線見右圖。,2.1 擬合建模方法,示例解決,,（相關系數R為0.99972，平均絕對偏差SD為0.00815 ）,DME飽和蒸氣壓和溫度之間的二次擬合,DME飽和蒸氣壓和溫度之間的一次線性擬合,（相關系數R為0.97296，平均絕對偏差SD為0.0707）,結論：通過比較左右圖以及各自的

6、相關系數和平均絕對偏差可知，發(fā)現對于DME飽和蒸氣壓和溫度之間的關系，用二次曲線擬合優(yōu)于線性擬合。,2.1 擬合建模方法,3、多變量的曲線擬合,前面介紹的曲線擬合方法只涉及單變量函數的曲線擬合，但實際在化工實驗數據處理及模型參數擬合時，通常會碰到多變量的參數擬合問題。一個典型的例子是傳熱實驗中努塞爾數、雷諾數及普朗特數之間的擬合問題： 根據若干組實驗測得的數據（Nu，Re，Pr），如何求出式中的參數c1、c2、c3，這

7、是一個有兩個變量的參數擬合問題。,2.1 擬合建模方法,3、多變量的曲線擬合,擬合式不是線性的。需要做變換，等式兩邊取對數得：,對于給定的序列x1i,x2i,yi,i=1,2,3…m,設擬合后的函數形式為： 均方誤差為： 由多元函數極值原理， Q(a0,a1,a2)最小條件為：,2.1 擬合建模方法,3、多變量的曲線擬合,整理得多變量一次多項式函數擬合的法方程 通過求解一元線

8、性方程組就可以得到多變量函數線性 擬合時的參數。,2.1 擬合建模方法,示例解決,,根據某傳熱實驗測得如下數據，請用下列方程的形式擬合實驗曲線。,擬合式的形式為：擬合式不是線性的。需要做變換，等式兩邊取對數得：,2.1 擬合建模方法,示例解決,,1）首先對Nu，Re，Pr取對數，得到,2）通過上面公式求解：a0 = 0.022992a1 = 0.8a2 = 0.3,2.1 擬合建模方法,多項式擬合有時實驗數據表現

9、為一曲線，相應的擬合函數未知，需要一種普適的函數擬合曲線。常用方法之一就是用多項式擬合。原則上任何連續(xù)函數均可用多項式展開：若將變量進行變換： 則多項式化為多元一次函數：,,2.1 擬合建模方法,,多項式擬合可用Excel中的LINEST()函數和“回歸”求多項式的參數b、a1、a2…an及其回歸統(tǒng)計。通常到三次方就有中等精度。在實際工作中，在滿足擬合精度的前提下多項式的階數要盡可能的低。對于N個數據點，用于擬合的多項

10、式最高階數為N-1,2.1 擬合建模方法,化工實際問題的提出,非線性方程問題無論是從理論上還是從計算公式上，都要比線性方程復雜的多。而對于具體的化工問題，初值和求解范圍常常可根據具體的化工知識來決定。常見的雷諾數和摩擦系數關系方程在雷諾數低于4000時有以下關系式：,(2-1),已知雷諾數Re，如何根據公式（2-1）求出摩擦系數λ，這是我們在管路設計中必須首先解決的問題。,2.2 迭代法求解非線性方程,2.2 迭代法求解非線性方程,化

11、工實際問題的提出,2.2 迭代法求解非線性方程,化工實際問題的提出,2.2 迭代法求解非線性方程,化工實際問題的提出,,簡單迭代法實例應用已知在一溫度下，碳酸鈣在純水中的溶度積：以及碳酸根離子水解平衡常數：求碳酸鈣溶解度，mol/L。（忽略水的解離）,2.2 迭代法求解非線性方程,,解：由溶度積和平衡常數定義及物料平衡關系，得解方程組得,2.2 迭代法求解非線性方程,,簡單迭代法進行迭代求解：

12、用變量x代表碳酸鈣在純水中的溶解度等價變換選取迭代式,2.2 迭代法求解非線性方程,,,,Excel迭代控制參數設置,2.2 迭代法求解非線性方程,,Excel輸入計算公式“=SQRT(7.8E-7*SQRT(B2)+2.9E-9)”,回車，輸出結果如下；,2.2 迭代法求解非線性方程,,直接迭代法迭代公式通過代數恒等變形，將方程f(x)=0化成與之等價的方程x=φ(x)。令xk+1=φ(xk)，此式即為

13、直接迭代法的迭代公式。給定初值x0,由迭代公式產生點列{xk}k=0,1,2...,若 則x*即是方程f(x)=0的根。,2.2 迭代法求解非線性方程,,收斂判定對于方程 f(x)=0 構造的多種迭代格式xk+1=φ(xk) ，怎樣判斷構造的迭代格式是否收斂？收斂是否與迭代的初值有關？根據數學知識，我們可以直接利用以下收斂條件：（1） 當x∈[a,b]有a≦φ(x)≦b（2） φ(x)在[a,b]上可導，并且存

14、在正數L<1，使任意的x∈[a,b ]，有|φ(x)|≦L。 若滿足上述收斂條件，則在[a,b]上有唯一的點x*滿足x*=φ(x*) ，此時稱x* 為φ(x) 的不動點。,2.2 迭代法求解非線性方程,,迭代格式構造迭代格式xk+1=φ(xk)對任意初值x0∈[a,b] ，均收斂于φ(x) 的不動點x*，并有下面誤差估計式：構造收斂迭代格式有兩個要素：（1）等價形式x=φ(x)應滿足 |φ’(x*)|<1；

15、（2）初值必須取自x* 的充分小鄰域，其大小決定于函數f(x)，及做出的等價形式x=φ(x) 。,2.2 迭代法求解非線性方程,,迭代控制 在一般的迭代計算中，都會給定精度控制量ε，當時，迭代終止。在Excel迭代計算中，所求解的方程一般迭代次數較少，很快就能達到要求精度。,2.2 迭代法求解非線性方程,,實例應用例.求方程f(x)=x3-3x+1=0的三個實根（ 精度控制在小數點后6位）。解：恒等變形，有如下三種形式

16、由此構成三個迭代公式,2.2 迭代法求解非線性方程,,解：已知原方程在(-∞,+∞)內的三個實根分別分布在(-2,-1), (0,1)和(1,2)三個區(qū)間內。（1）采用迭代公式(a)和(b)，取同一初值x0=0.5,計算原方程在(0,1)內的實根,2.2 迭代法求解非線性方程,,在B3和D3中分別輸入迭代公式“=（B2^3+1）/3”、”=1/(3-D2^2)”，填充完畢。拖拉填充柄的結果如下圖。,2.2 迭代法求解非線性

17、方程,,（2）采用迭代公式(c)，分別取初值x0=1.5和x0=-1.5，計算原方程在(1,2)及(-2,-1)內的實根。,2.2 迭代法求解非線性方程,,（3）如果用迭代公式(a)，取初值x0=1.5,計算原方程在(1,2)內的實根，結果會如何呢？,不在（1,2）內,2.2 迭代法求解非線性方程,,迭代公式與初值選用同學們可自行驗證以下結論：對于迭代公式(a)和(b)，如果在區(qū)間(-2,-1)或(1,2)區(qū)間內取初值，則這兩個迭代公

18、式產生的點列不收斂于相應區(qū)間。同樣迭代公式(c)在(0,1)內取初值也是如此。由此可見，方程f(x)=0構成的直接迭代公式xk+1=φ(xk)是否收斂于該方程的根，既于初值x0的選取有關，也與迭代函數φ(x)有關。,2.2 迭代法求解非線性方程,,直接迭代的幾何意義將方程f(x)=0化為等價方程x=φ(x)，f(x)=0的根x*即是直線y=x與曲線y=φ(x)交點的橫坐標。右圖中，曲線y=φ(x)位于直線y=x的下方，任給定初值x

19、0,迭代式xk+1=φ(xk)收斂于x*。,2.2 迭代法求解非線性方程,,而右圖曲線y=φ(x)位于直線y=x的上方，任給定初值x0,迭代式xk+1=φ(xk)不收斂。,2.2 迭代法求解非線性方程,,松弛迭代法當迭代公式收斂很慢時，我們可以從x=φ(x)出發(fā)構造新的迭代形式，以加快收斂速度，這里先介紹松弛迭代法。迭代公式 其中， 為第k次迭代的預報值；

20、 稱為松弛因子。,2.2 迭代法求解非線性方程,,實例應用用松弛迭代法求解例3，以直接迭代公式(a)出發(fā)構造松弛迭代公式三個區(qū)間(-2,-1),(0,1)和(1,2)所取初值仍為-1.5,0.5,1.5。,2.2 迭代法求解非線性方程,,,埃特金（Aitken）迭代法為了避免確定松弛因子 的麻煩，Aitken又對松弛迭代法做了改進，提出了兩次校正的

21、加速迭代公式，即為下式。,2.2 迭代法求解非線性方程,,實例應用同樣用埃特金（Aitken）迭代法來解例3，基礎迭代公式仍是公式（a），初值選取不變。,2.2 迭代法求解非線性方程,,,牛頓（Newton）迭代法迭代公式不同于直接迭代法是將非線性方程f(x)=0恒等變形得到迭代公式，牛頓法則是將非線性方程f(x)=0在x0點展開，即并令用線性方程p(x)=0近似代替非線性方程f(x)=0，從中解得,2.2 迭代法求解非線

22、性方程,,令 作為f(x)=0的根的第一級近似值。一般地，記作為方程f(x)=0的根的第k+1級近似值，此式即為牛頓迭代公式。,2.2 迭代法求解非線性方程,,實例應用用牛頓（Newton）迭代法解例3，計算迭代公式：初值選取不變。,2.2 迭代法求解非線性方程,,牛頓迭代幾何意義方程f(x)=0的解就是曲線y=f(x)與x軸交點的橫坐標x*。設xk為初值，過點(xk,f(xk

23、))作y=f(x)的切線，則切線方程為令與 x軸的交點橫坐標為,令 ，即是第K次迭代點。,2.2 迭代法求解非線性方程,,牛頓迭代法的應用求重根求f(x)=0的m重根的牛頓公式，收斂速度快于單根公式。求復根牛頓迭代公式 不僅可以求單實根，也可以用來求復根，但初值必須是復數。,2.2 迭代法求解非線性方程,,牛

24、頓下降法牛頓迭代法在單根附近具有較快的收斂速度（至少平方收斂）。因此在運用該方法時初始點x0應取在單根x*的附近。當x0偏離x*較遠時，在某些情形下就不能保證收斂。此種情形下，通常采用牛頓下降法，其迭代公式為式中，λ∈(0,1]為阻尼因子，應滿足|f(xk+1)|<|f(xk)|。,2.2 迭代法求解非線性方程,,割線法迭代公式在牛頓迭代法中，需要先求得非線性方程的導數，但有時導數并不好求。針對此問題，可用一階差

25、商代替牛頓迭代公式中的導數f’(x)，就得到割線法迭代公式割線法中需要給定兩個初值x0和x1。,2.2 迭代法求解非線性方程,,實例應用同樣用割線法解例3， 區(qū)間(-2,-1),(0,1)和(1,2)上的初值分別選為{-1.5，-1.51}，{0.5,0.49}和{1.5,1.49}。,2.2 迭代法求解非線性方程,2.3 求解線性方程組,化工實際問題的提出,化工實際問題的提出,2.3 求解線性方程組,,線性方程組線性

26、方程組是各個方程關于未知量均為一次的方程組。其一般形式如下：其中aij、bi為已知常數，xj為未知數。,2.3 求解線性方程組,用線性代數中的概念來表達，則線性方程組可以寫成：A 是m×n矩陣，x是含有n個元素列向量，b是含有m個元素列向量。,,2.3 求解線性方程組,行列式解法根據Cramer法則，線性聯立方程有唯一解的條件是其系數行列式為非零值：此時方程的解為：,,2.3 求解線性

27、方程組,,|Di|是方程組的系數行列式，但其中的第i列元素被常數列陣b所取代。例如|D1|是用常數b列陣取代系數矩陣的第一列元素所得的行列式：,2.3 求解線性方程組,,示例演示例.利用Excel求下面三元一次方程組的解。,2.3 求解線性方程組,,Excel操作步驟：1）輸入系數矩陣A，列向量b。利用函數MDETERM( )計算系數行列式|A|=-9，方程組有唯一解。2）在計算行列式|D1|、|D2|和|D3|的值。3）由

28、計算得方程組的解。,2.3 求解線性方程組,,2.3 求解線性方程組,矩陣解法線性方程組的矩陣形式：經變換可得：因此，若系數矩陣A可逆，則解線性一次方程組就變成簡單的矩陣乘法。,,2.3 求解線性方程組,,示例演示同樣用矩陣法求解例3.Excel造作步驟：1）輸入系數矩陣A，列向量b。2）利用函數MINVERSE( )求系數矩陣的逆矩陣。3）利用函數MMULT( )求得方程組的解。注意：三鍵確認。,2.3 求解線性

29、方程組,,2.3 求解線性方程組,,行列式法與矩陣法行列式法和矩陣法所能求解的方程組的特點？多元一次非齊次系數矩陣為方陣系數矩陣可逆或行列式不為0。,2.3 求解線性方程組,,Newton-Raphson迭代法N-R迭代法原理N-R迭代法屬于間接解非線性方程組的方法。原理類似于一元方程求解。然而因有不止一個自變量，導數變?yōu)槠珜?，迭代公式成為迭代增量的線性方程組。設有非線性方程組：,Excel解非線性方程組,Excel解非

30、線性方程組,對方程組的每一個方程在其近似解 處用Teller級數展開，只取線性項得：,,Excel解非線性方程組,其中 是方程 在 處的一階偏導數。因此原非線性方程組就轉化為△xi的線性方程組：,,Excel解非線性方程組,將初始值 代入函數fi和偏導數 ，然后用解線性方程組的方法計算得到 ，從而得到自變量x

31、i的第一次迭代值：繼續(xù)迭代，可得：,,Excel解非線性方程組,迭代上式，直至：δ 為一小正數。 即為滿足指定精度的原方程的解。,,Excel解非線性方程組,,示例演示例.用迭代法求解下面非線性方程組的解.(計算精度10-15),Excel解非線性方程組,,首先,每個方程對各自變量求偏導：,Excel解非線性方程組,,將求的偏導代入方程組(a),于是有線性方程組(b):,,Excel解非線

32、性方程組,Excel操作步驟:1）輸入x1、x2、x3的初始值0.5。2）計算線性方程組（b）的系數矩陣和常數列矩陣。3）用矩陣法解線性方程組（b），輸入公式“=MMULT(MINVERSE(D2:F4),G2:G4)”，三鍵確認，得第一次迭代改變量Δxi1。4）用矩陣加法得第一次迭代后的自變量xi1(=xi0+Δxi1)。5）重復步驟2)、3)和4)，直至兩次迭代改變量小于10-15 。,,Excel解非線性方程組,2.4

33、求微分方程數值解,化工實際問題的提出,例子1. 方形蓄水箱加水和排水問題要得到水位高度和流量的關系（方形水槽）,2.4 求微分方程數值解,化工實際問題的提出,例子2. 三角剖面水槽進水和排水問題 (去年考試題目)要得到水位高度和流量的關系（三角剖面水槽）,2.4 求微分方程數值解,1. 原理 化學反應動力學方程一般用一階微分方程表示，其通式為： [a, b]是自變量的定義域，f(x, y)為已知函數。,,,,,,

34、,,,,1. 原理 若一函數y=F(x)代入微分方程，使得式(1)成立，即: 則該函數就是微分方程的解。,,,,,,,,,,2.4 求微分方程數值解,1. 原理 解微分方程用到積分，得到的解析式是含一個任意常數C的通解。若有初始條件（初值）： y0＝F (x0) 則通解的常數C可以確定，得到特解。,,,,,,,,,,2.4 求微分方程數值解,1. 原理 例如對于一級反

35、應：A→ B，其速率方程為： [A]t是反應物A在時間t時的濃度，其通解為：,,,,,,,,,,,,2.4 求微分方程數值解,1. 原理 若已知t＝0時A的濃度為[A]0，則在初始條件下的特解為：,,,,,,,,,,,2.4 求微分方程數值解,1. 原理 并非所有微分方程都有解析解。事實上除了一些簡單的基元反應，大多數反應動力學難以得到解析解或解析式很復雜，甚至不存在解析解。此時必須求助于數值解。 另一

36、方面反應動力學關心的問題是在t時刻反應體系內各物質的濃度[A]t , [B]t , ……，有足夠精度的近似解即可。,,,,,,,,,,,2.4 求微分方程數值解,1. 原理 數值解在化工“三傳一反”的模型化中都廣泛應用。 從小試管放大到工業(yè)反應器：反應動力學模型+傳遞模型。 數值解有普適性，可用于復雜的微分方程體系及任何初始條件。,,,,,,,,,,,2.4 求微分方程數值解,2. Euler法 常微分方

37、程數值解采用離散方法，即找出一種有效的數值計算方法，計算自變量的離散點：x0，x1，x2，…，xn以及對應的y近似值y0，y1，y2，…，yn。最簡單的是Euler法，通常取等間距的x值： x1-x0 =x2-x1=…= xn-xn-1=h h稱為步長。,,,,,,,,,,,2.4 求微分方程數值解,2. Euler法 以初始條件(x0, y0)代入微分方程式(1)，得到初始點P0

38、(x0, y0)的斜率：以差商近似微商：,,,,,,,,,,,2.4 求微分方程數值解,2. Euler法 因此有： 即x向前移動h，在x1處得到y(tǒng)1。 同樣可以由x1, y1代入式(1)得到y(tǒng)2。依此類推，可得到所有x0，x1，x2，…，xn對應的y近似值y0，y1，y2，…，yn。這些數值點連成一個折線函數，用以代替原來的曲線函數。,,,,,,,,,,,2.4 求微分方程數值解,2. Euler法,,,

39、,,,,,,,,2.4 求微分方程數值解,2. Euler法 Euler法簡單易懂，幾何意義明確，但誤差太大。由圖可知，由于誤差的累加，隨著x向前推進，折線越來越偏離原來的曲線。,,,,,,,,,,,2.4 求微分方程數值解,3. Euler法示例例：對一級反應動力學方程 設定：初始濃度[A]0=0.2000mol/L，反應速率常數k=0.01s-1。由上式可推得其求濃度近似值的遞推公式為：,,,,,,,,,,,2

40、.4 求微分方程數值解,3. Euler法示例具體解法如下：1）A列輸入時間，間隔為20s，直到140s。2）B列用Euler法計算在時間t時的濃度[A]t。在B5單元格輸入初始濃度0.2000mol/L。B6輸入遞推公式:=B5-B5*$D$1*$D$2。3）選定B6單元格，向下拖拽到B12。由于B5單元格為相對引用，每下一格，濃度值更新一次，得到各個時刻由Euler法計算的濃度。,,,,,,,,,,,2.4 求微分方程數值解

41、,3. Euler法示例4）C列是根據 計算的解析濃度值。在C5單元格輸入公式： =$B$5*EXP(-$D$1*A5)。 然后自動填充得其余值。5）D列為Euler法計算值與解析值的相對誤差: =100×(解析值-Euler值)/解析值。,,,,,,,,,,,2.4 求微分方程數值解,3. Euler法示例6) 減小步長h，可改進精度。如步長分別為5和

42、1時，在140s的誤差分別降為3.56％和0.70％。,,,,,,,,,,2.4 求微分方程數值解,4. Runge-Kutta法 Euler法產生較大誤差的原因是f(x,y)為曲線，用Teller公式展開： Euler公式只取了線性項(斜率），忽略了高次項。,,,,,,,,,,,2.4 求微分方程數值解,4. Runge-Kutta法 R-K法則包括了Teller展開式的高次項，其中最常用的是四階R-K公

43、式。在遞推公式里x取值為：xi, xi +h/2, xi +h。則微分方程通式的四階R-K公式為：,,,,,,,,,,,2.4 求微分方程數值解,4. Runge-Kutta法微分方程通式的四階R-K公式為：,,,,,,,,,,,,,2.4 求微分方程數值解,5. Runge-Kutta法示例 例子同Euler法一樣。由于微分方程僅涉及因變量y，R-K法四項表達式簡化為：,,,,,,,,,,,,2.4 求微分方程數值解,5.

44、Runge-Kutta法示例具體步驟：1）在F5單元格輸入初始濃度0.2000mol/L。2）根據T1-T4的計算公式， 在B6輸入：=-$D$1*$D$2*F5 在C6輸入：=-$D$1*$D$2*(F5+B6/2) 在D6輸入：=-$D$1*$D$2*(F5+C6/2) 在E6輸入：=-$D$1*$D$2*(F5+D6),,,,,,,,,,,,2.4 求微分方程數值解,5. Runge-Kutt

45、a法示例具體步驟：3）根據遞推公式，在F6輸入：=F5+(B6+2*C6+2*D6+E6)/64）選定區(qū)域B6:F6，拖拽填充柄到第12行。5）通過與解析解的相對誤差可以發(fā)現，R-K法相當精確，完全可以滿足一般反應動力學研究的需要。,,,,,,,,,,,,2.4 求微分方程數值解,5. Runge-Kutta法示例,,,,,,,,,,,,2.4 求微分方程數值解,本章小結,2.1 擬合建模方法2.2 迭代法求解非線性方程2