等離子體物理中分數階模型數值解法研究進展

梅立泉, 郭士民

(西安交通大學數學與統計學院，西安 710049)

0 引言

等離子體是指由帶電粒子(電子、離子、帶電荷的塵埃等)和中性粒子組成的、具有集體行為的準中性氣體團。恒星、星際介質、核聚變裝置中的電離氣體、黑洞吸積盤、地球磁層等物質均屬于等離子體的范疇。宇宙中超過99%的可見物質由等離子體組成，包括作為天體物理和空間物理主要研究對象的各種尺度下的天體，如恒星、星際介質、星系、星系團等等。當前空間工程應用的重點區域–近地外太空都是等離子體(一般地球表層60 公里到1 000 公里的大氣層都是等離子體，像神九飛船飛行的高度為220 公里到330 公里)。所以，等離子體物理在空間科學研究和空間工程應用中具有非常重要的地位。天體物理以及空間物理中觀測到的大量現象都是由等離子體的原理和規律描述的。這些現象包括行星形成、太陽耀斑、恒星形成、太陽以及恒星星風、粒子加速、黑洞吸積與噴流、伽瑪射線暴、黑洞撕裂恒星事件、超新星爆發、星系形成與演化、宇宙大尺度結構等等。目前，在這些現象的數值模擬、半解析研究中，大多數工作都不得不只專注于大尺度的行為，但現象本身是多尺度的。因此，當前最急需的工作是仔細研究這些天文現象更具體的物理過程，研究其更精細的數學模型，并開展模型求解的數值方法研究和數值模擬工作。磁流體動力學模型是科學研究與工程應用中的一個基本數學模型，在航空航天、天體物理以及受控核聚變中有重要的應用前景。

目前等離子體物理中的數值模擬基本都是基于整數階的動力學方程[1–4]。但是，這些整數階方程只是等離子體物理中實際問題的一種近似模型，在考慮了反常擴散之后，這些問題完全用整數階方程刻畫是不夠準確的，需要考慮分數階模型。

近幾十年來，分數階微分方程逐漸引入并越來越多地應用于反常擴散、粘彈性材料、信號與圖像處理、系統識別、石油滲流、管道的邊界層效應、金融以及分形理論等應用領域。以統計物理學中的連續時間隨機行走模型為例，由于時間和空間的非局域性，模型中粒子束的傳播速度并不符合經典的Brown 運動理論，不滿足Fick 定律而表現為反常擴散行為，從而相應的數學模型也有別于經典的擴散方程。事實上，許多復雜動力系統也都包含反常擴散。由于分數階微積分具有的歷史依賴與非局部的特性比較適合刻畫反常擴散中的記憶性和非局部性質，因而分數階方程比整數階方程更能有效地描述這些復雜系統。分數階微分方程已經引起了人們廣泛的關注，逐漸成為一個新的活躍的研究領域。隨著涉及的應用領域越來越多，分數階微積分方程的研究在理論分析和數值模擬等多方面都迎來了很多新的挑戰。

