楊苗苗， 葛永斌

（寧夏大學(xué) 數(shù)學(xué)統(tǒng)計(jì)學(xué)院，寧夏 銀川750021）

對(duì)流擴(kuò)散反應(yīng)方程作為一種基本的運(yùn)動(dòng)方程，在很多科學(xué)和工程技術(shù)領(lǐng)域都有著非常廣泛地應(yīng)用，可用于描述流體的流動(dòng)與傳熱、燃燒與爆炸、核工業(yè)中反應(yīng)堆的冷卻及工業(yè)生產(chǎn)中化學(xué)堆的反應(yīng)等現(xiàn)象.因此對(duì)于該類方程精確數(shù)值解的研究具有很高的理論價(jià)值和實(shí)際作用.

求解該類方程的方法有很多種，最常用的是有限差分法.目前已經(jīng)發(fā)展了很多具有高精度并且緊致的差分格式，如針對(duì)對(duì)流擴(kuò)散方程：Sun等［1］利用理查森外推法與算子插值法改進(jìn)了四階緊致差分公式，使格式具有六階精度；田振夫［2］在迎風(fēng)變換的基礎(chǔ)上，將對(duì)流擴(kuò)散方程轉(zhuǎn)化成含源純擴(kuò)散方程，并利用Hermite法推導(dǎo)出指數(shù)型的差分格式，具有四階精度；梁昌弘等［3］利用3個(gè)點(diǎn)上的未知函數(shù)值及四階Padé格式將一階導(dǎo)數(shù)代換的方法，借助顯式緊致格式和隱式緊致格式的思想，得到了一種性能更優(yōu)的四階混合型格式；Tian等［4］建立了一種可用于求解定常對(duì)流擴(kuò)散方程的高精度指數(shù)型有限差分格式，適合大梯度問(wèn)題或邊界層問(wèn)題的求解；Chen等［5］在差分格式的構(gòu)造中將攝動(dòng)方法與一階迎風(fēng)格式相結(jié)合，得到了求解二維定常對(duì)流擴(kuò)散方程的一種高精度格式；Gupta等［6］采用泰勒級(jí)數(shù)展開法得到具有四階精度的有限差分格式；Spotz［7］則利用截?cái)嗾`差余項(xiàng)修正的方法給出了一種多項(xiàng)式型的四階緊致差分格式.盡管文獻(xiàn)［6-7］格式的推導(dǎo)方法不一樣，但本質(zhì)上是同一格式；Zhao等［8］將該方法推廣到變系數(shù)的對(duì)流擴(kuò)散反應(yīng)方程中，還給出了所得格式的收斂性分析；劉明會(huì)［9］基于原方程代入和四階緊致差分公式法，將二階導(dǎo)數(shù)代替后得到了一種新的格式.文獻(xiàn)［8-9］格式盡管推導(dǎo)過(guò)程及方法不一樣，但本質(zhì)上仍是同一格式.田振夫［10］用截?cái)嗾`差余項(xiàng)修正法，將余項(xiàng)中的三階導(dǎo)數(shù)和四階導(dǎo)數(shù)用一階和二階導(dǎo)數(shù)替換，代入原方程整理便得到四階精度的差分格式，將該方法推廣到二維后，則得文獻(xiàn)［11］中的格式；其他針對(duì)對(duì)流擴(kuò)散反應(yīng)方程的數(shù)值求解方法還有，魏劍英［12］首先由常系數(shù)變易法得到3個(gè)點(diǎn)上的指數(shù)型差分格式，再利用源項(xiàng)在中心點(diǎn)做二階泰勒級(jí)數(shù)展開，得到求解常系數(shù)的四階差分格式；田芳等［13］借助常系數(shù)指數(shù)型高精度緊致差分格式，采用殘量修正法得到變系數(shù)對(duì)流擴(kuò)散反應(yīng)方程的緊致差分格式；田芳等［14］還基于泰勒級(jí)數(shù)展開，給出了一、二階導(dǎo)數(shù)的表達(dá)式，構(gòu)造了非均勻網(wǎng)格上的對(duì)流擴(kuò)散反應(yīng)方程的多項(xiàng)式型的差分格式；蘭斌等［15］采用坐標(biāo)變換法將原方程由物理空間的非均勻網(wǎng)格變換為計(jì)算空間的均勻網(wǎng)格，再結(jié)合中心差分格式得到了求解對(duì)流擴(kuò)散反應(yīng)方程的緊致差分格式.

