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

第4章

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

4.1微分方程在化工中的應(yīng)用

4.2歐拉（Euler）公式

4.3龍格-庫(kù)塔方法

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

4.5程序示例

目錄4.1微分方程在化工中的應(yīng)用微分方程在化工中應(yīng)用的簡(jiǎn)單而又典型的例子是套管式換熱器的穩(wěn)態(tài)溫度分布。首先作以下假設(shè)：

1、套管內(nèi)側(cè)為液體，其溫度只隨套管的長(zhǎng)度改變而改變，忽略溫度的徑向變化；套管環(huán)隙為蒸汽，其溫度在任何位置均為恒定值，可認(rèn)為是飽和蒸汽的溫度。

2、忽略套管內(nèi)側(cè)流體的縱向熱傳導(dǎo)。

3、在整個(gè)套管長(zhǎng)度方向上，總傳熱系數(shù)K不變。據(jù)以上假設(shè)，可以得到圖中所示微元的能量平衡方程（示意圖見(jiàn)圖4-1）：

流入的熱量+傳入的熱量-流出的熱量=0（4-1）即：

蒸汽入口流體入口，u，t0冷凝液出口流體出口，u，tL4.14.44.34.2總目錄4.54.1微分方程在化工中的應(yīng)用化簡(jiǎn)上式得其溫度的微分方程：（4-3）微分方程中各變量的含義如下：通過(guò)求解微分方程（4-3），就可以得到管內(nèi)流體的溫度隨管子長(zhǎng)度而改變的曲線，為化工模擬和設(shè)計(jì)提供依據(jù)。如果方程（4-3）中傳熱系數(shù)、物流性質(zhì)不隨溫度或位置的變化，那么，方程（4-3）是可以解析求解的（數(shù)值求解當(dāng)然可以），得到我們常見(jiàn)的換熱器傳熱方程；如果上述性質(zhì)可能隨溫度或位置而改變，那么方程（4-3）就只能用數(shù)值的方法求解了。事實(shí)上傳熱系數(shù)也好，物流性質(zhì)也好都會(huì)隨著溫度的改變而改變，故在深入研究換熱器各點(diǎn)溫度分布時(shí)，擬采用微分方程的數(shù)值求解為好。

4.14.44.34.2總目錄4.54.1微分方程在化工中的應(yīng)用另一個(gè)在化工中常見(jiàn)的微分方程是物料冷卻過(guò)程的數(shù)學(xué)模型，其模型可用下式表示：

它含有自變量t（時(shí)間）、未知函數(shù)

T（隨時(shí)間變化的物料溫度）、T0（環(huán)境溫度）、k（降溫速率）以及溫度的一階導(dǎo)數(shù),是一個(gè)常微分方程。在微分方程中我們稱自變量函數(shù)只有一個(gè)的微分方程為常微分方程，自變量函數(shù)個(gè)數(shù)為兩個(gè)或兩個(gè)以上的微分方程為偏微分方程。給定微分方程及其初始條件，稱為初值問(wèn)題；給定微分方程及其邊界條件，稱為邊值問(wèn)題。在化工模擬中主要碰到的是常微分方程的初值問(wèn)題：或記為

只有一些特殊形式的,才能找到它的解析解；對(duì)于大多數(shù)常微分方程的初值問(wèn)題,只能計(jì)算它的數(shù)值解。常微分方程初值問(wèn)題的數(shù)值解就是求y(x)在求解區(qū)間[a,b]上各個(gè)分點(diǎn)序列xn，n=1,2,…,m的數(shù)值解yn。在計(jì)算中約定y(xn)表示常微分方程準(zhǔn)確解的值，yn表示y(xn)的近似值。下面我們將向大家介紹幾種常用的常微分方程數(shù)值求解法。

4.14.44.34.2總目錄4.54.2歐拉（Euler）公式

4.2.1向前歐拉公式

4.2.2向后歐拉公式

4.2.3中心歐拉公式

4.2.4梯形公式

4.14.44.34.2總目錄4.54.2.1向前歐拉公式

