非線性對流反應(yīng)擴散方程爆破解的B方法

李智杰，王 妍

（西安建筑科技大學(xué) 理學(xué)院，陜西 西安 710055）

0 引 言

考慮如下完全非線性反應(yīng)擴散方程的爆破問題，其中α,ε,λ ＞0 ，p＞m≥1，L為正常數(shù)，初值u0(x)滿足

完全非線性系統(tǒng)作為一種重要的物理模型，常用于描述復(fù)雜介質(zhì)的流動問題，如超變形核、液滴裂變以及慣性聚變等［1］。1950 年，Zeldovich 和Kompaneets 提出了一種非線性多孔介質(zhì)方程ut=α(um)xx，用于描述絕熱氣體在多孔介質(zhì)中的流動過程。例如，當m=2 時，對應(yīng)水文學(xué)中的Boussinesq方程的無量綱重構(gòu)［2］。此后，含有源項λup及對流項ε(un)x的多孔介質(zhì)方程相繼被提出。例如，液滴模型中的B(n,n)方程即為后者在n=m時的特殊形式［3-4］，而當n≠m時，則通常表示泡沫排水方程［5］或斜板上薄黏性液體運動的過程［6］等。

此后，由于研究實際物理模型的需要，一類能夠描述復(fù)雜多孔介質(zhì)間作用系統(tǒng)的方程

開始被關(guān)注［7-8］。該方程（2）的各個參數(shù)α、ε、λ、m、n以及p在不同組合下，可以涵蓋大量非線性拋物問題相關(guān)的數(shù)學(xué)模型，并且該方程的一些退化形式已有大量研究[9-14]。此外，由于當p＜m時方程（2）已被證明其所有解都一致有界［15］，因此，目前的研究大都集中于p＞m的情形。

處理非線性奇異問題的數(shù)值方法有局部加密、多步格式［16］、自適應(yīng)移動網(wǎng)格［17］等。由于爆破點附近解趨于無窮大，一些傳統(tǒng)的數(shù)值方法不能精確地模擬爆破點附近解的漸近行為，因而估計的爆破時間會偏離真正的爆破時間，尤其是考慮到問題（1）中還存在對流項對解產(chǎn)生的復(fù)雜影響［18］。因此，需要一種能夠精確模擬爆破點附近解的漸近行為，同時能降低運算復(fù)雜度的數(shù)值方法。B方法作為一種兼顧以上需求的時間半離散數(shù)值方法，其構(gòu)造思想由LE ROUX 于1992 年首先提出。在實際問題中，許多爆破模型通常由一個非線性源項和一個擴散項2部分組成，而且非線性項是造成爆破的主要原因，擴散項只是延遲了爆破發(fā)生的時間。前者對解的擾動較后者占優(yōu)，且越臨近爆破發(fā)生時越明顯。因此，為使數(shù)值解能夠在追蹤爆破時間時具有更高的精度，LE ROUX 等首先將方程變?yōu)闆]有空間算子的演化方程的形式，即僅考慮對爆破過程起重要作用的非線性常微分方程，并求出該方程的精確解。之后，再對這一非線性常微分方程進行離散并將空間算子回填，最終得到一種能夠保持解的幾何結(jié)構(gòu)的數(shù)值方法?？紤]到該方法對爆破解具有極好的逼近特性，故以爆破（blow-up）的英文首字母命名，稱作B方法。

LE ROUX等于1994年完善了B方法的相關(guān)理論之后，利用B方法在1998年和1999年研究等離子體物理學(xué)中一類帶有高階源項的非線性擴散方程的爆破問題［19，20］，給出了解在有限時間內(nèi)爆破的理論依據(jù)，證明了數(shù)值解的漸近性與收斂性，并于2000 年驗證了B 方法在該方程高維情況下的適用性［21］。2009 年，LE ROUX 又利用B 方法研究方程組解的漸近行為［22］，得到了爆破時間的估計以及數(shù)值爆破解的收斂性。2015年，BECK等進一步發(fā)展了LE ROUX 的理論，使用B 方法研究擬線性拋物方程（組）及半線性波動方程的爆破問題［23］。2018年，霍冠澤等將B方法用在處理一類拋物方程的猝滅問題上［24］，證明了數(shù)值解的存在性并通過算例驗證了B 方法在求解方程猝滅解時的有效性。2020年，霍冠澤等又使用B 方法研究了一類四階拋物方程的爆破問題［25］，展示B 方法在處理該類問題時的優(yōu)勢。

