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

李智杰，王 妍

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

0 引 言

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

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

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

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

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

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

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

1 方程解的漸近行為

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

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

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

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

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

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

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

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

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

4 算例

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 語(yǔ)

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