對流-擴散方程數(shù)值解的四次B樣條方法

齊梓萱，錢 江

(河海大學(xué)理學(xué)院，江蘇 南京 211100)

樣條函數(shù)是具有一定光滑性的分段或分片定義的多項式函數(shù)，研究多元樣條函數(shù)一般的方法是光滑余因子協(xié)調(diào)法[1-2]，其他經(jīng)典方法包括B網(wǎng)方法[3]與B樣條方法[4-6]。樣條函數(shù)方法廣泛應(yīng)用于數(shù)值逼近[7]、計算幾何[2]、計算機輔助幾何設(shè)計[6]、有限元、微分方程數(shù)值解等領(lǐng)域。

具體而言，文獻[8]利用光滑余因子協(xié)調(diào)法計算出2型三角剖分上的二元三次樣條基函數(shù)，構(gòu)造樣條擬插值算子[9]及分析擬插值算子導(dǎo)數(shù)逼近[10]廖肇源[11]與張勝剛[12]計算了三次樣條基函數(shù)并將其應(yīng)用于Poinsson方程的數(shù)值求解，取得了好的效果。由于計算的簡便性，三次樣條函數(shù)還被廣泛應(yīng)用于其它微分方程數(shù)值解。明祖芬和張大凱[13]給出了三次樣條函數(shù)及其導(dǎo)數(shù)的基本關(guān)系，并將其應(yīng)用于對流擴散問題，給出格式的相容性、收斂性和穩(wěn)定性條件。吳蓓蓓和徐麗[14]對三次B-樣條和三次三角B-樣條基函數(shù)引入權(quán)因子，給出了對流擴散方程的混合三次B-樣條配點法。董建強等[15]基于樣條插值和Pade逼近公式，構(gòu)造了一種求解一維對流擴散方程的高精度差分格式。CHANG等[16]提出一種基于樣條的運動方法，并用四次樣條曲線表示物體的外部輪廓，將其應(yīng)用于基于圖像的IBVS的輪廓跟蹤任務(wù)中。KERSEY[17]將樣條混合差值與樣條的光滑性推廣到實際問題。文獻[18]提出一種基于局部算法的B樣條函數(shù)擬合任意形式曲線的新方法，并通過實驗數(shù)據(jù)驗證其可用于自動逆向工程領(lǐng)域。

然而考慮到更高次樣條函數(shù)計算的復(fù)雜性與逼近具有更高的精度，高次樣條函數(shù)在微分方程數(shù)值解中的應(yīng)用值得進一步研究。本文將利用光滑余因子協(xié)調(diào)法求解一元四次B樣條基函數(shù)，并將其應(yīng)用于求解一類對流-擴散方程。

1 均勻四次B樣條基函數(shù)的計算

定理1[2]. 設(shè)一個一元n次樣條s(x)∈Sn(x0,x1,···,xN)的互異節(jié)點為{x0,x1,···,xN}，且具有局部支集(x0,xN)，則此樣條函數(shù)可表示為

并且相鄰兩區(qū)間的樣條函數(shù)滿足si(x)-si-1(x)=ci(x-xi)n，i=0,1,···,N且滿足協(xié)調(diào)方程。

利用定理1，可計算出任意次非均勻B樣條基函數(shù)，擬將樣條逼近方法應(yīng)用于微分方程數(shù)值解，本文計算出均勻節(jié)點上的四次B樣條基函數(shù)，且計算有限區(qū)間端點處的非均勻四次B樣條基函數(shù)。

節(jié)點為{xi,xi+1,xi+2,xi+3,xi+4,xi+5}均勻步長為h的四次B樣條函數(shù)。

令ci+5=λ，則ci+4=-5λ，ci+3=10λ，ci+2=-10λ，ci+1=5λ，ci=-λ，得到四次B樣條函數(shù)為

再利用B樣條基函數(shù)的單位分解性可得

定理2.步長為h的均勻節(jié)點{xi,xi+1,xi+2,xi+3,xi+4,xi+5}四次B樣條函數(shù)為且在子區(qū)間Ii:=[xi,xi+1]中存在唯一一組四次B樣條基函數(shù)

如圖1所示，所得的均勻四次B樣條基函數(shù)具有非負性、單位分解性。

圖1 四次B樣條基函數(shù)Fig. 1 Univariate quartic b-spline bases

2 區(qū)間端點帶重節(jié)點的樣條基函數(shù)

本文考慮在有界閉區(qū)間[x0,xn]內(nèi)均勻步長為h，將得到在x0，xn處具有重節(jié)點的四次B樣條的具體表示，從而得到每個支集上的四次樣條函數(shù)。

情形1.考慮x-1=x0＜x1＜x2＜x3＜x4，設(shè)在節(jié)點x-1=x0,x1,x2,x3,x4處的光滑余因子分別為。可得協(xié)調(diào)方程為，即

因此節(jié)點為x-1=x0＜x1＜x2＜x3＜x4時的四次B樣條函數(shù)為

情形2.x-2=x-1=x0＜x1＜x2＜x3，設(shè)節(jié)點x-2=x-1=x0，x1，x2，x3處的光滑余因子分別為。協(xié)調(diào)方程為

同理可得在節(jié)點為x-2=x-1=x0＜x1＜x2＜x3時的四次B樣條函數(shù)為

情形3.x-3=x-2=x-1=x0＜x1＜x2,記節(jié)點x-3=x-2=x-1=x0,x1,x2處的光滑余因子分別為。協(xié)調(diào)方程為

