拋物型方程的差分方法

微分方程數(shù)值解法第2章拋物型方程的差分方法

2.1差分格式建立的基礎(chǔ)

2.2顯式差分格式

2.3隱式差分格式

2.4解三對(duì)角形方程的追趕法

2.5差分格式的穩(wěn)定性和收斂性

2.6非線性拋物型方程的差分解法舉例

2.7二維拋物型方程的差分格式

2.8交替方向的隱式差分格式(ADI

格式)

本章，我們研究線性拋物型方程的差分解法，主要討論差分方程的構(gòu)造方法和有關(guān)的理論問(wèn)題以及研究方法等，重點(diǎn)在于一維線性拋物型方程的差分方法，對(duì)于非線性以及多維拋物型方程的差分解法也進(jìn)行了研究。其中,為平面上某一區(qū)域。(2.1)眾所周知，一維線性拋物型方程的一般形式為

(2)初邊值問(wèn)題(或稱混合問(wèn)題)通常考慮的定解問(wèn)題有：(1)初值問(wèn)題(或稱Cauchy問(wèn)題)在區(qū)域上求函數(shù)，使?jié)M足(2.2)為給定的初始函數(shù)。(2.3)(2.4)

在區(qū)域上求函數(shù)，使?jié)M足邊值條件初值條件為了構(gòu)造微分方程(2.1)的有限差分逼近，首先將求解區(qū)域用二組平行于軸和軸的直線構(gòu)成的網(wǎng)格覆蓋，網(wǎng)格邊長(zhǎng)在方向?yàn)椋诜较驗(yàn)?如圖2.1所示)。分別稱為空間方向和時(shí)間方向的步長(zhǎng)，網(wǎng)格線的交點(diǎn)稱為網(wǎng)格的結(jié)點(diǎn)。對(duì)初值問(wèn)題來(lái)說(shuō)，網(wǎng)格是2.1差分格式建立的基礎(chǔ)在上的結(jié)點(diǎn)稱為邊界結(jié)點(diǎn)，屬于內(nèi)的結(jié)點(diǎn)

稱為內(nèi)部結(jié)點(diǎn)。對(duì)于初邊值問(wèn)題，設(shè)，則網(wǎng)格是研究導(dǎo)數(shù)的差商近似表達(dá)式。為此對(duì)二元函數(shù)定義，且假定具有我們需要的有界偏導(dǎo)數(shù)。在上的結(jié)點(diǎn)稱為邊界結(jié)點(diǎn)，屬于內(nèi)的結(jié)點(diǎn)稱為內(nèi)部結(jié)點(diǎn)。差分方程就是在網(wǎng)格點(diǎn)上求出微分方程解的近似值的一種方法，因此又稱為網(wǎng)格法。構(gòu)造逼近微分方程的差分方程的方法。由Taylor展開(kāi)，有

則在處對(duì)的一階偏導(dǎo)數(shù)有三個(gè)可能的近似:(2.5)(2.6)(2.7)向前差商向后差商中心差商

顯然，用差商近似導(dǎo)數(shù)存在誤差，令(2.8)則

關(guān)于導(dǎo)數(shù)的近似差商表達(dá)式，也可以通過(guò)線性算子作為推導(dǎo)工具得到，定義：截?cái)嗾`差,階為用向后差商近似導(dǎo)數(shù)的截?cái)嗾`差階也為而中心差商近似導(dǎo)數(shù)的截?cái)嗾`差階為為方向偏導(dǎo)數(shù)算子為方向位移算子，為方向平均算子，其中:

方向的差分算子:(2.9)前差算子:,(2.10)后差算子:,中心差算子:(2.11),,

建立差分算子和導(dǎo)數(shù)算子之間的關(guān)系，由Talyor

展開(kāi)，有由得(2.12)或者(2.13)同理有因?yàn)楣?2.14)同理(2.15)因?yàn)?/p>

(2.16)則(2.17)

式(2.14)，(2.15)，(2.17)分別給出了偏導(dǎo)數(shù)算子關(guān)于前差、后差、中心差的級(jí)數(shù)表達(dá)式雙曲正弦3246(2.18.1)(2.18.2)(2.18.3)

利用這些關(guān)系式就可給出偏導(dǎo)數(shù)的差分表達(dá)式返回又由可得二階偏導(dǎo)數(shù)的差分表達(dá)式(2.19.1)(2.19.2)(2.19.3)返回返回4235(2.20.1)(2.20.2)(2.20.3)(2.21.1)(2.21.2)(2.21.3)對(duì)于三階、四階偏導(dǎo)數(shù)的差分表達(dá)式為

從以上這些偏導(dǎo)數(shù)的差分表達(dá)式，我們可以得到偏導(dǎo)數(shù)的各種精度的近似表達(dá)式。且

又由二階導(dǎo)數(shù)的前差表達(dá)式(2.19.1)，得因此在的前差表達(dá)式中取第一項(xiàng)，則有即截?cái)嗾`差階為?，F(xiàn)在研究構(gòu)造微分方程(2.1)的差分方程的方法，為此記微分方程(2.1)為(2.22)

