1908年法国物理学家Langevin在研究布朗运动时提出：溶液中介质分子对布朗粒子的作用力可分为两个部分,一部分是流体对布朗粒子的黏滞阻力,另一部分是溶液分子无规则热运动施加在布朗粒子上的一种涨落不定的力.由此, Langevin在牛顿第二定律的基础上提出Langevin方程^[1]

(1.1) $ \begin{equation} m\frac{{\rm d}^2x}{{\rm d}t^2} = -\lambda\frac{{\rm d}x}{{\rm d}t}+F(t)+\xi(t), \end{equation} $

其中, $ m $是粒子的质量, $ \lambda $是黏滞系数, $ F $是外部力, $ \xi(t) $是随机力,服从高斯分布. $ \xi(t) $一般也可用$ \frac{{\rm d}W(t)}{{\rm d}t} $来代替, $ W(t) $表示经典布朗运动,该式将由质量和加速度共同决定的粒子运动力等价表示为,与速度有关的阻力、外部力和随机力等三个力对粒子共同作用的结果.

进一步在整数阶Langevin方程中引入分数阶算子.与整数阶导数相比,分数阶导数具有记忆性和全局相关性,使变量对过去时刻的状态有依赖性. West(2002)观察到股票价格波动具有“肥尾性”和短暂的时间相关性,衡量市场波动的指标具有逆幂律衰减的时间相关性,于是提出用分数阶Langevin方程作为金融时间序列的动力学模型^[2-3].关于分数阶Langevin方程的研究还包括长期相关性^[4-5]、位移相关性等方面^[6-7].

关于数值求解分数阶Langevin方程,通常构造的数值解有显式算法和隐式算法.文献[8]在没有确定性场的情况下引入分数阶Langevin方程,考虑粒子的初始条件和初始速度,利用拉普拉斯变换的方法和三参数Mittag-Leffler函数得到方程的解.文献[9]提出Jacobi-Gauss-Lobatto方法来求解具有三点边界条件的非线性分数阶Langevin方程,并通过数值算例检验了该方法的准确性.

在各类数值解法中,显式算法计算过程较为简便,但精度比较差,稳定性不好,隐式算法虽然提高了精度,但算法复杂,计算量大.文献[10]预估校正算法在结合了两者的优点基础上,同时简化了算法的计算量.文献[11]提出的用预估校正Euler方法来求解随机微分方程,该方法能克服使用显式欧拉方法时经常遇到的一些数值不稳定性,数值实验证明了其提出的对称预估校正欧拉方法的渐近稳定性.文献[12]采用路径逼近的方法,建立了用马尔可夫切换求解随机微分方程的数值方法,并且用预测校正Euler-Maruyama方法来克服近似路径模拟过程中的误差传播,证明了数值解对精确解的强收敛性,并揭示了在一定条件下误差在系数函数上的有序性.文献[13]将预估校正算法运用于Itô型随机微分方程,推导出预估校正方法的稳定性条件,证明在适当的条件下,预估校正方法对于所有非常小的时间步长保持时间指数稳定性.文献[14]针对Itô型随机延迟微分方程,在Lipschitz和线性增长条件下,预估校正方法被证明是min$ (1/2, p) $阶均方收敛的,并证明了对于某些$ p $,预估校正方法的渐近MS稳定性界限范围远大于Euler-Maruyama方法.文献[15]用预估校正方法求解随机比例延迟方程,证明了预估校正方法是1/2阶强收敛且线性稳定的.文献[16]提出线性隐式预估校正方法用于求解具有非光滑初始数据的空间分数阶反应-扩散方程,并且证明该方法是无条件稳定的,时间步长中引起的振荡行为随着空间分数导数的阶数减小(慢扩散)而减小.

本文尝试用预估校正算法来求解分数阶Langevin方程.将(1.1)式中质量$ m $取值为1,将加速度$ a $定义为速度$ v $的分数阶导数,定义外部力$ F(t) $与初始速度$ v_0 $和时间$ t $有关,给速度的耗散参数$ \lambda $添加合适的标度单位,考虑如下分数阶Langevin方程(参见文献[1])

