带有Caputo-Fabrizio导数的分数阶微分方程的快速高阶算法的研究

带有Caputo-Fabrizio导数的分数阶微分方程的快速高阶算法的研究

傅博¹, 王世宇¹, 高婷婷¹, 吕学琴^,¹^,²^,^*

¹哈尔滨师范大学数学科学学院哈尔滨 150025

²天津中德应用技术大学天津 300350

A Fast High Order Method for Fractional Differential Equations with the Caputo-Fabrizio Derivative

Fu Bo¹, Wang Shiyu¹, Gao Tingting¹, Lv Xueqin^,¹^,²^,^*

¹School of Mathematical Sciences, Harbin Normal University, Harbin 150025

²Tianjin Sino-German University of Applied Sciences, Tianjin 300350

通讯作者: ^*吕学琴, E-mail: xueqinhsd@163.com

收稿日期: 2021-11-9 修回日期: 2022-10-8

基金资助:

天津市技术计划项目(21YDTPJC00120)
哈尔滨师范大学高等教育教学改革研究重点项目(XJGZ2022010)

Received: 2021-11-9 Revised: 2022-10-8

Fund supported:

Tianjin Technical Expert Project(21YDTPJC00120)
Key Research Project of Higher Education Teaching Reform in Harbin Normal University(XJGZ2022010)

摘要

对于数值求解常微分方程的传统方法而言, 其计算量和存储量一般都是比较大的, 为了解决这个问题, 该文基于L2方案和一个简单的递归关系提出了一个快速高阶数值方法求解带有Caputo-Fabrizio 导的分数阶微分方程, 该算法在很大程度上减少了存储量和计算量,进一步的, 该文分析了快速算法的可行性、误差估计和稳定性分析, 并通过数值算例论证的所提方法的正确性.

关键词： Caputo-Fabrizio导数; 快速高阶数值算法; L2方案

Abstract

The computational work and storage of numerically solving the ODEs are generally huge for the traditional direct methods. To overcome this difficulty, we presents a fast high order numerical method for fractional ordinary differential equations with the Caputo-Fabrizio derivative based on a L2 scheme and a simply recurrence relation. The algorithm greatly reduces the storage capacity and the total calculation cost. Furthermore, we also analyzes the feasibility of algorithm, error estimation and stability analysis of the fast scheme. Illustrative examples are included to demonstrate the performance of our technique.

Keywords： Caputo-Fabrizio derivative; Fast high order numerical method; L2 scheme

PDF (629KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

傅博, 王世宇, 高婷婷, 吕学琴. 带有Caputo-Fabrizio导数的分数阶微分方程的快速高阶算法的研究[J]. 数学物理学报, 2023, 43(3): 896-912

Fu Bo, Wang Shiyu, Gao Tingting, Lv Xueqin. A Fast High Order Method for Fractional Differential Equations with the Caputo-Fabrizio Derivative[J]. Acta Mathematica Scientia, 2023, 43(3): 896-912

1 引言

该文研究了一类带有非局部边界条件的微分方程数值求解,这类问题常应用于化学工程、热弹性学.在过去的几十年里, 分数微积分受到了物理学家和数学家的广泛关注, 因为它能准确地反映许多应用科学中的物理过程. 分数阶微分方程的数值分析已成为一个深入的研究课题. 分数阶微积分在近年来得到了广泛的研究^[1-2]. 它的研究结果可以说明若干物理问题^[3]. Caputo和Fabrizio^[4]提出了一种新的分数阶导数. 另外, 在文献[5]中给出了具有非局部非奇异核的分数阶导数. 对于一般的右端函数 $f$ , 通常很难找到它的解析解, 因此有必要研究新的数值方法来求解该方程.

Cao和Xu^[6-7]考虑了分数阶常微分方程, 给出了一种改进的逐块逼近分数阶时间导数的方法. Yang^[8]等致力于分数阶多步方法在分数阶扩散波方程中的应用. Sun 和Wu^[9]以及Lin和Xu^[10]分别提出了时间分数阶扩散方程时间离散的 $L1$ 格式, 并证明了该格式在时间上的收敛阶是 $2-\alpha$ . Jiang^[11]等给出了近似Caputo分数阶导数的快速逼近方法, 它是基于指数渐进法. Huang^[12]等证明了该方法的收敛性至少是3阶的. 在文献{[13 ⇓⇓⇓⇓⇓-19]} 中学者们研究了时间分数阶扩散方程的快速求解, 在文献[20-21] 中研究了解的起始时间奇点性的特别处理.

本文的主要目标是应用L2格式和简单递归关系来求解分数阶微分方程, 包括由Caputo和Fabrizio提出的光滑核分数阶导数. Caputo-Fabrizio分数阶导数有助于理解一些现象, 如二级流体的热分析、电化学现象、异常扩散等.

在该文中, 研究分数阶线性常微分方程

$\begin{equation}\label{1.1} {\cal D}_{t}^{\alpha}u(t)=f(t), \alpha\in (0,1), 0<t<T, \end{equation}$

(1.1)

初始条件 $u(0)=u_{0}$ 和右端项 $f(t)=\lambda [u(t)-u_{0}\exp(-\frac{\alpha}{1-\alpha}t)],$ 算子 ${\cal D}_{t}^{\alpha}$ 是 Caputo-Fabrizio 导数, 定义如下

$\begin{equation}\label{1.2} {\cal D}_{t}^{\alpha}u(t)=\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t}u'(\tau)\exp(-\alpha\frac{t-\tau}{1-\alpha}){\rm d}\tau, \end{equation}$

(1.2)

其中 $\alpha\in(0,1)$ , $M(\alpha)$ 满足 $M(0)=M(1)=1.$

在文献[7]中, 该方程(1.1)的唯一解如下

$u(t)=u_{0}\exp(\frac{\lambda \alpha t}{M(\alpha)-\lambda(1-\alpha)}),$

结果表明, 当 $u_{0}\lambda < 0$ 时, 解是递减的, 这将被用来分析算法的稳定性.

该文的主要内容和结构如下: 在第2节中, 介绍了一些基本的定义; 在第3节中, 详细描述了Caputo-Fabrizio导数的算法; 第4节给出了解的存在性、误差估计和稳定性分析. 验证本文的理论结果的数值算例在第5节给出.

2 主要定义和概念

${\bf 定义 2.1}$ $\alpha$ 阶Caputo分数阶时间导数(UFDt), 给出如下定义

$\begin{equation}\label{2.1} D_{t}^{\alpha}u(t)=\frac{ 1}{\Gamma(1-\alpha)}\int_{a}^{t}\frac{u'(\tau)}{(t-\tau)^{\alpha}}{\rm d}\tau. \end{equation}$

(2.1)

${\bf 定义 2.2}$ 设 $u\in H^{1}(a, b), b > a, \alpha\in[0,1]$ , Caputo-Fabrizio分数阶导数(NFDt)的定义如下

$\begin{equation}\label{2.2} {\cal D}_{t}^{\alpha}u(t)=\frac{ M(\alpha)}{1-\alpha}\int_{a}^{t}u'(\tau)\exp(-\alpha\frac{t-\tau}{1-\alpha}){\rm d}\tau, \end{equation}$

(2.2)

其中 $\alpha \in(0,1)$ , $M(\alpha)$ 满足 $M(0)=M(1)=1$ .

3 快速有限差分算子的构造

在该节中, 考虑方程(1.2)的高阶快速算法, 在区间 $[T]$ 上取一组网格点 $\Omega_{t}:=\{t_{n},$ $n=0,1,\cdots, N\}$ , 其中 $t_{0}=0, t_{N}=T,$ 并且 $h=\frac{T}{N}, t_{k-\frac{1}{2}}=t_{k-1}+\frac{1}{2}h, I_{k}=[t_{k-1}, t_{k}],$ $k=1,2,\cdots,N.$

L2 格式如下

$\begin{eqnarray*} \prod_{2,j+1}u(\tau)&=&u(t_{j-1})\frac{(\tau-t_{j})(\tau-t_{j+1})}{2h^2}+u(t_{j})\frac{(\tau-t_{j-1})(\tau-t_{j+1})}{-h^2} \\ && +u(t_{j+1})\frac{(\tau-t_{j-1})(\tau-t_{j})}{2h^2}. \end{eqnarray*}$

它是一个二阶多项式, 在点 $\{t_{j-1},t_{j},t_{j+1}\}$ 内插值 $u(t)$ . 首先将(1.2)式中的卷积积分分为 $D_{h}^{I}(t_{k})$ 和 $D_{h}^{II}(t_{k})$ , 即

$\begin{eqnarray*} {\cal D}_{t}^{\alpha}u(t)|_{t=t_{k}} &=&\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{k}}u'(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau\\ &=&\frac{ M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}u'(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau+\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{k-1}}u'(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau\\ &=&D_{h}^{I}(t_{k})+D_{h}^{II}(t_{k}). \end{eqnarray*}$

