分数阶扩散方程的无网格数值模拟

Mesh-Less Numerical Simulation of the Fractional Diffusion Equation

摘要: 文章借助径向基无网格方法数值求解分数阶扩散方程，时间上使用有限差分方法离散时间导数，空间上分别选取Multi Quadrics (MQ)，Thin Plate Spline (TPS)和Cubic三种径向基函数近似未知函数，比较得出三种径向基函数的逼近精度类似，但Cubic径向基函数无须选择形状参数，数值结果验证了该方法的可行性和有效性。

Abstract: In this paper, the mesh-less method based on radial basis functions is used to solve the fractional diffusion equation numerically. In terms of time, finite difference method is used to discrete time derivatives. In terms of space, three approximate unknown functions of radial basis functions of Multi Quadrics (MQ), Thin Plate Spline (TPS) and Cubic are selected respectively. The approxima-tion accuracy of three kinds of radial basis function is similar, but Cubic radial basis function does not need to select shape parameter. Numerical results demonstrate the feasibility and effectiveness of the proposed method.

1. 引言

在分形介质中分子扩散现象不能用标准的扩散方程进行描述，称之为反常扩散 [1] 。近年来，该现象引起了人们的不断关注，在半导体、核磁共振、多孔介质、湍流、固体表面扩散及经济金融等研究中有着广泛应用。与正扩散相比，反常扩散表现出长期相互作用和历史依赖性的显著特征，可以借助分数阶扩散方程进行刻画。但分数阶导数具有非局域性，使求解变得更加困难，并且在绝大多数情况下分数阶方程无法得到解析解 [2] [3] ，这使得反常扩散问题的数值算法开发显得尤为重要 [4] 。

由文献可知，许多学者正致力于分数阶方程的数值方法研究。例如白鹭等 [5] 提出了一种求解时间分数阶偏微分方程的数值算法。余梓航 [6] 对不同阶数分数阶导数构造算子，考虑了一维及二维常系数空间分数阶扩散方程的有限差分方法。陈景华 [7] ，沈淑君 [8] ，任晶 [9] 对相关数值研究的收敛性、准确性和稳定性进行了广泛的讨论。宋灵宇 [10] 等采用Modified Kansa法计算椭圆型偏微分方程—双调和方程的数值解。陈文 [11] 等首次应用Kansa方法求解分数阶扩散方程，选用Multi Quadrics(MQ)作为基函数成功求解了一维分数阶扩散方程，但求解精度极大依赖于形状参数的选取，为改变这一不足，考虑重新选取不同类型的径向基函数进行比较。

在给出定解问题的基础上，使用有限差分方法离散Caputo时间分数阶导数，分别选取MQ，Thin Plate Spline (TPS)和Cubic三种径向基函数用于近似未知函数，比较不同基函数的优缺点，数值结果表明，Cubic基函数在达到相同精度的情况下可以克服形状参数选取的限制。

2. 定解问题

我们考虑时间分数阶扩散方程

${\partial }^{\alpha }u/\partial t^{\alpha }+u=k{\nabla }^{2}u+f\left(x,t\right),0<\alpha <1,x\in \Omega ,t>0$ (1)

边界条件为

$u\left(x,t\right)=g\left(x,t\right),x\in \partial \Omega ,t>0$ (2)

初始条件为

$u\left(x,0\right)=w\left(x\right),t=0$ (3)

其中 $u\left(x,t\right)$ 表示溶质浓度， $f\left(x,t\right)$ 表示原象， $g\left(x,t\right)$ 表示边界溶质浓度， $w\left(x\right)$ 表示初始溶质浓度， $\alpha$ 为时间导数的阶数， ${\partial }^{\alpha }u\left(x,t\right)/\partial t^{\alpha }$ 卡普托分数阶导数，k为扩散系数， $x=\left[x,y,z\right]$ ， ${\nabla }^2$ 表示拉普拉斯算子。 $\Omega$ 表示有界域， $\partial \Omega$ 是它的边界， $\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}$ 是如下形式的Caputo分数阶导数

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\{\begin{array}{l}\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{\partial u\left(x,\xi \right)}{\partial \xi }\frac{1}{{\left(t-\xi \right)}^{\alpha }}\text{d}\xi ,}0<\alpha <1\\ \frac{\partial u\left(x,\xi \right)}{\partial \xi },\alpha =1\end{array}$ (4)

3. 数值方法

在离散化时间的基础上，借助有限差分和径向基函数得到求解问题的数值格式。

3.1. 时间分数阶导数的离散化

假设求解时间为T，记时间步长为 $\tau$ ，则求解时间离散为 $T=N\tau$ ，第 $n+1$ 步 $t_{n+1}=\left(n+1\right)\tau$ 时，离散积分区间得到

