第4章常微分方程數(shù)值解

第4章

常微分方程數(shù)值解4.1微分方程在化工中的應(yīng)用

4.2歐拉（Euler）公式

4.3龍格-庫塔方法

4.4常微分方程組的數(shù)值解法

4.5程序示例及應(yīng)用

總目錄4.1微分方程在化工中的應(yīng)用微分方程在化工中應(yīng)用的簡單而又典型的例子是套管式換熱器的穩(wěn)態(tài)溫度分布。首先作以下假設(shè)：1、套管內(nèi)側(cè)為液體，其溫度只隨套管的長度改變而改變，忽略溫度的徑向變化；套管環(huán)隙為蒸汽，其溫度在任何位置均為恒定值，可認為是飽和蒸汽的溫度。2、忽略套管內(nèi)側(cè)流體的縱向熱傳導(dǎo)。3、在整個套管長度方向上，總傳熱系數(shù)K不變?？偰夸洷菊履夸?.14.24.34.44.54.1微分方程在化工中的應(yīng)用蒸汽入口流體入口，u，t0冷凝液出口流體出口，u，tL圖4-1套管式換熱器溫度分布示意圖流入的熱量+傳入的熱量-流出的熱量=0（4-1）（4-2）（4-3）總目錄本章目錄4.14.24.34.44.54.1微分方程在化工中的應(yīng)用另一個在化工中常見的微分方程是物料冷卻過程的數(shù)學模型，其模型可用下式表示：

（4-3）在微分方程中我們稱自變量函數(shù)只有一個的微分方程為常微分方程，自變量函數(shù)個數(shù)為兩個或兩個以上的微分方程為偏微分方程。給定微分方程及其初始條件，稱為初值問題；給定微分方程及其邊界條件，稱為邊值問題。總目錄本章目錄4.14.24.34.44.54.1微分方程在化工中的應(yīng)用在化工模擬中主要碰到的是常微分方程的初值問題：（4-5）或?qū)τ诖蠖鄶?shù)常微分方程的初值問題,只能計算它的數(shù)值解。常微分方程初值問題的數(shù)值解就是求y(x)在求解區(qū)間[a,b]上各個分點序列xn，n=1,2,…,m的數(shù)值解yn。在計算中約定y(xn)表示常微分方程準確解的值，yn表示y(xn)的近似值?？偰夸洷菊履夸?.14.24.34.44.54.2歐拉（Euler）公式4.2.1向前歐拉公式4.2.2向后歐拉公式

4.2.3中心歐拉公式4.2.4梯形公式

總目錄本章目錄4.14.24.34.44.54.2.1向前歐拉公式（4-6）下式為計算近似值的向前歐拉公式：圖4-2歐拉折線法幾何示意圖總目錄本章目錄4.14.24.34.44.54.2.1向前歐拉公式實例例4.1：假定某物體的溫度w因自熱而產(chǎn)生的熱量可以使物體在每秒鐘內(nèi)以4%的速度增長，同時該物體由于散熱可使其溫度在每秒種內(nèi)下降100k，則物體溫度隨時間變化的微分方程：（t以秒為單位)

分別以初始溫x(0)=1500k，y(0)=2500k，z(0)=3500k用歐拉公式預(yù)測24秒后的物體溫度趨勢?？偰夸洷菊履夸?.14.24.34.44.54.2.1向前歐拉公式實例解：

w0分別以x0=1500，y0=2500，z0=3500代入。計算結(jié)果見表4-1。圖4-3三種初始值的溫度變化曲線

表4-1總目錄本章目錄4.14.24.34.44.54.2.1向前歐拉公式實例從表4-1可以看到當自熱引起物體溫度升高的速度小于散熱引起溫度下降的速度，物體的溫度隨時間而逐漸減少：當自熱引起物體溫度升高的速度與散熱引起溫度下降的速度平衡時，物體的溫度保持不變；當自熱引起物體溫度升高的速度大于散熱引起溫度下降的速度，物體的溫度隨時間而增長。在圖4-3中L1，L2，L3分別表示初始值3500，2500和1500的三條溫度變化趨勢曲線?？偰夸洷菊履夸?.14.24.34.44.5h充分小時，以上迭代收斂。記，則

h充分小時，可保證，其中L為李普希茲條件。4.2.2向后歐拉公式向后歐拉公式：（4-6）式(4-7)是yn+1的非線性方程，即隱式歐拉公式，用迭代法求得yn+1。初始值由向前歐拉公式提供。最簡單的迭代公式為：總目錄本章目錄4.14.24.34.44.54.2.3中心歐拉公式y(tǒng)(x)的在x=x1處的中心差商式：又，可得到y(tǒng)(x2)的近似值y2計算公式：類似地，可得到計算y(xn+1)近似值yn+1的計算公式：