令

$\begin{eqnarray*} &&D_{h}^{I}(t_{k})=\frac{ M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}u'(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau, \\ && D_{h}^{II}(t_{k})=\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{k-1}}u'(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau, \end{eqnarray*}$

$D_{h}^{I}(t_{k})$ 被近似如下

$\begin{eqnarray*}\label{3.1} D_{h}^{I}(t_{k}) &\approx&\frac{M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}[\prod_{2,k}u(\tau)]' \exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau\\ &=&a_{k,k-2}u(t_{k-2})+a_{k,k-1}u(t_{k-1})+a_{k,k}u(t_{k}), \end{eqnarray*}$

(3.1)

其中

$\begin{eqnarray*} &&a_{k,k-2}=\frac{M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}{\varphi'_{k,k-2}(\tau)} \exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau, \\ && a_{k,k-1}=\frac{M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}{\varphi'_{k,k-1}(\tau)} \exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau, \\ && a_{k,k}=\frac{M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}{\varphi'_{k,k}(\tau)} \exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau, \\ && \varphi_{k,k-2}(\tau)=\frac{(\tau-t_{k-1})(\tau-t_{k})}{2h^2},\\ &&\varphi_{k,k-1}(\tau)=\frac{(\tau-t_{k-2})(\tau-t_{k})}{-h^2},\\ &&\varphi_{k,k}(\tau)=\frac{(\tau-t_{k-2})(\tau-t_{k-1})}{2h^2}. \end{eqnarray*}$

对于 $D_{h}^{II}(t_{k})$ 的部分, 本文应用分部积分来消除 $u'(\tau)$ 并得到

$\begin{eqnarray*}\label{3.2} D_{h}^{II}(t_{k}) &=&\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{k-1}}u'(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau\\ &=&\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{k-1}}\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}u(\tau)\\ &=&\frac{ M(\alpha)}{1-\alpha}\bigg[u(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha})\bigg|_{0}^{t_{k-1}} -\int_{0}^{t_{k-1}}\frac{\alpha}{1-\alpha}u(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau\bigg]\\ &=&\frac{ M(\alpha)}{1-\alpha}\bigg[u(t_{k-1})\exp(\frac{-\alpha h}{1-\alpha})-u(0)\exp(\frac{-\alpha t_{k}}{1-\alpha})\\ &&-\frac{\alpha}{1-\alpha}\int_{0}^{t_{k-1}}u(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau\bigg]. \end{eqnarray*}$

(3.2)

令

$H(t_{k})=\int_{0}^{t_{k-1}}\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha})u(\tau){\rm d}\tau,$

对于 $k=1,2,\cdots,N$ , 计算 $H(t_{k})$ , 该文得到如下的递归关系

$\begin{equation}\label{3.3} H(t_{k})=\exp(\frac{-\alpha h}{1-\alpha})H(t_{k-1})+ \int_{t_{k-2}}^{t_{k-1}}\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha})u(\tau){\rm d}\tau. \end{equation}$

(3.3)

在每一步, 只需要计算 $H(t_{k})$ , 因为 $H(t_{k-1})$ 是已知的.该文可以通过L2格式对 $u(t)$ 进行插值来计算(3.3)式右边的积分, 然后对其进行近似计算, 有

$\begin{eqnarray*} \int_{t_{k-2}}^{t_{k-1}}\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha})u(\tau){\rm d}\tau &\approx&\int_{t_{k-2}}^{t_{k-1}}\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha})\prod_{2,k}u(\tau){\rm d}\tau\\ &=&\hat{a}_{k,k-2}u(t_{k-2})+\hat{a}_{k,k-1}u(t_{k-1})+\hat{a}_{k,k}u(t_{k}), \end{eqnarray*}$

其中

$\hat{a}_{k,i}=\int_{t_{k-2}}^{t_{k-1}}\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha})\varphi_{k,i}{\rm d}\tau, i=k-2,k-1,k,$

因此

$\begin{equation}\label{3.4} D_{h}^{II}(t_{k})\approx\frac{ M(\alpha)}{1-\alpha} \bigg[u(t_{k-1})\exp(\frac{-\alpha h}{1-\alpha})-u(0)\exp(\frac{-\alpha t_{k}}{1-\alpha})-\frac{\alpha}{1-\alpha}H(t_{k})\bigg], \end{equation}$

(3.4)

其中

$H(t_{k})=\exp(\frac{-\alpha h}{1-\alpha})H(t_{k-1})+\hat{a}_{k,k-2}u(t_{k-2})+ \hat{a}_{k,k-1}u(t_{k-1})+\hat{a}_{k,k}u(t_{k}), k\geq2.$

下面本文引入有限差分算子 $L_{h}^{\alpha}$ 来离散函数 $\{u^{k}\}_{k=0}^{N}$ .

当 $2\leq k\leq N$ , 定义

$\begin{eqnarray*}\label{3.5} L_{h}^{\alpha}u_{k}&=&(a_{k,k-2}u_{k-2}+a_{k,k-1}u_{k-1}+a_{k,k}u_{k})\\ &&+\frac{ M(\alpha)}{1-\alpha}\bigg[\exp(\frac{-\alpha h}{1-\alpha})u_{k-1}-\exp(\frac{-\alpha t_{k}}{1-\alpha})u_{0}-\frac{\alpha}{1-\alpha}H(t_{k})\bigg]. \end{eqnarray*}$

(3.5)

当 $k=1$ 时, 在区间 $I_{1}=[t_{0},t_{1}]$ 上用二次插值逼近 $u(t)$ , 即

$\begin{equation}\label{3.6} u(\tau)\approx\varphi_{1,0}(\tau)u_{0}+\varphi_{1,1}(\tau)u_{\frac{1}{2}}+\varphi_{1,2}(\tau)u_{1}, \end{equation}$

(3.6)

其中

$\begin{eqnarray*} \varphi_{1,0}(\tau)=2\frac{(\tau-t_{1})(\tau-t_{\frac{1}{2}})}{h^2}, \varphi_{1,1}(\tau)=-4\frac{(\tau-t_{0})(\tau-t_{1})}{h^2}, \varphi_{1,2}(\tau)=2\frac{(\tau-t_{0})(\tau-t_{\frac{1}{2}})}{h^2}, \end{eqnarray*}$

根据泰勒公式, 可以得到如下插值关系

$\begin{equation} \label{3.7} u_{\frac{1}{2}}\approx\frac{3}{8}u_{0}+\frac{3}{4}u_{1}-\frac{1}{8}u_{2}, \end{equation}$

(3.7)

把上面的公式代入(3.6)式, 得到一个近似值

$\begin{eqnarray*}\label{3.8} u(\tau)&\approx&\bigg[\varphi_{1,0}(\tau)+\frac{3}{8}\varphi_{1,1}(\tau)\bigg]u_{0} +\bigg[\varphi_{1,2}(\tau)+\frac{3}{4}\varphi_{1,1}(\tau)\bigg]u_{1}-\frac{1}{8}\varphi_{1,1}(\tau)u_{2}\\ &\doteq&\omega_{0}(\tau)u_{0}+\omega_{1}(\tau)u_{1}+\omega_{2}(\tau)u_{2}, \end{eqnarray*}$

(3.8)

将(3.8)式带入(1.2)式得

