于 倩 楊 青
( 山東師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,250358,濟南 )
有限體積法[1]又稱控制體積方法,它是從方程組的積分形式出發(fā),對積分形式進行離散,得到有限體積格式,該方法構(gòu)造的離散格式簡單,滿足局部守恒.近年來廣泛應(yīng)用于空氣動力學(xué).緊致有限體積法[2,3]是將緊致差分[4]的思想與有限體積法結(jié)合所得到的一種高精度算法,已被用于多種類型的方程[5-7].
下邊我們考慮如下形式的一維拋物界面問題
(1)
(2)
(3)
本節(jié)針對問題(1)-(3)建立緊致有限體積格式.
下面根據(jù)界點γ的位置對Ih進行調(diào)整. 如果存在xi, 使得xi=γ, 則Ih不變.如果γ∈(xi-1,xi), 當(dāng)|γ-xi-1|>|γ-xi|時, 我們記xi=γ, 否則xi-1=γ; 也就是說通過節(jié)點轉(zhuǎn)移的方法, 使界點成為網(wǎng)格中的一個節(jié)點. 為了敘述方便, 記xk=γ,h1=xk-xk-1,h2=xk+1-xk.
其次, 我們對于網(wǎng)格Ih進行對偶剖分,xi-1/2=(xi+xi-1)/2,i=1,2,…,N. 對偶剖分單元記為
(4)
下面我們對T1,T2,T3進行離散,分四種情況進行討論.
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的積分可以使用兩點以上的求積公式: Gauss求積公式或者是Simpson求積公式, 格式形式與緊差分格式近似, 而這種格式是由有限體積方法導(dǎo)出的緊格式, 因而得到高階緊致有限體積格式.
取相應(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ù)都達到了四階或四階以上. 因此, 時間的誤差階數(shù)都達到了二階或二階以上, 說明了本文所建立格式的有效性.
圖1 N=20, M=400, 間斷點為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