對(duì)于常微分方程初值問(wèn)題［式（4-5）］,在求解區(qū)間[a,b]上作等距分割，步長(zhǎng)，記xn

=xn-1+h,

n=1,2…,m。用差商近似導(dǎo)數(shù)計(jì)算常微分方程。

做的在x=x0處的一階向前差商得：又，于是得到：故y(x1)的近似值y1可按

求得。類似地，由

得到計(jì)算近似值的向前歐拉公式：

由yn直接算出yn+1值的計(jì)算格式稱為顯式格式，向前歐拉公式是顯式格式。而歐拉方法的幾何意義是以y1作為斜率，通過(guò)點(diǎn)(x0,y0)做一條直線，它與直線x=x1的交點(diǎn)就是y1。依此類推，是以yn+1作為斜率，經(jīng)過(guò)點(diǎn)的直線與直線x=xn+1的交點(diǎn)。故歐拉法也稱為歐拉折線法，如右圖所示。4.14.44.34.2總目錄4.54.2.1向前歐拉公式例4.1：假定某物體的溫度w因自熱而產(chǎn)生的熱量可以使物體在每秒鐘內(nèi)以4%的速度增長(zhǎng)，同時(shí)該物體由于散熱可使其溫度在每秒種內(nèi)下降100k，則物體溫度隨時(shí)間變化的微分方程：（t以秒為單位)

分別以初始溫度x(0)=1500k,y(0)=2500k,z(0)=3500k，用歐拉公式預(yù)測(cè)24秒后的物體溫度趨勢(shì)。解：

w0分別以x0=1500,y0=2500,z0=3500代入，計(jì)算結(jié)果見(jiàn)表4-1(下頁(yè))。

從表4-1可以看到當(dāng)自熱引起物體溫度升高的速度小于散熱引起溫度下降的速度，物體的溫度隨時(shí)間而逐漸減少：當(dāng)自熱引起物體溫度升高的速度與散熱引起溫度下降的速度平衡時(shí)，物體的溫度保持不變；當(dāng)自熱引起物體溫度升高的速度大于散熱引起溫度下降的速度，物體的溫度隨時(shí)間而增長(zhǎng)。在圖4-3（下頁(yè)）中L1,L2,L3分別表示初始值3500,2500和1500的三條溫度變化趨勢(shì)曲線。

4.14.44.34.2總目錄4.54.2.1向前歐拉公式圖4-3三種初始值的溫度變化曲線

表4-1

4.14.44.34.2總目錄4.54.2.2向后歐拉公式

做出的在x1處的一階向后差商式：

而，得到的近似值y1的計(jì)算公式：

類似地，可得到計(jì)算y(xn+1)近似值yn+1的計(jì)算公式：

公式(4-7)稱為向后歐拉公式。通常為非線性函數(shù)，因此式(4-7)是關(guān)于yn+1的非線性方程，稱為隱式歐拉公式，需要通過(guò)迭代法求得yn+1。其中初始值可由向前歐拉公式提供。最簡(jiǎn)單的迭代公式為：可以證明，h充分小時(shí)，以上迭代收斂。事實(shí)上，記，則

h充分小時(shí)，可以保證，其中L為李普希茲條件。

4.14.44.34.2總目錄4.54.2.3中心歐拉公式

做出y(x)的在x=x1處的中心差商式：又，可得到y(tǒng)(x2)的近似值y2的計(jì)算公式：類似地，可得到計(jì)算y(xn+1)近似值yn+1的計(jì)算公式：

(4-8)公式(4-8)稱為中心格式。按公式(4-8)，需要知道yn-1,yn的值才能求得yn+1的值。因此，要先用其它公式計(jì)算出y1，再用中心格式算出y2,y3,…。y1可用向前歐拉公式計(jì)算，為提高精度，也可用向后歐拉公式計(jì)算。

4.14.44.34.2總目錄4.54.2.4梯形公式