$\begin{eqnarray*}\label{3.9} {\cal D}_{t}^{\alpha}u(t)|_{t=t_{1}} &=&\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}u'(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\\ &\approx&\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}[\omega_{0}'(\tau)u_{0}+\omega_{1}'(\tau)u_{1}+\omega_{2}'(\tau)u_{2}]\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\\ &=&a_{1,0}u_{0}+a_{1,1}u_{1}+a_{1,2}u_{2}\\ &\doteq& L_{h}^{\alpha}u_{1}, k=1, \end{eqnarray*}$

(3.9)

其中

$a_{1,i}=\frac{M(\alpha)}{1-\alpha}\int_{t_{0}}^{t_{1}}{\omega_{i}'(\tau)} \exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau, i=0,1,2.$

(3.10)

由(3.5)和(3.9)式, 本文得到如下形式

$\begin{equation}\label{3.10} L_{h}^{\alpha}u(t_{k})=f(t_{k}), ~k=1,2,\cdots,N. \end{equation}$

4 理论分析

当 $k=1$ , 由(3.9)和(3.10)式可得

$a_{1,0}u_{0}+a_{1,1}u_{1}+a_{1,2}u_{2}=f(t_{1})=\lambda (u_{1}-u_{0}\exp(\frac{ -\alpha t_{1}}{1-\alpha})),$

不难得到

$(a_{1,1}-\lambda)u_{1}+a_{1,2}u_{2}=-(a_{1,0}+\lambda\exp(\frac{ -\alpha t_{1}}{1-\alpha}))u_{0}.$

为了便于书写, 引入一些符号表示, 令

$\sigma=\exp(\frac{-\alpha}{1-\alpha}),b_{10}=a_{1,0}+\lambda\sigma^{t_{1}},b_{11}=a_{1,1}-\lambda,b_{12}=a_{1,2},$

于是得到

$\begin{equation}\label{4.1} b_{11}u_{1}+b_{12}u_{2}=-b_{10}u_{0}. \end{equation}$

(4.1)

当 $k=2$ , 由(3.5)和(3.10)式可得

$\begin{eqnarray*} &&a_{2,0}u_{0}+a_{2,1}u_{1}+a_{2,2}u_{2}+\frac{M(\alpha)}{1-\alpha} \bigg[\exp(\frac{ -\alpha h}{1-\alpha})u_{1}-\exp(\frac{ -\alpha t_{2}}{1-\alpha})u_{0}-\frac{\alpha}{1-\alpha}H(t_{2})\bigg]\\ &=&f(t_{2}) =\lambda(u_{2}-\exp(\frac{ -\alpha t_{2}}{1-\alpha})u_{0}), \end{eqnarray*}$

其中 $H(t_{2})=\hat{a}_{2,0}u_{0}+\hat{a}_{2,1}u_{1}+\hat{a}_{2,2}u_{2},$ 因此, 通过合并同类项得

$\begin{eqnarray*} &&(a_{2,1}+\frac{M(\alpha)}{1-\alpha}\exp(\frac{ -\alpha h}{1-\alpha})- \frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{a}_{2,1})u_{1}+(a_{2,2}-\frac{\alpha M(\alpha)} {(1-\alpha)^{2}}\hat{a}_{2,2}-\lambda)u_{2}\\ &=&-(a_{2,0}-\frac{M(\alpha)}{1-\alpha}\exp(\frac{ -\alpha t_{2}}{1-\alpha})-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{a}_{2,0}+\lambda\exp(\frac{ -\alpha t_{2}}{1-\alpha}))u_{0}, \end{eqnarray*}$

令

$b_{20}=a_{2,0}-\frac{ M(\alpha)}{(1-\alpha)}(-\sigma^{ t_2})-\frac{\alpha M(\alpha)}{(1-\alpha)^2}\hat{a}_{2,0} +\lambda\sigma^{t_{2}}, b_{21}=a_{2,1}-\frac{M(\alpha)}{1-\alpha}\sigma^{h}-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{a}_{2,1},$

$\begin{equation} b_{22}=a_{2,2}-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{a}_{2,2}-\lambda, \label{4.2} b_{21}u_{1}+b_{22}u_{2}=-b_{20}u_{0}. \end{equation}$

(4.2)

当 $k=3$ , 结果如下

$\begin{eqnarray*} &&(a_{3,1}u_{1}+a_{3,2}u_{2}+a_{3,3}u_{3})+\frac{ M(\alpha)}{1-\alpha} \bigg[\exp(\frac{-\alpha h}{1-\alpha})u_{2}-\exp(\frac{-\alpha t_{3}}{1-\alpha})u_{0}-\frac{\alpha}{1-\alpha}H(t_{3})\bigg]\\ &=&f(t_{3})=\lambda(u_{3}-u_{0}\exp(\frac{-\alpha t_{3}}{1-\alpha})), \end{eqnarray*}$

其中

$\begin{eqnarray*} H(t_{3})=\exp(\frac{-\alpha h}{1-\alpha})H(t_{2})+\hat{a}_{3,1}u_{1}+ \hat{a}_{3,2}u_{2}+\hat{a}_{3,3}u_{3}, \end{eqnarray*}$

因此

$\begin{eqnarray*} &&\bigg[a_{3,1}-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}(\hat{a}_{2,1}\exp(\frac{ -\alpha h}{1-\alpha}) +\hat{a}_{3,1})\bigg]u_{1}+\bigg[a_{3,2}+\frac{M(\alpha)}{1-\alpha}\exp(\frac{ -\alpha h}{1-\alpha})\\ &&-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}(\hat{a}_{2,2}\exp(\frac{ -\alpha h}{1-\alpha})+{a}_{3,2})\bigg]u_{2} +(a_{3,3}-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{a}_{3,3}-\lambda)u_{3}\\ &=&-\bigg[\frac{ M(\alpha)}{1-\alpha}(-\exp(\frac{-\alpha t_{3}}{1-\alpha}))-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{a}_{2,0}\exp(\frac{ -\alpha h}{1-\alpha})+\exp(\frac{-\alpha t_{3}}{1-\alpha})\bigg]u_{0}, \end{eqnarray*}$

令

$\begin{eqnarray*} &&b_{30}=\frac{ M(\alpha)}{1-\alpha}(-\sigma^{t_{3}})-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{a}_{2,0}\sigma^{h}+\lambda\sigma^{t_{3}},b_{31}=a_{3,1}-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}(\hat{a}_{2,1}\sigma^{h}+\hat{a}_{3,1}), \\ &&b_{32}=a_{3,2}+\frac{M(\alpha)}{1-\alpha}\sigma^{h}-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}(\hat{a}_{2,2}\sigma^{h}+\hat{a}_{3,2}),b_{33}=a_{3,3}-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{a}_{3,3}-\lambda, \end{eqnarray*}$

于是得到

$\begin{equation}\label{4.3} b_{31}u_{1}+b_{32}u_{2}+b_{33}u_{3}=-b_{30}u_{0}. \end{equation}$

(4.3)

同理, 可以得到 $4\leq k\leq N$ 时的公式, 表示如下

$\begin{eqnarray*} && (a_{N,N-2}u_{N-2}+a_{N,N-1}u_{N-1}+a_{N,N}u_{N})\\ &&+\frac{ M(\alpha)}{1-\alpha}\bigg[\exp(\frac{ -\alpha h}{1-\alpha}) u_{N-1}-\exp(\frac{ -\alpha t_{N}}{1-\alpha})u_{0} -\frac{(\alpha)}{1-\alpha}H(t_{N})\bigg]\\ &=&f(t_{N})=\lambda(u_{N})-u_{0}\exp(\frac{ -\alpha t_{N}}{1-\alpha}), \end{eqnarray*}$

其中

