1 问题的提出

分数阶微分方程最重要的是非局部性质, 它能有效地避免整数阶导数的局部性, 所以，分数阶导数可以更加精确地描述对历史有依赖的问题。这种非局部性质使得分数阶微分方程在许多工程和科学领域，特别是在描述具有记忆过程、遗传特性^[1-3]和异质材料等方面越来越有吸引力。分数阶微分方程已被广泛应用于各种复杂的黏性流体流动模型，如异常扩散过程^[4]、分数动力学^[5]、粘弹性材料^[6]、输送管道边界层效应^[7]、氧气通过毛细血管输送到组织^[8]等。同时，利用分数阶导数描述问题时，使用少量的参数就能获得与实际相吻合的结果，所以，与整数阶微分方程相比较，在处理复杂系统时，使用分数阶微分方程建模将更加简单、参数意义更加清楚。一般来说，对分数阶微分方程的求解是非常困难的。近几年，学者们对分数阶积分微分方程的基本理论及其数值求解进行了大量的研究。因为，一方面分数阶微分方程的求解可以转化为等价的积分方程进行求解；另一方面为使得数学模型更加符合客观实际，对分数阶积分微分方程的求解已经不可避免，所以，对分数阶积分微分方程的理论及其数值方法的研究越来越受到学者们的关注。本文提出了一种求解以下非线性分数阶初值问题的有效数值方法。

