一維拋物方程界面問題的緊致有限體積格式

于 倩 楊 青

( 山東師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，250358，濟(jì)南 )

1 引 言

有限體積法[1]又稱控制體積方法,它是從方程組的積分形式出發(fā),對積分形式進(jìn)行離散,得到有限體積格式,該方法構(gòu)造的離散格式簡單,滿足局部守恒.近年來廣泛應(yīng)用于空氣動力學(xué).緊致有限體積法[2,3]是將緊致差分[4]的思想與有限體積法結(jié)合所得到的一種高精度算法,已被用于多種類型的方程[5-7].

下邊我們考慮如下形式的一維拋物界面問題

(1)

(2)

(3)

本節(jié)針對問題(1)-(3)建立緊致有限體積格式.

2 緊致有限體積格式的建立

下面根據(jù)界點(diǎn)γ的位置對Ih進(jìn)行調(diào)整. 如果存在xi, 使得xi=γ, 則Ih不變.如果γ∈(xi-1,xi), 當(dāng)|γ-xi-1|>|γ-xi|時, 我們記xi=γ, 否則xi-1=γ; 也就是說通過節(jié)點(diǎn)轉(zhuǎn)移的方法, 使界點(diǎn)成為網(wǎng)格中的一個節(jié)點(diǎn). 為了敘述方便, 記xk=γ,h1=xk-xk-1,h2=xk+1-xk.

其次, 我們對于網(wǎng)格Ih進(jìn)行對偶剖分,xi-1/2=(xi+xi-1)/2,i=1,2,…,N. 對偶剖分單元記為

(4)

下面我們對T1,T2,T3進(jìn)行離散,分四種情況進(jìn)行討論.

1) 情況1:|i-k|≥2.

首先離散T1, 在區(qū)間[xi-1,xi+1]上, 使用插值條件(xi-1,ut(xi-1,tj-1/2)),(xi,ut(xi,tj-1/2)),

(5)

其中ξ=(x-xi)/h. 于是

其中i=1，2,…,k-2,k+2,…,N-1,j=2,3,…,M+1.

(6)

然后在時間方向上用中心差商離散, 得到

其中i=1，2,…,k-2,k+2,…,N-1,j=2,3,…,M+1.

(7)

類似地, 對于T3離散, 得到

其中i=1，2,…,k-2,k+2,…,N-1,j=2,3,…,M+1.

(8)

最后考慮T2的離散,

其中i=1，2,…,k-2,k+2,…,N-1,j=2,3,…,M+1.

(9)

又因為由原方程得:

其中i=1，2,…,k-2,k+2,…,N-1,j=2,3,…,M+1.

代入方程得

(10)

(11)

其中i=1，2，…，k-2,k+2，…，N-1,j=2,3，…，M+1.

(12)

2) 情況2:i=k-1.

由插值多項式

得

(13)

(14)

其中

對于T1的離散:

(15)

T2的離散形式類似于|i-k|≥2的情況, 即把式子(11)、(14)、(15)代入(5), 整理可得

(16)

3) 情況3:i=k.

因為

(17)

利用插值多項式得

(18)

(19)

由(17),(18),(19)整理得

(20)

又因為

(21)

利用方程(4)得

(22)

將(22)代入(21)整理得

化簡得

(23)

又因為

易得

(24)

同理，

(25)

(26)

根據(jù)(26), 類似可得T3的離散格式.由(11)得T2的離散格式，將T1、T2、T3的離散代入(5)中得

(27)

4) 情況4：i=k+1. 計算方法同i=k-1時, 利用插值多項式計算系數(shù)e,m,n代入(5)得

(28)

其中

由(12)、(16)、(27)、(28)組成了一個線性代數(shù)方程組, 該系數(shù)矩陣具有很好的三對角性質(zhì), 可以用追趕法求解. 具體計算的時候, 右端項中的關(guān)于f的積分可以使用兩點(diǎn)以上的求積公式: Gauss求積公式或者是Simpson求積公式, 格式形式與緊差分格式近似, 而這種格式是由有限體積方法導(dǎo)出的緊格式, 因而得到高階緊致有限體積格式.

3 數(shù)值算例

取相應(yīng)的精確解為

則函數(shù)f(x,t)和邊值條件φ1(x),g1(x),φ2(x)可由u(x,t)相應(yīng)得到， 應(yīng)用上述介紹的有限體積格式來計算該定解問題, 其中h=1/M,τ=1/N.取N=M2,下面圖1給出了N=M2時差分格式解U的誤差和收斂階的計算結(jié)果. 可以看出, 在不同的范數(shù)意義下, 空間的誤差階數(shù)都達(dá)到了四階或四階以上. 因此, 時間的誤差階數(shù)都達(dá)到了二階或二階以上, 說明了本文所建立格式的有效性.

圖1 N=20, M=400, 間斷點(diǎn)為r=3/8時, 數(shù)值解與真解的圖形

h最大模誤差收斂階L2模誤差收斂階120 3.53e-06— 2.08e-06—1306.68e-07≈4 3.84e-07≈41401.73e-07≈59.81e-08≈51504.62e-08≈62.59e-08≈61601.08e-08≈85.93e-08≈8