$\begin{eqnarray*} H(t_{N})&=&\exp(\frac{ -\alpha h}{1-\alpha})H(t_{N-1})+\hat{a}_{N,N-2}u_{N-2}+\hat{a}_{N,N-1}u_{N-1}+\hat{a}_{N,N}u_{N}\\ &=&\exp(\frac{ -(N-2)\alpha h}{1-\alpha})(\hat{a}_{2,0}u_{0}+\hat{a}_{2,1}u_{1}+\hat{a}_{2,2}u_{2})\\ &&+\exp(\frac{-(N-3)\alpha h}{1-\alpha})(\hat{a}_{3,1}u_{1}+\hat{a}_{3,2}u_{2}+\hat{a}_{3,3}u_{3})\\ &&+\exp(\frac{ -(N-4)\alpha h}{1-\alpha})(\hat{a}_{4,2}u_{2}+\hat{a}_{4,3}u_{3}+\hat{a}_{4,4}u_{4}) +\cdots\\ &&+\exp(\frac{-2\alpha h}{1-\alpha})(\hat{a}_{N-2,N-4}u_{N-4}+\hat{a}_{N-2,N-3}u_{N-3}+\hat{a}_{N-2,N-2}u_{N-2})\\ &&+\exp(\frac{\alpha h}{1-\alpha})(\hat{a}_{N-1,N-3}u_{N-3}+\hat{a}_{N-1,N-2}u_{N-2}+\hat{a}_{N-1,N-1}u_{N-1})\\ &&+(\hat{a}_{N,N-2}u_{N-2}+\hat{a}_{N,N-1}u_{N-1}+\hat{a}_{N,N}u_{N}), \end{eqnarray*}$

因此, 不难得到下列公式

$\begin{eqnarray*} &&\bigg[-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\exp(\frac{ -(N-2)\alpha h}{1-\alpha})\hat{a}_{2,1}+\exp(\frac{ -(N-3)\alpha h}{1-\alpha})\hat{a}_{3,1})\bigg]u_{1}\\ &&+ \bigg[-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\exp(\frac{ -(N-2)\alpha h}{1-\alpha})\hat{a}_{2,2}+\exp(\frac{ -(N-3)\alpha h}{1-\alpha})\hat{a}_{3,2}+\exp(\frac{ -(N-4)\alpha h}{1-\alpha})\hat{a}_{4,2})\bigg]u_{2}\\ &&+ \bigg[-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\exp(\frac{ -(N-3)\alpha h}{1-\alpha})\hat{a}_{3,3}+\exp(\frac{ -(N-4)\alpha h}{1-\alpha})\hat{a}_{4,3})\bigg]u_{3}+\cdots \\ &&+ \bigg[a_{N,N-2}-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\exp(\frac{ -2\alpha h}{1-\alpha})\hat{a}_{N-2,N-2}+\exp(\frac{ -\alpha h}{1-\alpha})\hat{a}_{N-1,N-2}+\hat{a}_{N,N-2})\bigg]u_{N-2}\\ &&+ \bigg[a_{N,N-1}+\frac{M(\alpha)}{1-\alpha}\exp(\frac{ -\alpha h}{1-\alpha}-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\exp(\frac{ -\alpha h}{1-\alpha}\hat{a}_{N-1,N-1}+\hat{a}_{N,N-1})\bigg]u_{N-1}\\ &&+ \bigg[a_{N,N}-\frac{\alpha M(\alpha)}{(1-\alpha)^2}\hat{a}_{N,N}-\lambda\bigg]u_{N}\\ &=&-\bigg[\frac{ M(\alpha)}{(1-\alpha)}(-\exp(\frac{ -\alpha t_{N}}{1-\alpha}))-\frac{\alpha M(\alpha)}{(1-\alpha)^2}\exp (\frac{ -(N-2)\alpha h}{1-\alpha})\hat{a}_{2,0}+\exp(\frac{ -\alpha t_{N}}{1-\alpha})\bigg]u_{0},\label{4.4} \end{eqnarray*}$

(4.4)

令

$\begin{eqnarray*} &&b_{i0}=\frac{ M(\alpha)}{(1-\alpha)}(-\sigma^{ t_i})-\frac{\alpha M(\alpha)}{(1-\alpha)^2}\hat{a}_{2,0}\sigma^{h i}+\lambda\sigma^{t_{i}}, i=4,\cdots, N,\\ &&b_{i1}=-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\hat{a}_{2,1}\sigma^{(2-i)h}+\hat{a}_{3,1}\sigma^{(3-i)h}), i=4,\cdots,N,\\ &&b_{i2}=-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\hat{a}_{2,2}\sigma^{(2-i)h}+\hat{a}_{3,2}\sigma^{(3-i)h}+\hat{a}_{4,2}\sigma^{(4-i)h}), i=5,\cdots,N,\\ &&b_{ii-2}=a_{i,i-2}-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\sigma^{-2h}\hat{a}_{i-2,i-2}+\sigma^{-h}\hat{a}_{i-1,i-2}+\hat{a}_{i,i-2}), i=4,\cdots, N, \\ &&b_{ii-1}=a_{i,i-1}+\frac{M(\alpha)\sigma^{h}}{1-\alpha}-\frac{\alpha M(\alpha)}{(1-\alpha)^2}(\sigma^{-h}\hat{a}_{i-1,i-1}+\hat{a}_{i,i-1}), i=4,\cdots, N, \\ &&b_{ii}=a_{i,i}-\frac{\alpha M(\alpha)}{(1-\alpha)^2}\hat{a}_{i,i}-\lambda, i=4,\cdots, N, \\ &&\tilde{f}_{i}=-b_{i0}u_{0},i=1,\cdots, N. \end{eqnarray*}$

那么方程组(3.10)可以改写为以下形式

$\begin{equation}\label{4.5} \begin{array}{ll} b_{11}u_{1}+b_{12}u_{2}=\tilde{f}_{1}, b_{21}u_{1}+b_{22}u_{2}=\tilde{f}_{2},\\ b_{31}u_{1}+b_{32}u_{2}+b_{33}u_{3}=\tilde{f}_{3},\\ b_{41}u_{1}+b_{42}u_{2}+b_{43}u_{3}+b_{44}u_{4}=\tilde{f}_{4},\\ \cdots~~~~~~~\cdots~~~~~~~\cdots~~~~~~~\cdots~~~~~~~\cdots\\ b_{N-1 1}u_{1}+\cdots+b_{N-1 N-1}u_{N-1}=\tilde{f}_{N-1},\\ b_{N 1}u_{1}+b_{N 2}u_{2}+\cdots+b_{N N-1}u_{N-1}+b_{N N}u_{N}=\tilde{f}_{N}. \end{array} \end{equation}$

(4.5)

4.1 方程组(4.5)的可解性分析

接下来本文给出一些引理, 这些引理将用于可解性分析和稳定性分析.

${\bf 引理 4.1}$ 所有符号如上所述, 令 $h<\frac{2(1-\alpha)}{\alpha}$ , 由此可见以下式子成立

$\begin{equation}\label{4.6} b_{11}b_{22}-b_{21}b_{12}> 0,~~\mbox{和} ~~b_{ii}>0, i=1,2,\cdots,N. \end{equation}$

(4.6)

${\bf证}$ 首先由上面公式, $b_{11}$ 很容易获得, 类似的, 计算 $b_{kk}, 2\leq k\leq N,$ 有

