基于非均匀网格的预测校正法求解非线性分数阶常微分方程初值问题

A Predictor-Corrector Method for Nonlinear Fractional Differential Equation Initial Value Problems on Graded Grids

摘要: 非线性分数阶微分方程初值问题的典型解在初始时刻具有弱奇异性，为了更好地逼近它的典型解，本文将Nguyen和Jang在均匀网格上的三阶预测校正方法推广到非均匀网格上。与文献中现有的预测校正方法相比，这种新方法在降低计算成本的同时，显著提高了数值精度。另外，通过数值算例得到的结果显示，当选择合适的网格参数时，误差的收敛阶大于3。

Abstract: In order to better approximate the classical solution of nonlinear fractional differential equation in-itial value problems, which has weak singularity at the initial time, the third order predictor- cor-rector method of Nguyen and Jang on uniform grids is extended to non-uniform grids. Compared with the existing predictor-corrector methods in the literature, this new method reduces the com-putation cost and significantly improves the numerical accuracy. In addition, numerical examples show that the convergence order of errors is greater than the third order when appropriate mesh parameters are selected.

1. 引言

近年来，由于分数阶微分方程在物理学和工程学中的各种应用，引起了研究人员的广泛关注。例如，学者普遍感兴趣的分数阶模型包括分数阶反常扩散模型 [1] - [6] 、分数阶水波模型 [7] 、分数阶Cable模型 [8] 、分数阶Allen-Cahn方程 [9] 、分数阶Bloch-Torrey方程 [10] [11] 、分数阶扩散波动方程 [12] [13] [14] [15] [16] 等。由于分数阶导数本身具有的历史性和非局部特性，所以寻找其精度较高的数值算法变得尤为重要。

Adams方法具有稳定性好、计算量合理、易于实现等优点，是常微分方程数值解的一种常用的离散方法。为了将这些方法推广到分数阶微分方程，已经做了大量的工作。Diethelm和他的合著者 [17] [18] [19] [20] 研究了一种使用Adam公式的预测–校正方法，并对均匀网格进行了详细的误差分析。随后，Liu，Roberts和Yan将这种预测校正方法推广到非均匀网格 [21] 。

对于非线性分数阶常微分方程初值问题，Nguyen和Jang [22] 提出了均匀网格上的预测–校正法，该方法以较低的计算成本获得了较高的精度。但是该类问题的主要困难之处是典型解在初始时具有弱奇异性，为了解决数值解中的这一困难，Zhou和Li [23] 将 [22] 的Adams型预测–校正方法推广到非均匀网格上，在近似方法中只采用了线性插值。我们在这些基础上，将采用二次插值的Adams型预测校正方法推广到了非均匀网格上，通过数值算例得到的结果显示，当选择合适的网格参数时，误差的收敛阶大于3。

2. 主要内容

在本篇文章中，我们要研究的分数阶微分方程如下

$\{\begin{array}{l}{\text{D}}_0^{\alpha }y\left(t\right)=f\left(t,y\left(t\right)\right)\text{,}0<t\le T,\hfill \\ y\left(0\right)=y_0.\hfill \end{array}$ (1)

其中 $\alpha \in \left(0,1\right)$ ， ${\text{D}}_0^{\alpha }y\left(t\right)$ 是 $\alpha$ 阶Caputo导数，即

${\text{D}}_0^{\alpha }y\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{-\alpha }{y}^{\prime }\left(\tau \right)\text{d}\tau }$ 。(2)

众所周知，一个连续函数 $y\left(t\right)$ 是方程(1)的解当且仅当它是Volterra积分方程 [24] 的解，所以我们有

$y\left(t\right)=y_0+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}f\left(\tau ,y\left(\tau \right)\right)\text{d}\tau }$ 。 (3)

这里的Volterra积分方程的理论分析和数值方法可以参考 [24] 和 [25] 等文献。

在本研究中，我们采取的是以下网格：N是任意正整数，将 $\left[0,T\right]$ 区间划分为N个子区间，网格点

$t_n:=T{\left(\frac{n}{N}\right)}^r$ ， $0\le n\le N$ ，r是待选择的网格分级参数。