$\begin{array}{c}\frac{{\partial }^{\alpha }u\left(x,t_{n+1}\right)}{\partial t^{\alpha }}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle \underset{0}{\overset{t_{n+1}}{\int }}\frac{\partial u\left(x,\xi \right)}{\partial \xi }\frac{1}{{\left(t_{n+1}-\xi \right)}^{\alpha }}\text{d}\xi }\\ =\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle \underset{k=0}{\overset{n}{\sum }}{\displaystyle \underset{k\tau }{\overset{\left(k+1\right)\tau }{\int }}\frac{\partial u\left(x,\xi \right)}{\partial \xi }}}\frac{\text{d}\xi }{{\left(t_{k+1}-\xi \right)}^{\alpha }}\\ \approx \frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle \underset{k=0}{\overset{n}{\sum }}{\displaystyle \underset{k\tau }{\overset{\left(k+1\right)\tau }{\int }}\frac{\partial u\left(x,{\xi }_k\right)}{\partial \xi }}}\cdot \frac{\text{d}\xi }{{\left(t_{k+1}-\xi \right)}^{\alpha }}\end{array}$ (5)

借助有限差分方法离散时间导数得

$\frac{\partial u\left(x,{\xi }_k\right)}{\partial \xi }=\frac{u\left(x,{\xi }_{k+1}\right)-u\left(x,{\xi }_k\right)}{\xi }+o\left(\tau \right)$ (6)

代入式(5)得到

$\begin{array}{c}\frac{{\partial }^{\alpha }u\left(x,t_{n+1}\right)}{\partial t^{\alpha }}\approx \frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle \underset{k=0}{\overset{n}{\sum }}\frac{u\left(x,t_{k+1}\right)-u\left(x,t_k\right)}{\tau }}{\displaystyle \underset{k\tau }{\overset{\left(k+1\right)\tau }{\int }}\frac{\text{d}\xi }{{\left(t_{k+1}-\xi \right)}^{\alpha }}}\\ =\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle \underset{k=0}{\overset{n}{\sum }}\frac{u\left(x,t_{n+1-k}\right)-u\left(x,t_{n-k}\right)}{\tau }}{\displaystyle \underset{k\tau }{\overset{\left(k+1\right)\tau }{\int }}\frac{\text{d}r}{r^{\alpha }}}\\ =\{\begin{array}{l}\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left(u^{n+1}-u^n\right)+\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{k=1}{\overset{n}{\sum }}\left(u^{n+1-k}-u^{n-k}\right)\left[{\left(k+1\right)}^{1-\alpha }-k^{1-\alpha }\right],n\ge 1}\\ \frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left(u^1-u^0\right)\end{array}\\ =\{\begin{array}{l}{a}_0\left(u^{n+1}-u^n\right)+a_0{\displaystyle \underset{k=1}{\overset{n}{\sum }}{b}_k}\left(u^{n+1-k}-u^{n-k}\right),n\ge 1\\ a_0\left(u^1-u^0\right),n=0\end{array}\end{array}$ (7)

其中 $a_0=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)},b_k={\left(k+1\right)}^{1-\alpha }-k^{1-\alpha },\left(k=0,1,2,\cdots ,n\right),u^0=u\left(x,t=0\right)=w\left(x\right)$ ，将式(7)代入式(1)得时间离散化方程