$\begin{eqnarray*}\label{4.7} b_{kk} &=&a_{k,k}-\frac{\alpha M(\alpha)}{(1-\alpha)^2}\hat{a}_{k,k}-\lambda\\ &=&\frac{ M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\omega_{1}(\tau)-\lambda\\ &=&\frac{ M(\alpha)}{1-\alpha} \bigg[\omega_{1}(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha})\bigg|_{t_{k-1}}^{t_{k}} -\int_{t_{k-2}}^{t_{k-1}}\frac{\alpha}{1-\alpha}\varphi_{k,k}(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\bigg]-\lambda\\ &\approx&\frac{ M(\alpha)}{1-\alpha}\bigg[1-\frac{\alpha h}{2(1-\alpha)}\bigg]-\lambda, \end{eqnarray*}$

(4.7)

由 $\lambda<0$ , $h<\frac{2(1-\alpha)}{\alpha}$ , 可得 $b_{ii}>0, i=1,2,\cdots,N.$

根据 $b_{12}$ 和 $b_{21}$ 的定义, 经过计算可以得到

$\begin{eqnarray*}\label{4.8} b_{12}&=&\frac{ M(\alpha)}{1-\alpha}\int_{t_0}^{t_{1}}\omega_{2}'(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\\ &=&\frac{ M(\alpha)}{1-\alpha}\int_{t_0}^{t_{1}}\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\omega_{2}(\tau)\\ &=&\frac{ M(\alpha)}{1-\alpha}\bigg[\omega_{2}(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha})\bigg|_{t_0}^{t_{1}} -\int_{t_0}^{t_{1}}\frac{\alpha}{1-\alpha}\omega_{2}(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\bigg]\\ &=&0,\\ \end{eqnarray*}$

(4.8)

$\begin{eqnarray*} \label{4.9} b_{21}&=&\frac{ M(\alpha)}{1-\alpha}\int_{t_1}^{t_{2}}\varphi_{2,1}'(\tau)\exp(-\alpha\frac{t_{2}-\tau}{1-\alpha}){\rm d}\tau+\frac{ch M(\alpha)}{1-\alpha}\\ &&-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}\int_{t_0}^{t_{1}}\frac{\alpha}{1-\alpha}\varphi_{2,1}(\tau)\exp(-\alpha\frac{t_{2}-\tau}{1-\alpha}){\rm d}\tau\\ &=&\frac{ \alpha hM(\alpha)}{(1-\alpha)^{2}}\exp(-\frac{\alpha h}{1-\alpha}), \end{eqnarray*}$

(4.9)

结合(4.6)-(4.9)式, $0<\alpha<1$ , $h<\frac{2(1-\alpha)}{\alpha}$ , 因此

$\begin{equation}\label{4.10} b_{11}b_{22}-b_{21}b_{12} =\bigg[\frac{ M(\alpha)}{1-\alpha}(1-\frac{\alpha h}{2(1-\alpha)})-\lambda\bigg]^2 > 0. \end{equation}$

(4.10)

所以引理4.1得证.

${\bf 定理 4.1}$ 假设 $h<\frac{2(1-\alpha)}{\alpha}$ , 对于任意给定的初始条件 $u_{0}$ 和 $f(t)$ , 线性方程组(3.10)有唯一解 $U=(u_{1},u_{2},\cdots,u_{N})^{T}$ .

${\bf证}$ 首先, 来解二阶线性方程组

$\begin{equation}\label{4.11} B_{2}U_{2}=F_{2}, \end{equation}$

(4.11)

其中 $B_{2}=\left( \begin{array}{cc} b_{11}~ & b_{12} \\ b_{21}~& b_{22} \\ \end{array} \right),$ $U_{2}=(u_{1},u_{2})^{\top},$ $F_{2}=(\tilde{f}_{1}, \tilde{f}_{2})^{\top}.$

根据引理4.1的(4.10)式, $| B_2| \neq 0$ , 因此方程组(4.11)有唯一解 $U_2$ . 进一步, 该文求解 $N-2$ 阶线性方程组

$\begin{equation}\label{4.12} \begin{array}{lll} \begin{array}{lll}B_{N-2}U_{N-2}=F_{N-2},\end{array} \end{array} \end{equation}$

(4.12)

其中 $U_{N-2}=(u_3,u_4,\cdots, u_N)^{\top}, F_{N-2}=(\tilde{f}_{3},\tilde{f}_{4},\cdots,\tilde{f}_{N})^{\top}, B_{N-2}=(b_{ij}), j\leq i, i,j=3,\cdots,N$ 是一个下三角矩阵, 由引理4.1可知, $B_{N-2}$ 的所有对角线上的项都是正的, 因此, 线性方程组(4.12)有唯一解 $U_{N-2}$ . 证毕.

4.2 误差估计

${\bf 定理 4.2}$ 假设 $u(t)\in C^{4}[T]$ 并且令

$R_{k}=D_{h}^{\alpha}u(t_{k})-L_{h}^{\alpha}u(t_{k}), k=1,2,\cdots,N,$

那么下面式子成立

$| R_{k}| =O(h^4), h=\frac{T}{N}.$

${\bf证}$ 首先考虑 $k=1$ 的情况, 应用二次插值误差估计, 该文有

$u(\tau)-I_{[t_{0},t_{1}]}u(\tau)=\frac{u'''(\xi)}{3!}(\tau-t_{0})(\tau-t_{\frac{1}{2}})(\tau-t_{1}), \xi\in (t_{0},t_{1}),$

由(3.5)-(3.7)式可得

$\begin{eqnarray*}\label{4.13} | R_{1}| &=&\bigg| \frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}u'(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau -(a_{1,0}u_{0}+a_{1,1}u_{1}+a_{1,2}u_{2})\bigg| \\ &=&\bigg| \frac{M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}u'(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau \\ &&-\frac{M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}[\omega_{0}'(\tau)u_{0} +\omega_{1}'(\tau)u_{1}+\omega_{2}'(\tau)u_{2}]\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\bigg| \\ &=&\bigg| \frac{M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}u'(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau-\frac{M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}[(\varphi_{1,0}(\tau)+\frac{3}{8}\varphi_{1,1}(\tau))'u_{0}\\ &&+(\varphi_{1,2}(\tau)+\frac{3}{4}\varphi_{1,1}(\tau))'u_{1}-\frac{1}{8}\varphi_{1,1}(\tau)'u_{2}]\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\bigg| \\ &=&\bigg| \frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}u'(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau-\frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}(I_{1}u(\tau))'\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\bigg| \\ &\leq&\frac{ M(\alpha)}{1-\alpha}\bigg| \int_{0}^{t_{1}}[u(\tau)-I_{1}u(\tau)]'\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\bigg| \\ &=&\bigg| \frac{ M(\alpha)}{1-\alpha}\int_{0}^{t_{1}}\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}[u(\tau)-I_{1}u(\tau)]\bigg| \\ &\leq&\frac{ M(\alpha)}{1-\alpha}\bigg| \frac{u'''(\xi)}{3!}(\tau-t_{0})(\tau-t_{\frac{1}{2}})(\tau-t_{1})\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha})\Big| _{t_{0}}^{{t_1}}\\ &&-\int_{0}^{t_{1}}\frac{u'''(\xi)}{3!}(\tau-t_{0})(\tau-t_{\frac{1}{2}})(\tau-t_{1}){\rm d} \exp(-\alpha\frac{t_{1}-\tau}{1-\alpha})\bigg| \\ &\leq&\frac{\alpha M(\alpha)}{(1-\alpha)^2}\bigg| \frac{u'''(\xi)}{3!}h^{3}\int_{0}^{t_{1}}{\rm d}\tau \bigg| =O(h^4). \end{eqnarray*}$

(4.13)

其次, 考虑误差 $R_{k}, 2\leq k\leq N$ , $R_{k}$ 来自局部项 $D_{h}^{I}(t_{k})$ 和局部项 $D_{h}^{II}(t_{k})$ 这两部分的误差, 分别记为 $R_{k}^{1}$ 和 $R_{k}^{2}$ . 应用 $u(\tau)$ 的二次插值 $\prod_{2,k}$ , 得到了以下误差估计公式