公式(4-8)稱為中心格式。按公式(4-8)，需要知道yn-1,yn的值才能求得yn+1的值。因此，要先用其它公式計算出y1，再用中心格式算出y2,y3,…。y1可用向前歐拉公式計算，為提高精度，也可用向后歐拉公式計算。(4-8)總目錄本章目錄4.14.24.34.44.54.2.4梯形公式梯形公式：(4-9)梯形公式也是隱式格式。用顯式的歐拉公式和隱式的梯形公式給出的一次預(yù)估-校正公式：(4-10)上式也稱為改進的歐拉公式，它可合并成：總目錄本章目錄4.14.24.34.44.54.2.4梯形公式實例例4.2：請用預(yù)估-校正公式（改進的歐拉公式）解右面初值問題：解：用下面的迭代公式，對每個點迭代4次，k=1,2,3,4。

總目錄本章目錄4.14.24.34.44.54.2.4梯形公式實例該方程的精確解是

計算結(jié)果如表4-2所示。表4-2計算結(jié)果總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法本課程主要研究其實際應(yīng)用,，故直接給出各類龍格-庫塔公式。1、二階龍格-庫塔其中c1=0,c2=1,a=1/2,b=1/2。其中c1=1/2,c2=1/2,a=1,b=1或(4-11)(4-12)總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法2、三階龍格-庫塔公式

(4-15)(4-14)(4-13)總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法3、四階龍格—庫塔公式

(4-17)(4-16)總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法實例例4.3：用四階龍格—庫塔公式(4-16)求解下面初值問題解：取步長h=0.2，計算公式為：表4-3計算結(jié)果總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法步長的選擇下面以四階龍格-庫塔方法為例，說明如何自動選擇步長，使計算結(jié)果滿足給定精度的要求。設(shè)從節(jié)點xn出發(fā)，先以h為步長，利用四階龍格-庫塔公式方法經(jīng)過一步計算得y(xn+1)的近似值，記為，由于公式的局部截斷誤差是y(h5)，故有當h不大時，c可近似地看作常數(shù)。然后將步長h對折，即取h/2為步長，從出發(fā)經(jīng)過兩步計算求y(xn+1)的近似值，記為，每一步計算的局部截斷誤差為c(h/2)5，于是就有(4-18)總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法步長的選擇把它與（4-18）式相比，可得：經(jīng)整理可得：總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法步長的選擇這表明以作為y(xn+1)的近似值，其誤差可用先后兩次計算結(jié)果之差來表示，因而，只需考察是否成立。若成立，則可將作為y(xn+1)的近似值；若不成立，則將步長再次對折進行計算，直到不等式成立為止，并取最后的作為計算結(jié)果。以上方法就是計算過程中自動選擇步長的方法，也稱為變步長方法?？偰夸洷菊履夸?.14.24.34.44.54.3龍格-庫塔方法步長的選擇用一個例子說明：由VB選用標準四階龍格-庫塔方法計算得

因此合理選擇步長既能保證精度又能減少計算量?？偰夸洷菊履夸?.14.24.34.44.54.4常微分方程組的數(shù)值解法4.4.1一階常微分方程組的數(shù)值解法4.4.2高階常微分方程數(shù)值方法

總目錄本章目錄4.14.24.34.44.54.4.1一階常微分方程組的數(shù)值解法將由m個一階方程組成的常微分方程初值問題：向量形式：其中：(4-19)(4-20)總目錄本章目錄4.14.24.34.44.54.4.1一階常微分方程組的數(shù)值解法下面以兩個方程組為例，給出相應(yīng)的計算公式。常微分方程組：歐拉公式：預(yù)估—校正公式：

(4-21)(4-22)總目錄本章目錄4.14.24.34.44.54.4.1一階常微分方程組的數(shù)值解法四階龍格—庫塔公式：(4-23)總目錄本章目錄4.14.24.34.44.54.4.1一階常微分方程組的數(shù)值解法

實例例4.4：兩種微生物，其數(shù)量分別是u=u(t)，v=v(t)，t的單位為分，其中一種微生物以吃另一種微生為生，兩種微生物的增長函數(shù)如下列常微分方程組所示，預(yù)測3分鐘后這一對微生物的數(shù)量。總目錄本章目錄4.14.24.34.44.54.4.1一階常微分方程組的數(shù)值解法

實例解：記用歐拉預(yù)估—校正公式(4-22)表4-4計算結(jié)果總目錄本章目錄4.14.24.34.44.54.4.2高階常微分方程數(shù)值方法以三階常微分方程為例說明高階常微分方程的數(shù)值計算步驟。

(4-24)令得到一階方程組：總目錄本章目錄4.14.24.34.44.54.5程序示例及應(yīng)用任務(wù)1用改進的歐拉公式求解常微分方程初值問題：

算法描述對給定的F(x,y)，用改進的歐拉公式求解常微分方程初值問題的解。

(計算實例VB程序見課本)總目錄本章目錄4.14.24.34.44.54.5程序示例及應(yīng)用任務(wù)2：用四階龍格—庫塔方法求解常微分方程初值問題：

算法描述 對給定的F（x,y），用四階龍格—庫塔方法求解常微分方程初值問題。(計算實例VB程序見