(1.2) $ \begin{equation} \left\{\begin{array}{ll} { } {}_{0}D_{t}^{\alpha}[v(t)] = -\lambda^\alpha v(t)+v_0\frac{t^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t), &0<t\leq T, 0<\alpha\leq1, \\ v^{(k)}(0) = v_0^{(k)}. \end{array}\right. \end{equation} $

考虑(1.2)式的分数阶Langevin方程, $ _{0}D_{t}^{\alpha} $表示Caputo定义的分数阶导数,取均匀网格节点$ t_{j} = jh(h = T/N) $, $ v(t_{j}) = v_{j} $, $ j = 0, 1, \cdots, N $, (1.2)式等价于(参见文献[17-18])

(2.1) $ \begin{equation} v(t) = \frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-\tau)^{\alpha-1} \bigg[-\lambda^\alpha v(\tau)+v_0\frac{\tau^{-\alpha}}{\Gamma(1-\alpha)}+\xi(\tau)\bigg]\, {\rm d}\tau \end{equation} $

或

(2.2) $ \begin{equation} v(t) = J^{\alpha}\bigg[-\lambda^\alpha v(t)+v_0\frac{t^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t)\bigg]. \end{equation} $

利用R2算法将Caputo积分离散化,其中, R2算法(参见文献[18])

$ D^{\alpha}f(t_n)\approx h^{-\alpha}\sum\limits_{j = 0}^{N}C_{n, j}^{(\alpha)}f_j, $

