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

拋物型方程的差分解法

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

很多實際的物理問題都可以用這類方程描述：熱傳導方程：

現(xiàn)以熱傳導方程為例，介紹拋物型方程的有限差分格式。設熱傳導方程：定解條件(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,…)構成的矩形網(wǎng)覆蓋了xt平面，網(wǎng)格點（xj,tk）稱為結點，簡記為（j,k），h、

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

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

內的結點稱為內部結點。txoh

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

（10.3.3）

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

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

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

*

***

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

+h2)。

為了提高截斷誤差的階，可以利用中心差商：j=1,2,…,N–1；k=0,1,2,…（10.3.9）得到

Richardson格式，其結點圖為：

*

***

（j,k）*截斷誤差階為o(

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

k=1,2,…r=a

/h2稱為網(wǎng)比，（10.3.10）式中出現(xiàn)了三層網(wǎng)格上的值，（10.3.10）才能逐層計算。故需要事先求得第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,…古典隱格式，其結點圖為：

（j,k）****截斷誤差為o(

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

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

整理,得

此即六點對稱格式，也稱為Crank-Nicolson格式，所用結點圖為：

***k+1 ***k

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

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

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

的傳播如下表：

r=?時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時古典顯格式的誤差傳播

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=?

時的古典顯格式，誤差方程為：

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

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

2)得到六點對稱格式：

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

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

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

式中容積含水量

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

值。v=q/

，q為通量。

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

2、有限元法設有微分方程

定義在由邊界

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

設{

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

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

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

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

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

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

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

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

從

i的表達式中可以看出,

j

僅在[xj-1,xj+1]上不為零，于是當j

0，j

N時，有積分（3）式，由分部積公式：

(3)(4)分別計算（4）式的各項積分,得：

（5）對（5）中的項進行離散：

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

則有

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