同理可得節(jié)點為x-3=x-2=x-1=x0＜x1＜x2時的四次B樣條函數(shù)為

情形4.x-4=x-3=x-2=x-1=x0＜x1，設(shè) 節(jié) 點x-4=x-3=x-2=x-1=x0,x1的光滑余因子為。協(xié)調(diào)方程為

同理可得在節(jié)點為x-4=x-3=x-2=x-1=x0＜x1時的四次B樣條函數(shù)為

定理3.節(jié)點為x-4=x-3=x-2=x-1=x0＜x1＜x2＜x3＜x4＜x5且xi+1-xi=h i=0,1,···,5時，在每個子區(qū)間內(nèi)存在唯一一組四次B樣條基函數(shù)為

證明：在單區(qū)間I0=[x0,x1]內(nèi)存在的B樣條函數(shù)有

由B樣條基函數(shù)的單位分解性知

由此得到式(15)。同理可證得式(16)，(17)，(18)。

證畢

類似于對左端點x0處帶重節(jié)點的四次B樣條的分析，同樣可以得到右端點xn處帶重節(jié)點的四次B樣條的結(jié)論。

定 理4.設(shè)xn-5＜xn-4＜xn-3＜xn-2＜xn-1＜xn=xn+1=xn+2=xn+3=xn+4且xi+1-xi=h,i=n-5,n-4,···,n-1，則在每個子區(qū)間內(nèi)存在唯一一組的四次B樣條基函數(shù)

證明：在單區(qū)間In-1=[xn-1,xn]中，存在的B樣條函數(shù)有

由B樣條基函數(shù)的單位分解性知

同理可證得式(22)，(23)，(24)。

證畢

3 對流-擴散方程的樣條逼近離散格式

考慮一維對流-擴散方程

其中，θ是關(guān)于空間x與時間t的二元函數(shù)，其為在時刻T點x處的體積濃度；u為對流速度；Dx為擴散系數(shù)，函數(shù)f(x)，g0(x)，g1(x)，g2(x)是充分光滑的。

表1為函數(shù)si(x)及其導(dǎo)數(shù)在各節(jié)點處的值。

表1 si(x)及其導(dǎo)數(shù)在各節(jié)點處的值Table 1 The value of si(x) and its derivative at each node

記點(xi,tk)處的逼近解為在區(qū)間[xi,xi+1]上為。

點(xi,tk)處時間變量采用一階向前差分，引入?yún)?shù)δ(0≤δ≤1)，式(27)可離散化為

將式(30)帶入式(31)中，生成具有n+4個未知量的n+1個方程。為了使方程組具有唯一解，還需要另外添加3個方程。由邊界條件可得

n+4個方程組的矩陣表示為

求出C0則第k+1時間層的逼近解便可由式(33)計算得到。

另外，若采用向前差分法，右邊值采用向后差分離散方程，則得到的矩陣表示為

由于所采用的四次樣條基函數(shù)具有非負性與單位分解性，因此求解過程是穩(wěn)定的，且精度高于基于三次樣條得到的求解結(jié)果。限于篇幅，后續(xù)將進一步研究基于四次B樣條的樣條擬插值算子，并將其應(yīng)用于求解偏微分方程，從理論上分析四次樣條求解微分方程的穩(wěn)定性與收斂性。

4 數(shù)值算例

考慮對流-擴散方程

其準確解為

方法1.樣條逼近方法。根據(jù)式(31)并取參數(shù)δ=0.5，得到

取時間步長Δt=0.04，空間步長h=0.02，得Ck+1滿足式(33)，其中，。C0滿足式(36)，可求得C0=(1.0304,1.0100,0.9900,0.9704,0.9512,0.9323,0.9139,···,0.4025,0.3945,0.3867,0.3791,0.3715,0.3642,0.3570)T。由此可通過式(33)解得相應(yīng)的Ck+1。

方法2.有限差分法。用差分方法離散方程，同樣取δ=0.5，離散化得

取Δt=0.04，h=0.02，可得式(34)，其中，,,。

表2取t=0.2，對x取不同值與文獻[13]中的三次樣條解做比較。

表2 四次樣條解與三次樣條解對比Table 2 The quartic spline solution is compared with the cubic spline solution

表3為表2的2種方法所求解與精確解的數(shù)值誤差對比。

表3 數(shù)值誤差對比Table 3 Numerical error comparison

表4取t=0.04，在[0.9,1.0]區(qū)間內(nèi)與差分方法的解進行比較。

表4 四次樣條解與差分方法解對比Table 4 The quartic spline solution is compared with the difference method solution

表5為2種方法所求解與精確解的數(shù)值誤差對比。

表5 數(shù)值誤差對比Table 5 Numerical error comparison

對于式(42)的離散，方法1和方法2均對時間變量在(xi,tk)處取一階向前差分，而對空間變量的離散，方法1引入?yún)?shù)δ并在(xi,tk)，(xi,tk+1)處取樣條函數(shù)的一階、二階導(dǎo)數(shù)形式；方法2需要考慮2種差分方式，同樣引入?yún)?shù)δ對方程進行一階及二階向前差分，對于右端點的處理采用一階及二階向后差分。進一步，本文算例中通過對2種方法進行數(shù)值比較可知，樣條方法的逼近效果更好，且Matlab計算表明，采用樣條逼近方法得到解的時間為5.939 365 s，而有限差分法運行的時間為10.668 284 s。因此，選擇樣條逼近方法更簡便實用。并且通過對于t=0.2時與三次樣條逼近解的對比可知，四次樣條的逼近比三次樣條逼近效果好。