分數階微積分在天體物理領域的物理意義由Podlubny 教授2008 年在牛津大學做訪問學者時給出，其中分數階微積分中所涉及的積分變量(時間)與霍金“時間簡史”中的時間理論有著非常重要的聯系(宇宙的時間是非均勻的，即宇宙大爆炸開始逐漸變慢的)。在等離子體中，帶電粒子的擴散過程直接作用于流體或磁流體的應力張量和粘滯效應。因此，帶電粒子的擴散過程能夠對等離子體流動行為產生十分重要的影響。在高溫、高壓、磁約束、強耦合等復雜等離子體物理環境中，由于狀態的非平衡性、運動的各向異性、時空分布的非局部性、速度空間的不均勻性、相互作用的長程性等諸多因素的影響，帶電粒子不再滿足布朗運動規律而是呈現出反常擴散行為，即帶電粒子的均方位移與時間之間的線性關系不再成立。反常擴散的一個重要特征是某時刻通過某點的流量不僅與其鄰域內的物理量有關，而且與整個空間中其它點處的物理量以及變化的歷史有關，從而呈現出非局部性質。此時，帶電粒子具有L′evy 飛行特性[5]，即帶電粒子的運動軌跡服從L′evy 分布，等離子體中的能量耗散具有冪律拖尾現象[6]。此類等離子體流體動力學演化過程具有明顯的歷史依賴性與長程相關性，主要表現為一定程度上的遲滯和松弛現象[7]。針對此種情況，經典的整數階梯度型定律不再成立，所以整數階的方程不能準確地刻畫基于反常擴散過程的流體動力學行為。而分數階微積分具有非局部性質，這使得其成為描述反常擴散過程強有力的數學工具，分數階方程可以精細地刻畫等離子體中基于反常擴散過程的流體動力學行為的歷史依賴性與長程相關性(非局部性質)。因此，將分數階模型應用到等離子體物理領域將有重要的意義和前景。

1 分數階微分方程數值解法主要研究進展

分數階微積分作為數學的一個分支，它和經典微積分一樣歷史悠久。早在300 多年前，數學家Leibniz 和L’Hospital 就以書信的形式研究過分數階導數。在被提出至今的三百多年里，最初在物理領域并未獲得廣泛關注與應用，發展非常緩慢，僅僅作為數學領域中的純理論問題被諸多學者研究。歷史上對分數階微分方程和分數階微積分理論做出過重要貢獻的數學家還包括Riemann、Euler、Laplace、Liouville、Abel、Fourier 等。近年來，隨著對物理現象認識程度的加深以及現代計算機運算能力的提高，物理、化學、生物等多個學科領域的分數階微積分建模越來越引起人們的重視[8–13]。

與整數階微分方程相比，分數階微分方程具有獨特的性質和優勢。這主要體現在以下兩個方面。

1) 從理論角度而言，由于分數階微分算子是通過積分的形式定義的，這表明函數在某點處的分數階導數不但與該點處的性質有關，而且與某個區域上的整體性質有關。因此，相對于整數階微分算子僅具有的局部性質，分數階微分算子能夠體現函數變化的全局相關性。當然，分數階微分算子是具有弱奇異核的擬微分算子且不滿足半群性質、交換律等性質。因此，分數階微分方程絕非整數階微分方程的簡單推廣。

2) 從應用角度而言，經典的整數階微分方程僅僅能夠體現物理過程在某個時刻或者空間某點處的變化和性質，而分數階微分方程所刻畫的性質則與該物理過程所依賴的整個發展歷史或者所涉及的整個空間有關。因此，在實際應用中，分數階微分方程可以更加精確地描述具有歷史依賴性(時間分數階問題)和長程相關性(空間分數階問題)的物理過程，能夠克服許多整數階模型與實驗結果不吻合的缺點。

當從實際問題中抽象出分數階模型后，一個亟待解決的問題是如何對此類數學模型進行求解。然而，由于分數階微分算子具有非局部性質和弱奇異性質，分數階微分方程的求解具有很大的難度。一般來說，當前主要有兩類方法來求解分數階微分方程。第一類是解析求解方法，如Fourier 變換法、Laplace 變換法、Mellin 變換法、Green 函數法、分離變量法、算子法、Adomain 分解法、變分迭代法以及同倫分析法等。但是，由于分數階算子的非局部性質，可以求出的解析解中會包含如Fox 函數、超幾何函數、Mittag-Leffler函數、Wright 函數等形式復雜的特殊函數，此時解析解形式將過于復雜而難以開展實際應用。另一方面，對于非線性、高維、變系數等情形比較復雜的分數階微分方程，人們很難構造其解析解。第二類是數值求解方法，如有限差分法、有限元方法、譜方法、有限體積法等。