L

是關(guān)于的線性算子，。包括二個(gè)相鄰時(shí)間層的網(wǎng)格結(jié)點(diǎn)的差分方程可以從Talor

展開(kāi)式推出返回設(shè)，于是(2.23)如果算子L不依賴于t，即，則(2.25)將式(2.17)，,代入算子L中，即在L中用中心差分算子代替了微分算子，于是有(2.24)返回3835目前通常用于解方程(2.1)的各種差分方程，都是方程(2.25)的近似表達(dá)式。下面各節(jié)，我們將以式(2.25)為基礎(chǔ)，對(duì)簡(jiǎn)單的拋物型方程，推導(dǎo)一些常用差分格式。對(duì)于用差分方法求偏導(dǎo)數(shù)方程的數(shù)值解來(lái)說(shuō)，設(shè)計(jì)差分方程，用之作為微分方程的近似，僅僅是第一步。本章除致力于這一研究外，特別著重討論了諸如差分格式的穩(wěn)定性、收斂性等基本問(wèn)題，它們也是本書(shū)研究的主要內(nèi)容之一。2.2顯式差分格式

現(xiàn)在，對(duì)拋物型方程(2.1)的幾種特殊情況，從方程(2.25)出發(fā)，構(gòu)造微分方程的有限差分近似。2.2.1一維常系數(shù)熱傳導(dǎo)方程的古典顯示格式

首先考慮一維熱傳導(dǎo)方程(2.26)的差分近似。差分方程的構(gòu)造由，方程(2.24)為代入式(2.19.3)，得