$ \left\{ {\begin{array}{*{20}{l}} {{{\rm{D}}^{\alpha + 1}}u = f(t,u(t)),\quad t \in [0,T]} \\ {u(0) = 0,\quad {u'}(0) = 0} \end{array}} \right. $

(1)

其中， $ \alpha \in (0,1),f \in C(0,T) $ , 并且 ${{\rm{D}}^\alpha }$ 是根据下面定义的 $\alpha $ 阶Caputo导数。

$ \begin{split}& {{\rm{D}}^\alpha }u = \frac{1}{{\Gamma (n - \alpha )}}\int_0^t {{{(t - \tau )}^{n - \alpha - 1}}} {u^{(n)}}(\tau ){\text{d}}\tau, \\& \quad {{t}} > 0,\quad {{n}} - 1 < \alpha < {{n}}{} \end{split}$

Caputo导数 ${{\rm{D}}^\alpha }$ 不同于传统的微分算子，它是一种全局算子。目前分数阶导数主要分为3种类型：Riesz 意义下的分数阶导数、Riemann-Liouville 分数阶导数和Caputo 型分数阶导数。1832年Liouville采用级数的形式给出了分数阶导数，1853 年Riemann 采用定积分的形式给出了另一种分数阶微分，Grünwald 和Krug 将Riemann 和Liouville 的结论进行了统一，得到了Riemann-Liouville 分数阶导数的定义。1967 年Caputo 提出了Caputo 型分数阶导数。Caputo 型分数阶导数因其广泛的应用背景和相对简便的运算被广泛地应用于各个领域。目前在分数阶微分方程理论的研究中，关于其数值解进行了大量的研究。在Caputo 型分数阶定义下的初值条件具有明确的物理意义，因而在工程物理学领域该方程具有重要的应用价值。但由于解的复杂性，不能直接应用于实际工程中，所以，求解Caputo 型分数阶微分方程的数值解受到广泛关注。

谱方法作为一种全球性的高精度数值方法，在求解微分方程过程中扮演着重要的角色，非常适用于求解分数阶微分方程。然而，现有的谱方法主要是一步法。这种方法要求解具有高度正则性。然而，通常对于分数阶问题，解不可能是任意光滑的，即使源项非常光滑。因此，一步法缺乏有效处理局部弱奇异的能力。最近，Guo等^[9]提出了一种求解单分数阶边值问题的hp型Chebyshev谱配置方法，并从理论上和数值上得到了hp收敛性。在他们工作的基础上，本文提出了一种求解非线性分数阶初值问题的hp型Legendre谱配置法，该方法灵活地采用时间步长和局部逼近阶数，可以实现任意精度的数值解。并且，本文所提出的数值格式对于数值实现和误差分析来说更加简单易行。

2 预备知识

设 ${I_h}$ 是区间 $I = \left[ {0,T} \right]$ 上的一个网格, ${I_h}: = \{ {t_n}: 0 = {t_0} < {t_1} < \cdots < {t_N} = T \}$ ，并且设 $ {h_n} = {t_n} - {t_{n - 1}} $ ， ${I_n} = ({t_{n - 1}},{t_n}]$ 。通常情况下，对于一个给定区间 $\varLambda $ 和某一权函数 $\chi (x)$ ，定义

$ {L}_{\chi }^{2}(\varLambda )=\left\{v|v是可测的并且{\Vert v\Vert }_{{L}_{\chi }^{2}(\varLambda )} < \infty \right\} $

其中， ${\left\| v \right\|_{L_\chi ^2(\varLambda )}} = {\left( { \displaystyle\int_\varLambda {{{\left| {v(x)} \right|}^2}\chi (x){\rm{d}}x} } \right)^{\frac{1}{2}}},\,H_\chi ^r(\varLambda)$ 是带权Sobolev空间。

设 ${L_k}(x),x \in ( - 1,1)$ 是标准的阶数为 $k$ 的Legendre多项式。Legendre 多项式集合构成一个完备的 ${L^2}( - 1,1)$ 正交系统，即

$ \int_{ - 1}^1 {{L_k}} (x){L_j}(x){\rm{d}}x = \frac{2}{{2k + 1}}{\delta _{k,j}} $

(2)

其中， $ {\delta _{k,j}} $ 是Kronecker符号。

定义 $k$ 阶变换的Legendre多项式

$ {L_{n,k}}(t) = {L_k}\left( {\frac{{2t - {t_{n - 1}} - {t_n}}}{{{h_n}}}} \right),\quad t \in {I_n},\quad k \geqslant 0 $

(3)

${L_{n,k}}(t),k \geqslant 0$ 的集合构成一个完备的 ${L^2}\left( {{I_n}} \right)$ 正交系统，即

$ \int_{{I_n}} {{L_{n,k}}} (t){L_{n,j}}(t){\rm{d}}t = \frac{{{h_n}}}{{2k + 1}}{\delta _{k,j}} $

(4)

现在来看Legendre-Gauss插值。对于任意整数 ${M_n} \geqslant 0$ ，采用 $\left\{ {{x_{n,j}},{\omega _{n,j}}} \right\}_{j = 0}^{{M_n}}$ 表示在区间 ${I_n}$ 上标准的Legendre-Gauss 插值的节点和对应的权函数值。设 ${\mathcal{P}_{{M_n}}}\left( {{I_n}} \right)$ 是区间 ${I_n}$ 上阶数至多为 ${M_n}$ 的多项式集合， $ {\left\{{t}_{n,j}\right\}}_{j=0}^{{M}_{n}}  $ 是区间 ${I_n}$ 上变换的Legendre-Gauss 求积节点，

$ {t_{n,j}} = \frac{1}{2}\left( {{h_n}{x_{n,j}} + {t_{n - 1}} + {t_n}} \right),\quad 1 \leqslant n \leqslant N,\quad 0 \leqslant j \leqslant {M_n} $

根据标准Legendre-Gauss 求积的性质，可得

$ \int_{{I_n}} \phi (t){\rm{d}}t = \frac{{{h_n}}}{2}\sum\limits_{j = 0}^{{M_n}} \phi \left( {{t_{n,j}}} \right){\omega _{n,j}},\quad \phi \in {\mathcal{P}_{2{M_n} + 1}}\left( {{I_n}} \right) $

(5)

和

$\begin{split}& \sum\limits_{j = 0}^{{M_n}} {{L_{n,p}}} \left( {{t_{n,j}}} \right){L_{n,q}}\left( {{t_{n,j}}} \right){\omega _{n,j}} = \frac{2}{{2p + 1}}{\delta _{p,q}},\\& \qquad 0 \leqslant p + q \leqslant 2{M_n} + 1 \end{split}$

(6)

采用 ${\mathcal{I}_{s,{M_n}}}:C\left( {{t_{n - 1}},{t_n}} \right) \to {\mathcal{P}_{{M_n}}}\left( {{t_{n - 1}},{t_n}} \right)$ 表示对 ${{s}}$ 的变换的Legendre-Gauss插值算子, 使得

$ {\mathcal{I}_{s,{M_n}}}v\left( {{t_{n,j}}} \right) = v\left( {{t_{n,j}}} \right),\quad 0 \leqslant j \leqslant {M_n} $

(7)

为便于使用, 引入另一个插值算子 $\mathcal{I}_{\xi ,{M_n}}^t: C\left( {{t_{n - 1}},t} \right) \to {\mathcal{P}_{{M_n}}}\left( {{t_{n - 1}},t} \right)$ ，定义为

$ \mathcal{I}_{\xi ,{M_n}}^tv\left( {{\xi _{n,j}}} \right) = v\left( {{\xi _{n,j}}} \right),\quad 0 \leqslant j \leqslant {M_n} $

(8)

其中，

$ {\xi _{n,j}} = \sigma \left( {{t_{n,j}},t} \right): = {t_{n - 1}} + \frac{{\left( {t - {t_{n - 1}}} \right)\left( {{t_{n,j}} - {t_{n - 1}}} \right)}}{{{h_n}}} $

注意到

$ \mathcal{I}_{\xi ,{M_n}}^tv(\xi ) = {\mathcal{I}_{\lambda ,{M_n}}}v(\sigma (\lambda ,t)) $

正如文献[10]中提到的，由于弱奇点的存在 $({\rm{e.g.}},\;\nu \in (0,1))$ ，考虑到加权插值求积公式的权值取决于核中的弱奇异因子 ${(t - s)^\nu }$ 是很自然的, 对于任意的 $\phi (s) \in {\mathcal{P}_{{M_k}}}\left( {{I_k}} \right)$ , 引入加权插值求积公式

$ \begin{split} & \int_{{I_k}} {{{(t - s)}^\nu }} \phi (s){\rm{d}}s = \sum\limits_{j = 0}^{{M_k}} \phi \left( {{t_{k,j}}} \right)\tilde \omega _{k,j}^v(t),\\&\qquad t \in {I_n},\quad k < n \end{split} $

(9)

其中，

$ \tilde \omega _{k,j}^v(t) = \int_{{I_k}} {{{(t - s)}^\nu }} {l_{k,j}}(s){\rm{d}}s $

(10)

这里 $\left\{ {{l_{k,j}}(s)} \right\}_{j = 0}^{{M_k}}$ 是对应插值点 $\left\{ {{t_{k,j}}} \right\}_{j = 0}^{{M_k}}$ 的 Lagrange 基多项式。函数 $ \tilde \omega _{k,j}^v(t) $ 可以根据Lagrange 多项式的性质精确计算^[11]。

对于任意的 $\psi (\xi ) \in {\mathcal{P}_{{M_n}}}\left( {\left[ {{t_{n - 1}},t} \right]} \right)$ ，引入加权插值求积公式

$ \int_{{t_{n - 1}}}^t {{{(t - \xi )}^\nu }} \psi (\xi )d\xi = \sum\limits_{j = 0}^{{M_n}} \psi \left( {{\xi _{n,j}}} \right)\hat \omega _{n,j}^\nu (t),\quad t \in {I_n} $

(11)

其中，

$ \hat \omega _{n,j}^\nu (t) = \int_{{t_{n - 1}}}^t {{{(t - \xi )}^\nu }} l_{{{n}},j}^t(\xi ){\rm{d}}\xi $

(12)

这里 $\left\{ {l_{n,j}^t(\xi )} \right\}_{j = 0}^{{M_n}}$ 是对应插值点 $\left\{ {{\xi _{n,j}}} \right\}_{j = 0}^{{M_n}}$ 的 Lagrange 基多项式。式(11)中的函数 $\hat \omega _{n,j}^\nu (t)$ 可以根据带权 ${(t - \xi )^\nu }$ Jacobi-Gauss求积公式精确计算。

3 hp型谱配置法

方程(1)可以转换为以下等价的Volterra积分方程^[12]：

$ u(t) = \frac{1}{{\Gamma (1 + \alpha )}}\int_0^t {{{(t - s)}^\alpha }} f(s,u(s)){\rm{d}}s $

(13)

本文的主要目的是基于等价式(13) 提出并分析一种有效的求解方程(1)的 ${\rm{hp}}$ 型Legendre-Gauss谱配置法。

设 ${u^n}(t)$ 是第 $n$ 个区间上式(13)的解，即

$ {u^n}(t): = u(t),\quad t \in {I_n},\quad 1 \leqslant n \leqslant N $

由式(13)可得，对于任意 $t \in {I_n}$ ，

$\begin{split}& {u^n}(t) = \frac{1}{{\Gamma (1 + \alpha )}}\Bigg( \sum\limits_{k = 1}^{n - 1} {\int_{{I_k}} {{{(t - s)}^\alpha }} } f\left( {s,{u^k}(s)} \right){\rm{d}}s +\Bigg.\\&\quad \Bigg. \int_{{t_{n - 1}}}^t {{{(t - \xi )}^\alpha }} f\left( {\xi ,{u^n}(\xi )} \right){\rm{d}}\xi \Bigg) \end{split} $

(14)

求解式(14)的hp型谱配置法的目的是寻找 ${U^n}(t) \in {\mathcal{P}_{{M_n}}}\left( {{I_n}} \right)$ ，使得

$\begin{split}& {U^n}(t) =\\& \quad{{\cal I}_{t,{M_n}}}\left[ {\frac{1}{{\Gamma (\alpha + 1)}}\left( {\sum\limits_{k = 1}^{n - 1} {\int_{{I_k}} {{{(t - s)}^\alpha }} } {{\cal I}_{s,{M_k}}}f\left( {s,{U^k}(s)} \right){\rm{d}}s} \right.} \right. +\\& \quad \left. {\left. { \int_{{t_{n - 1}}}^t {{{(t - \xi )}^\alpha }} {\cal I}_{\xi ,{M_n}}^tf\left( {\xi ,{U^n}(\xi )} \right){\rm{d}}\xi } \right)} \right]\\[-12pt] \end{split} $

(15)

其中， ${U^k}(t) \in {\mathcal{P}_{{M_k}}}\left( {{I_k}} \right)$ 是区间 ${I_k}$ 上 ${u^k}(t)$ 的数值解。用 $U$ 表示全局数值解, 即

$ U(t) = {U^k}(t),\quad t \in {I_k},\quad 1 \leqslant k \leqslant N $

现在描述数值实现, 并给出了式(15)的算法。设

$ {U^n}(t) = \sum\limits_{p = 0}^{{M_n}} {u_p^n} {L_{n,p}}(t) $

(16)

由式(6)和式(9)可得，对于 $1 \leqslant k \leqslant n - 1$ ,

$\begin{split}& {\mathcal{I}_{t,{M_n}}}\int_{{I_k}} {{{(t - s)}^\alpha }} {\mathcal{I}_{s,{M_k}}}f\left( {s,{U^k}(s)} \right){\rm{d}}s = \\& \quad {\mathcal{I}_{t,{M_n}}}\sum\limits_{{p'} = 0}^{{M_k}} f \left( {{t_{k,{p'}}},{U^k}\left( {{t_{k,{p'}}}} \right)} \right)\tilde \omega _{k,{{\text{p}}'}}^\alpha (t) = \\& \quad \sum\limits_{{p'} = 0}^{{M_k}} f \left( {{t_{k,{p'}}},{U^k}\left( {{t_{k,{p'}}}} \right)} \right){\mathcal{I}_{t,{M_n}}}\tilde \omega _{k,{{\text{p}}'}}^\alpha (t) = \\& \quad \sum\limits_{{p'} = 0}^{{M_k}} f \left( {{t_{k,{p'}}},{U^k}\left( {{t_{k,{p'}}}} \right)} \right)\sum\limits_{p = 0}^{{M_n}} {\omega _{k,{p'},p}^{n,\alpha }} {L_{n,p}}(t) = \\& \quad \sum\limits_{p = 0}^{{M_n}} {a_{k,p}^n} {L_{n,p}}(t) \\[-8pt] \end{split} $

(17)

其中，

$ \begin{split}& a_{k,p}^n = \sum\limits_{{p'} = 0}^{{M_k}} f \left( {{t_{k,{p'}}},{U^k}\left( {{t_{k,{p'}}}} \right)} \right)\omega _{k,{p'},p}^{n,\alpha } , \omega _{k,{p^\prime },p}^{n,\alpha } = \\& \qquad \frac{{2p + 1}}{2}\sum\limits_{j = 0}^{{M_k}} {\tilde \omega _{k,{{\text{p}}'}}^\alpha } \left( {{t_{n,j}}} \right){L_{n,p}}\left( {{t_{n,j}}} \right){\omega _{n,j}} \end{split} $

(18)

由式(6)，(8)和(11)可得

$ \begin{split}& {\mathcal{I}_{t,{M_n}}}\left( {\int_{{t_{n - 1}}}^t {{{(t - \xi )}^\alpha }} \mathcal{I}_{\xi ,{M_n}}^tf\left( {\xi ,{U^n}(\xi )} \right){\rm{d}}\xi } \right) = \\& \qquad {\mathcal{I}_{t,{M_n}}}\sum\limits_{j = 0}^{{M_n}} f \left( {{\xi _{n,j}},{U^n}\left( {{\xi _{n,j}}} \right)} \right)\hat \omega _{n,j}^\alpha (t) =\\& \qquad \sum\limits_{p = 0}^{{M_n}} {b_p^n} {L_{n,p}}(t) \\[-8pt] \end{split} $

(19)

其中，

$ \begin{split} b_p^n =&\frac{{2p + 1}}{2}\cdot\\& \sum\limits_{i,j = 0}^{{M_n}} f \left( {{\xi _{n,j}},{U^n}\left( {{\xi _{n,j}}} \right)} \right)\hat \omega _{n,j}^\alpha \left( {{t_{n,i}}} \right){L_{n,p}}\left( {{t_{n,i}}} \right){\omega _{n,i}} \end{split}$

(20)

由式(15)～(20)可得

$ \begin{split}& \sum\limits_{p = 0}^{{M_n}} {u_p^n} {L_{n,p}}(t) = \frac{1}{{\Gamma (\alpha + 1)}}\sum\limits_{p = 0}^{{M_n}} {\sum\limits_{k = 1}^{n - 1} {a_{k,p}^n} } {L_{n,p}}(t) +\\& \qquad \quad \frac{1}{{\Gamma (\alpha + 1)}}\sum\limits_{p = 0}^{{M_n}} {b_p^n} {L_{n,p}}(t) \end{split} $

(21)

根据 Legendre 多项式的正交性，比较式(11)的展开系数得到，对于 $0 \leqslant p \leqslant {M_n}$ ，

$ u_p^n = \frac{1}{{\Gamma (\alpha + 1)}}\left( {\sum\limits_{k = 1}^{n - 1} {a_{k,p}^n} + b_p^n} \right) $

(22)

可采用迭代法计算系统式（22）的展开系数（例如，Newton-Raphson 迭代）。

4 数值结果

现提供一系列数值实验来验证所提出算法的有效性。首先，定义在网格点上离散的 ${L^2}$ 误差和 ${L^\infty }$ 误差如下：

$ \begin{gathered} {E_{{L^2}(0,T)}}: = {\left( {\sum\limits_{k = 1}^N {\frac{{{h_k}}}{2}} \sum\limits_{j = 0}^{{M_k}} {{{\left( {u\left( {{t_{k,j}}} \right) - U\left( {{t_{k,j}}} \right)} \right)}^2}} {\omega _{k,j}}} \right)^{\frac{1}{2}}} \\ {E_{{L^\infty }(0,T)}}: = \mathop {\max }\limits_{_{1 \leqslant k \leqslant N}} \left( {\mathop {\max }\limits_{_{0 \leqslant j \leqslant {M_k}}} \left| {u\left( {{t_{k,j}}} \right) - U\left( {{t_{k,j}}} \right)} \right|} \right) \\ \end{gathered} $

例1　考虑问题

$ \left\{ {\begin{array}{*{20}{l}} {{{\rm{D}}^{1.5}}y(t) + y(t) = {t^\alpha } + \dfrac{{\Gamma (\alpha + 1)}}{{\Gamma (\alpha - 1/2)}}{t^{\alpha - 3/2}}} \\ {u(0) = 0,\quad {u^\prime }(0) = 0} \end{array}} \right. $

(23)

问题(23)的真解是 $u(t) = {t^\alpha }$ 。

使用系统(15)求解问题(23), 取 $\alpha = \dfrac{{111}}{{17}}$ 。

在图1（a）和1（b）的第1个图中，分别列出了问题（23）在相同的 $M$ （多项式阶数）以及变化的 $h$ （网格长度）下的数值误差（ ${\rm{h}}$ 型） $ {E_{{L^2}（0，T）}} $ 和 $ {E_{{L^\infty }（0，T）}} $ ，其中，选取 $M$ = $1，2，4，6$ 。由此可以观察到h型的快速收敛速度。

在图1(a)和1(b)的第2个图中，分别列出了问题(23)在相同的 $h$ （网格长度）以及变化的 $M$ （多项式阶数）下的数值误差( ${\text{p}}$ 型) $ {E_{{L^2}(0，T)}} $ 和 $ {E_{{L^\infty }(0，T)}} $ ，其中，选取 $h$ = $1，\dfrac{1}{2}，\dfrac{1}{4}，\dfrac{1}{6}$ $\Big(\Big.$ 或 $1，\dfrac{1}{2}，\dfrac{1}{4}，\dfrac{1}{8}$ $\Big.\Big)$ ，由此可以观察到数值误差的收敛速度随着M的减小按代数衰减。