本文針對(duì)一維、二維定常對(duì)流擴(kuò)散反應(yīng)方程，基于截?cái)嗾`差余項(xiàng)修正的方法推導(dǎo)建立了一種新的緊致差分格式，具有四階精度.盡管本文方法與文獻(xiàn)［7-11］一樣采用到了截?cái)嗾`差余項(xiàng)修正法，但本文僅用到了一階導(dǎo)數(shù)的表達(dá)式來(lái)修正未知函數(shù)中的三階和四階導(dǎo)數(shù)項(xiàng)，而沒(méi)有用到二階導(dǎo)數(shù)項(xiàng)，這樣做的目的是使格式具有更小的耗散誤差.

1 一維問(wèn)題差分格式的建立

首先，考慮如下一維對(duì)流擴(kuò)散反應(yīng)方程：

其中a＞0為擴(kuò)散項(xiàng)系數(shù)，一般是常數(shù)；u（x）是未知函數(shù)，p（x）是對(duì)流項(xiàng)系數(shù)，s（x）是反應(yīng)項(xiàng)系數(shù)，f（x）為源項(xiàng)，并且要求u（x）、p（x）、s（x）、f（x）具有充分的光滑性.當(dāng)s（x）≡0時(shí)，模型方程（1）為對(duì)流擴(kuò)散方程.

將區(qū)間［c，d］分為N等份：x i＝c＋ih，i＝0，1，…，N，定義網(wǎng)格步長(zhǎng)h＝（d－c）／N；中心差分算子：

由泰勒級(jí)數(shù)展開可得：

將（3）和（4）式代入（1）式得

由（1）式得

對(duì)（6）式兩邊關(guān)于x求導(dǎo)可得

則（7）式為文獻(xiàn)［7-10］中處理三階導(dǎo)數(shù)項(xiàng)的方法，即三階導(dǎo)數(shù)由一階和二階導(dǎo)數(shù)來(lái)表示.本文為了減小耗散誤差，將（6）式代入（7）式并整理可得

（8）式為本文處理三階導(dǎo)數(shù)項(xiàng)的方法，即u xxx i僅由一階導(dǎo)數(shù)來(lái)表示.這是本文格式區(qū)別于文獻(xiàn)［7-10］的地方.同理，對(duì)四階導(dǎo)數(shù)的處理如下：

即u xxxx i也僅用到一階導(dǎo)數(shù)u x i，沒(méi)有如文獻(xiàn)［7-10］一樣用到二階導(dǎo)數(shù)u xx i.

將（8）和（9）式代入（5）式并舍去 Ο（h4）得