在時間分數階偏微分方程的數值求解方面，Sanz-Serna[14]考慮一類非線性偏積分微分方程，得到了關于時間的半離散格式，進而證明了格式的收斂精度為1 階。Gorenflo 等[15]利用離散隨機游走模型對時間分數階擴散方程進行了研究，但是缺乏相關的理論分析。隨后，Liu 等[16]在Gorenflo 的研究基礎上給出了離散非馬爾可夫模型的穩定性和收斂性分析。Sun 和Wu[17]采用有限差分法構造了兩類含有不同階Caputo 分數階導數的分數階次擴散方程(0α＜ 1)和分數階擴散波動方程(1α＜ 2)的離散格式，并證明兩種格式在時間方向上的精度分別為2?α和3?α階，國際上將其稱之為L1 公式。Lin 和Xu[18]利用有限差分/譜方法數值求解了Caputo 分數階導數意義下的時間分數階擴散方程，證明了數值格式的無條件穩定性，并且給出了數值格式的收斂階。Gao 等[19]提出了求解時間分數階次擴散方程的有限差分格式，被稱為L1?2 公式，收斂精度為3?α。隨后，Alikhanov[20]使用超收斂點構造了時間方向收斂階為3?α的有限差分格式，并給出了詳細的穩定性和收斂性證明，稱為L2?1σ格式。

在空間分數階偏微分方程的數值求解方面，Liu 等[21]提出了求解空間分數階Fokker-Plank 方程的分數階行方法，該方法將分數階偏微分方程轉化為常微分方程系統。然后，利用BDF 方法進行求解。Ervin 和Roop[22]定義了一類分數階函數空間，證明了這種分數階空間與分數階Sobolev 空間在一定條件下是等價的，建立了空間分數階對流–擴散方程的Galerkin 變分形式，利用有限元方法對該方程進行了求解。2014 年，Xu 和Hesthaven[23]利用間斷有限元方法求解了空間分數階擴散方程和空間分數階對流–擴散方程，建立了數值格式的誤差估計，并且通過數值算例驗證了理論分析。Wang 等[24]利用有限元方法數值求解了帶有非齊次Dirichlet 邊界條件的空間分數階擴散方程，他們將標準有限元方法與半解析方法相結合，使得提出的算法對計算機存儲需求由O(N2)降低到O(N)，其中N為離散系統的階數。Sousa[25]采用有限差分方法求解了帶有源項的一維空間分數階對流擴散方程，設計了具有二階收斂性質的有限差分格式，并且分析了Riemann-Liouville 分數階導數的階數對數值格式穩定性的影響。對非規則區域，Yang 等[26]考慮任意非規則區域，給出了非線性Riesz 空間分數階擴散方程的有限元格式，討論了格式的穩定性及收斂性。Lee[27]提出了基于算子分裂的空間分數階反應–擴散方程的Fourier 譜方法，并對Allen-Cahn 模型、FitzHugh-Nagumo 模型、Gray-Scott 模型開展了數值模擬工作。2019 年在文獻[28]中，我們建立了求解空間分數階Cahn-Hilliard 和Allen-Cahn 方程能量穩定2 階時間精度的Fourier 譜方法。