在网格点 $t=t_{n+1}$ 处，可改写(3)为

$\begin{array}{c}y\left(t_{n+1}\right)=y_0+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{j=0}{\overset{n-1}{\sum }}{\displaystyle {\int }_{t_j}^{t_{j+1}}{\left(t-\tau \right)}^{\alpha -1}f\left(\tau ,y\left(\tau \right)\right)\text{d}\tau }}+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_{t_n}^{t_{n+1}}{\left(t-\tau \right)}^{\alpha -1}f\left(\tau ,y\left(\tau \right)\right)\text{d}\tau }\\ \triangleq y_0+F_h\left(t_{n+1}\right)+F_l\left(t_{n+1}\right)\end{array}$ (4)

这里的 $F_h\left(t_{n+1}\right)$ 和 $F_l\left(t_{n+1}\right)$ 分别表示Volterra积分算子的历史部分和局部部分。

本篇文章中将在分级网格上来讨论二次插值的预测校正方法，通过与均匀网格的对比体现了本方案的优越性，具体操作在第3节中会作详细介绍。

3. 二次插值的预测校正方法

为方便起见，定义 $y_j=y\left(t_j\right)$ 为解在 $t_j$ 处的精确解， $j=0,\cdots ,N$ 。定义 ${\tilde{y}}_{j+1}^P$ ， ${\tilde{y}}_{j+1}\approx y_{j+1}$ 分别是 $y_{j+1}$ 的预测校正值和近似校正值。类似地可定义 $f_j=f\left(t_j,y_j\right)$ ， ${\tilde{f}}_j=\tilde{f}\left(t_j,y_j\right)$ 以及 ${\tilde{f}}^P{}_j=f\left(t_j,{\tilde{y}}^P{}_j\right)$ 。Volterra积分算子的历史部分 $F_h\left(t_{n+1}\right)$ 和局部部分 $F_l\left(t_{n+1}\right)$ 分别用 ${\tilde{F}}_h\left(t_{n+1}\right)$ 与 ${\tilde{F}}_l\left(t_{n+1}\right)$ 来近似。

在每个区间 $\left[t_{j-1},t_{j+1}\right]$ 上，这里 $j\ge 1$ ，我们用二次拉格朗日插值来逼近 $f\left(t\right)$ ，即：

$f\left(t\right)\approx f_{j-1}\frac{\left(t-t_j\right)\left(t-t_{j+1}\right)}{\left(t_{j-1}-t_j\right)\left(t_{j-1}-t_{j+1}\right)}+f_j\frac{\left(t-t_{j-1}\right)\left(t-t_{j+1}\right)}{\left(t_j-t_{j-1}\right)\left(t_j-t_{j+1}\right)}+f_{j+1}\frac{\left(t-t_{j-1}\right)\left(t-t_j\right)}{\left(t_{j+1}-t_{j-1}\right)\left(t_{j+1}-t_j\right)}$ (5)

再来处理区间 $\left[t_0,t_1\right]$ ， $f\left(t\right)$ 在该区间上可用

$f\left(t\right)\approx f_0\frac{\left(t-t_{1/2}\right)\left(t-t_1\right)}{\left(t_0-t_{1/2}\right)\left(t_0-t_1\right)}+f_{1/2}\frac{\left(t-t_0\right)\left(t-t_1\right)}{\left(t_{1/2}-t_0\right)\left(t_{1/2}-t_1\right)}+f_1\frac{\left(t-t_0\right)\left(t-t_{1/2}\right)}{\left(t_1-t_0\right)\left(t_1-t_{1/2}\right)}$ (6)

来近似。将上述两式代入(4)式中，我们得到

${\tilde{y}}_{n+1}=y_0+{\tilde{F}}_h\left(t_{n+1}\right)+{\tilde{F}}_l\left(t_{n+1}\right)$ 。 (7)

这里的记忆部分

