化工應(yīng)用數(shù)學(xué)-2第二章-數(shù)學(xué)模型基礎(chǔ)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

19、0,迭代式xk+1=φ(xk)收斂于x*。,2.2 迭代法求解非線性方程,,而右圖曲線y=φ(x)位于直線y=x的上方，任給定初值x0,迭代式xk+1=φ(xk)不收斂。,2.2 迭代法求解非線性方程,,松弛迭代法當(dāng)?shù)绞諗亢苈龝r(shí)，我們可以從x=φ(x)出發(fā)構(gòu)造新的迭代形式，以加快收斂速度，這里先介紹松弛迭代法。迭代公式 其中， 為第k次迭代的預(yù)報(bào)值；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

44、Runge-Kutta法示例具體步驟：1）在F5單元格輸入初始濃度0.2000mol/L。2）根據(jù)T1-T4的計(jì)算公式， 在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 求微分方程數(shù)值解,5. Runge-Kutt

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