在時間–空間分數階偏微分方程的數值求解方面，Deng[29]構造了時空分數階Fokker-Planck 方程的半離散及全離散有限元格式，并詳細分析了這些格式的適定性、穩定性及收斂性。Li 等[30]考慮了非線性時空分數階亞擴散及超擴散方程的有限元格式，并給出了誤差分析及其數值模擬。Yang 等[31]研究了二維時間–空間分數階擴散方程的數值算法，該方法首先對方程進行矩陣變換，得到一組分數階微分系統，然后利用有限差分方法、有限元方法對該系統進行數值求解。Li 等[32]分析了時間分數階擴散–波動方程的Crank-Nicolson ADI 有限元格式。Zayernour 和Karniadakis[33]提出了求解時間–空間分數階對流方程的間斷譜元方法，通過數值算例驗證了算法的譜收斂性。Hanert 和Piret[34]通過Chebyshev 擬譜方法求解了含有Caputo 分數階導數的時間–空間分數階擴散方程，數值算例表明：當模型方程的解具有光滑性質時，該數值方法能夠達到指數收斂。為了數值模擬二維問題，Li 等[35]研究了時間–空間分數階相場模型(時間–空間分數階Allen-Cahn 方程)，并且建立了方程的時間–空間全離散格式：在時間方向上，利用有限差分格式對Caputo 分數階導數進行離散；在空間方向上，利用配置方法對L′evy 過程的無窮小生成元進行離散。文獻[36]構造了求解時間–空間分數階次擴散和超擴散方程的局部間斷Galerkin 方法，證明了方法的穩定性和收斂性。文獻[37]建立了求解二維分布階時間–空間分數階反應–擴散方程的有限差分/Legendre-Galerkin 譜方法。文獻[38]構造了求解三維時間–空間分數階擴散方程的ADI 譜Galerkin 方法，證明了方法的無條件穩定性和空間最優誤差估計。

2 分數階色散方程數值解法研究進展

色散現象是指波包中不同頻率的波以不同的速度傳播，即波的傳播速度依賴于頻率，隨著時間的演化這些波通常會分散開。具有這樣性質的未知函數的方程，稱為色散方程。非線性色散偏微分方程是一類非常重要的非線性發展方程。它在應用數學和物理領域中占有重要的地位，如流體力學、固體物理、量子力學、等離子體物理、非線性光學等，在化學和生物中也有廣泛的應用。

數學物理中有許多非線性色散和波動方程，主要有非線性波動方程、Maxwell-Klein-Gordon 方程、Yang-Mills 方程、正則長波方程、非線性Schr¨odinger 方程、Korteweg-de Vries(KdV)類型的方程(修正的KdV、廣義的KdV 方程等)和一些系統(Kadomtsev-Petviashvili 方程、Davey-Stewartson 方程等)。這類方程一般具有一些共同特征，如可用散射反演化方法求解、存在達布變換、具有多個守恒律和延長結構。

非線性色散偏微分方程中Schr¨odinger 方程作為量子力學中的最基本的物理方程描述了物理系統的量子態怎樣隨時間演化的偏微分方程，奠定了近代量子力學的基礎。它用于描述諸如Bose-Einstein 凝聚的多體理論和凝聚態物理學。Schr¨odinger 波動方程在整個量子力學張占有重要的地位。它在流體力學、光學、化學、電磁學尤其光纖通信中有廣泛的應用。耦合的Schr¨odinger 方程是Schr¨odinger 方程的向量形式，可以用于描述許多物理現象，如沿著正交偏振軸的脈沖傳播、兩組分的玻色愛因斯坦和怪波等。在非相對論極限狀態下，耦合的Schr¨odinger 方程可用于描述在具有低頻離子響應的高頻電子等離子體波中傳播的Langmuir 包絡孤子的動力學行為。在非線性光學中，耦合的Schr¨odinger 方程可以用于描述至少在兩個通道中同時模擬多模孤子脈沖的傳播。

經典的Schr¨odinger 方程是基于自由粒子的Feynman 路徑積分滿足布朗運動的假設下得到的，但是在現實中，很多物理現象并不滿足該假設。Laskin[39]提出了用L′evy 路徑積分來替換Feynman 路徑積分，從而得到帶有Riesz 空間分數階導數的Schr¨odinger 方程并開啟了分數階量子力學的大門。隨后，Hu 和Kallianpur[40]提出了帶有分數階拉普拉斯算子的Schr¨odinger 方程并給出其概率形式的解。Muslih 等[41]利用分數變分原理導出了時間分數階Schr¨odinger 方程，Wang 和Xu[42]在分數路徑積分和分數布朗運動的基礎上發展了一些廣義分數階Schr¨odinger 方程。

