刘 鑫,张建影
(长春工业大学 数学与统计学院,长春 130012)
分数阶微积分方程是在传统微积分方程的基础上,用分数阶导数(积分)项代替传统微积分方程的时间或空间导数(积分)项.与整数阶方程相比,分数阶微积分方程在描述一些特殊问题时更有意义.例如: 多孔介质渗流问题[1];非牛顿流体的非局域性问题[2];描述地下水在非饱和土中的渗流过程[3];将分数阶方程与最优控制结合起来开发高效算法[4]等.但由于分数阶算子的特殊性,使得寻找分数阶微积分方程的精确解析解较困难.相比于求解精确解析解,数值解法更灵活,也适用于更多情况.目前关于分数阶微积分方程数值解的研究已有一定进展,例如: Mukhtar等[5]用Adomian分解变换方法和q-同伦分析变换方法,得到了分数阶方程多维模型的数值结果;Ali等[6]用最优同伦渐近方法推导出Caputo型分数阶导数的半解析解;Burqan等[7]利用Laplace变换和残差幂级数方法求解了分数阶方程的级数解.
时间分数阶方程是分数阶方程的重要分支,如何开发其数值算法目前已得到广泛关注.相比于宏观的有限差分方法,在求解时间分数阶微分方程时使用的复杂格式降低了算法的有效性.介于统计力学与流体力学之间的人工系统——格子Boltzmann模型,其更简单并适用于数值模拟.本文用格子Boltzmann方法(LBM)求解修正的时间序列分数阶方程.
考虑一个修正的时间分数阶扩散方程[8]:
(1)
(2)
整理可得
(4)
其中
下面用格子Boltzmann方法恢复宏观方程.定义
(6)
这里fα(x,t)表示在时间为t、位置为x、速度为eα的粒子的分布函数,α为粒子速度方向.格子Boltzmann方程为
(7)
(8)
引入小Knudsen数ε,定义ε=l/L=Δt,其中l为平均自由程,L为宏观特征长度.在引入ε的情况下,对分布函数fα(x,t)进行Chapman-Enskog展开,有
(9)
(10)
进行时间多尺度展开可得
t=t0+εt1+ε2t2+ε3t3+ε4t4+O(ε5),
(11)
从而有
(12)
源项有如下表达形式:
Ωα=εΩ1+ε2Ω2+ε3Ω3+ε4Ω4+O(ε5).
(13)
将式(8),(9),(11)~(13)代入式(7),并比较ε的各阶系数有
(16)
设源项满足如下条件:
平衡态分布函数满足:
(19)
(20)
(21)
(22)
根据式(17)~(21)对式(22)进行整理,即可恢复出宏观方程
(23)
图1 D1Q3的格子模型Fig.1 Lattice model of D1Q3
这里=-Dλεc2
下面推导平衡态分布函数.D1Q3的格子模型如图1所示.
图2 D2Q9的格子模型Fig.2 Lattice model of D2Q9
根据式(19)~(21)可得平衡态分布函数为
(24)
(25)
其中c=|eα|是粒子沿各方向的移动速度.D2Q9的格子模型如图2所示.
同理可得平衡态分布函数为
下面计算一个一维实例,因此所涉及的矢量x这里均为标量x.定义uN(xk,t)为数值解,uE(xk,t)为精确解.选择相对误差、最大绝对误差和全局相对误差进行验证,各误差表达式分别如下:
(29)
(30)
(31)
例1考虑带有如下初边值条件的时间分数阶方程:
(32)
方程(32)的源项为
(33)
其精确解为
u(x,t)=t2ex.
(34)
图3 例1的数值计算结果(A)、误差分析(B)和收敛阶(C)Fig.3 Numerical calculation results (A),error analysis (B) and convergence order (C) of example 1
图3给出了例1的数值计算结果、误差分析和收敛阶.由图3(A)可见,数值解和精确解曲线吻合度较高,几乎完全重叠;由图3(B)可见,最大绝对误差为0.009 10,最大相对误差为0.004 52,两种误差都较小;由图3(C)可见,在一维情形下格子Boltzmann模型具有空间一阶收敛阶.
(35)
相对误差、最大绝对误差和全局相对误差的表达式如下:
例2考虑带有如下初边值条件的时间分数阶方程:
(39)
方程(39)的源项为
(40)
精确解为
u(x,y,t)=t2ex+y.
(41)
数值计算的参数选取: 时间步长和空间步长分别为
格子数为65×65,时间t=1.0,松弛因子τ=1.25,分数阶取值为β=0.001,α=0.001,扩散系数为D=0.001.
图4 例2的数值计算结果(A),(B)、误差分析(C)和收敛阶(D)Fig.4 Numerical calculation results (A),(B),error analysis (C) and convergence order (D) of example 2
图4给出了例2的二维数值计算结果、误差分析和收敛阶.由图4(A),(B)可见,方程(39)的数值解和精确解基本一致;由图4(C)可见,最大相对误差是0.012 2,最大绝对误差是0.009 06;由图4(D)可见,在二维情形下格子Boltzmann模型具有空间一阶精度.
综上所述,本文基于格子Boltzmann方法,针对修正的时间分数阶方程进行了讨论.首先将方程离散化处理,得到了标准的反应扩散方程.其次,根据Taylor展开式和Chapman-Enskog多尺度展开等技术,用格子Boltzmann方法准确地恢复了宏观方程,并推导出各方向的平衡态分布函数表达式.最后,用一个一维实例和一个二维实例进行数值计算,所得数值解与精确解吻合较好,误差在合理范围内.因此,格子Boltzmann模型在求解修正的时间分数阶方程的数值解上效果较好.