在兩點(diǎn)之間進(jìn)行梯形近似計(jì)算有：則得梯形公式：梯形公式也是隱式格式，計(jì)算中為了保證一定的精確度，又避免用迭代過(guò)程不菲的計(jì)算量，可先用顯式公式算出初始值，再用隱式公式進(jìn)行一次修正。稱為預(yù)估-校正過(guò)程。例如，下面是用顯式的歐拉公式和隱式的梯形公式給出的一次預(yù)估-校正公式：上式也稱為改進(jìn)的歐拉公式，它可合并成

如果想要獲得較高的計(jì)算精度，可進(jìn)行多次迭代計(jì)算，也就是進(jìn)行多次校正計(jì)算。下面的例子對(duì)每一個(gè)點(diǎn)進(jìn)行了4次迭代計(jì)算。

4.14.44.34.2總目錄4.54.2.4梯形公式例4.2：請(qǐng)用預(yù)估-校正公式（改進(jìn)的歐拉公式）解下面初值問(wèn)題：解：。用下面的迭代公式，對(duì)每個(gè)點(diǎn)迭代4次，k=1,2,3,4。

該方程的精確解是，計(jì)算結(jié)果如表4-2所示。

4.14.44.34.2總目錄4.5VB程序比較4.14.44.34.2總目錄4.54.3龍格-庫(kù)塔方法

龍格-庫(kù)塔法是求解常微分方程較常用的一種方法，它通過(guò)巧妙的線性組合，在顯式格式的情況下獲得理想的計(jì)算精度，大大提高了計(jì)算速度。該方法的推導(dǎo)過(guò)程不是本課程要研究的對(duì)象，本課程主要研究其實(shí)際應(yīng)用,，故直接給出各類龍格-庫(kù)塔公式。

1、二階龍格-庫(kù)塔其中c1=1/2,c2=1/2,a=1,b=1或

其中c1=0,c2=1,a=1/2,b=1/2。

二階龍格-庫(kù)塔公式的局部截?cái)嗾`差為O(h3),是二階精度的計(jì)算公式。也可建立高階的龍格-庫(kù)塔公式，如三階龍格-庫(kù)塔公式、四階龍格-庫(kù)塔公式。較常用的是四階龍格-庫(kù)塔公式，四階龍格-庫(kù)塔公式的局部截?cái)嗾`差界為O(h5),是四階精度的計(jì)算公式。

4.14.44.34.2總目錄4.54.3龍格-庫(kù)塔方法2、三階龍格-庫(kù)塔公式

4.14.44.34.2總目錄4.54.3龍格-庫(kù)塔方法3、四階龍格—庫(kù)塔公式

4.14.44.34.2總目錄4.54.3龍格-庫(kù)塔方法例4.3：用四階龍格—庫(kù)塔公式(4-16)求解下面初值問(wèn)題：解：取步長(zhǎng)h=0.2，計(jì)算公式為：

計(jì)算結(jié)果列表4-3中。

4.14.44.34.2總目錄4.54.3龍格-庫(kù)塔方法——步長(zhǎng)的選擇

怎樣選取合適的步長(zhǎng)，這在實(shí)際計(jì)算中是很重要的。由于步長(zhǎng)愈小，每步計(jì)算的截?cái)嗾`差就愈小；但在一定的求解范圍內(nèi)，需要完成的步數(shù)就愈多，這不但引起計(jì)算量的增加，而且計(jì)算步數(shù)的增加往往造成舍入誤差的嚴(yán)重積累，反而降低了計(jì)算精度。上面介紹的龍格-庫(kù)塔方法是對(duì)定步長(zhǎng)（即步長(zhǎng)h為常數(shù)）而言的，但為了保證精確度，一種有效的措施是在計(jì)算過(guò)程中自動(dòng)進(jìn)行步長(zhǎng)調(diào)整，即變步長(zhǎng)技巧。下面以四階龍格-庫(kù)塔方法為例，說(shuō)明如何自動(dòng)選擇步長(zhǎng)，使計(jì)算結(jié)果滿足給定精度的要求。設(shè)從節(jié)點(diǎn)xn出發(fā)，先以h為步長(zhǎng)，利用四階龍格-庫(kù)塔公式方法經(jīng)過(guò)一步計(jì)算得y(xn+1)的近似值，記為，由于公式的局部截?cái)嗾`差是y(h5)，故有