Zhao 等[43]針對Riesz 分數階導數建立了一類具有四階精度的緊致差分算法，并將所建立的算法應用到二維空間分數階非線性Schr¨odinger 方程的數值求解，證明了數值格式的穩定性和收斂性。同年，Wang 等[44]利用有限差分方法數值求解了空間分數階耦合非線性Schr¨odinger 方程，證明了離散系統解的存在性和唯一性，并且分析了隱式格式的收斂性。Duo 和Zhang[45]結合時間分裂方法、Crank-Nicolson 方法和松弛方法構造了三種質量守恒的Fourier 譜方法來求解空間分數階的Schr¨odinger 方程。文獻[46]利用有限差分/Legendre-Galerkin 譜方法對二維時間分數階非線性擴散–波動方程進行了數值求解，證明了數值格式的穩定性和收斂性，并且對Sine-Gordon 方程的環形孤立子開展數值模擬。文獻[47]利用有限差分/Hermite-Galerkin 譜方法對無界區域上的多維時間分數階非線性反應擴散方程進行了數值求解，證明了數值格式的無條件穩定性，并將方法應用到了分數階Allen-Cahn 和Gray-Scott 模型的數值模擬。文獻[48–49]分別構造了求解二維非線性空間分數階Schr¨oinger 方程和耦合非線性空間分數階Schr¨odinger 方程的譜Galerkin 方法，證明了方法的穩定性和空間最優的誤差估計。文獻[50]構造了空間分數階Klein-Gordon-Schr¨odinger 方程IEQ-Crank-Nicolson-Fourier 譜方法，該算法最大的特點是所有引入的輔助變量都是半顯處理的，算法能嚴格保證能量守恒性質，并且是對稱正定的線性格式，可以用預處理的共軛梯度法來加速求解。

3 分數階磁流體方程數值解法研究進展

當從宏觀角度研究等離子體在電磁場(如外磁場、外電流磁場、電磁擾動等)中的運動時，必須研究電流體運動的流體力學方程和刻畫電磁場運動的Maxwell 方程耦合在一起的磁流體動力學(Magnetohydrodynamics, MHD)方程。

在許多等離子體環境中，磁約束、強耦合、超高溫、超高壓等物理條件會產生相互作用的長程性、時空分布的非局部性、粒子運動的各向異性、狀態的非平衡性、速度空間的不均勻性等諸多復雜因素，這些復雜因素會使得帶電粒子不再滿足布朗運動規律而是具備反常擴散的特性[51]。在此類等離子體中，磁流體動力學演化過程具有明顯的非局部性質[52]，即歷史依賴性與長程相關性。為了彌補經典的MHD 方程的不足之處，人們將非局部算子引入到磁流體動力學并建立了非局部MHD 方程(或稱為廣義MHD 方程)[53]。非局部MHD 方程可以精細地刻畫等離子體中基于反常擴散過程的磁流體動力學行為的非局部性質，具有十分重要的研究意義。

