研究了求解非线性分数阶微分方程的hp 型Legendre 谱配置法。首先提出将多分数阶微分方程转化成等价的Volterra 积分方程, 其次构造了近似求解原方程的数值方法, 最后通过数值实验说明了该算法的理论正确性以及所构造数值方法的有效性。 关键词 : 非线性分数阶微分方程 Legendre谱配置法 hp型误差界 Abstract : The hp-version Legendre spectral collocation method for solving nonlinear fractional differential equations was studied. At first, the multifractional differential equation was transformed into an equivalent Volterra integral equation. Then a numerical method to approximate the original equation was constructed. Finally, the correctness of the algorithm theory and the effectiveness of the proposed numerical method were demonstrated by numerical experiments. Key words : nonlinear fractional differential equations Legendre spectral collocation method hp-version error bounds 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. $ ${{\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)$ ${\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)$ 是标准的阶数为 的Legendre多项式。Legendre 多项式集合构成一个完备的 ${L^2}( - 1,1)$ 正交系统，即

$ {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 $ $ \int_{{I_n}} {{L_{n,k}}} (t){L_{n,j}}(t){\rm{d}}t = \frac{{{h_n}}}{{2k + 1}}{\delta _{k,j}} $ $ \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) $ $\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}$ ${\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} $

为便于使用, 引入另一个插值算子 $\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} $ $ \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} $

对于任意的 $\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} $ $ \hat \omega _{n,j}^\nu (t) = \int_{{t_{n - 1}}}^t {{{(t - \xi )}^\nu }} l_{{{n}},j}^t(\xi ){\rm{d}}\xi $ $ u(t) = \frac{1}{{\Gamma (1 + \alpha )}}\int_0^t {{{(t - s)}^\alpha }} f(s,u(s)){\rm{d}}s $ $\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} $ $\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} $ $\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} $ $ \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} $ $ \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} $ $ \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}$ $ \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} $ $ u_p^n = \frac{1}{{\Gamma (\alpha + 1)}}\left( {\sum\limits_{k = 1}^{n - 1} {a_{k,p}^n} + b_p^n} \right) $ $ \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} $ $ \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. $ CARRERAS B A, LYNCH V E, ZASLAVSKY G M. Anomalous diffusion and exit time distribution of particle tracers in plasma turbulence model[J]. Physics of Plasmas, 2001, 8(12): 5096-5103. DOI:10.1063/1.1416180 TARASOV V E. Fractional dynamics: applications of fractional calculus to dynamics of particles, fields and media[M]. New York: Springer, 2011. MAINARDI F. Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models[M]. London: Imperial College Press, 2010. SUGIMOTO N. Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[J]. Journal of Fluid Mechanics, 1991, 225: 631-653. DOI:10.1017/S0022112091002203 SRIVASTAVA V, RAI K N. A multi-term fractional diffusion equation for oxygen delivery through a capillary to tissues[J]. Mathematical and Computer Modelling, 2010, 51(5/6): 616-624. GUO Y L, WANG Z Q. An hp-version Chebyshev collocation method for nonlinear fractional differential equations[J]. Applied Numerical Mathematics, 2020, 158: 194-211. DOI:10.1016/j.apnum.2020.08.003 BRUNNER H. Collocation methods for Volterra integral and related functional differential equations[M]. New York: Cambridge University Press, 2004. WANG C L, WANG Z Q, JIA H L. An hp-version spectral collocation method for nonlinear Volterra integro-differential equation with weakly singular kernels[J]. Journal of Scientific Computing, 2017, 72(2): 647-678. DOI:10.1007/s10915-017-0373-3 ZAKY M A, AMEEN I G. On the rate of convergence of spectral collocation methods for nonlinear multi-order fractional initial value problems[J]. Computational and Applied Mathematics, 2019, 38(3): 144. DOI:10.1007/s40314-019-0922-5