由此可見，利用B 方法研究微分方程的奇異問題，尤其是爆破問題時確實能夠得到非常好的效果。不過，以上問題都圍繞只包含拉普拉斯算子及非線性反應(yīng)項的方程（組）展開，對含有對流項的問題研究較少；而且含有對流項的方程沒有變分結(jié)構(gòu)，故之前文獻所使用的能量泛函方法也將不再適用。本文通過構(gòu)造輔助函數(shù)方法證明問題（1）的解在有限時間內(nèi)爆破理論，即通過對應(yīng)橢圓方程的特征函數(shù)定義新的檢驗函數(shù)，并證明其滿足一個常微分方程不等式。本文中還給出了問題（1）B方法數(shù)值格式的具體構(gòu)造過程及其收斂性、數(shù)值解存在性的證明和爆破時間估計。通過數(shù)值算例驗證相關(guān)結(jié)果的正確性，并將B 方法與3 種傳統(tǒng)數(shù)值格式在不同區(qū)間長度下的模擬結(jié)果相比較，展示了B方法在處理相關(guān)問題時的優(yōu)勢。同時，算例也刻畫了對流項及對流系數(shù)對于數(shù)值解爆破過程的影響。

1 方程解的漸近行為

首先給出問題（1）的解能夠發(fā)生爆破的充分條件。在區(qū)域[-L,L]引入Laplace 算子的線性特征值問題

成立，并由此可知問題（1）的解會在有限時間內(nèi)爆破。

2 B方法數(shù)值格式及其數(shù)值解性質(zhì)

定理2［27］若帶有Dirichlet邊值條件的橢圓問題Lv=f(x,v,?v)滿足：

（i）非線性項f(x,v,?v)關(guān)于?v至多是二次增長的，即存在函數(shù)h(x,v)關(guān)于x∈[-L,L]可測，關(guān)于v連續(xù)，以及函數(shù)H(x,v,vx)關(guān)于其自變量均連續(xù)，使得f能夠被表示為

3 構(gòu)造全離散數(shù)值格式

由于隱格式（7）的全離散格式難以直接計算，因此需要尋找它的一個逼近序列，并證明在第n層給定迭代首項的情況下，該序列能夠經(jīng)過有限次迭代之后收斂到第n+1層。根據(jù)文獻［28］的結(jié)論，若

4）生成采樣坐標對，繪制函數(shù)圖像。

利用B方法及有限差分方法構(gòu)造了全離散數(shù)值格式。值得注意的是，在選取空間差分方式時有多種選擇，此處選擇有限差分方法，是由于其所對應(yīng)的迭代矩陣相對簡潔，在利用計算機進行數(shù)值模擬時可以大幅縮減計算時間。另外，可以根據(jù)不同的需求對空間離散方法進行調(diào)整，以滿足不同類型的實際模型。

4 算例

以一維非線性對流反應(yīng)擴散方程為例，采用所構(gòu)造的格式進行模擬，并與3 種傳統(tǒng)數(shù)值格式的計算結(jié)果比較，驗證本文數(shù)值格式在處理爆破問題時的精確性及可靠性。

固 定 參 數(shù)α=0.1 ，ε=0.3 ，λ=2 ，m=3 ，p=6，討論隨著區(qū)域范圍的變化，4種格式相應(yīng)的數(shù)值模擬情況。令u(x,t)為待求未知量，得到方程及其初邊值條件如下：

如前所述，隨著時間推移，方程（17）的解將在有限時間內(nèi)出現(xiàn)爆破現(xiàn)象。選擇了4 種數(shù)值格式模擬問題（17）解的爆破過程，分別為B 方法、切比雪夫譜方法、有限差分法及4-5 階龍格庫塔法?？臻g步長統(tǒng)一取作Δx=0.01，并選取T=0.23 作為數(shù)值爆破時間進行觀察。這樣既能夠保證刻畫解在t時刻附近爆破現(xiàn)象，又不會因為過于靠近實際爆破時間而使得模擬失真。首先以L=1/2 為例，展示爆破發(fā)生前，4 種格式數(shù)值解隨時間變化的過程。其中有限差分法對應(yīng)數(shù)值解的爆破現(xiàn)象發(fā)生過早，此時將爆破時間取為T=4×10-4，即解出現(xiàn)奇性的前一刻。給定參數(shù)下不同算法的模擬結(jié)果如圖1所示。

圖1 給定參數(shù)下不同算法的模擬結(jié)果Fig.1 Simulation results of different algorithms under given parameters

從圖1 可以看到：在臨近爆破點時，B 方法依舊能夠保持解的光滑性，而切比雪夫譜方法對應(yīng)解則產(chǎn)生較大的震蕩；有限差分法在T時刻之后產(chǎn)生了奇性，以至于計算機在模擬時會由于部分結(jié)點震蕩過大而出現(xiàn)嚴重的網(wǎng)格間斷。4～5 階龍格庫塔法雖然能夠保證解的光滑性，但爆破過程中數(shù)值解出現(xiàn)的形變較為嚴重，爆破點位置也有較大偏差。可見，B方法在模擬數(shù)值解爆破過程的精確度方面更具優(yōu)勢。