對(duì)上式中的對(duì)流項(xiàng)和反應(yīng)項(xiàng)系數(shù)的導(dǎo)數(shù)利用中心差分公式（2）離散可得由該過(guò)程可知，本文格式的截?cái)嗾`差為Ο（h4），即格式（10）在空間上可達(dá)四階精度.對(duì)于f（x）的導(dǎo)數(shù)項(xiàng)的計(jì)算，可直接代入精確解，也可用中心差分離散，并不會(huì)影響格式的整體精度.另外，注意到本文一維格式只用到3個(gè)網(wǎng)格點(diǎn)，形成三對(duì)角型的代數(shù)方程組，可采用追趕法進(jìn)行求解.

2 誤差分析

下面分別對(duì)格式（10）的耗散誤差和色散誤差進(jìn)行分析.

為了簡(jiǎn)化計(jì)算，將本文格式（10）化為常系數(shù)

則有

同理

將（12）～（14）式代入（11）式得

化簡(jiǎn)得

因此，本文差分格式的特征函數(shù)可表示為

其中

另外，為了進(jìn)行對(duì)比分析，通過(guò)計(jì)算可得文獻(xiàn)［3］中MHOC格式的特征函數(shù)為

其中

文獻(xiàn)［8］中FOC格式的特征函數(shù)為

其中

圖1～4分別給出了當(dāng)η＝1，Pe取不同數(shù)值時(shí)，MHOC格式［3］、FOC格式［8］和本文格式的耗散性及色散性的誤差分析.不難發(fā)現(xiàn)，對(duì)耗散誤差而言，文本格式比MHOC格式和FOC格式的耗散誤差均??；但就色散性而言，本文格式比FOC格式要好，但不及MHOC格式.

圖1 當(dāng)Pe＝10時(shí)耗散誤差Fig.1 The dissipation error when Pe＝10

圖2 當(dāng)Pe＝1 000時(shí)耗散誤差Fig.2 The dissipation error when Pe＝1 000

圖3 當(dāng)Pe＝10時(shí)色散誤差Fig.3 The dispersion error when Pe＝10

圖4 當(dāng)Pe＝1 000時(shí)色散誤差Fig.4 The dispersion error when Pe＝1 000

3 二維問(wèn)題差分格式的建立

考慮如下二維對(duì)流擴(kuò)散反應(yīng)方程

其中，a，b＞0為擴(kuò)散項(xiàng)系數(shù)，一般是常數(shù)；u（x，y）是未知函數(shù)，p（x，y）和q（x，y）是對(duì)流項(xiàng)系數(shù)，s（x，y）是反應(yīng)項(xiàng)系數(shù)，f（x，y）為源項(xiàng)，并且要求u（x，y）、p（x，y）、q（x，y）、s（x，y）、f（x，y）具有充分的光滑性.

將區(qū)間［c，d］分為N等份：x i＝c＋ih，y j＝c＋jh，i，j＝0，1，…，N，定義網(wǎng)格步長(zhǎng)h＝（d－c）／N.

將（15）式可轉(zhuǎn)化成如下2個(gè)一維形式的方程：

顯然方程（16）和（17）從形式上看可認(rèn)為是一維的.因此，可直接利用一維對(duì)流擴(kuò)散反應(yīng)方程的四階精度差分方法對(duì)其進(jìn)行離散，即可得在x方向有

其中

對(duì)（16）式中f1（x，y）兩邊關(guān)于x求偏導(dǎo)，可得：

將（20）和（21）式代入（19）式整理得

同理可得在y方向有

其中

于是有

注意（15）式為一種輔助關(guān)系，由（18）和（23）式的離散整理可得

其中

定義中心差分算子：

將（26）式代入（25）式進(jìn)行離散：

整理可得（15）式的離散格式，即為本文格式：

其中

其中，A0，A1，…，A8中的函數(shù)p、q、s及其導(dǎo)數(shù)均在（i，j）點(diǎn)處取值.關(guān)于其一階和二階導(dǎo)數(shù)項(xiàng)的計(jì)算，均采用中心差分公式進(jìn)行計(jì)算.

由推導(dǎo)過(guò)程可知，該格式的截?cái)嗾`差為O（h4），即格式（27）在空間上可達(dá)四階精度.另外，注意到本文二維格式只用到9個(gè)網(wǎng)格點(diǎn)，故所得格式為緊致格式.

4 數(shù)值算例

為了驗(yàn)證本文一維和二維格式的精確性和有效性，現(xiàn)將以下有精確解的數(shù)值算例分別與MHOC格式［3］、Chen格式［5］、HOC格式［7］和FOC格式［8］的最大絕對(duì)誤差的數(shù)值結(jié)果進(jìn)行比較.其中，將一維和二維問(wèn)題的最大絕對(duì)誤差和收斂階定義如下：

