求解拋物型方程的兩層高精度差分格式

陳佳欣，李旭生

（遼寧石油化工大學(xué)理學(xué)院，遼寧撫順113001）

拋物型方程常常應(yīng)用于擴(kuò)散、熱傳導(dǎo)等領(lǐng)域，而且應(yīng)用的范圍很廣，研究方向也比較多[1?2]。學(xué)者們對(duì)拋物型方程進(jìn)行了大量分析，進(jìn)而豐富了拋物型方程理論，也加深了拋物型方程與其他自然學(xué)科領(lǐng)域的聯(lián)系，但是拋物型方程形式復(fù)雜，解析形式的解很難求出，所以討論方程的數(shù)值解法十分有意義。隨著對(duì)拋物型方程研究的不斷深入，對(duì)拋物型方程數(shù)值解法的研究也越來越廣泛，通常要面臨如何解決實(shí)際問題，提出一個(gè)好的方法顯得尤為重要。在構(gòu)建不同的差分方程過程中，主要的目的是得到一個(gè)具有高精度的截?cái)嗾`差且結(jié)果更為精確和穩(wěn)定的差分格式。就目前來說，構(gòu)造差分格式的方法有很多種，如有限差分法、有限元法、有限體積法、譜方法、pade逼近法等[3?4]。最常用的是有限差分法，這種方法的基本思想是將連續(xù)的問題離散化，對(duì)求解的區(qū)間用網(wǎng)格進(jìn)行剖分，然后將微分算子離散化，從而把微分方程的定解問題轉(zhuǎn)化為線性代數(shù)方程的求解問題[5]。

姜蘊(yùn)芝等[6]根據(jù)Saul′ev算法對(duì)古典隱格式進(jìn)行了一定改進(jìn)，得到求解一維拋物方程的新格式，一種無條件且穩(wěn)定的半隱格式，該格式的截?cái)嗾`差為并從半隱格式出發(fā)構(gòu)造出一組顯格式。半隱格式的穩(wěn)定性與Saul′ev算法的穩(wěn)定性相近，顯格式的差分格式精度要高于Saul′ev算法。吐克孜·艾肯等[7]用C?N算法將偏微分方程轉(zhuǎn)化為常微分方程組，通過使用指數(shù)函數(shù)的Trotter積分公式將該方程組系數(shù)矩陣分塊成多個(gè)小矩陣，并通過C?N算法求得結(jié)果。該方法計(jì)算量少且恒穩(wěn)定。張莉等[8]基于組合差商法提出一個(gè)絕對(duì)穩(wěn)定的隱式格式，并將空間上的6個(gè)點(diǎn)進(jìn)行taylor展開，使其達(dá)到一定精度。曾文平[9]構(gòu)造一組3層高精度的雙參數(shù)隱式差分格式，當(dāng)參數(shù)α=0.5、β=0時(shí)，得到一個(gè)2層格式，當(dāng)參數(shù)為非負(fù)數(shù)時(shí)，該格式是絕對(duì)穩(wěn)定的，其截?cái)嗾`差在空間和時(shí)間步長(zhǎng)上的精度分別為6階精度和2階精度。W.R.Hao等[10]通過區(qū)域分裂的方法得到一個(gè)求解拋物型方程的高精度顯式格式，該格式的截?cái)嗾`差為O(τ3+h3)。詹涌強(qiáng)等[11]通過使用待定系數(shù)法構(gòu)造一個(gè)求解拋物方程的高精度隱式算法，該隱式差分格式的截?cái)嗾`差為O(τ3+h4),當(dāng)r＞時(shí)，該差分格式是穩(wěn)定的。周敏等[12]首先進(jìn)行網(wǎng)格剖分，然后選取線性無關(guān)的節(jié)點(diǎn)，利用這些節(jié)點(diǎn)構(gòu)造出一個(gè)含有多個(gè)參數(shù)的差分格式，并對(duì)該格式的某一點(diǎn)進(jìn)行泰勒展開，使用待定系數(shù)法使截?cái)嗾`差達(dá)到一定的精度，然后求解參數(shù)。得到一個(gè)2層8個(gè)點(diǎn)的隱式差分格式，此格式的截?cái)嗾`差為O(τ3+h5)，在0.001＜r＜0.231或0.236＜r＜0.772具有穩(wěn)定性。M.S.Ma等[13]構(gòu)造了一個(gè)時(shí)間步長(zhǎng)和空間步長(zhǎng)均為4階的高精度隱式格式，該隱格式的穩(wěn)定性條件為0＜r≤1。文獻(xiàn)[14]構(gòu)造出一個(gè)求解拋物型方程的3層高精度隱式格式，數(shù)值算例證明了格式的精確性和可行性。

本文利用待定系數(shù)法和泰勒展開法構(gòu)造出新的隱式差分格式，該差分格式精度為時(shí)間上的3階精度，空間上的5階精度，比文獻(xiàn)[8—11]精度高，該格式數(shù)值解比文獻(xiàn)[12—13]的數(shù)值解更接近精確解，文獻(xiàn)[14]的精確性與本文精確性相當(dāng)，但是文獻(xiàn)[14]的算法為3層差分格式，3層差分格式在迭代計(jì)算時(shí)不是自啟動(dòng)的。提出的差分格式為2層格式，相比3層差分格式，在保持較好的精度的同時(shí)，差分格式也較為簡(jiǎn)單，然后對(duì)格式的穩(wěn)定性進(jìn)行了理論分析，最后得出一系列數(shù)值解驗(yàn)證此格式的可行性。

1 差分格式的構(gòu)造

考慮如下拋物型方程的初邊值問題：

式中，L為空間長(zhǎng)度；T為時(shí)間長(zhǎng)度；u(x,0)為初始值；u(0,t)、u(L,t)為邊界值；a為正常數(shù)。

首先選取空間步長(zhǎng)為h=LK，然后選取時(shí)間步長(zhǎng)為τ=TM，其中K、M都是正整數(shù)，然后用兩簇平行線將矩形區(qū)域分割成矩形網(wǎng)格，再用適當(dāng)?shù)牟钌檀娣匠讨邢鄳?yīng)的偏微商，從而求得所要的差分格式。

對(duì)[0,L]×[0,T]這個(gè)矩形區(qū)域進(jìn)行網(wǎng)格剖分，在網(wǎng)格點(diǎn)上選取局部節(jié)點(diǎn)集,(xj+1,tk+1),(xj-1,tk),(xj,tk)(xj+1,tk)，其中xj=jh,tk=kτ。將6個(gè)節(jié)點(diǎn)的u值點(diǎn)分別在節(jié)點(diǎn)（xj,tk）進(jìn)行泰勒展開，得：

式中，ukj=u(xj,tk)。

然后分別求出差商的漸近表達(dá)式，其中Δtunj=：

當(dāng)問題(1)的解是充分光滑且連續(xù)時(shí)，通過aτ=h2r得：

使用式(12)整理得：

利用所求得的差商建立一個(gè)含有多個(gè)未知參數(shù)的差分算子：

式中，ci為未知參數(shù)。

目的是用來逼近本文求解的微分算子：

而在本文中所建立的差分格式中的未知數(shù)c1、c2、c3、c4、c5，即本文需要求解的系數(shù)。

為了使式(20)差分格式的截?cái)嗾`差精度達(dá)到O(τ3+h5)，則需要滿足：