${\tilde{F}}_h\left(t_{n+1}\right)=\frac{1}{\Gamma \left(\alpha \right)}\left[A_{n+1}^{0,0}{\tilde{f}}_0+A_{n+1}^{1,0}{\tilde{f}}_{1/2}+A_{n+1}^{2,0}{\tilde{f}}_1\right]+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}\left[A_{n+1}^{0,j}{\tilde{f}}_{j-1}+A_{n+1}^{1,j}{\tilde{f}}_j+A_{n+1}^{2,0}{\tilde{f}}_{j+1}\right]}$ (8)

以及局部部分

${\tilde{F}}_l\left(t_{n+1}\right)=\frac{1}{\Gamma \left(\alpha \right)}\left[A_{n+1}^{0,n}{\tilde{f}}_{n-1}+A_{n+1}^{1,n}{\tilde{f}}_n+A_{n+1}^{2,n}{\tilde{f}}_{n+1}^P\right]$ 。 (9)

其中的预测近似值 ${\tilde{y}}_{n+1}^P$ 用(10)式来近似，即在区间 $\left[t_n,t_{n+1}\right]$ 上用 ${\tilde{f}}_{n-2}$ ， ${\tilde{f}}_{n-1}$ 以及 ${\tilde{f}}_n$ 的二次插值来逼近 $f\left(t\right)$ 。

${\tilde{y}}_{n+1}^P=y_0+{\tilde{F}}_h\left(t_{n+1}\right)+\frac{1}{\Gamma \left(\alpha \right)}\left[b_{n+1}^{0,n}{\tilde{f}}_{n-2}+b_{n+1}^{1,n}{\tilde{f}}_{n-1}+b_{n+1}^{2,n}{\tilde{f}}_n\right]$ (10)

下面是需要计算的系数部分，计算公式如下

$A_{n+1}^{0,0}={\displaystyle {\int }_{t_0}^{t_1}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_{1/2}\right)\left(t-t_1\right)}{\left(t_0-t_{1/2}\right)\left(t_0-t_1\right)}\text{d}t}$ ， $A_{n+1}^{1,0}={\displaystyle {\int }_{t_0}^{t_1}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_0\right)\left(t-t_1\right)}{\left(t_{1/2}-t_0\right)\left(t_{1/2}-t_1\right)}\text{d}t}$ ，

$A_{n+1}^{2,0}={\displaystyle {\int }_{t_0}^{t_1}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_0\right)\left(t-t_{1/2}\right)}{\left(t_1-t_0\right)\left(t_1-t_{1/2}\right)}\text{d}t}$ ， $b_{n+1}^{0,n}={\displaystyle {\int }_{t_n}^{t_{n+1}}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_{n-1}\right)\left(t-t_n\right)}{\left(t_{n-2}-t_{n-1}\right)\left(t_{n-2}-t_n\right)}\text{d}t}$ ，

$b_{n+1}^{1,n}={\displaystyle {\int }_{t_n}^{t_{n+1}}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_{n-2}\right)\left(t-t_n\right)}{\left(t_{n-1}-t_{n-2}\right)\left(t_{n-1}-t_n\right)}\text{d}t}$ ， $b_{n+1}^{0,n}={\displaystyle {\int }_{t_n}^{t_{n+1}}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_{n-1}\right)\left(t-t_n\right)}{\left(t_{n-2}-t_{n-1}\right)\left(t_{n-2}-t_n\right)}\text{d}t}$ 。

其余需计算的系数如下(其中 $j=1,\cdots ,n$ )

$A_{n+1}^{0,j}={\displaystyle {\int }_{t_j}^{t_{j+1}}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_j\right)\left(t-t_{j+1}\right)}{\left(t_{j-1}-t_j\right)\left(t_{j-1}-t_{j+1}\right)}\text{d}t}$ ， $A_{n+1}^{0,j}={\displaystyle {\int }_{t_j}^{t_{j+1}}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_{j-1}\right)\left(t-t_{j+1}\right)}{\left(t_j-t_{j-1}\right)\left(t_j-t_{j+1}\right)}\text{d}t}$ ，

$A_{n+1}^{2,j}={\displaystyle {\int }_{t_j}^{t_{j+1}}{\left(t_{n+1}-t\right)}^{\alpha -1}\frac{\left(t-t_{j-1}\right)\left(t-t_j\right)}{\left(t_{j+1}-t_{j-1}\right)\left(t_{j+1}-t_j\right)}\text{d}t}$ 。