對問題所在區(qū)域進行調(diào)整，分別取L=1/2 、π/4、π/2、2，并將4種算法在對應(yīng)爆破時刻數(shù)值解的形狀通過等比例壓縮至相同的坐標區(qū)間，以便更為直觀的比較區(qū)間長度的變化對各算法數(shù)值解帶來的影響，結(jié)果如圖2 所示。由圖2 可以看到：隨著區(qū)間長度的增加，B方法的空間結(jié)點分布依舊均勻，曲線平滑，同時方程爆破點位置的偏移現(xiàn)象也刻畫的較為準確；切比雪夫譜方法雖然獲得了與B 方法相近的爆破點，但由于數(shù)值解本身的劇烈震蕩，且在爆破發(fā)生時數(shù)值解不恒為正，使得數(shù)值解與理論結(jié)果不符；有限差分法由于誤差較大，導(dǎo)致僅僅在非常短的時間之內(nèi)，數(shù)值解便出現(xiàn)奇性且爆破位置偏差嚴重；4～5階龍格庫塔法數(shù)值解雖然光滑，但行波迎風面結(jié)點分布不均，爆破位置同樣偏差較大，刻畫不準確，并且這一現(xiàn)象隨著L的增大逐漸明顯。

圖2 不同區(qū)間長度下4種格式模擬結(jié)果Fig.2 Simulation results of four formats under different interval lengths

以B 方法在L=1/2 ，T=0.2 時刻的模擬結(jié)果為例。圖3 展示了取不同的對流項系數(shù)時，隨著ε增大，對流項所產(chǎn)生的阻尼效果越來越強的物理現(xiàn)象。當ε逐漸變大，對流效果愈加明顯的同時，同一時刻數(shù)值解的最大值也在相應(yīng)地減小?？梢?，ε對該問題解的影響與其物理意義保持一致。

圖3 不同對流系數(shù)對數(shù)值解的影響Fig.3 The influence of different convection coefficients on numerical solution

此外，B方法在利用計算機進行模擬時，相較于其他算法具有更高的性價比。表1 給出了上述4 種格式相應(yīng)程序的運行情況。表1中：運行時間為同一設(shè)備在相同狀態(tài)下，各程序運行10次所用的平均時間；函數(shù)調(diào)用次數(shù)為包括主函數(shù)、子函數(shù)及嵌套函數(shù)等在內(nèi)的總調(diào)用次數(shù)；網(wǎng)格密度為單位面積內(nèi)的平均網(wǎng)格數(shù)。實驗設(shè)備為一臺搭載2.6 GHz Intel Core i7-10750H處理器，內(nèi)存為16 GiB的筆記本電腦。

表1 各數(shù)值格式運行情況Tab.1 Operating information in various formats

通過表1可以看出：與B方法相比，切比雪夫譜方法雖然計算時間較短，但調(diào)用函數(shù)次數(shù)較多，網(wǎng)格密度大且產(chǎn)生震蕩現(xiàn)象，因此可以考慮用來確定爆破點的位置，但不適合模擬爆破問題的數(shù)值結(jié)果。有限差分法雖然計算時間最短，函數(shù)調(diào)用次數(shù)較少且網(wǎng)格密度小，但結(jié)合圖1、2 來看，其精度最低，誤差最大；4～5階龍格庫塔法雖然未產(chǎn)生震蕩且網(wǎng)格密度適中，但函數(shù)總調(diào)用次數(shù)過多，程序遍歷時間過長，運算效率低下；B 方法雖然計算時間略長，但函數(shù)總調(diào)用次數(shù)最少，即使在低性能計算機上運行也能夠保證計算效率，因此更適用于模擬此類含奇異解的非線性方程。

5 結(jié) 語

得益于B 方法對數(shù)值解奇異行為的控制，該方法在處理非線性對流擴散方程的爆破問題時有非常好的表現(xiàn)，使得其能夠在較短時間內(nèi)，準確地模擬解的爆破行為。同時，該方法用空間換時間，提升了對運行內(nèi)存的利用率，是目前利用計算機模擬非線性微分方程爆破問題有力的工具之一。但B方法作為一種時間上半離散格式，在構(gòu)造全離散格式的過程中，有時仍需要依賴誤差較大、精度較低的傳統(tǒng)差分方法離散空間算子，甚至在處理高維問題時需要進行較多的預(yù)處理工作。因此，以B 方法為基礎(chǔ)的數(shù)值格式在未來仍有諸多亟待解決的問題，以及較大的提升空間。