第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)隙為蒸汽，其溫度在任何位置均為恒定值，可認(rèn)為是飽和蒸汽的溫度。2、忽略套管內(nèi)側(cè)流體的縱向熱傳導(dǎo)。3、在整個(gè)套管長度方向上，總傳熱系數(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)用另一個(gè)在化工中常見的微分方程是物料冷卻過程的數(shù)學(xué)模型，其模型可用下式表示：

（4-3）在微分方程中我們稱自變量函數(shù)只有一個(gè)的微分方程為常微分方程，自變量函數(shù)個(gè)數(shù)為兩個(gè)或兩個(gè)以上的微分方程為偏微分方程。給定微分方程及其初始條件，稱為初值問題；給定微分方程及其邊界條件，稱為邊值問題。總目錄本章目錄4.14.24.34.44.54.1微分方程在化工中的應(yīng)用在化工模擬中主要碰到的是常微分方程的初值問題：（4-5）或?qū)τ诖蠖鄶?shù)常微分方程的初值問題,只能計(jì)算它的數(shù)值解。常微分方程初值問題的數(shù)值解就是求y(x)在求解區(qū)間[a,b]上各個(gè)分點(diǎn)序列xn，n=1,2,…,m的數(shù)值解yn。在計(jì)算中約定y(xn)表示常微分方程準(zhǔn)確解的值，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）下式為計(jì)算近似值的向前歐拉公式：圖4-2歐拉折線法幾何示意圖總目錄本章目錄4.14.24.34.44.54.2.1向前歐拉公式實(shí)例例4.1：假定某物體的溫度w因自熱而產(chǎn)生的熱量可以使物體在每秒鐘內(nèi)以4%的速度增長，同時(shí)該物體由于散熱可使其溫度在每秒種內(nèi)下降100k，則物體溫度隨時(shí)間變化的微分方程：（t以秒為單位)

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

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

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

h充分小時(shí)，可保證，其中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計(jì)算公式：類似地，可得到計(jì)算y(xn+1)近似值yn+1的計(jì)算公式：

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

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

計(jì)算結(jié)果如表4-2所示。表4-2計(jì)算結(jié)果總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法本課程主要研究其實(shí)際應(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龍格-庫塔方法實(shí)例例4.3：用四階龍格—庫塔公式(4-16)求解下面初值問題解：取步長h=0.2，計(jì)算公式為：表4-3計(jì)算結(jié)果總目錄本章目錄4.14.24.34.44.54.3龍格-庫塔方法步長的選擇下面以四階龍格-庫塔方法為例，說明如何自動(dòng)選擇步長，使計(jì)算結(jié)果滿足給定精度的要求。設(shè)從節(jié)點(diǎn)xn出發(fā)，先以h為步長，利用四階龍格-庫塔公式方法經(jīng)過一步計(jì)算得y(xn+1)的近似值，記為，由于公式的局部截?cái)嗾`差是y(h5)，故有當(dāng)h不大時(shí)，c可近似地看作常數(shù)。然后將步長h對(duì)折，即取h/2為步長，從出發(fā)經(jīng)過兩步計(jì)算求y(xn+1)的近似值，記為，每一步計(jì)算的局部截?cái)嗾`差為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)的近似值，其誤差可用先后兩次計(jì)算結(jié)果之差來表示，因而，只需考察是否成立。若成立，則可將作為y(xn+1)的近似值；若不成立，則將步長再次對(duì)折進(jìn)行計(jì)算，直到不等式成立為止，并取最后的作為計(jì)算結(jié)果。以上方法就是計(jì)算過程中自動(dòng)選擇步長的方法，也稱為變步長方法?？偰夸洷菊履夸?.14.24.34.44.54.3龍格-庫塔方法步長的選擇用一個(gè)例子說明：由VB選用標(biāo)準(zhǔn)四階龍格-庫塔方法計(jì)算得

因此合理選擇步長既能保證精度又能減少計(jì)算量?？偰夸洷菊履夸?.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個(gè)一階方程組成的常微分方程初值問題：向量形式：其中：(4-19)(4-20)總目錄本章目錄4.14.24.34.44.54.4.1一階常微分方程組的數(shù)值解法下面以兩個(gè)方程組為例，給出相應(yīng)的計(jì)算公式。常微分方程組：歐拉公式：預(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ù)值解法

實(shí)例例4.4：兩種微生物，其數(shù)量分別是u=u(t)，v=v(t)，t的單位為分，其中一種微生物以吃另一種微生為生，兩種微生物的增長函數(shù)如下列常微分方程組所示，預(yù)測3分鐘后這一對(duì)微生物的數(shù)量?？偰夸洷菊履夸?.14.24.34.44.54.4.1一階常微分方程組的數(shù)值解法

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

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

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

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

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