基于五次B 樣條的對(duì)流-擴(kuò)散方程數(shù)值解法

錢 江，王永杰

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

1975 年,王仁宏教授利用函數(shù)論與代數(shù)幾何的方法開創(chuàng)性地建立了任意剖分下多元樣條函數(shù)的理論框架,提出了研究多元樣條最一般的方法,即光滑余因子協(xié)調(diào)法[1-3],并由此，取得了豐碩的多元樣條理論與應(yīng)用研究成果。具體而言,利用光滑余因子協(xié)調(diào)法可以計(jì)算出2 型三角剖分上的二元三次樣條基函數(shù)[4]；進(jìn)一步,利用保多項(xiàng)式性,構(gòu)造出基于線性泛函的樣條擬插值算子[5]，推導(dǎo)出樣條擬插值的導(dǎo)數(shù)逼近方法[6]。此外,多元樣條在其他領(lǐng)域也獲得了不錯(cuò)的成果。文獻(xiàn)[7]利用三次樣條方法給出了Poisson 方程的數(shù)值求解。文獻(xiàn)[8]建立了三次樣條函數(shù)與導(dǎo)數(shù)的計(jì)算公式,并應(yīng)用于求解對(duì)流擴(kuò)散問題。文獻(xiàn)[9]提出了一種樣條運(yùn)動(dòng)方法,將四次樣條曲線表示物體的外部輪廓,從而探索基于圖像的IBVS 的輪廓跟蹤任務(wù)。白根柱[10]構(gòu)造了一段特定條件下的三次有理B 樣條曲線,即曲線通過給定的控制多邊形始末端點(diǎn),且與始末控制多邊形相切,后以該曲線為軸線的管道光滑拼接出軸線異面管道。王慧等人[11]為了高效地獲得擬合曲線,在均勻三次B 樣條曲線的基礎(chǔ)上,提出一種雙層最小二乘漸進(jìn)迭代逼近算法。文獻(xiàn)[12]在混合B 樣條的基礎(chǔ)上構(gòu)造新的技術(shù),給出一種非結(jié)構(gòu)化樣條實(shí)體等幾何拓?fù)鋬?yōu)化方法。沈菀薔等人[13]在B 樣條基礎(chǔ)上給出一種變次數(shù)樣條曲線的細(xì)分算法,該算法可在細(xì)分前指定每段的次數(shù)和兩段曲線之間的連續(xù)性。文獻(xiàn)[14]為了提高三維模型之間拼接曲面的精度和效率,提出一種基于三次均勻B 樣條曲線曲面的網(wǎng)格融合方法。文獻(xiàn)[15]為了證明半離散格式解的有界性與收斂性,用三次B樣條有限元法給出一類四階主項(xiàng)帶有變系數(shù)的拋物方程。文獻(xiàn)[16]為了獲得的平滑航跡整體過渡自然,采用B 樣條曲線插補(bǔ)法求解航跡的平滑問題。齊遠(yuǎn)節(jié)等[17]提出了兩種新的求解對(duì)流擴(kuò)散方程的三次樣條差分格式。王森[18]把牛頓插值公式拓展到埃爾米特插值公式和樣條函數(shù)插值的使用領(lǐng)域，給出了解決插值問題簡(jiǎn)便高效的方法。

考慮到更高次的樣條函數(shù)在逼近理論中應(yīng)具有更高的精度[19],而計(jì)算高次的B 樣條需要較復(fù)雜的方法與技巧,由此將高次樣條函數(shù)應(yīng)用于微分方程數(shù)值解值得進(jìn)一步研究。作者擬考慮建立五次B 樣條方法,同時(shí)注意到在比較五次B 樣條微分正交法時(shí)發(fā)現(xiàn),文章[20]僅僅給出無重節(jié)點(diǎn)情況下的五次B 樣條函數(shù),且該函數(shù)并不具有單位分解性，在逼近曲線的過程中優(yōu)勢(shì)不明顯。而利用de Boor-Cox 公式計(jì)算均勻節(jié)點(diǎn)下的五次B 樣條需要558 次運(yùn)算,同時(shí)計(jì)算帶有重節(jié)點(diǎn)的五次B樣條運(yùn)算步驟更多,更加繁瑣。在多種情況考慮下,本文將利用光滑余因子協(xié)調(diào)法,從最小支集的角度計(jì)算在均勻剖分下的五次B 樣條基函數(shù)[21],并將其運(yùn)用在求解一類對(duì)流擴(kuò)散方程中。

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