算子之間的關(guān)系則(2.27)其中為步長(zhǎng)比。返回在上式中，如果僅僅保留二階中心差分，且設(shè)為相應(yīng)差分方程解在結(jié)點(diǎn)(mh,nk)上的值，則(2.28)代入的表達(dá)式，則得差分方程(2.29)將格式(2.29)應(yīng)用于解初值問(wèn)題(初邊值問(wèn)題)古典顯式差分格式圖2.2差分格式(2.29)也可簡(jiǎn)單地由導(dǎo)數(shù)的差商近似表達(dá)式得到代入微分方程(2.26)，并令差分方程解為即可。雖然在邊界結(jié)點(diǎn)上，差分方程和微分方程具有相同的初值或者初邊值條件，但是，一般而言，結(jié)點(diǎn)上微分方程的精確解和古典顯式差分格式(2.29)的精確解不相等。(2.30)記假定具有下面推導(dǎo)中所需要的有界偏導(dǎo)數(shù)，則由展開(kāi)，有截?cái)嗾`差42(2.31)則由式(2.26)，(2.29)，(2.30)，(2.31)得(2.32)從式(2.31)有或(2.33)從而，上式右邊量描寫(xiě)了古典顯式差分格式(2.29)在點(diǎn)對(duì)微分方程的近似程度，將其定義為差分格式在點(diǎn)的截?cái)嗾`差，記為，即(2.34)假定在所考慮的區(qū)域保持有界，則古典顯式差分格式的截?cái)嗾`差階為。從式(2.33)又可見(jiàn)到，如令，因?yàn)楣式財(cái)嗾`差的階可以提高，這時(shí)。(2.35.1)或者(2.35.2)相應(yīng)的截?cái)嗾`差階為。通常，格式可用圖2.3表示。為了提高截?cái)嗾`差的階，我們也可用在式(2.27)中保留四階中心差分項(xiàng)的辦法達(dá)到，這時(shí)有差分格式(2.27)m,n+1m-2,nm-1,nm,nm+1,nm+2,n圖2.3m,n+1m-1,nm,nm+1,n圖2.2返回2.2.2系數(shù)依賴于的一維熱傳導(dǎo)方程的顯式格式(2.36)這時(shí)，。L保留右邊前二項(xiàng)，由，則有差分方程(2.37)則這一差分格式可用圖2.4表示，其中，這是一個(gè)顯式差分格式，其截?cái)嗾`差階為。m,n+1m-1,nm,nm+1,n圖2.4由方程右邊進(jìn)一步，考慮熱傳導(dǎo)方程(2.38)的差分近似。12在上式中保留前二項(xiàng)，并且和分別用和代替，則得差分方程(2.39)也可通過(guò)直接用中心差分算子代替微分算子的辦法獲得方程(2.38)的差分近似(2.40)這也是一個(gè)顯式差分格式。格式(2.39)和(2.40)的截?cái)嗾`差階都是。易見(jiàn)，由注：均在處計(jì)算。Delta顯然，微分方程(2.36)，(2.38)中的如果為，即其自變量包括空間變量和時(shí)間變量，這時(shí)差分格式(2.37)，(2.39)，(2.40)同樣是微分方程的具有截?cái)嗾`差階的差分近似，這時(shí)格式(2.37)，(2.39)中和,格式(2.40)中和分別換成

,。代入格式(2.40)即為格式(2.39)，差分格式(2.40)的推導(dǎo)方法,即在微分方程中直接用差分算子代替正如前面已經(jīng)指出的是推導(dǎo)差分格式的一個(gè)常用方法。2.3隱式差分格式

隱式差分格式特點(diǎn)：

1.具有二個(gè)或二個(gè)以上結(jié)點(diǎn)處的值未知；

2.計(jì)算工作量較大；

3.穩(wěn)定性較好。得由推導(dǎo)其最簡(jiǎn)單的隱式差分逼近─古典隱式格式?，F(xiàn)在對(duì)熱傳導(dǎo)方程2.3.1古典隱式格式1715格式用圖2.5表示，其截?cái)嗾`差階為，與古典顯式差分格式相同。或者(2.41)保留二階導(dǎo)數(shù)項(xiàng)，且以替代，則得差分格式我們也可通過(guò)直接用差分算子代替的方法，即代入微分方程，得到格式(2.41)。古典隱式差分格式m,n+1m-1,n+1m+1,n+1m,n圖2.5隱式差分格式是解熱傳導(dǎo)方程(2.26)的常用的差分格式，由式(2.24)，有2.3.2隱式格式由得(2.42)42兩邊僅保留二項(xiàng)，用代替，則得差分格式(2.43)這是一個(gè)隱式差分格式，稱為差分格式，截?cái)嗾`差階為。(2.44)由于格式(2.44)中包括六個(gè)結(jié)點(diǎn)，故也稱為六點(diǎn)格式(如圖2.6所示)。m,n+1m-1,n+1m+1,n+1m,n圖2.6m-1,nm+1,n44也可將代入微分方程(2.26)，得到格式。由式(2.19.3)，可令則可得另一精度較高的六點(diǎn)差分格式，如前在式(2.42)中僅保留直到的項(xiàng)，即有13代入上式，則有如下差分格式：(2.45)稱為差分格式。38截?cái)嗾`差階

2323因?yàn)?8前面，我們已經(jīng)推導(dǎo)了熱傳導(dǎo)方程(2.26)的古典顯式格式。古典隱式格式及格式等。實(shí)際上，它們都可以作為本節(jié)推導(dǎo)的加權(quán)六點(diǎn)隱式格式的特殊情形。2.3.3加權(quán)六點(diǎn)隱式格式由得到即用代替，則得差分格式或者(2.46)這是一個(gè)六點(diǎn)差分格式(如圖2.7所示)，稱為加權(quán)六點(diǎn)差分格式。40時(shí),為古典顯式格式；時(shí),為古典隱式格式；時(shí),為格式；加權(quán)六點(diǎn)格式亦可直接由差商代替導(dǎo)數(shù)得到2.3.4系數(shù)依賴于的一維熱傳導(dǎo)方程的一個(gè)隱式格式的推導(dǎo)