$\begin{eqnarray*} \label{4.14} &&u(\tau)-\prod_{2,k}u(\tau)=\frac{u'''(\xi)}{3!}(\tau-t_{k-1})(\tau-t_{k-2})(\tau-t_{k}), \xi\in (t_{k-2},t_{k}), \end{eqnarray*}$

(4.14)

$\begin{eqnarray*} \label{4.15} | R_{k}^{1}| &=&\bigg| \frac{ M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}u'(\tau)\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau -(a_{k,k-2}u_{k-2}+a_{k,k-1}u_{k-1}+a_{k,k}u_{k})\bigg| \\ &=&\bigg| \frac{M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}u'(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau -\frac{M(\alpha)}{1-\alpha}\int_{t_{k-1}}^{t_{k}}[\prod_{2,k}u(\tau)]'\exp(-\alpha\frac{t_{k}-\tau}{1-\alpha}){\rm d}\tau\bigg| \\ &\leq& \frac{M(\alpha)}{1-\alpha}\bigg| \int_{t_{k-1}}^{t_{k}}[u(\tau)-\prod_{2,k}u(\tau)]'\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau\bigg| \\ &=&\frac{M(\alpha)}{1-\alpha}\bigg| \int_{t_{k-1}}^{t_{k}}\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}[u(\tau)-\prod_{2,k}u(\tau)]\bigg| \\ &=&\frac{M(\alpha)}{1-\alpha}\bigg| \exp(-\alpha\frac{t_{1}-\tau}{1-\alpha})[u(\tau)-\prod_{2,k}u(\tau)]\Big| _{t_{k-1}}^{t_{k}}\\ &&-\frac{\alpha}{1-\alpha}\int_{t_{k-1}}^{t_{k}}\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha})[u(\tau)-\prod_{2,k}u(\tau)]{\rm d}\tau \bigg| \\ &=&\frac{\alpha M(\alpha)}{1-\alpha}\bigg| \int_{t_{k-1}}^{t_{k}}\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha})[u(\tau)-\prod_{2,k}u(\tau)]{\rm d}\tau\bigg | \\ &=&\frac{\alpha M(\alpha)}{1-\alpha}\bigg| \int_{t_{k-1}}^{t_{k}}\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}) \frac{u'''(\xi)}{3!}(\tau-t_{k-1})(\tau-t_{k-2})(\tau-t_{k}){\rm d}\tau\bigg| \\ &\leq& \max\limits_{\xi\in(t_{k-2},t_{k})}\frac{| u'''(\xi)| }{3!}h^4\\ &=&O(h^4). \end{eqnarray*}$

(4.15)

$R_{k}^{2}$ 的误差来自 $H(t_{k})$ 右侧的第二项, 与(4.15)式类似, 很容易计算出来 $| R_{k}^{2}| =O(h^4)$ . 综上所述, 此定理成立.

4.3 稳定性分析

${\bf引理 4.2}$ 所有符号如前面所示, 令 $h<\frac{2(1-\alpha)}{\alpha}$ 并且

$\begin{eqnarray*} &&\beta_{1}:=-b_{10}-\lambda \exp(\frac{-\alpha t_{1}}{1-\alpha}), \beta_{2}:=-b_{20}-\lambda \exp(\frac{-\alpha t_{2}}{1-\alpha}), \\ &&\beta_{i}:=-b_{i0}-\lambda \exp(\frac{-\alpha t_{i}}{1-\alpha})-\sum\limits_{j=1}^{i-1}b_{ij}, \tilde{\beta}_{i}=b_{ii}-\lambda, i=3,4,\cdots,N, \\ && D=(b_{11}-\lambda)(b_{22}-\lambda)-b_{12}b_{21}, D_{1}=\beta_{1}(b_{22}-\lambda)-\beta_{2}b_{12}, D_{2}=\beta_{2}(b_{11}-\lambda)-\beta_{1}b_{21}, \end{eqnarray*}$

因此下列结论成立

(1) $| D| \neq 0, D_{1}<D, D_{2}<D,$

(2) $b_{3i}<0, i=0,1,2,$

(3) $\beta_{i}>0, \tilde{\beta}_{i}>0, i=3,4,\cdots,N.$

${\bf证}$ (1)因为 $\lambda<0, h<\frac{2(1-\alpha)}{\alpha}$ , 由此可见

$\begin{eqnarray*} D&=&(b_{11}-\lambda)(b_{22}-\lambda)-b_{12}b_{21}\\ &=&\bigg(\frac{M(\alpha)}{1-\alpha}\int_{t_{0}}^{t_{1}}\omega_{1}'(\tau)\exp(-\alpha\frac{t_{1}-\tau}{1-\alpha}){\rm d}\tau-2\lambda \bigg) (a_{22}-\frac{\alpha M(\alpha)}{(1-\alpha)^{2}}\hat{\alpha}_{22}-2\lambda)-b_{12}b_{21}\\ &=&\bigg[\frac{M(\alpha)}{1-\alpha}(1-\frac{\alpha h}{2(1-\alpha)})-2\lambda\bigg]^2 \neq0, \end{eqnarray*}$

由

$\exp(\frac{-\alpha t_{1}}{1-\alpha})<1, \lambda<0, -(b_{10}+b_{22})=\frac{\alpha M(\alpha)}{(1-\alpha)}(\frac{\alpha h}{2(1-\alpha)}-1)(1+\sigma^{t_1})+\lambda-\sigma^{t_{1}}<0,$

因此

$D_{1}-D=(b_{22}-\lambda)\bigg[-(b_{10}+b_{22})+\lambda(1-\exp(\frac{-\alpha t_{1}}{1-\alpha}))\bigg]<0,$

同样, 由

$\exp(\frac{-\alpha t_{1}}{1-\alpha})<1, b_{10}<0, b_{21}>0, \lambda<0,$

$-(b_{20}+b_{11})=\frac{\alpha M(\alpha)}{(1-\alpha)}(\frac{\alpha h}{2(1-\alpha)}-1)(1+\sigma^{t_2})+\lambda-\sigma^{t_{2}}<0,$

因此

$D_{2}-D=(b_{11}-\lambda)\bigg[-(b_{20}+b_{11})+\lambda(1-\exp(\frac{-\alpha t_{2}}{1-\alpha}))\bigg]-b_{21} \bigg[-b_{10}-\lambda \exp(\frac{-\alpha t_{1}}{1-\alpha})\bigg]<0.$

(2)根据 $b_{3i}, i=0,1,2$ 的定义, 分别计算

$\begin{eqnarray*} b_{30} &=&-\frac{M(\alpha)}{1-\alpha)}\sigma^{t_3}-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}\hat{a}_{2,0}\sigma^{3h}+\sigma^{t_{3}}\\ &=&-\frac{M(\alpha)}{1-\alpha}\sigma^{t_3}-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}\int_{t_0}^{t_{1}}\exp(-\alpha\frac{t_{2}-\tau}{1-\alpha})\varphi_{2,0}(\tau){\rm d}\tau \sigma^{3h}+\sigma^{t_{3}}\\ &\approx&-\frac{M(\alpha)}{1-\alpha}\sigma^{t_3}-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}\frac{h}{2}\exp(\frac{-2\alpha h}{1-\alpha})\sigma^{3h}+\sigma^{t_{3}}\\ &=&-\bigg[\frac{M(\alpha)}{1-\alpha}\sigma^{t_3}+\frac{ \alpha h M(\alpha)}{2(1-\alpha)^2}\exp(\frac{-2\alpha h}{1-\alpha})\sigma^{3h}-\sigma^{t_{3}}\bigg] <0, \end{eqnarray*}$

同样, 省略计算过程, 该文可以得到

$\begin{eqnarray*} b_{31} &=&a_{3,1}-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}(\hat{a}_{2,1}\sigma^{-h}+\hat{a}_{3,1})\\ &=&-\frac{\alpha M(\alpha)}{(1-\alpha)^2}\bigg[\frac{h}{2}\exp(\frac{-\alpha h}{1-\alpha})\sigma^{-h}+\frac{h}{2}\exp(\frac{-2\alpha h}{1-\alpha})\bigg]<0, \\ b_{32} &=&a_{3,2}-\frac{M(\alpha)}{1-\alpha}\sigma^{h}-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}(\hat{a}_{2,2}\sigma^{-h}+\hat{a}_{3,2})\\ &=&-\frac{ M(\alpha)}{1-\alpha}\exp(\frac{-\alpha h}{1-\alpha})\bigg[2+\frac{\alpha h}{2(1-\alpha)}+\frac{h}{2}\exp(\frac{-\alpha h}{1-\alpha})\bigg]<0. \end{eqnarray*}$

(3) 本文把 $\beta_{i}, (i=3,4,\cdots,N)$ 分成两部分, 并写成

$\begin{eqnarray*} \beta_{i1}=-b_{i0}-\lambda \exp(\frac{-\alpha t_{i}}{1-\alpha}), \beta_{i2}=-\sum\limits_{j=1}^{i-1}b_{ij}, \end{eqnarray*}$

现在, 证明 $\beta_{i1}>0, \beta_{i2}>0.$ 因为 $\hat{a}_{2,0}>0, \lambda<0$ 和 $0<\alpha<1$ ,因此

$\beta_{i1} =-b_{i0}-\lambda \exp(\frac{-\alpha t_{i}}{1-\alpha})) =\frac{M(\alpha)}{1-\alpha}\sigma^{t_i}+\frac{ \alpha M(\alpha)}{(1-\alpha)^2}\hat{a}_{2,0}\sigma^{ih}-\lambda \exp(\frac{-\alpha t_{i}}{1-\alpha})>0.$

下一步, 对 $j$ 证明 $b_{ij}<0, j=1,2,\cdots,i-1.$

如果 $j=1$ , 由 $\hat{a}_{2,1}>0, \hat{a}_{3,1}>0$ , 得到

$b_{i1}=-\frac{ \alpha M(\alpha)}{1-\alpha}(\hat{a}_{2,1}\sigma^{(2-i)h}+\hat{a}_{3,1}\sigma^{(3-i)h})<0.$

如果 $j=i-2$ , 经过计算

$a_{i,i-2}=0, \hat{a}_{i-2,i-2}=0, \hat{a}_{i,i-2}=\hat{a}_{i-1,i-2}=\frac{h}{2}\exp(\frac{-2\alpha h}{1-\alpha})>0,$

因此

$b_{i i-2}=a_{i,i-2}-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}(\hat{a}_{i-2,i-2}\sigma^{2h}+\hat{a}_{i-1,i-2}\sigma^{h}+\hat{a}_{i,i-2})<0.$

如果 $j=i-1$ , 由

$a_{i,i-1}=-\frac{ M(\alpha)}{1-\alpha}\exp(\frac{-\alpha h}{1-\alpha})(1+\frac{\alpha h}{2(1-\alpha)}), \hat{a}_{i-1,i-1}>0, \hat{a}_{i,i-1}>0,$

因此

$\begin{eqnarray*} b_{i,i-1} &=&a_{i,i-1}+\frac{ M(\alpha)}{1-\alpha}\sigma^{h}-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}(\hat{a}_{i-1,i-1}\sigma^{h}+\hat{a}_{i,i-1})\\ &=&-\frac{ M(\alpha)}{1-\alpha}\exp(\frac{-\alpha h}{1-\alpha})\frac{\alpha h}{2(1-\alpha)}-\frac{ \alpha M(\alpha)}{(1-\alpha)^2}(\hat{a}_{i-1,i-1}\sigma^{h}+\hat{a}_{i,i-1})\\ &<&0. \end{eqnarray*}$

所以 $\beta_{i2}=-\sum\limits_{j=1}^{i-1}b_{ij}>0, \beta_{i}>0.$ 从而 $\tilde{\beta}_{i}=b_{ii}-\lambda=\frac{ M(\alpha)}{1-\alpha}(1-\frac{\alpha h}{2(1-\alpha)})-2\lambda>0.$ 证毕.

${\bf 定理 4.3}$ 假设 $u_{0}>0, \lambda<0$ , 方程组(4.5)在 $h<\frac{2(1-\alpha)}{\alpha}$ 的条件下是稳定的, 即数值解 $u_{i}$ 满足

$\begin{equation}\label{4.16} 0<u_{i}<u_{0}, i=1,2,\cdots,N. \end{equation}$

(4.16)

${\bf证}$ 首先, 研究以下方程

$\begin{eqnarray*} &&b_{10}u_{0}+b_{11}u_{1}+b_{12}u_{2}=\lambda\bigg[u_{1}-u_{0}\exp(-\frac{\alpha t_{1}}{1-\alpha})\bigg], \\ && b_{20}u_{0}+b_{21}u_{1}+b_{22}u_{2}=\lambda\bigg[u_{2}-u_{0}\exp(-\frac{\alpha t_{2}}{1-\alpha})\bigg], \end{eqnarray*}$

这个方程的解可表示为

$u_{1}=\frac{D_{1}}{D}u_{0},~~~u_{2}=\frac{D_{2}}{D}u_{0},$

其中

$\begin{eqnarray*} &&D=(b_{11}-\lambda)(b_{22}-\lambda)-b_{12}b_{21}, \\ &&D_{1}=(b_{22}-\lambda)(-b_{10}-\lambda \exp(\frac{-\alpha t_{1}}{1-\alpha}))-b_{12}(-b_{20}-\lambda \exp(\frac{-\alpha t_{2}}{1-\alpha})), \\ &&D_{2}=(b_{11}-\lambda)(-b_{20}-\lambda \exp(\frac{-\alpha t_{2}}{1-\alpha}))-b_{21}(-b_{10}-\lambda \exp(\frac{-\alpha t_{1}}{1-\alpha})). \end{eqnarray*}$

由引理4.2的(1)可知, $| D| \neq 0, D_{1}-D<0, D_{2}-D<0$ , 因此, 得到

$\begin{equation}\label{4.17} 0<u_{i}<u_{0}, i=1,2. \end{equation}$

(4.17)

对于 $3\leq i\leq N$ , 有以下公式

$\begin{equation}\label{4.18} u_{i}=\bigg[(-b_{i0}-\lambda \exp(\frac{-\alpha t_i}{1-\alpha}))u_{0}-\sum\limits_{j=1}^{i-1}b_{ij}u_{j}\bigg]/(b_{ii}-\lambda),\end{equation}$

(4.18)

由上式可知

$u_{3}=\bigg[(-b_{30}-\lambda \exp(\frac{-\alpha t_3}{1-\alpha}))u_{0}-b_{31}u_{1}-b_{32}u_{2}\bigg]/(b_{33}-\lambda).$

根据引理4.1和引理4.2的(2)可知, $b_{33}>0, b_{3i}<0, i=0,1,2$ , 并且由(4.17)式得出如下形式

$u_{3}<\bigg[(-b_{30}-\lambda \exp(\frac{-\alpha t_3}{1-\alpha}))-b_{31}-b_{32}\bigg]u_{0}/(b_{33}-\lambda),$

进一步, 可以直接验证 $0<-b_{30}-\lambda \exp(\frac{-\alpha t_3}{1-\alpha})-b_{31}-b_{32}<(b_{33}-\lambda),$ 所以

$\begin{equation}\label{4.19} 0<u_{3}<u_{0}. \end{equation}$

(4.19)

接下来, 本文用数学归纳法来证明剩下的结果. 假设(4.12)式成立 $i=k-1\leq N,$ 该文将证明 $i=k.$ 时的情况, 由(4.14)式得

$\begin{equation}\label{4.20} u_{k}=\bigg[(-b_{k0}-\lambda \exp(\frac{-\alpha t_{k}}{1-\alpha}))u_{0}-\sum\limits_{j=1}^{k-1}b_{kj}u_{j}\bigg]/(b_{k k}-\lambda),\end{equation}$

(4.20)

不难证明, (4.16)式右边的系数都是正的. 该文从归纳假设中推导出 $u_{k}<\frac{\beta_{k}}{\tilde{\beta}_{k}}u_{0},$ 从引理4.2的(3)得到

$\beta_{k} =-b_{k0}-\lambda \exp(\frac{-\alpha t_{k}}{1-\alpha})-\sum\limits_{j=1}^{k-1}b_{kj}>0, \tilde{\beta}_{k}=b_{k k}-\lambda>0,$

由于方程组(4.5)对于常数解是精确的, 因此得到 $b_{k0}+\sum\limits_{j=1}^{k-1}b_{kj}=0,$ 因此 $\beta_{k}-\tilde{\beta}_{k}<0,$ 所以 $0<u_{k}<u_0.$ 证毕.

5 数值算例

该文下面将给出一些例子来说明本方案的可行性, 该算例所有的计算均是通过软件Mathematica12 实现得到的. 些例子的目的是为了证明所提出的快速高阶方案在收敛速度和CPU 时间方面的准确性和效率, 最后通过算例验证了算法的稳定性和正则性.

${\bf 例1}$ 该文考虑文献[7]中带有下面几个右端项 $f(t,u(t))$ 的初值问题

$\begin{array}{lll} (1) f(t,u(t))=G(t),\\ (2) f(t,u(t))=G(t)-t^3+u(t),\\ (3) f(t,u(t))=G(t)+t^6-u^2(t),\\ \end{array}\\$

其中 $u(t)=t^3, M(\alpha)=1$ , $T=1$ 并且

$G(t)=M(\alpha)\bigg[\frac{3}{\alpha}t^2-\frac{6(1-\alpha)}{\alpha^2}t+\frac{6(1-\alpha)^2}{\alpha^3}-\frac{6(1-\alpha)^2}{\alpha^3}\exp(-\frac{\alpha}{1-\alpha}t)\bigg].$

表1 -表3展示了 $\alpha= 0.3,~0.7$ 时的误差和收敛阶, 并与文献[7]中的算例进行了比较. CPU时间的叠加如图1 -图3所示, 可以看出随着节点数量的增加, 快速方案比直接方案快得多.

表1 当右端项为(1)时的误差和收敛阶比较

新窗口打开| 下载CSV

表2 当右端项为(2)时的误差和收敛阶比较

新窗口打开| 下载CSV

表3 当右端项为(3)时的误差和收敛阶比较

新窗口打开| 下载CSV

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 (1)的CPU时间比较

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 (2)的CPU时间比较

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 (3)的CPU时间比较

${\bf 例 2}$ 下面将对稳定性进行的分析, 为此, 该文将方程(2-1)的右端项改成如下所示.

$f(t)=\frac{M(\alpha)}{\alpha^2+(1-\alpha)^2} \bigg[(1-\alpha)\sin t+\alpha \cos t-\alpha \exp(-\frac{\alpha}{1-\alpha}t)\bigg],$

其中解析解是 $u(t)=\sin t,$ $T = 100$ . 从表4的结果可以看出, 经过长时间的计算, 数值解可以很好地近似解析解. 算例表明, 该方法具有良好的稳定性.

表4 $\alpha=0.3,0.5$ 和 $0.7$ 的最大误差

新窗口打开| 下载CSV

${\bf 例3}$ 为了验证解的正则性, 该文考虑不同的右端项 $f$ , 计算随步长变化的收敛阶数. 考虑计算次数 $T$ 与步长 $h$ 的误差

$E(h)=|x_{h}(T)-x_{\frac{h}{2}}(T)|$

并且收敛阶被定义为

$r=\log_{2}\frac{E(h)}{E(\frac{h}{2})}$

计算结果如表5所示.

表5 $\alpha=0.3$ 时的收敛阶数

新窗口打开| 下载CSV

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Oldham

K B

, Spanier

. The Fractional Calculus. New York: Academic Press, 1974

[本文引用: 1]

[2]

Samko

, Kilbas

A A

, Marichev

. Fractional Integrals and Derivatives: Theory and Applications. Oxfordshire: Taylor & Francis, 1993

[本文引用: 1]

[3]

, Wang

, Hu

Measuring memory with the order of fractional derivative

Sci Rep, 2013, 3: Article number 3431

[本文引用: 1]

[4]

Michele

, Mauro

A new definition of fractional derivative without singular kernel

Progr Fract Differ Appl, 2015, 1(2): 73-85

[本文引用: 1]

[5]

Atangana

, Baleanu

New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model

arXiv preprint arXiv:1602.03408, 2016

[本文引用: 1]

[6]

Cao

, Xu

A high order schema for the numerical solution of the fractional ordinary differential equations

J Comput Phys, 2013, 38: 154-168

[本文引用: 1]

[7]

Cao

, Wang

, Xu

A high-order scheme for fractional ordinary differential equations with the Caputo Fabrizio derivative

Communications on Applied Mathematics and Computation, 2019, 1: 1-21

DOI:10.1007/s42967-019-0010-2 [本文引用: 4]

[8]

Yang

, Huang

, Liang

, Tang

Numerical solution of fractional diffusion-wave equation based on fractional multistep method

Appl Math Model, 2014, 38(14): 3652-3661

DOI:10.1016/j.apm.2013.11.069 URL [本文引用: 1]

[9]

Sun

, Wu

A fully discrete difference scheme for a diffusion-wave system

Appl Numer Math, 2006, 56(2): 193-209

DOI:10.1016/j.apnum.2005.03.003 URL [本文引用: 1]

[10]

Lin

, Xu

Finite difference/spectral approximations for the time-fractional diffusion equation

J Comput Phys, 2007, 225(2): 1533-1552

DOI:10.1016/j.jcp.2007.02.001 URL [本文引用: 1]

[11]

Jiang

, Zhang

Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations

Commun Comput Phys, 2017, 21(3): 650-678

DOI:10.4208/cicp.OA-2016-0136 URL [本文引用: 1]

The computational work and storage of numerically solving the time fractional PDEs are generally huge for the traditional direct methods since they require totalmemory andwork, whereNTandNSrepresent the total number of time steps and grid points in space, respectively. To overcome this difficulty, we present an efficient algorithm for the evaluation of the Caputo fractional derivativeof orderα∈(0,1). The algorithm is based on an efficient sum-of-exponentials (SOE) approximation for the kernelt–1–αon the interval [Δt,T] with a uniform absolute errorε. We give the theoretical analysis to show that the number of exponentialsNexpneeded is of orderforT≫1 orforTH1 for fixed accuracyε. The resulting algorithm requires onlystorage andwork when numerically solving the time fractional PDEs. Furthermore, we also give the stability and error analysis of the new scheme, and present several numerical examples to demonstrate the performance of our scheme.

[12]

Huang

, Tang

, vazquez

Convergence analysis of a block-by-block method for fractional differential equations

Numer Math Theor Methods Appl, 2012, 5(2): 229-241

DOI:10.4208/nmtma URL [本文引用: 1]

[13]

Alikhanov

A A

A new difference scheme for the time fractional diffusion equation

J Comput Phys, 2015, 280: 424-438

DOI:10.1016/j.jcp.2014.09.031 URL [本文引用: 1]

[14]

Baffet

, Hesthaven

A kernel compression scheme for fractional differential equations

SIAM J Numer Anal, 2017, 55(2): 496-520

DOI:10.1137/15M1043960 URL [本文引用: 1]

[15]

Baffet

, Hesthaven

High-order accurate adaptive kernel compression time-stepping schemes for fractional differential equaitons

J Sci Comput, 2017, 72(3): 1169-1195

DOI:10.1007/s10915-017-0393-z URL [本文引用: 1]

[16]

Deng

W H

Finite element method for the space and time fractional Fokker-Planck equation

SIAM J Numer Anal, 2008, 47(1): 204-226

DOI:10.1137/080714130 URL [本文引用: 1]

[17]

R H

, Ng

M K

, Sun

H W

A fast direct method for block triangular Toeplitz-like with tri-diagonal block systems from time-fractional partial differential equations

J Comput Phys, 2015, 303: 203-211

DOI:10.1016/j.jcp.2015.09.042 URL [本文引用: 1]

[18]

Liu

, Shen

, Anh

, Turner

Analysis of a discrete non-Markovian random walk approximation for the time fractional diffusion equation

ANZIAM J, 2005, 46: C488-C504

[本文引用: 1]

[19]

McLean

Fast summation by interval clustering for an evolution equation with memory

J Sci Comput, 2012, 34(6): A3039-A3056

[本文引用: 1]

[20]

Hou

, Xu

. A fractional spectral method with applications to some singular problems Adv Comput Math, 2017, 343(5): 911-944

[本文引用: 1]

[21]

Stynes

, O'Riordan

, Gracia

Error analysis of a finite difference method of graded meses for a time-fractional Caputo-Fabrizio derivatives

Eur Phys J C, 2016, 76: 362

DOI:10.1140/epjc/s10052-016-4209-3 URL [本文引用: 1]

1974

... 该文研究了一类带有非局部边界条件的微分方程数值求解,这类问题常应用于化学工程、热弹性学.在过去的几十年里, 分数微积分受到了物理学家和数学家的广泛关注, 因为它能准确地反映许多应用科学中的物理过程. 分数阶微分方程的数值分析已成为一个深入的研究课题. 分数阶微积分在近年来得到了广泛的研究^[1-2]. 它的研究结果可以说明若干物理问题^[3]. Caputo和Fabrizio^[4]提出了一种新的分数阶导数. 另外, 在文献[5]中给出了具有非局部非奇异核的分数阶导数. 对于一般的右端函数

$f$