偏微分方程數(shù)值解2

拋物型方程的差分解法

拋物型方程是指如下形式的方程：

很多實(shí)際的物理問題都可以用這類方程描述：熱傳導(dǎo)方程：

現(xiàn)以熱傳導(dǎo)方程為例，介紹拋物型方程的有限差分格式。設(shè)熱傳導(dǎo)方程：定解條件(10.3.1)(10.3.2)求（10.3.1）滿足（10.3.2）的解。10.3.1矩形網(wǎng)格用兩組平行直線族xj=jh,

tk=k

(j=0,

1,…，k=0，

1,…)構(gòu)成的矩形網(wǎng)覆蓋了xt平面，網(wǎng)格點(diǎn)（xj,tk）稱為結(jié)點(diǎn)，簡記為（j,k），h、

為常數(shù)，分別稱為空間步長及時(shí)間步長，或稱h為沿x方向的步長，稱

為沿t方向的步長，，N為正整數(shù)。在t=0上的結(jié)點(diǎn)稱為邊界結(jié)點(diǎn)，其余所有屬于

內(nèi)的結(jié)點(diǎn)稱為內(nèi)部結(jié)點(diǎn)。txoh

(xj,tk)10.3.2．古典差分格式于平面區(qū)域上考慮傳導(dǎo)方程：a為正常數(shù)

（10.3.3）

(10.3.4)于結(jié)點(diǎn)（j,k）處偏導(dǎo)數(shù)與差商之間有如下近似的關(guān)系：利用上述表達(dá)式得到LU在(j,k)處的關(guān)系式：

(10.3.5)視為u(xj,tk)的近似值。

令，j=1,2,…,N–1；k=0,1,2,…則有：（10.3.6）差分方程（10.3.6）稱為解熱傳導(dǎo)方程（10.3.3）的古典顯格式，它所用到的結(jié)點(diǎn)如下圖：

*

***

（j,k）將（10.3.6）寫成便于計(jì)算的格式：（10.3.7）稱為網(wǎng)比，利用（10.3.7）及初邊值條件（4）在網(wǎng)格上的值（10.3.8）即可算出k=1,2,…，各層上的值。截?cái)嗾`差階為0(

+h2)。

為了提高截?cái)嗾`差的階，可以利用中心差商：j=1,2,…,N–1；k=0,1,2,…（10.3.9）得到

Richardson格式，其結(jié)點(diǎn)圖為：

*

***

（j,k）*截?cái)嗾`差階為o(

2+h2)，較古典顯格式高。將（10.3.9）式改寫成適于計(jì)算的形式：j=1,2,…,N–1；

k=1,2,…r=a

/h2稱為網(wǎng)比，（10.3.10）式中出現(xiàn)了三層網(wǎng)格上的值，（10.3.10）才能逐層計(jì)算。故需要事先求得第k－1層的值

和第k層的值,如果利用向后差商

j=1,2,…,N–1；k=0,1,2,…(10.3.11)(10.3.12)j=1,2,…,N–1；k=0,1,2,…古典隱格式，其結(jié)點(diǎn)圖為：

（j,k）****截?cái)嗾`差為o(

+h2)，與古典顯格式相同。

10.3.3．六點(diǎn)對稱格式取該點(diǎn)的中心差商，從而對于方程（10.3.3）式，在點(diǎn)列方程，，將以上各式代入（10.3.3）式得到差分方程：

整理,得

此即六點(diǎn)對稱格式，也稱為Crank-Nicolson格式，所用結(jié)點(diǎn)圖為：

***k+1 ***k

j+1jj–1（10.3.13）10.3.4．穩(wěn)定性（1）當(dāng)步長無限縮小時(shí)，差分方程的解是否逼近于微分方程（2）計(jì)算過程中產(chǎn)生的誤差在以后的計(jì)算中是無限增加，還是可以控制？（穩(wěn)定性）的解？（收斂性）穩(wěn)定性問題是研究拋物型差分方程的一個(gè)中心課題！

考察Richardson格式的穩(wěn)定性。

用表示計(jì)算所產(chǎn)生的誤差，如果右端無誤差存在，則滿足：取（10.3.14）假設(shè)k-1層之前無誤差存在。即，而在第k層產(chǎn)生了誤差。，這一層其它點(diǎn)也無誤差，而且在計(jì)算過程中不再產(chǎn)生新的誤差，利用（10.3.14）式算出誤差

的傳播如下表：

r=?時(shí)Richardson格式的誤差傳播

j

j0–4j0–3j0–2j0–1j0

j0+1j0+2j0+3j0+4k

-2

-4

7

4

-6

17

-24

17

-6

-8

31

-68

89

-68

31

-8

-10

49

-144

277

-388

277

-144

49

-10

71

-260