$a_0{u}^{n+1}+u^{n+1}-k{\nabla }^2{u}^{n+1}=\{\begin{array}{l}{a}_0\left\{u^n-{\displaystyle \underset{k=1}{\overset{n}{\sum }}{b}_k\left(u^{n+1-k}-u^{n-k}\right)}\right\}+f^{n+1},n\ge 1\\ a_0{u}_0+f^1,n=0\end{array}$ (8)

其中， $f^{n+1}=f\left(x,t_{n+1}\right),n=0,1,\cdots ,N$ 。

3.2. 空间上的函数逼近

在固定时间层 $t_{n+1}=\left(n+1\right)\tau$ 上，借助径向基函数逼近未知函数得

$u\left(x_i,t_{n+1}\right)={\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j^{n+1}}\phi \left(r_{ij}\right)+{\lambda }_{M+1}^{n+1}{x}_i+{\lambda }_{M+2}^{n+1}$ (9)

其中 $\left\{{\lambda }_j^{n+1}\right\}$ 是该层的未知系数， $\phi \left(r_{ij}\right)$ 是径向基函数， $r_{ij}={\Vert x_i-x_j\Vert }_2$ 。

除了式(9)中的M个点的公式外，还需要以下2个方程的正则化条件

${\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j^{n+1}}={\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j^{n+1}}{x}_j=0$ (10)

将式(9)代入式(8)，考虑式(10)和边界条件，可以得到以下矩阵形式的离散方程

$B{\left\{\lambda \right\}}^{n+1}=b^{n+1}$ (11)

其中

$B={\left[\begin{array}{ccccccc}L\left({\phi }_{11}\right)& \cdots & L\left({\phi }_{1j}\right)& \cdots & L\left({\phi }_{1M}\right)& L\left(x_1\right)& L\left(1\right)\\ ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮\\ L\left({\phi }_{i1}\right)& \cdots & L\left({\phi }_{1j}\right)& \cdots & L\left({\phi }_{iM}\right)& L\left(x_i\right)& L\left(1\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ L\left({\phi }_{M1}\right)& \cdots & L\left({\phi }_{Mj}\right)& \cdots & L\left({\phi }_{MM}\right)& L\left(x_M\right)& L\left(1\right)\\ x_1& \cdots & x_j& \cdots & x_M& 0& 0\\ 1& \cdots & 1& \cdots & 1& 0& 0\end{array}\right]}_{\left(M+2\right)\times \left(M+2\right)}$ (12)

其中 ${\phi }_{ij}=\phi \left(r_{ij}\right)$ ，L代表算子

$L(\ast )=\{\begin{array}{l}\left(a_0+1-\kappa {\nabla }^2\right)(\ast ),1<i<M\\ (\ast ),i=1\text{or}1=M,\end{array}$ (13)

且 $b^{n+1}={\left[b_1^{n+1},\cdots ,b_i^{n+1},\cdots ,b_M^{n+1},0,0,0,0\right]}^{\text{T}}$ ，

$\begin{array}{l}{b}_i^{n+1}=\{\begin{array}{l}{a}_0{u}_i^0+f_i^1,n=0,1<i<M\\ a_0\left\{u_i^n-{\displaystyle \underset{k=1}{\overset{n}{\sum }}{b}_k\left(u_i^{n+1-k}-u_i^{n-k}\right)}\right\}+f_i^{n+1},n\ge 1,1<i<M\\ g\left(x_i,t_{n+1}\right),i=1\text{or}i=M\end{array}\\ u_i^{n+1}=u\left(x_i,t_{n+1}\right),f_i^{n+1}=f\left(x_i,t_{n+1}\right),并且n=0,1,\cdots ,N\end{array}$ (14)

极坐标下的拉普拉斯算子表示为

${\nabla }^2\left(\phi \left(r\right)\right)=\frac{\partial \phi \left(r\right)}{\partial r}\left(\frac{{\partial }^{2}r}{\partial x^2}+\frac{{\partial }^{2}r}{\partial y^2}+\frac{{\partial }^{2}r}{\partial z^2}\right)+\frac{{\partial }^2\phi \left(r\right)}{\partial r^2}\left[{\left(\frac{\partial r}{\partial x}\right)}^2+{\left(\frac{\partial r}{\partial y}\right)}^2+{\left(\frac{\partial r}{\partial z}\right)}^2\right]$ (15)

特别地，当径向基函数分别取MQ，TPS和Cubic时，可以分别求出其一二阶导数的具体形式，详见表1。

Table 1. The first and second derivatives of three basis functions

表1. 三个基函数的一阶、二阶导数

4. 数值实验

本部分通过数值算例考察使用三种不同基函数求解一维分数阶扩散方程的优缺点，使用如下三类误差衡量求解精度

$\begin{array}{l}{L}_1=\underset{j}{\max }\left|u_j-U_j\right|,\\ L_2=\sqrt{{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left|u_j-U_j\right|}^2}},\\ \text{RMS}=\sqrt{\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left|u_j-U_j\right|}^2}}\end{array}$ (16)

其中 $u_j$ 是精确解， $U_j$ 是数值解。

下面分别考虑相同计算时间内不同时间步长和不同计算时间内相同时间步长两类情况下，三个基函数的逼近误差情况。

Table 2. Error of different time steps at T = 1

表2. T = 1时不同时间步长的误差

表2为T = 1时不同时间步长的计算误差结果，由表中数值可以看出，MQ基函数中形状参数取不同值(c = 0.1和c = 0.15)，TPS基函数中形状参数取不同值(β = 3和β = 15)时，计算结果有显著差异，说明MQ和TPS对应结果依赖于形状参数这一局限性。随着时间步长τ的减小，三种基函数对应的求解误差均在减小，且具有相同精度。

Table 3. Error at different times τ = 0.01

表3. τ = 0.01时不同时间的误差

表3为τ = 0.01时不同计算时间的误差结果，由此表类似得出MQ和TPS基函数的计算结果依赖于形状参数的选择，此外，在相同的时间步长下Cubic基函数对应求解结果和TPS，MQ对应结果具有相同精度。

同时，我们分别做出了数值解和对应L₁误差的图形。图1从左依次为T = 1，τ = 0.1时MQ、TPS、Cubic的数值解，图2从左依次为T = 1，τ = 0.1时MQ、TPS、Cubic的误差。

从图1和图2中可以观察到，在恰当选取形状参数(MQ (c = 0.1)、TPS (β = 3))时，三种径向基函数具有类似的数值求解结果和求解精度，而此时Cubic径向基函数具有不依赖于形状参数的显著优势。

Figure 1. Numerical solution comparison

图1. 数值解比较

Figure 2. Error comparison

图2. 误差比较

5. 结论

文章应用无网格方法求解一维时间分数阶扩散方程。在时间上用有限差分法离散，空间上分别选取MQ，TPS和Cubic三种径向基函数近似未知函数，通过数值实验比较得出三种径向基函数具有类似的求解结果和逼近精度，但Cubic径向基函数无须选择形状参数，研究结果有助于无网格数值方法在微分方程数值求解相关领域的推广和应用。

基金项目

大学生创新创业训练计划项目(X202110452136)。