由其展開(kāi)式可得(2.47)的差分逼近?？紤]方程已知11令代入式(2.48)，則(2.48)因此43格式(2.49.1)具有截?cái)嗾`差階。這是一個(gè)隱式差分格式(如圖2.8所示)。(2.49.1)因此得差分方程(2.49.2)可寫(xiě)成形式m,n+1m-1,n+1m+1,n+1m,n圖2.8m-1,nm+1,n前節(jié)引進(jìn)的隱式差分方程，在要求解未知函數(shù)值的時(shí)間層上包括三個(gè)未知函數(shù)值。因此，這些隱式差分格式僅僅適合于解如圖中所示的邊值問(wèn)題。在每一時(shí)間層，需要求解的隱式差分方程形成了一個(gè)線性代數(shù)方程組，它的系數(shù)矩陣是三對(duì)角形矩陣，即僅在主對(duì)角線及其相鄰二條對(duì)角線上有非零元素。方程組寫(xiě)成一般形式是2.4解三對(duì)角形方程的追趕法m,n+1m-1,n+1m+1,n+1m,nm,n+1m-1,n+1m+1,n+1m,nm-1,nm+1,n(2.50)這一類方程可用追趕法求解。由方程組(2.50)中的第一個(gè)方程解出，得將此式代入方程組(2.50)的第二個(gè)方程，得到即令，則上式可寫(xiě)為其中完全類似地，可以推出下面的公式(2.51)其中注意當(dāng)時(shí),。即將關(guān)系式代入式(2.50)中最后一個(gè)方程，得到若令則有。如果已經(jīng)算出，那么解向量的最后一個(gè)分量就已求得，為了求得的所有分量，只有利用方程(2.51)即可逐步求出,因此，整個(gè)求解過(guò)程分為兩大步：

,第一步依次確定計(jì)算公式可歸結(jié)為第二步依相反次序確定通常，第1步稱為“追”的過(guò)程，第2步稱為“趕”的過(guò)程，整個(gè)求解過(guò)程稱為追趕法。(2)(3)則上述追趕法過(guò)程是穩(wěn)定的。(1)可以論證，如果例2.2

說(shuō)明用方法數(shù)值解如下定解問(wèn)題的過(guò)程：由前已知格式為如果選擇，則，要解的方程組寫(xiě)成矩陣形式是(2.52)相應(yīng)于上述定解問(wèn)題的差分方程組為其中，為七階方陣，為列向量，它們的表達(dá)式從式(2.52)可知。因?yàn)樵谇蟮趯訒r(shí)，已計(jì)算得，(它們?cè)谥谐霈F(xiàn))由邊值條件已知，故方程組右邊已知，且又因此可用追趕法求解方程組(2.52)，由方程組右邊值及可求出，然后順次，可求出。我們先看一個(gè)數(shù)值例子，考慮初邊值問(wèn)題(2.53)其中2.5差分格式的穩(wěn)定性和收斂性2.5.1問(wèn)題的提出利用顯式差分格式(2.29)，即式中。連同初值條件邊值條件逐層解出結(jié)點(diǎn)處的值?，F(xiàn)在對(duì)，取二種，使與。圖2.9和圖2.10中的曲線表示不同時(shí)刻微分方程的精確解，圖中“·”表示差分方程的解。圖2.9所示時(shí)的計(jì)算結(jié)果是曲線自上而下依次為微分方程的精確解。黑點(diǎn)是用差分格式在時(shí)算出的相應(yīng)各層上的近似值。二者符合得很好，由于對(duì)稱性我們只給出一半圖形。圖2.10是當(dāng)時(shí)差分方程解和微分方程精確解的圖示，黑點(diǎn)仍表示差分方程解，其中分別為在時(shí)的計(jì)算結(jié)果。從圖中看出，隨著的增大，差分方程的解越來(lái)越遠(yuǎn)離微分方程的解。由此可見(jiàn)，值的不同，得出的結(jié)果有很大的差別，如的結(jié)果是可用的，但是時(shí)的結(jié)果就完全沒(méi)有用。當(dāng)然上面各種情況所得的差分方程解是由計(jì)算機(jī)得到的，不可能是差分方程理論上的準(zhǔn)確解，而是差分方程的近似解，我們用表示。顯然與之間存在著差別，差分方程的準(zhǔn)確解與微分方程的解之間，如前所述,也是有差別的。因而從計(jì)算機(jī)上解得的差分方程近似解與微分方程解之間的差別實(shí)質(zhì)上包括兩方面的差別，即(2.54)下面我們先研究上式右邊第二項(xiàng)，即差分方程的理論解與計(jì)算機(jī)上解得的近似解之間的差別是隨著的增大而無(wú)限增加還是有所控制。如果這種差別是無(wú)限增加，則稱差分格式不穩(wěn)定，顯然不穩(wěn)定的格式是不能使用的,因?yàn)檎`差的無(wú)限增加淹沒(méi)了真解。上例中時(shí)就是差分方程不穩(wěn)定的情況。從差分方程比如格式(2.29)可知，在求

