等离子体物理中分数阶模型数值解法研究进展

梅立泉, 郭士民

(西安交通大学数学与统计学院，西安 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) 未形成成熟的数值计算软件，严重滞后于应用的需求。

目前已有的工作主要针对一维或二维标量分数阶微分方程开展数值求解，但能更好描述真实世界问题的三维分数阶非线性耦合方程组的数值求解却比较少。另外，分数阶方程的高效算法设计、算法的误差估计、带有非光滑初值的分数阶问题、分布阶问题等都是值得研究的方向。相比于大型设备的高投入，数值模拟是性价比非常高的研究手段。因此，进一步完善和发展等离子体物理中分数阶微分方程数值解法是非常适合我们这样一个发展中的新兴国家着力发展的方向，未来前景可期。