當(dāng)h不大時(shí)，c可近似地看作常數(shù)。然后將步長(zhǎng)h對(duì)折，即取h/2為步長(zhǎng)，從出發(fā)經(jīng)過(guò)兩步計(jì)算求y(xn+1)的近似值，記為，每一步計(jì)算的局部截?cái)嗾`差為c(h/2)5，于是就有

4.14.44.34.2總目錄4.54.3龍格-庫(kù)塔方法——步長(zhǎng)的選擇把它與（4-18）式相比，可得

經(jīng)整理可得

這表明以作為y(xn+1)的近似值，其誤差可用先后兩次計(jì)算結(jié)果之差來(lái)表示，因而，只需考察是否成立。若成立，則可將作為y(xn+1)的近似值；若不成立，則將步長(zhǎng)再次對(duì)折進(jìn)行計(jì)算，直到不等式成立為止，并取最后的作為計(jì)算結(jié)果。以上方法就是計(jì)算過(guò)程中自動(dòng)選擇步長(zhǎng)的方法，也稱為變步長(zhǎng)方法。從表面上看，為了選擇適當(dāng)?shù)牟介L(zhǎng)，每一步的計(jì)算量增加了，但從總體考慮，工作量往往還是減少的。龍格-庫(kù)塔方法的主要優(yōu)點(diǎn)是計(jì)算精確度較高，能滿足通常的計(jì)算要求；容易編制程序；每次計(jì)算y(xn+1)，只用到前一步的計(jì)算結(jié)果yn，因此，在已知初始值y0的條件下，就可自動(dòng)地進(jìn)行計(jì)算，是單步法，而且計(jì)算過(guò)程中可隨時(shí)改變步長(zhǎng)。它的缺點(diǎn)是每前進(jìn)一步需要計(jì)算多次f(x,y)的值，因此，計(jì)算工作量較大，且其截?cái)嗾`差難以估計(jì)。在實(shí)際應(yīng)用上，一般當(dāng)要求更高精確度時(shí)，采用的辦法是縮小步長(zhǎng)，而不是采用更高階的公式，因?yàn)楦唠A公式的計(jì)算太復(fù)雜，一般選用標(biāo)準(zhǔn)四階龍格-庫(kù)塔方法即可。

4.14.44.34.2總目錄4.54.3龍格-庫(kù)塔方法——步長(zhǎng)的選擇下面用一個(gè)例子說(shuō)明：由VB選用標(biāo)準(zhǔn)四階龍格-庫(kù)塔方法計(jì)算得由積分法算得y(2)=2.23607,當(dāng)h=0.5時(shí)相差0.00013，而h=0.25誤差小于0.000001，但當(dāng)h=0.125時(shí)雖然足夠精確但計(jì)算次數(shù)卻比h=0.25多了一倍。因此合理選擇步長(zhǎng)既能保證精度又能減少計(jì)算量。4.14.44.34.2總目錄4.54.4常微分方程組的數(shù)值解法

4.4.1一階常微分方程組的數(shù)值解法

4.4.2高階常微分方程數(shù)值方法

4.14.44.34.2總目錄4.54.4.1一階常微分方程組的數(shù)值解法將由m個(gè)一階方程組成的常微分方程初值問(wèn)題：

寫(xiě)成向量形式：

其中：

4.14.44.34.2總目錄4.54.4.1一階常微分方程組的數(shù)值解法前面對(duì)常微分方程所用的各種方法，都可以平行地應(yīng)用到常微分方程組的數(shù)值解中。下面以兩個(gè)方程組為例，給出相應(yīng)的計(jì)算公式。

常微分方程組：歐拉公式：預(yù)估—校正公式：

4.14.44.34.2總目錄4.54.4.1一階常微分方程組的數(shù)值解法四階龍格—庫(kù)塔公式：

4.14.44.34.2總目錄4.54.4.1一階常微分方程組的數(shù)值解法例4.4：兩種微生物，其數(shù)量分別是u=u(t)，v=v(t)，t的單位為分，其