式中的uc和ue分別代表數(shù)值解和精確解，h1和h2代表不同空間步長(zhǎng)，Error1和Error2分別對(duì)應(yīng)步長(zhǎng)為h1和h2時(shí)的最大絕對(duì)誤差.

問(wèn)題1［3］

該問(wèn)題的精確解為

問(wèn)題2［8］

該問(wèn)題的精確解為

表2 問(wèn)題2的最大絕對(duì)誤差和收斂階Tab.2 The maximum absolute error and convergence rate for Problem 2

問(wèn)題3［7］

該問(wèn)題的精確解為

表3 問(wèn)題3的最大絕對(duì)誤差和收斂階Tab.3 The maximum absolute error and convergence rate for Problem 3

問(wèn)題4

該問(wèn)題的精確解為

表4 問(wèn)題4的最大絕對(duì)誤差和收斂階Tab.4 The maximum absolute error and convergence rate for Problem 4

表1～4列出了針對(duì)一維問(wèn)題和二維問(wèn)題1～4的4種格式在不同h下的最大絕對(duì)誤差和收斂階.不難發(fā)現(xiàn)，盡管以上4種格式的精度基本上都達(dá)到了四階，但是本文計(jì)算所得的最大絕對(duì)誤差要更小些，即本文格式的計(jì)算結(jié)果更接近于精確解，也具有更高的計(jì)算精度.需要說(shuō)明的是，對(duì)于MHOC格式［3］、Chen格式［5］、HOC格式［7］和FOC格式［8］的計(jì)算結(jié)果，采用了Phoebe Solver軟件進(jìn)行計(jì)算得到.Phoebe Solver軟件為Web版，網(wǎng)址為www.phoebesolver.com.對(duì)于問(wèn)題4沒(méi)有采用Chen格式［5］和HOC格式［7］計(jì)算是因?yàn)樵搯?wèn)題是對(duì)流擴(kuò)散反應(yīng)方程，而這2篇文獻(xiàn)中的方程模型為對(duì)流擴(kuò)散方程，不包括反應(yīng)項(xiàng).因此，這2篇文獻(xiàn)所提格式不適用于該問(wèn)題的求解.

表1 問(wèn)題1的最大絕對(duì)誤差和收斂階Tab.1 The maximum absolute error and convergence rate for Problem 1

5 結(jié)論

本文首先針對(duì)方程（1）提出了一種新的緊致差分格式，具有四階精度，即在空間上采取泰勒級(jí)數(shù)展開法和截?cái)嗾`差余項(xiàng)修正法，先利用一階導(dǎo)數(shù)項(xiàng)來(lái)表示由原方程得到的未知函數(shù)的二階導(dǎo)數(shù)項(xiàng)，再表示出誤差余項(xiàng)中的三階和四階導(dǎo)數(shù)項(xiàng)，最終得到了一個(gè)三點(diǎn)四階的緊致差分格式，格式對(duì)應(yīng)的三對(duì)角線型方程組可采用追趕法進(jìn)行求解.接下來(lái)，對(duì)本文格式的耗散誤差和色散誤差進(jìn)行了分析，與文獻(xiàn)中格式的對(duì)比結(jié)果表明本文格式具有更小的耗散誤差.然后將該方法推廣到二維，得到了求解二維對(duì)流擴(kuò)散反應(yīng)方程的一種具有四階精度的緊致差分格式.最后給出了數(shù)值實(shí)驗(yàn)，通過(guò)與文獻(xiàn)中格式數(shù)值結(jié)果的比較發(fā)現(xiàn)本文格式具有更小的計(jì)算誤差和更高的計(jì)算精度.

致謝寧夏自治區(qū)重點(diǎn)研發(fā)項(xiàng)目（2018BEE03007）和寧夏大學(xué)研究生創(chuàng)新項(xiàng)目（GIP2019010）對(duì)本文給予了資助，謹(jǐn)致謝意.