首先給出光滑余因子協(xié)調(diào)法[1-3]的定理作為計(jì)算五次B 樣條基函數(shù)的準(zhǔn)備工作。

定理 1設(shè)一個(gè)一元 n 次樣條s(x)∈S(x0,x1,x2,…,xN) 的互異節(jié)點(diǎn)為{x0,x1,…,xN},且具有局部支集(x0,xN),則此樣條函數(shù)可表示為

利用定理1,可以計(jì)算出任意次均勻B 樣條基函數(shù),擬將樣條逼近方法應(yīng)用于微分方程數(shù)值解,本文將給出在均勻節(jié)點(diǎn)上不同情況下的五次B 樣條基函數(shù),嘗試得到運(yùn)用五次樣條基函數(shù)擬合對(duì)流-擴(kuò)散方程的方法。

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

這樣,給出了步長(zhǎng)為h,均勻節(jié)點(diǎn)為{xi,xi+1,xi+2,xi+3,xi+4,xi+5,xi+6} 的五次B 樣條函數(shù),即

且在子區(qū)間Ii:=[xi,xi+1] 中存在唯一一組五次B 樣條基函數(shù)

如圖1 所示,通過光滑余因子協(xié)調(diào)法[2-4]得到的均勻五次B 樣條基函數(shù)具有非負(fù)性、單位分解性等優(yōu)秀的性質(zhì)。

圖1 均勻五次B 樣條基函數(shù)

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

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

情形1.考慮x-1=x0＜x1＜x2＜x3＜x4＜x5,設(shè)在節(jié)點(diǎn)x-1=x0,x1,x2,x3,x4,x5處的光滑余因子分別為?？傻脜f(xié)調(diào)方程為

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

情形4.在x-4=x-3=x-2=x-1=x0＜x1＜x2的情況下,我們?cè)O(shè)節(jié)點(diǎn)x-4=x-3=x-2=x-1=x0,x1,x2的光滑余因子為

同理,在節(jié)點(diǎn)為x-4=x-3=x-2=x-1=x0＜x1＜x2時(shí)的五次B 樣條函數(shù)為

情形5.當(dāng)x-5=x-4=x-3=x-2=x-1=x0＜x1時(shí),我們?cè)O(shè)節(jié)點(diǎn)x-4=x-3=x-2=x-1=x0,x1處的光滑余因子分別為

同理可得節(jié)點(diǎn)為x-5=x-4=x-3=x-2=x-1=x0＜x1時(shí)的五次B 樣條函數(shù)為

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

在I0=[x0,x1] 上,基函數(shù)為

圖2 I0=[x 0,x1] 上的五次B 樣條基函數(shù)圖像

在I1=[x1,x2] 上,基函數(shù)為

圖3 I1=[x1,x2]上的五次B 樣條基函數(shù)圖像

在I2=[x2,x3]上,基函數(shù)為

圖4 I2=[x2,x3]上的五次B 樣條基函數(shù)圖像

在I3=[x3,x4]上,基函數(shù)為

圖5 I3=[x3,x4]上的五次B 樣條基函數(shù)圖像

在I4=[x4,x5]上,基函數(shù)為

圖6 I4=[x4,x5]上的五次B 樣條基函數(shù)圖像

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

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

由此得到式(18)。同理可證得式(19),(20),(21),(22)。

證畢

類似于對(duì)左端點(diǎn)處帶重節(jié)點(diǎn)的五次B 樣條的分析,同樣可以得到右端點(diǎn)處帶重節(jié)點(diǎn)的五次B樣條的結(jié)論。

定 理 3設(shè)xn-6＜xn-5＜xn-4＜xn-3＜xn-2＜xn-1＜xn=xn+1=xn+2=xn+3=xn+4=xn+5且xi+1-xi=h,i=n-6,n-5,…,n-1,則在每個(gè)子區(qū)間內(nèi)存在唯一一組的五次B 樣條基函數(shù),即在In-1=[xn-1,xn]上,基函數(shù)為

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

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

3 對(duì)流-擴(kuò)散方程的樣條逼近格式

令u(x,t)為如下初邊值問題的真解[22]：

其中,u為在時(shí)刻T點(diǎn)x處的體積濃度,v是對(duì)流速度,Dx為擴(kuò)散系數(shù)。

對(duì)方程(32)在點(diǎn)(xi,tk)按時(shí)間步長(zhǎng)進(jìn)行離散,引入?yún)?shù)δ(0 ≤δ≤1),我們得到