$ C_{n, j}^{(\alpha)} = \frac{1}{\Gamma(2-\alpha)} \left\{\begin{array}{ll}(1-\alpha)n^{-\alpha}-b_{n-1}^{(\alpha)}, & j = 0, \\ b_{n-j}^{(\alpha)}-b_{n-j-1}^{(\alpha)}, & 1\leq j\leq n-1, \\ 1, & j = n, \\ 0, & \mbox{其他}. \end{array}\right. $

$ b_j^{(\alpha)} = (j+1)^{(1-\alpha)}-j^{(1-\alpha)}, $

得到

(2.3) $ \begin{equation} v_h(t_{n+1}) = h^{\alpha}\sum\limits_{j = 0}^{n+1}a_{j, n+1}\bigg[-\lambda^\alpha v_h(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg], \end{equation} $

其中,系数

$ a_{j, n} = \frac{1}{\Gamma(2+\alpha)}\left\{\begin{array}{ll}(1+\alpha)n^{\alpha}-n^{1+\alpha}-(n-1)^{1+\alpha}, & j = 0, \\ (n-j+1)^{(1+\alpha)}-2(n-j)^{(1+\alpha)}+(n-j-1)^{(1+\alpha)}, & 1\leq j\leq n-1, \\ 1, & j = n. \end{array}\right. $

差分逼近格式(2.3)式称为分数阶Adams-Moulton方法,等号右边的$ v_h(t_{n+1}) $项是未知的,用R0算法先算出预估值(参见文献[18])

(2.4) $ \begin{equation} v_h^p(t_{n+1}) = h^{\alpha}\sum\limits_{j = 0}^{n}b_{j, n+1}\bigg[-\lambda^\alpha v_h(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg], \end{equation} $

其中$ b_{j, n} = \frac{(n-j)^{\alpha}-(n-j-1)^{\alpha}}{\Gamma(1+\alpha)} $.

将(2.4)式代入(2.3)式右边,得到(1.2)式分数阶Langevin方程的预估校正算法

(2.5) $ \begin{eqnarray} v_h(t_{n+1})& = &h^{\alpha}\sum\limits_{j = 0}^{n}a_{j, n+1} \bigg[-\lambda^\alpha v_h(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg] {}\\ &&+h^{\alpha}\bigg[-\lambda^\alpha v_h^p(t_{n+1})+v_0\frac{t_{n+1}^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_{n+1})\bigg]. \end{eqnarray} $

引理3.1  假设$ v(t)\in C^1[0, T] $,那么

(3.1) $ \begin{equation} \bigg|J^\alpha v(t_n)-h^\alpha\sum\limits_{j = 0}^{n-1}b_{j, n}v(t_j)\bigg|\leq\frac{1}{\Gamma(1+\alpha)}\| v'\| _{\infty}t_n^\alpha h. \end{equation} $

证  因为

$ J^{\alpha}v(t_n) = \frac{1}{\Gamma(\alpha)}\int_{0}^{t_n}(t_n-\tau)^{\alpha-1}v(\tau)\, {\rm d}\tau = \frac{1}{\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h}(t_n-\tau)^{\alpha-1}v(\tau)\, {\rm d}\tau, $

$ \begin{eqnarray*} h^\alpha\sum\limits_{j = 0}^{n-1}b_{j, n}v(t_j) & = &h^\alpha\frac{1}{\Gamma(1+\alpha)}\sum\limits_{j = 0}^{n-1}[(n-j)^\alpha-(n-j-1)^\alpha]v(t_j) {\nonumber} \\ & = &\frac{1}{\alpha\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}[(t_n-jh)^\alpha-(t_n-(j+1)h)^\alpha]v(t_j) {\nonumber} \\ & = &\frac{1}{\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h}(t_n-\tau)^{\alpha-1}v(t_j)\, {\rm d}\tau, {\nonumber} \end{eqnarray*} $

所以有

$ \begin{eqnarray*} J^{\alpha}v(t_n)-h^\alpha\sum\limits_{j = 0}^{n-1}b_{j, n}v(t_j) & = &\frac{1}{\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h}(t_n-\tau)^{\alpha-1}v(\tau)\, {\rm d}\tau\\ && -\frac{1}{\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h}(t_n-\tau)^{\alpha-1}v(t_j)\, {\rm d}\tau {\nonumber} \\ & = &\frac{1}{\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h}[(t_n-\tau)^{\alpha-1}(v(\tau)-v(t_j))]\, {\rm d}\tau. {\nonumber} \end{eqnarray*} $

由中值定理及$ t_n = nh $,令$ \xi\in [t_j, \tau] $,有

$ \begin{eqnarray*} &&J^{\alpha}v(t_n)-h^\alpha\sum\limits_{j = 0}^{n-1}b_{j, n}v(t_j) {\nonumber} \\ & = &\frac{1}{\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h}[(t_n-\tau)^{\alpha-1}(v (\tau)-v(t_j))]\, {\rm d}\tau {\nonumber} \\ & = &\frac{1}{\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h}[(t_n-\tau)^{\alpha-1}v'(\xi) (\tau-jh)]\, {\rm d}\tau{\nonumber} \\ &\leq&\| v'\| _{\infty} \frac{1}{\Gamma(\alpha)}\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h}[(t_n-\tau)^{\alpha-1} (\tau-jh)]\, {\rm d}\tau {\nonumber} \\ & = &\| v'\| _{\infty} \frac{1}{\Gamma(\alpha)}(-\frac{1}{\alpha})\sum\limits_{j = 0}^{n-1}\int_{jh}^{(j+1)h} (\tau-jh)\, {\rm d}(t_n-\tau)^{\alpha} {\nonumber} \\ & = &\| v'\| _{\infty} \frac{1}{\Gamma(\alpha)}(-\frac{1}{\alpha})\sum\limits_{j = 0}^{n-1} \bigg[(\tau-jh)(t_n-\tau)^\alpha|_{jh}^{(j+1)h}-\int_{jh}^{(j+1)h} (t_n-\tau)^{\alpha}\, {\rm d}(\tau-jh)\bigg] {\nonumber} \\ & = &\| v'\| _{\infty} \frac{1}{\Gamma(\alpha)}(-\frac{1}{\alpha})\sum\limits_{j = 0}^{n-1} \bigg[h(nh-(j+1)h)^\alpha-(-\frac{1}{\alpha+1})\int_{jh}^{(j+1)h} \, {\rm d}(t_n-\tau)^{\alpha+1}\bigg] \\ & = &\| v'\| _{\infty} \frac{1}{\Gamma(\alpha)}(-\frac{1}{\alpha})\sum\limits_{j = 0}^{n-1} \bigg[h^{\alpha+1}(n-j-1)^\alpha+\frac{1}{\alpha+1}[(nh-(j+1)h)^{\alpha+1}-(nh-jh)^{\alpha+1}]\bigg] {\nonumber} \\ & = &\| v'\| _{\infty} \frac{1}{\Gamma(\alpha)}(-\frac{1}{\alpha})\sum\limits_{j = 0}^{n-1} \bigg[h^{\alpha+1}(n-j-1)^\alpha+\frac{1}{\alpha+1}h^{\alpha+1}[(n-j-1)^{\alpha+1}-(n-j)^{\alpha+1}]\bigg] {\nonumber} \\ & = &\| v'\| _{\infty} \frac{1}{\Gamma(\alpha)}\frac{h^{\alpha+1}}{\alpha}\sum\limits_{j = 0}^{n-1} \bigg[\frac{1}{\alpha+1}[(n-j)^{\alpha+1}-(n-j-1)^{\alpha+1}]-(n-j-1)^{\alpha}\bigg] {\nonumber} \\ & = &\| v'\| _{\infty} \frac{1}{\Gamma(\alpha)}\frac{h^{\alpha+1}}{\alpha} \bigg[\frac{1}{\alpha+1}n^{\alpha+1}-\sum\limits_{j = 0}^{n-1}j^\alpha\bigg] {\nonumber} \\ & = &\| v'\| _{\infty} \frac{h^{\alpha+1}}{\Gamma(\alpha+1)} \bigg[\int_{0}^{n}\tau^\alpha \, {\rm d}\tau-\sum\limits_{j = 0}^{n-1}j^\alpha\bigg]. \end{eqnarray*} $

由于

$ \int_{0}^{n}\tau^\alpha\, {\rm d}\tau\leq 1^\alpha(1-0)+2^\alpha(2-1)+\cdots+n^\alpha(n-(n-1)) = \sum\limits_{j = 0}^{n}j^\alpha, $

所以

$ \begin{eqnarray*} \bigg|J^{\alpha}v(t_n)-h^\alpha\sum\limits_{j = 0}^{n-1}b_{j, n}v(t_j)\bigg| &\leq&\| v'\| _{\infty} \frac{h^{\alpha+1}}{\Gamma(\alpha+1)} \bigg[\int_{0}^{n}\tau^\alpha \, {\rm d}\tau-\sum\limits_{j = 0}^{n-1}j^\alpha\bigg] {\nonumber} \\ &\leq&\| v'\| _{\infty} \frac{h^{\alpha+1}}{\Gamma(\alpha+1)} \bigg[\sum\limits_{j = 0}^{n}j^\alpha-\sum\limits_{j = 0}^{n-1}j^\alpha\bigg]\\ & = &\| v'\| _{\infty} \frac{h^{\alpha+1}}{\Gamma(\alpha+1)}n^\alpha = \| v'\| _{\infty} \frac{1}{\Gamma(\alpha+1)}t_n^\alpha h. {\nonumber} \end{eqnarray*} $

引理3.1得证.

引理3.2  假设$ v(t)\in C^2[0, T] $,那么存在依赖于$ \alpha $的常数$ C_\alpha $,使得

(3.2) $ \begin{equation} \bigg|J^\alpha v(t_n)-h^\alpha\sum\limits_{j = 0}^{n}a_{j, n}v(t_j)\bigg|\leq C_\alpha\| v'\| _{\infty}t_n^\alpha h^2. \end{equation} $

证  过程与引理3.1证明类似.

定理3.1  当$ 0<\alpha<1 $时,假设$ _{0}^{c}D_{t}^{\alpha}v(t)\in C^2[0, T] $,那么

(3.3) $ \begin{equation} \max\limits_{0\leq j\leq N}|v(t_j)-v_h(t_j)| = o(h^{1+\alpha}), \end{equation} $

其中$ N = [\frac{T}{h}] $.

证  令$ q = 1+\alpha $,利用数学归纳法(参见文献[19]). $ j = 0 $时,成立; $ j = 1 $时,成立.假设$ j = 2, 3, \cdots, k $时, (3.3)式成立,要证明$ j = k+1 $时成立.

先考察预估值$ v_h^p(t_{k+1}) $的误差

$ \begin{eqnarray*} \big|v(t_{k+1})-v_h^p(t_{k+1})\big| & = &\bigg|J^{\alpha} \bigg[-\lambda^\alpha v(t_{k+1})+v_0\frac{t_{k+1}^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_{k+1})\bigg] \\ &&-h^{\alpha}\sum\limits_{j = 0}^{k}b_{j, k+1} \bigg[-\lambda^\alpha v_h(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg] \bigg| {\nonumber} \\ & = &\bigg|J^{\alpha}\bigg[-\lambda^\alpha v(t_{k+1})+v_0\frac{t_{k+1}^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_{k+1})\bigg] \\ &&-h^{\alpha}\sum\limits_{j = 0}^{k}b_{j, k+1} \bigg[\bigg[-\lambda^\alpha v(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg] {\nonumber} \\ && -\bigg[-\lambda^\alpha v(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg] \\ &&+\bigg[-\lambda^\alpha v_h(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg]\bigg]\bigg| {\nonumber} \\ & = &\bigg|J^{\alpha}\bigg[-\lambda^\alpha v(t_{k+1})+v_0\frac{t_{k+1}^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_{k+1})\bigg] \\ &&-h^{\alpha}\sum\limits_{j = 0}^{k}b_{j, k+1}\bigg[-\lambda^\alpha v(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg] {\nonumber} \\ &&+ h^{\alpha}\sum\limits_{j = 0}^{k}b_{j, k+1} \bigg[\bigg[-\lambda^\alpha v(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg] \\ &&-\bigg[-\lambda^\alpha v_h(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg]\bigg]\bigg| {\nonumber} \\ &\leq&\bigg|J^\alpha{}_0^c D_{t}^{\alpha}v(t_{k+1})-h^{\alpha}\sum\limits_{j = 0}^{k}b_{j, k+1}{}_0^c D_{t} ^{\alpha}v(t_{k+1})\bigg|\\ &&+h^{\alpha}\sum\limits_{j = 0}^{k}b_{j, k+1} \bigg|[-\lambda^{\alpha}v(t_j)]-[-\lambda^{\alpha}v_h(t_j)]\bigg|. \end{eqnarray*} $

由$ j = k $时, $ \max\limits_{0\leq j\leq N}|v(t_j)-v_h(t_j)| = o(h^{q}) $成立,结合引理3.1,可得

$ \begin{eqnarray*} \big|v(t_{k+1})-v_h^p(t_{k+1})\big|&\leq &C_1t_{k+1}^\alpha h+h^\alpha\lambda^\alpha\sum\limits_{j = 0}^{k}b_{j, k+1}h^q {\nonumber} \\ & = &C_1t_{k+1}^\alpha h+\lambda^\alpha h^\alpha(k+1)^\alpha h^q\frac{1}{\Gamma(1-\alpha)} {\nonumber} \\ &\leq&C_1t_{k+1}^\alpha h+C_2t_{k+1}^\alpha h^q. {\nonumber} \end{eqnarray*} $

对校正误差进行分析

$ \begin{eqnarray*} \big|v(t_{k+1})-v_h(t_{k+1})\big| & = &\bigg| J^\alpha{}_0^c D_{t}^{\alpha}v(t_{k+1})-h^{\alpha}\sum\limits_{j = 0}^{k+1}a_{j, k+1}{}_0^c D_{t}^{\alpha}v(t_{k+1}) {\nonumber} \\ &&+\bigg(h^{\alpha}\sum\limits_{j = 0}^{k}a_{j, k+1} \bigg[-\lambda^\alpha v(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg]\\ && +h^{\alpha}\bigg[-\lambda^\alpha v(t_{k+1})+v_0\frac{t_{k+1}^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_{k+1})\bigg]\bigg) {\nonumber} \\ &&-\bigg(h^{\alpha}\sum\limits_{j = 0}^{k}a_{j, k+1} \bigg[-\lambda^\alpha v_h(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg] \\ && +h^{\alpha}\bigg[-\lambda^\alpha v_h^p(t_{k+1})+v_0\frac{t_{k+1}^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_{k+1})\bigg]\bigg)\bigg| \\ &\leq&\bigg| J^\alpha{}_0^c D_{t}^{\alpha}v(t_{k+1})-h^{\alpha}\sum\limits_{j = 0}^{k+1}a_{j, k+1}{}_0^c D_{t}^{\alpha}v(t_{k+1})\bigg| {\nonumber} \\ &&+h^{\alpha}\sum\limits_{j = 0}^{k}a_{j, k+1}\bigg| \bigg[-\lambda^\alpha v(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg]\\ &&-\bigg[-\lambda^\alpha v_h(t_j)+v_0\frac{t_j^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_j)\bigg]\bigg| {\nonumber} \\ &&+h^\alpha\bigg|\bigg[-\lambda^\alpha v(t_{k+1})+v_0\frac{t_{k+1}^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_{k+1})\bigg]\\ &&-\bigg[-\lambda^\alpha v_h^p(t_{k+1})+v_0\frac{t_{k+1}^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t_{k+1})\bigg]\bigg|. \end{eqnarray*} $

由预估值$ v_h^p(t_{k+1}) $误差和引理3.2,进一步得到

$ \begin{eqnarray*} \big|v(t_{k+1})-v_h(t_{k+1})\big| &\leq &c_3t_{k+1}^\alpha h^2+ch^\alpha h^q\sum\limits_{j = 0}^{k}a_{j, k+1}+c'h^\alpha(h+h^q) {\nonumber} \\ &\leq &c_3t_{k+1}^\alpha h^2+c_4h^{1+\alpha}+c_5t_{k+1}^\alpha h^{\alpha+q} \leq ch^q, {\nonumber} \end{eqnarray*} $

因此,定理3.1得证.

考虑如下分数阶Langevin方程^[2]

(4.1) $ \begin{equation} \left\{\begin{array}{ll} { } {}_{0}D_{t}^{\alpha}[v(t)] = -\lambda^\alpha v(t)+v_0\frac{t^{-\alpha}}{\Gamma(1-\alpha)}+\xi(t), &0<t\leq T, 0<\alpha\leq1, \\ v^{(k)}(0) = v_0^{(k)}. \end{array}\right. \end{equation} $

用方程的逆拉普拉斯变换,得到分数阶Langevin方程的一般解为^[1]

(4.2) $ \begin{equation} v(t) = E_\alpha[-(t)^{\alpha}]+\int_{0}^{t}\frac{E_{\alpha, \alpha}[-(t-t')^\alpha]} {(t-t')^{1-\alpha}}\xi(t')\, {\rm d}t', \end{equation} $

其中, Mittag-Leffler函数表示为

(4.3) $ \begin{equation} E_{\alpha, \beta}(z) = \sum\limits_{k = 0}^{\infty}\frac{z^k}{\Gamma(k\alpha+\beta)}, \alpha>0, \beta>0, \end{equation} $

$ E_\alpha(z) $是$ \beta = 1 $时的特殊情况.

为了验证用预估校正算法(2.5)式的方法来数值求解分数阶Langevin方程的有效性,将其与方程的真解(4.2)式进行对比(参见文献[20-22]),得到在不同$ \alpha $和$ h $取值下,分数阶Langevin方程预估校正算法数值解和真解的图形,其中,令$ T = 20, \lambda = 1, v_0 = 10 $.

为了削弱随机项对数值解的影响,得到的数值解是由50条样本轨迹取均值来表现的,而真值则由随机得到的某一条轨迹来表示.通过得到的真值和数值解的比较图形可以看出,数值解与真值趋势相同,且始终在真值附近波动.

进一步,用数值算例结果证明预估校正算法数值解误差的收敛阶.设收敛阶为$ \gamma $,对充分小的步长$ h $和任意固定一点$ \tau = nh $,有

$ E|v(t_n)-v_h(t_n)|\leq ch^\gamma. $

记误差

$ e(h) = E|v(t_n)-v_h(t_n)|, $

就有

$ e(h)\leq ch^\gamma. $

在两边同时取对数变形为

$ \log e(h)\approx \log c+\gamma \log h. $

数值算例已经得到不同$ \alpha $, $ h $取值下的真解与预估校正算法数值解的误差,考虑每一个$ \alpha $取值下,用最小二乘法拟合线性方程的$ \log e(h) = \log c+\gamma \log h $的系数$ \gamma $,得到对应的收敛阶.误差及收敛阶的值具体如表 1.

在每一个$ \alpha $取值下,可以在$ \log e(h)-\log h $坐标轴中,画出对应的点,并画出用最小二乘法拟合得到的直线,同时给出斜率为$ 1+\alpha $的直线,得到的结果如图 6所示.

固定$ h = 0.2 $,作出不同$ \alpha $下的数值解,得到如图 7所示图形,可以看到,当$ \alpha $比较小时,数值解一开始波动比较大,后面慢慢趋于稳定,而当$ \alpha $接近1时,数值解从一开始就比较稳定.

本文针对一类分数阶Langevin方程,利用预估校正算法对其进行数值求解.利用R2算法,求解$ v(t_{n+1}) $值.表达式中未知的$ v(t_{n+1}) $值用R0算法通过前$ t $个时刻的$ v $值来预估,最终得到预估校正算法的数值解.误差分析证明在$ 0<\alpha<1 $条件下,预估校正算法是$ (1+\alpha) $阶收敛的.数值试验部分将预估校正算法的数值解与通过Laplace变化求得的解析解进行误差分析和收敛误差阶验证.在不同$ \alpha $, $ h $取值下,预估校正算法的数值解的误差都很小且稳定.