641

-109

1311

-109

641

-260

71

r≤1/2時(shí)古典顯格式的誤差傳播

j

j0–4j0–3j0–2j0–1j0

j0+1j0+2j0+3j0+4k

0.500.50.2500.5

00.25

0.125

00.375

00.375

00.125

0.0625

00.25

00.375

00.25

00.0625

如果選用

r=?

時(shí)的古典顯格式，誤差方程為：

差分格式關(guān)于初值穩(wěn)定的實(shí)際含義是：如果其解在某一層存在誤差，則由它引起的以后各層上的誤差不超過原始誤差的M倍（M為與

無關(guān)的常數(shù)）。因此，在穩(wěn)定的條件下，只要初始誤差足夠小，以后各層的誤差也能足夠小。以上構(gòu)造的幾種差分格式中，古典顯格式：r≤1/2時(shí)穩(wěn)定古典隱格式：絕對穩(wěn)定Richardson格式：絕對不穩(wěn)定六點(diǎn)對稱格式：絕對穩(wěn)定。穩(wěn)定性概念：初邊值問題:有限差分法求解的拋物型方程PROJECT求t=0.1時(shí)刻的u值（解析解為：）分別用向前差分顯格式、隱格式、Richardson格式、和六點(diǎn)對稱格式并做各種方法的誤差分析，你能得出什么結(jié)論？（1）水流為穩(wěn)態(tài)和無其反應(yīng)情況下的溶質(zhì)運(yùn)移方程考慮六點(diǎn)對稱格式，在點(diǎn)列方程，令，，有：舍去o(h2)，o(

2)得到六點(diǎn)對稱格式：

令上式變?yōu)榱钌鲜娇梢詫懗桑?/p>

j=1,2,…,N–1對時(shí)間變量用向后差分，對空間變量用中心差分，可得到隱格式：令整理得：

（2）．非穩(wěn)態(tài)方程非穩(wěn)態(tài)一維垂直流情況下，土壤溶質(zhì)的基本方程為：

式中容積含水量

的求法如下：式中，h為負(fù)壓水頭；t為時(shí)間；z為到原點(diǎn)的距離（cm），先解非飽和垂直水流方程向下為正；C(h)為容水度，，K(h)為土壤導(dǎo)水率，可由一些常用的經(jīng)驗(yàn)公式算出。求出h后，用水分特征曲線換算成相應(yīng)的

值。v=q/

，q為通量。

式中：將以上方程整理后可寫成：

2、有限元法設(shè)有微分方程

定義在由邊界

圍成的區(qū)域以上，L為微分算子。

設(shè){

j}（j=1,2,…，n,…），是一完備的函數(shù)系。伽遼金方法是求形如的近似解，其中aj(j=1,2,…,n)為待定常數(shù)。un稱為試探函數(shù)，

j稱為形狀函數(shù)（或基函數(shù)，插值函數(shù)，為{

j}(j=1,2,…,n,…)中前n個(gè)線性無關(guān)的函數(shù)）。（1）(2)若u是方程（1）的精確解，則必有在Lu和f是連續(xù)函數(shù)的條件下，就等價(jià)于但在（2）中，只有n個(gè)待定常數(shù)，所以只需n

正交條件即可。這是一個(gè)關(guān)于a1，a2，…，an的線性方程組，即為所求的近似解。

稱為伽遼金方程組，解之，得到一組，于是，對應(yīng)的函數(shù)

一維溶質(zhì)運(yùn)稱模型的有限元法(1)．考慮水流為穩(wěn)態(tài)和無其它反應(yīng)情況下的溶質(zhì)運(yùn)移方程對區(qū)間[0,L]進(jìn)行剖分其結(jié)點(diǎn)為x0,…,xN。令0=a=x0，L=b=xN，伽遼金方法即是求形如。

的解，使其系數(shù)滿足方程式中N為結(jié)點(diǎn)總數(shù)；Cj為結(jié)點(diǎn)j在t時(shí)刻的濃度；

j為線性插值基函數(shù)，其表達(dá)式為：

從

i的表達(dá)式中可以看出,

j

僅在[xj-1,xj+1]上不為零，于是當(dāng)j

0，j

N時(shí)，有積分（3）式，由分部積公式：

(3)(4)分別計(jì)算（4）式的各項(xiàng)積分,得：

（5）對（5）中的項(xiàng)進(jìn)行離散：

（6）再將（5）式代入（6）式，并取時(shí)刻的時(shí)間水平濃度，得：整理得：令

則有

（2）．非穩(wěn)態(tài)方程非穩(wěn)態(tài)一維垂直水流情況下，土壤溶質(zhì)運(yùn)移的基本方程為：對區(qū)間[0,L]進(jìn)行剖分，其結(jié)節(jié)為