其中,是u(x,t)在點(diǎn)(xi,tk),xi=ih,tk=kτ處的逼近解。h=1∕I,τ=T∕M,I＞1 和K＞1 是兩個(gè)整數(shù)。

接下來,利用均勻節(jié)點(diǎn)下的五次B 樣條基函數(shù)來擬合對(duì)流-擴(kuò)散方程。記點(diǎn)()xi,tk處的逼近解為

首先給出均勻節(jié)點(diǎn)下的五次B 樣條函數(shù)si(x)及其導(dǎo)數(shù)在各節(jié)點(diǎn)處的值,如表1。

表1 及其導(dǎo)數(shù)在各節(jié)點(diǎn)處的值

根據(jù)表1 可知,逼近解及其關(guān)于x的一階和二階導(dǎo)函數(shù)分別為

將式(37)代入式(35)中,生成具有n+5 個(gè)未知量的n+1 個(gè)方程。為了使方程組具有唯一解,還需要另外添加4 個(gè)方程。由邊界條件可得

n+5 個(gè)方程組的矩陣矩陣表示為

其中,P,R∈R(n+5)×(n+5);Ck+1,Ck,b∈R(n+5),

求出C0則第k+1 時(shí)間層的逼近解就可由(39)式計(jì)算得到。

4 數(shù)值算例

考慮模型問題(32)-(34),選擇精確解為

其中x∈[0,1],t∈[0,1]。相應(yīng)的初邊值可以由它得到。

選擇用五次B 樣條基函數(shù)逼近該模型。取時(shí)間步長(zhǎng)τ=0.04,空間步長(zhǎng)h=0.02,δ=0.5。在表2 中,我們?nèi)=0.2,對(duì)取不同值與精確解作對(duì)比。

表2 五次樣條解、四次樣條解與精確解對(duì)比

表3 為表2 的兩種方法得到的解與精確解的數(shù)值誤差對(duì)比。

表3 與精確解的數(shù)值誤差對(duì)比

從表2 和表3 可以看到,用五次B 樣條逼近得到的解準(zhǔn)確率高,而且得到的解與精確解的誤差遠(yuǎn)小于比它低一次的四次B 樣條解誤差,逼近效果更好。用對(duì)流-擴(kuò)散方程舉例,雖然常用三次樣條逼近,但在相同步長(zhǎng)的情況下,在相同的有限閉區(qū)間上需要更少的五次B 樣條，取得更好的逼近效果。

需要指出的是,為了更好地分析五次樣條的逼近效果,首先需要建立具有保高次多項(xiàng)式性的五次樣條擬插值算子,這樣可以建立樣條解與精確解之間的誤差分析,目前由于這方面工作計(jì)算量很大，主要限于三次樣條情形。限于篇幅,只計(jì)算出滿足單位分解性,即保常數(shù)的五次B 樣條(包括均勻節(jié)點(diǎn)與帶重節(jié)點(diǎn)切線),再結(jié)合樣條的分段光滑性配置系數(shù),從而遞推計(jì)算出這些系數(shù),避免了Gauss 消去法求解線性方程組的繁瑣過程(注：遞推計(jì)算復(fù)雜度往往是ο(n2),而消去法計(jì)算復(fù)雜度是ο(n3),與差分法的不連續(xù)性。數(shù)值算例也表明,采用五次B 樣條計(jì)算出的數(shù)值解逼近效果優(yōu)于采用四次B 樣條的計(jì)算結(jié)果。另外,本文所提的五次B 樣條遞推算法不僅應(yīng)用于求解對(duì)流-擴(kuò)散方程,而且適用于一般的線性雙曲型方程或方程組。

5 結(jié)語(yǔ)

由于所采用的的五次樣條函數(shù)具有非負(fù)性與單位分解性,因此具有穩(wěn)定的求解過程,精度高于基于三次和四次樣條函數(shù)得到的結(jié)果,并且用光滑余因子協(xié)調(diào)法計(jì)算出的樣條函數(shù)在區(qū)間端點(diǎn)處的連接處光滑性更高。限于篇幅,后續(xù)將進(jìn)一步研究基于五次B 樣條的樣條擬插值算子,將其應(yīng)用于求解偏微分方程,從理論上分析五次樣條求解微分方程的穩(wěn)定性與收斂性。

在以后的工作中,作者將建立保高次多項(xiàng)式性的五次樣條擬插值算子,并分析逼近誤差,將其應(yīng)用于求解微分方程或方程組。