求解各個(gè)系數(shù)分別為：

將c1、c2、c3、c4、c5分別代入差分算子中求得差分格式為：

2 穩(wěn)定性和收斂性

令差分格式中的系數(shù)分別對(duì)應(yīng)A、B、C、D、E、F。

再利用Fourier分析法將差分格式化為矩陣形式，通過令ukj=vkeijθ進(jìn)而可以將差分格式轉(zhuǎn)化為矩陣的形式。

通過對(duì)矩陣進(jìn)行消元，其目的是消去eijβ，得到新的格式為：

然后分別將A、B、C、D、E、F代入式(28)得：

通過計(jì)算得出所需要的增長(zhǎng)因子為：

為了保證建立的差分格式的穩(wěn)定性，則需滿足|G|≤1的充分條件為：

3 數(shù)值算例

考慮拋物型方程中一維常系數(shù)擴(kuò)散方程為：

表1 當(dāng)r=0.1、n=100時(shí)不同格式的數(shù)值解

表2 格式1數(shù)值解與精確解的比較結(jié)果

4 結(jié) 論

通過待定系數(shù)法構(gòu)造出一個(gè)求解拋物型方程的兩層差分格式，格式的截?cái)嗾`差為O(τ3+h5)，并對(duì)格式的穩(wěn)定性進(jìn)行分析，即，該格式是穩(wěn)定的。從表1可以看出，文獻(xiàn)[12—13]的絕對(duì)誤差比格式1絕對(duì)誤差大4個(gè)數(shù)量級(jí)，文獻(xiàn)[14]與格式1精確性相當(dāng)，但是文獻(xiàn)[14]為3層差分格式。從表2可以看出，格式1數(shù)值解與精確解基本吻合，此格式具有精確性和可行性。