上述过程可归结为以下：

$\{\begin{array}{l}{\tilde{y}}_{n+1}=y_0+{\tilde{F}}_h\left(t_{n+1}\right)+{\tilde{F}}_l\left(t_{n+1}\right)\hfill \\ {\tilde{y}}_{n+1}^P=y_0+{\tilde{F}}_h\left(t_{n+1}\right)+\frac{1}{\Gamma \left(\alpha \right)}\left[b_{n+1}^{0,n}{\tilde{f}}_{n-2}+b_{n+1}^{1,n}{\tilde{f}}_{n-1}+b_{n+1}^{2,n}{\tilde{f}}_n\right]\hfill \\ {\tilde{F}}_l\left(t_{n+1}\right)=\frac{1}{\Gamma \left(\alpha \right)}\left[A_{n+1}^{0,n}{\tilde{f}}_{n-1}+A_{n+1}^{1,n}{\tilde{f}}_n+A_{n+1}^{2,n}{\tilde{f}}_{n+1}^P\right]\hfill \\ {\tilde{F}}_h\left(t_{n+1}\right)=\frac{1}{\Gamma \left(\alpha \right)}\left[A_{n+1}^{0,0}{\tilde{f}}_0+A_{n+1}^{1,0}{\tilde{f}}_{1/2}+A_{n+1}^{2,0}{\tilde{f}}_1\right]+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}\left[A_{n+1}^{0,j}{\tilde{f}}_{j-1}+A_{n+1}^{1,j}{\tilde{f}}_j+A_{n+1}^{2,0}{\tilde{f}}_{j+1}\right]}\hfill \end{array}$ (11)

另外，在计算的循环过程启动之前，我们还须计算 ${\tilde{y}}_1$ ， ${\tilde{y}}_2$ 这两个值，具体的计算过程如下：

1) ${\mathcal{L}}^0\left({\tilde{f}}_0\right)\to {\tilde{y}}_{1/4}^P\to {\mathcal{L}}^1\left({\tilde{f}}_0,{\tilde{f}}_{1/4}^P\right)\to {\tilde{y}}_{1/4}$

2) ${\mathcal{L}}^1\left({\tilde{f}}_0,{\tilde{f}}_{1/4}\right)\to {\tilde{F}}_h\left(t_{1/2}\right);{\mathcal{L}}^0\left({\tilde{f}}_{1/4}\right)\to {\tilde{y}}_{1/2}^{P_1}\to {\mathcal{L}}^1\left({\tilde{f}}_{1/4},{\tilde{f}}_{1/2}^{P_1}\right)\to {\tilde{y}}_{1/2}^{P_2}\to {\mathcal{L}}^2\left({\tilde{f}}_0,{\tilde{f}}_{1/4},{\tilde{f}}_{1/2}^{P_2}\right)\to {\tilde{y}}_{1/2}$

3) ${\mathcal{L}}^1\left({\tilde{f}}_0,{\tilde{f}}_{1/2}\right)\to {\tilde{F}}_h\left(t_1\right);{\mathcal{L}}^0\left({\tilde{f}}_{1/2}\right)\to {\tilde{y}}_1^{P_1}\to {\mathcal{L}}^1\left({\tilde{f}}_{1/2},{\tilde{f}}_1^{P_1}\right)\to {\tilde{y}}_1^{P_2}\to {\mathcal{L}}^2\left({\tilde{f}}_0,{\tilde{f}}_{1/2},{\tilde{f}}_1^{P_2}\right)\to {\tilde{y}}_1$

4) ${\mathcal{L}}^2\left({\tilde{f}}_0,{\tilde{f}}_{1/2},{\tilde{f}}_1\right)\to {\tilde{F}}_h\left(t_2\right);{\mathcal{L}}^0\left({\tilde{f}}_1\right)\to {\tilde{y}}_2^{P_1}\to {\mathcal{L}}^1\left({\tilde{f}}_1,{\tilde{f}}_2^{P_1}\right)\to {\tilde{y}}_2^{P_2}\to {\mathcal{L}}^2\left({\tilde{f}}_0,{\tilde{f}}_1,{\tilde{f}}_2^{P_2}\right)\to {\tilde{y}}_2$