對分數階MHD 方程，Zhang 等[54]利用有限差分方法求解了含有Caputo 分數階導數的2 維非定常時間分數階MHD 方程，證明了數值格式的無條件穩定性，并且建立了數值格式的誤差估計；在數值模擬部分，討論了分數階導數、Hartmann 數、壓力梯度參數等物理量與速度場之間的關系。Zhao 等[55]求解了含有齊次Dirichlet 邊界條件的二維時間–空間非局部MHD 方程，該論文利用L1 插值方法對時間方向上的非局部算子進行離散，空間變量利用有限差分方法進行離散，非線性項通過隱式格式進行處理，通過對非局部流動現象開展了數值模擬。對2 維不可壓縮非定常時間分數階MHD 方程，Bai 等[56]分別通過L1 插值逼近和有限差分格式離散時間方向的Caputo 分數階導數和空間變量建立數值格式，通過數值算例驗證了數值格式的收斂性，并討論了分數階導數的階數、Hartmann 數、雷諾數等物理參數對速度場、壓力的作用規律。Rasheed 與Anwar[57]構造了求解時間分數階MHD 方程的有限差分/有限元格式，時間方向的Caputo 分數階導數采用有限差分進行離散，空間變量采用有限元進行逼近；對方程中的非線性項采用隱式格式，在每個時間步上問題被離散為一個非線性代數方程組。針對此代數方程組，文中通過Newton 迭代法進行求解。文獻[58]建立了基于時間分裂方法的2 維不可壓縮分數階MHD 方程的全離散數值格式，該數值格式保證磁場在離散層面上也是無散度的，并對頂蓋驅動方腔流開展數值模擬。該文獻還研究了分數階Laplace 算子的階數對速度場、磁場等物理量的影響。

4 總結與展望

等離子體物理在空間科學研究和空間工程應用中具有非常重要的地位。目前等離子體物理中的數值模擬絕大多數都是基于整數階的動力學方程，經典的整數階微分方程不能精確地描述在可控核聚變裝置、大型等離子體研究裝置、強湍流等離子體、磁約束等離子體、空間等離子體中發現的基于反常擴散的非局部磁流體動力學行為。針對具有非局部性質的等離子體，分數階微分方程可以更加精確地描述具有歷史依賴性質(時間分數階問題)和長程相關性質(空間分數階問題)的物理過程。

本文綜述了作者所熟悉的等離子體物理中分數階微分方程數值模擬方面的研究進展。等離子體物理中分數階微分方程數值方法的研究已逐步展開，但相對于整數階模型，對分數階模型的研究仍處于初級階段。分數階偏微分方程的數值求解存在以下三個難點。

1) 非局部性：在分數階偏微分方程中，非局部性質是通過非局部算子的積分定義來體現的。在數值求解時間非局部模型時，非局部性質將導致當前時刻的數值解依賴于前面所有時刻的數值結果；在數值求解空間非局部模型時，非局部性質會使得離散后的系數矩陣為稠密矩陣。因此，求解非局部MHD 方程的數值方法具有計算復雜度大的難點。

2) 非線性：非局部MHD 方程具有內稟的非線性特性，主要體現在流體力學方程中的非線性對流項和洛倫茲力項，以及磁感應方程中的非線性對流項。因此，非局部MHD 方程具有較多的非線性項。在數值求解中，對這些非線性項的處理是一個研究難點。

3) 多物理場耦合性：無論是分數階色散方程，還是非局部MHD 方程，往往都是多個方程耦合在一起而構成的方程組，在計算過程中涉及到多個不同的物理效應和物理量同時求解。由于這些不同的物理效應和物理量所對應的算子各具特點，在數值求解這些物理效應和物理量時所采取的方法和策略也是不同的。這對構造等離子體物理中分數階微分方程高效穩定的數值算法帶來了很大的困難。

總的來講，對等離子體物理中分數階微分方程數值解法的研究仍不成熟，主要有：

1) 長時間歷程的計算和大空間域的計算等挑戰性難題仍未很好的解決；

2) 未形成成熟的數值計算軟件，嚴重滯后于應用的需求。

目前已有的工作主要針對一維或二維標量分數階微分方程開展數值求解，但能更好描述真實世界問題的三維分數階非線性耦合方程組的數值求解卻比較少。另外，分數階方程的高效算法設計、算法的誤差估計、帶有非光滑初值的分數階問題、分布階問題等都是值得研究的方向。相比于大型設備的高投入，數值模擬是性價比非常高的研究手段。因此，進一步完善和發展等離子體物理中分數階微分方程數值解法是非常適合我們這樣一個發展中的新興國家著力發展的方向，未來前景可期。