第一層的差分方程解時(shí)，用到第0層上的值，也就是初始值。由于計(jì)算機(jī)存儲(chǔ)數(shù)據(jù)為二進(jìn)制數(shù)位的限制,不可能完全精確地存儲(chǔ)在機(jī)器中，也就是計(jì)算用到的是帶有誤差的初始值。一般來(lái)說(shuō)，在計(jì)算時(shí)又出現(xiàn)了誤差，因此中包括了由于參加運(yùn)算而出現(xiàn)的誤差，即初始誤差的傳遞，以及本身計(jì)算過(guò)程中出現(xiàn)的誤差。這樣，在第時(shí)間層計(jì)算時(shí)得到的是由于前面的誤差傳遞和本身計(jì)算中出現(xiàn)的誤差引起的。下面我們給出研究差分格式穩(wěn)定性的最直接的方法，就是在第0層的一個(gè)結(jié)點(diǎn)上給出一個(gè)誤差，然后研究這個(gè)誤差的發(fā)展情況，即圖方法。假定在固定的某個(gè)結(jié)點(diǎn)引入一個(gè)誤差，即把改成了，而在這一層的其他結(jié)點(diǎn)上的初值還是，假定用帶有初始誤差的初值按差分格式去計(jì)算以后各排結(jié)點(diǎn)上的值，且假定計(jì)算時(shí)沒(méi)有引入其他誤差，我們把得到的值記做，這樣滿足原來(lái)的差分格式。假如我們使用差分格式(2.29)，于是2.5.2圖方法顯然兩解之差滿足(2.55)(2.56)以下分析當(dāng)和時(shí)，隨著增加而變化的情況。先看的情況，由式(2.55)得由此利用條件(2.56)即可算出的值(見(jiàn)表2.3)。表2.3由表2.3可知，用顯式差分格式(2.29)()計(jì)算時(shí)，由初始數(shù)據(jù)的誤差，在以后各層所引起的誤差是逐層減小的，這說(shuō)明差分格式(2.29)當(dāng)時(shí)是穩(wěn)定的。再看的情形，由(2.55)得由此利用條件(2.56)即可得出的值(見(jiàn)表2.4)。表2.4可見(jiàn)，用顯式差分格式(2.29)()計(jì)算時(shí)，由初始數(shù)據(jù)的誤差所引起的誤差在以后各層的計(jì)算中逐層迅速增大，以致不能控制，因此差分格式(2.29)在時(shí)是不穩(wěn)定的。用圖方法討論格式的穩(wěn)定性能直觀地看到差分格式是穩(wěn)定性，缺點(diǎn)是必先固定。2.5.3穩(wěn)定性定義、穩(wěn)定性分析的矩陣方法

以下討論求初邊值問(wèn)題差分方程寫(xiě)成矩陣形式為：(2.57)其中為維列向量；；為已知向量；為包括邊值條件的向量；為階方陣，可以隨而改變。如果差分方程為顯式，則對(duì)所有的(2.58)如果，則隱式格式可寫(xiě)成顯式形式設(shè)是初始值引進(jìn)的誤差向量，而在邊值以及其他各層計(jì)算中未引入其它任何誤差。由于的引入，差分方程的解為。則我們說(shuō)差分格式是穩(wěn)定的，其中是某一向量范數(shù)。穩(wěn)定性的定義：時(shí)，對(duì)于任何的，差分格式得到的解滿足不等式對(duì)于任意給定的，存在與無(wú)關(guān)且依賴于的正數(shù)，使當(dāng)設(shè)向量，則常用的向量范數(shù)有：(1)(2)(3)它們分別稱為2-范數(shù)，1-范數(shù)和無(wú)窮范數(shù)，其中2-范數(shù)亦稱為歐氏范數(shù)。

(1)，其中，為的共軛轉(zhuǎn)置矩陣，為的最大特征值；

(