其中 ${\mathcal{L}}^0\left(f_a\right)$ 为f在 $t=a$ 处的常数插值即 ${\mathcal{L}}^0\left(f_a\right)=f\left(a\right)$ ， ${\mathcal{L}}^1\left(f_a,f_b\right)$ 为f在网格 $\left(a,b\right)$ 处的线性插值， ${\mathcal{L}}^2\left(f_a,f_b,f_c\right)$ 为f在网格 $\left(a,b,c\right)$ 处的二次插值，此启动方案参照文献 [22] 。

4. 数值实验

在本小节中，我们给出了一个数值例子来说明均匀网格上与非均匀网格上的二次插值预测校正方法的收敛阶数。

例1考虑如下分数阶微分方程

$\{\begin{array}{l}\begin{array}{cc}{D}_0^{\alpha }y\left(t\right)=\frac{\Gamma \left(4+\alpha \right)}{6}{t}^3+t^{3+\alpha }-y\left(t\right),& 0<t\le T\end{array}\hfill \\ y\left(0\right)=0\hfill \end{array}$

该问题的精确解为： $y\left(t\right)=t^{3+\alpha }$ 。

图1是当 $\alpha =0.4,r=\frac{1+\alpha }{2\alpha }$ 时，例1的精确解与数值解的图像。由图可以看出，我们采取的方法所算

出的值几乎与精确解的函数值吻合，证明了该方法的有效性。

Figure 1. Image of the exact and numerical solution of example 1 when N = 20

图1. 当N = 20时例1的精确解与数值解的图像

表1表示用均匀网格上与分级网格上的二次插值预测校正方法来计算例1的最大节点误差及收敛阶

(其中 $\alpha =0.4$ )。另外 $er{r}_N:=\underset{0\le j\le N}{\max }\left|y\left(t_j\right)-y_j\right|$ ， $p_N:={\log }_2\left(er{r}_N/er{r}_{2N}\right)$ 。

Table 1. The maximum node error and convergence order of example 1

表1. 例1的最大节点误差及收敛阶

例2考虑如下非线性分数阶微分方程

$\{\begin{array}{l}\begin{array}{cc}{D}_0^{\alpha }y\left(t\right)=\frac{\Gamma \left(4+\alpha \right)}{6}{t}^3+t^{6+2\alpha }-y^2\left(t\right),& 0<t\le T\end{array}\hfill \\ y\left(0\right)=0\hfill \end{array}$

该问题的精确解为： $y\left(t\right)=t^{3+\alpha }$ 。

同样地，我们画出例2的精确解和数值解图像如图2。可以看出图像基本吻合，这表明我们的方法可以去计算这样的非线性分数阶微分方程初值问题。

由表1、表2，当 $r=1$ 时即取均匀网格时， $p_N$ 随着N的增大却在减小，对比 $r=\frac{1+\alpha }{2\alpha }$ 时， $p_N$ 随着N

的增大而增大，并且数值略大于均匀网格时的情形。因此，非均匀网格的收敛阶优于均匀网格的收敛阶，这个数值结果证明了方法的优越性。

Figure 2. Image of the exact and numerical solution of example 2 when N = 20

图2. 当N = 20时例2的精确解与数值解的图像

Table 2. The maximum node error and convergence order of example 2

表2. 例2的最大节点误差及收敛阶

5. 结论

分数阶常微分方程在数学模型中的应用受到越来越广泛的关注。不同的分数阶常微分方程模型已被应用到越来越多的领域中。并且发现它在研究一些具有记忆过程、遗传性质以及异质材料时比整数阶方程模型更有优势。本文研究的问题是非线性分数阶初值问题的数值解，我们将Nguyen和Jang的三阶预测校正方法推广到非均匀网格上，主要是在非均匀网格上进行二次插值，然后用预测校正法进行求解，数值实验表明误差阶超过三阶，我们接下来的工作将对本文所提出的数值方法进行误差分析并用于求解分数阶偏微分方程。