High-Order Schemes for Nonlinear Fractional Differential Equations

Abstract

We propose high-order schemes for nonlinear fractional initial value problems. We split the fractional integral into a history term and a local term. We take advantage of the sum of exponentials (SOE) scheme in order to approximate the history term. We also use a low-order quadrature scheme to approximate the fractional integral appearing in the local term and then apply a spectral deferred correction (SDC) method for the approximation of the local term. The resulting one-step time-stepping methods have high orders of convergence, which make adaptive implementation and accuracy control relatively simple. We prove the convergence and stability of the proposed schemes. Finally, we provide numerical examples to demonstrate the high-order convergence and adaptive implementation.

1. Introduction

Fractional differential equations (FDEs) have gained mathematical significance in recent decades due to the growing interest in dynamics of anomalous diffusion processes in amorphous materials. They have governed significant models in many fields of applied sciences and engineering, such as viscoelastic materials, image processing, option pricing models, control theory, etc. Surveys or collections of applications can be found in [1,2,3,4,5].

This paper is concerned with approximating the solution $y\left(t\right)$ to the fractional initial-value problem

$$\begin{array}{cc}\hfill D_t^{\alpha }y\left(t\right)& =f(t,y(t\left)\right),0\le t\le T,\hfill \\ \hfill y\left(0\right)& =y_0,\hfill \end{array}$$

where $y_0,y\left(t\right)\in {\mathbb{R}}^d,f:[0,T]\times {\mathbb{R}}^d\to {\mathbb{R}}^d$, and $D_t^{\alpha }y$ denotes the Caputo fractional derivative of order $\alpha$ ($0<\alpha <1$). Throughout this paper, we assume without loss of generality that $d=1$.

A number of numerical methods have been proposed for solving FDEs [6,7,8,9,10,11]. Two common approaches have been considered by many authors. The first approach is based on a direct discretization of the fractional derivative operators in the considered FDEs [6,7,12]. The second approach to solving FDEs is based on discretizing the integral in the equivalent forms [13,14,15,16]. There are also several numerical methods for FDEs, for example see [17,18,19,20].

Low-order approaches are the most frequently employed methods for FDEs, and it has been proven to be difficult to construct high-order and adaptive schemes [13,15,21,22]. In this paper, we present high-order computational schemes for solving nonlinear FDEs. The proposed schemes deal with integral forms of the FDEs and then split the fractional integral into a history term and a local term. We employ a sum of exponentials scheme (SOE) for the approximation of the history term and low-order quadrature approximations of the fractional integral. We further modify a spectral deferred correction (SDC) scheme to increase the order of the local approximation. In contrast to the methods inspired by [13], here we use an accurate sum of exponentials method as a kernel compression scheme and a completely different deferred correction scheme as a high-order adaptive method for FDEs.

The remainder of the paper is organized as follows. In Section 2, we provide an overview of the proposed methods. Section 3 discusses the sum of exponentials method. Section 4 introduces high-order time-stepping techniques for solving FDEs and discuses the relevant error analysis. Some numerical results are presented in Section 5 before the conclusion is given in Section 6.

2. Overview

We discuss a numerical method for the fractional initial value problem

$$\left\{\begin{array}{c}{D}_t^{\alpha }y\left(t\right)=f(t,y\left(t\right)),0\le t\le T,\\ y\left(0\right)=y_0,\end{array}\right.$$

(1)

where $y_0$ may be an arbitrary real number, and where $0<\alpha <1$. In (1), $D_t^{\alpha }$ denotes the Caputo differential operator, defined by

$$D_t^{\alpha }u\left(t\right)=J^{n-\alpha }{D}^{n}u\left(t\right),$$

where $n:=⌈\alpha ⌉$ is the smallest integer $\ge \alpha$. Here, $D^n$ is the usual differential operator of (integer) order n, and for $\mu >0$, $J^{\mu }$ is the integral operator of order $\mu$ in the sense of Riemann–Liouville, defined by

$$J^{\mu }u\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{t}{\int }}}{\left(t-s\right)}^{\mu -1}u\left(s\right)ds.$$

We assume the function f to be such that a unique solution to (1) exists on some interval $[0,T]$. The initial-value problem (1) is known to be equivalent to the Volterra integral equation

$$y\left(t\right)=y\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{t}{\int }}}{\left(t-s\right)}^{\alpha -1}f(s,y\left(s\right))ds.$$

(2)

Consider the fractional integral

$$I^{\alpha }f\left(t\right)={\displaystyle \underset{0}{\overset{t}{\int }}}{\left(t-s\right)}^{\alpha -1}f\left(s\right)ds.$$

Fix $t\ge 0$ and $h>0$, and let $\gamma \in (0,h]$. It is clear that

$$I^{\alpha }f(t+\gamma )={\displaystyle \underset{0}{\overset{t}{\int }}}{\left(t-s+\gamma \right)}^{\alpha -1}f\left(s\right)ds+{\displaystyle \underset{0}{\overset{\gamma }{\int }}}{\left(\gamma -s\right)}^{\alpha -1}f(t+s)ds.$$

(3)

Owing to (2) and (3), we obtain

$$y(t+\gamma )=H(t,\gamma )+L(t,\gamma ),$$

(4)

where

$$H(t,\gamma )=y\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{t}{\int }}}{\left(t-s+\gamma \right)}^{\alpha -1}f(s,y\left(s\right))ds$$

(5)

and

$$L(t,\gamma )=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{\gamma }{\int }}}{\left(\gamma -s\right)}^{\alpha -1}f(t+s,y(t+s))ds.$$

(6)

We call H and L the history and local parts, respectively.

We can observe that Equation (4) can be written as a Volterra integral equation (VIE),

$$Y\left(\gamma \right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{\gamma }{\int }}}{\left(\gamma -s\right)}^{\alpha -1}F(s,Y\left(s\right))ds+H\left(\gamma \right),$$

(7)

where $Y\left(\gamma \right)=y(t+\gamma )$, $F(s,y)=f(t+s,y)$ and $H\left(\gamma \right)=H(t,\gamma )$.

3. Approximation of the History Term

In this section, we describe a kernel compression scheme for the approximation of the history term. This will reduce the cost of evaluating the history term.

It is proved in [6] that for any specified error $ϵ>0$, we are able to find positive real numbers $s_i^{\alpha }$ and $w_i^{\alpha }$, $i=1,\dots ,N_{ϵ}$, such that

$$\left|\frac{1}{t^{1-\alpha }}-{\displaystyle \sum _{i=1}^{N_{ϵ}}}{w}_i^{\alpha }{e}^{-s_i^{\alpha }t}\right|<ϵ,\forall t\in [h,T],$$

where

$$N_{ϵ}=\mathcal{O}\left(log\frac{1}{ϵ}\left(loglog\frac{1}{ϵ}+log\frac{T}{h}\right)+log\frac{1}{h}\left(loglog\frac{1}{ϵ}+log\frac{1}{h}\right)\right).$$

Roughly, for fixed $ϵ$, the number of exponentials $N_{ϵ}$ needed is of order $N_{ϵ}=\mathcal{O}\left(log{N}_T\right)$ if $T\gg 1$ or $\mathcal{O}\left({log}^2{N}_T\right)$ if $T\approx 1$, where $N_T=\frac{T}{h}$ is the number of time-steps.

Approximating the kernel ${\left(t-s+\gamma \right)}^{\alpha -1}$ in $H(t,\gamma )$ by the sum of exponentials (SOE), we can write the history term (5) as

$$\begin{array}{cc}H(t,\gamma )\hfill & \approx y\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{t}{\int }}}\left({\displaystyle \sum _{i=1}^{N_{ϵ}}}{w}_i^{\alpha }{e}^{-s_i^{\alpha }\left(t-s+\gamma \right)}\right)f(s,y\left(s\right))ds\hfill \\ \hfill & =y_0+{\displaystyle \sum _{i=1}^{N_{ϵ}}}{\sigma }_i^{\alpha ,\gamma }{U}_i^{\alpha }\left(t\right),\hfill \end{array}$$

(8)

where

$${\lambda }_i=-s_i^{\alpha },{\sigma }_i^{\alpha ,\gamma }=w_i^{\alpha }{e}^{{\lambda }_i\gamma }/\Gamma \left(\alpha \right)$$

and

$$U_i^{\alpha }\left(t\right)={\displaystyle \underset{0}{\overset{t}{\int }}}{e}^{{\lambda }_i\left(t-s\right)}f(s,y\left(s\right))ds.$$

Observe that each $U_i^{\alpha }\left(t\right)$ is the solution to the IVP

$$\begin{array}{cc}\hfill \frac{d}{dt}{U}_i^{\alpha }\left(t\right)& ={\lambda }_i{U}_i^{\alpha }\left(t\right)+f,\hfill \\ \hfill U_i^{\alpha }\left(0\right)& =0.\hfill \end{array}$$

To simplify the notation, we consider $U_1^{\alpha }\left(t\right),\dots ,U_{N_{ϵ}}^{\alpha }\left(t\right)$ to be components of a vector $\mathsf{\Phi }=\mathsf{\Phi }\left(t\right)$; similarly, $\Lambda =diag({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{N_{ϵ}})$, $\mathbf{0}$${}_{N_{ϵ}}=\left(0,\dots ,0\right)\in {\mathbb{R}}^{N_{ϵ}}$, and $\mathbf{1}$${}_{N_{ϵ}}=\left(1,\dots ,1\right)\in {\mathbb{R}}^{N_{ϵ}}$. Thus we recover

$${\mathsf{\Phi }}^{\prime }=\mathsf{\Phi }\Lambda +f{\mathbf{1}}_{N_{ϵ}},\mathsf{\Phi }\left(0\right)={\mathbf{0}}_{N_{ϵ}}.$$

(9)

A stable Runge-Kutta (RK) methods may be used for approximating (9).

4. Approximation of the History Term

Let N be a positive integer and T be the given time. We divide the interval $[0,T]$ into

$$0=t_0<t_1<\dots <t_N=T,$$

with step sizes ${\tau }_i=$$t_i-t_{i-1}$, $1\le i\le N$. Let $y_k$ be the approximate solution of $y\left(t_k\right)$, $k=0,1,$ …$,N$.

At $t=t_n$, we consider (7),

$$Y\left(\gamma \right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{\gamma }{\int }}}{\left(\gamma -s\right)}^{\alpha -1}F(s,Y\left(s\right))ds+H\left(\gamma \right),$$

(10)

in $(0,{\tau }_n]$, where $H\left(\gamma \right)$ is given. Divide the interval $[0,{\tau }_n]$ into $0={\gamma }_0\le {\gamma }_1<\dots <{\gamma }_p\le {\tau }_n$, and let V be the numerical approximation to Y. Suppose that we know the approximate solutions $V_k$ at ${\gamma }_k$, for each $k=0,\dots ,j-1$, and let ${\mathcal{P}}_j=\left({\gamma }_1,\dots ,{\gamma }_j\right)$. We compute $V_j$ by solving

$$V_j={\displaystyle \sum _{k=0}^j}{\omega }_{jk}F({\gamma }_k,V_k)+H_j,$$

(11)

for $V_j$, where $H_j=H\left({\gamma }_j\right)$, and ${\omega }_{jk}={\omega }_{jk}\left({\mathcal{P}}_j\right)$ are the weights characterizing the scheme.

4.1. Quadrature Rules on Non-Uniform Meshes

We define

$$t_{n,j}=t_n+{\gamma }_j,Y_{n,j}=Y_n\left({\gamma }_j\right)=y_{n,j}=y(t_n+{\gamma }_j),H_n\left({\gamma }_j\right)=H(t_n,{\gamma }_j).$$

To obtain

$$y\left(t_{n,j}\right)=H(t_n,{\gamma }_j)+L(t_n,{\gamma }_j),$$

we need to approximate the fractional integral

$$L(t_n,{\gamma }_j)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{{\gamma }_j}{\int }}}{\left({\gamma }_j-s\right)}^{\alpha -1}f(t_n+s,y(t_n+s))ds.$$

To do so, we use the following approach

$$L(t_n,{\gamma }_j)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \sum _{k=0}^{j-1}}{\displaystyle \underset{{\gamma }_k}{\overset{{\gamma }_{k+1}}{\int }}}{\left({\gamma }_j-s\right)}_k^{\alpha -1}{\tilde{f}}_k(t_n+s,y(t_n+s))ds,$$

where ${\tilde{f}}_k(t_n+s,y(t_n+s))$, $k=0,1,$ …$,j-1$ is an approximation of $f(t_n+s,y(t_n+s))$, where $s\in [{\gamma }_k,{\gamma }_{k+1})$. Notice that different choices of ${\tilde{f}}_k$ lead to different schemes.

(i) By choosing ${\tilde{f}}_k$ as

$${\tilde{f}}_k(t_n+s,y(t_n+s))=f(t_n+{\gamma }_k,y(t_n+{\gamma }_k))=f(t_{n,k},y\left(t_{n,k}\right)),$$

the fractional rectangle method is derived as

$$L(t_n,{\gamma }_j)={\displaystyle \sum _{k=0}^{j-1}}{\omega }_{jk}f(t_{n,k},y\left(t_{n,k}\right)),$$

where

$$\begin{array}{cc}\hfill {\omega }_{jk}& =\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{{\gamma }_k}{\overset{{\gamma }_{k+1}}{\int }}}{\left({\gamma }_j-s\right)}_k^{\alpha -1}ds\hfill \\ \hfill & =\frac{{\left({\gamma }_j-{\gamma }_k\right)}^{\alpha }-{\left({\gamma }_j-{\gamma }_{k+1}\right)}^{\alpha }}{\Gamma (1+\alpha )},k=0,1,\dots ,j-1.\hfill \end{array}$$

(ii) If ${\tilde{f}}_k$ is selected as

$${\tilde{f}}_k(t_n+s,y(t_n+s))=\frac{t_n+s-{\gamma }_{k+1}}{{\gamma }_k-{\gamma }_{k+1}}f(t_{n,k},y\left(t_{n,k}\right))+\frac{t_n+s-{\gamma }_k}{{\gamma }_{k+1}-{\gamma }_k}f(t_{n,k+1},y\left(t_{n,k+1}\right)),$$

then the fractional trapezoid method is given by

$$L(t_n,{\gamma }_j)={\displaystyle \sum _{k=0}^j}{\tilde{\omega }}_{jk}f(t_{n,k},y\left(t_{n,k}\right))$$

in which

$$\begin{array}{cc}\hfill {\tilde{\omega }}_{jk}& =\frac{1}{\Gamma (2+\alpha )}\left\{\begin{array}{cc}\frac{1}{{\gamma }_1}{A}_0,& \mathrm{if}k=0,\\ \frac{1}{{\gamma }_{k+1}-{\gamma }_k}{A}_k+\frac{1}{{\gamma }_{k-1}-{\gamma }_k}{B}_k,& \mathrm{if}k=1,2,\dots ,j-1,\\ {\left({\gamma }_j-{\gamma }_{j-1}\right)}^{\alpha },& \mathrm{if}k=j,\end{array}\right.\hfill \\ \hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{A}_0={\left({\gamma }_j-{\gamma }_1\right)}^{\alpha +1}-{\gamma }_j^{\alpha +1}+\left(\alpha +1\right){\gamma }_1{\gamma }_j^{\alpha },\\ A_k={\left({\gamma }_j-{\gamma }_{k+1}\right)}^{\alpha +1}-{\left({\gamma }_j-{\gamma }_k\right)}^{\alpha +1}+\left(\alpha +1\right)\left({\gamma }_{k+1}-{\gamma }_k\right){\left({\gamma }_j-{\gamma }_k\right)}^{\alpha },\\ B_k={\left({\gamma }_j-{\gamma }_k\right)}^{\alpha +1}-{\left({\gamma }_j-{\gamma }_{k-1}\right)}^{\alpha +1}+\left(\alpha +1\right)\left({\gamma }_k-{\gamma }_{k-1}\right){\left({\gamma }_j-{\gamma }_k\right)}^{\alpha }.\end{array}\right.\hfill \\ \hfill \end{array}$$

It is proved in [8] that for $g\in C^1[0,T]$, we have that

$$\left|\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{{\gamma }_j}{\int }}}{\left({\gamma }_j-s\right)}^{\alpha -1}g\left(s\right)ds-{\displaystyle \sum _{k=0}^{j-1}}{\omega }_{jk}g\left({\gamma }_k\right)\right|\le \frac{{∥g^{\prime }∥}_{\infty }}{\Gamma (1+\alpha )}{T}^{\alpha }\underset{k}{max}{\tau }_k$$

and for $g\in C^2[0,T]$,

$$\left|\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{{\gamma }_j}{\int }}}{\left({\gamma }_j-s\right)}^{\alpha -1}g\left(s\right)ds-{\displaystyle \sum _{k=0}^j}{\tilde{\omega }}_{jk}g\left({\gamma }_k\right)\right|\le \frac{{∥g^{″}∥}_{\infty }}{2\Gamma (1+\alpha )}{T}^{\alpha }{\left(\underset{k}{max}{\tau }_k\right)}^2.$$

4.2. The Spectral Deferred Correction Framework

Let ${\tau }_n$ denote the desired time-step size. Then, (4) can be written as

$$y(t+\gamma )=H(t,\gamma )+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{\gamma }{\int }}}{\left(\gamma -s\right)}^{\alpha -1}f(t+s,y(t+s))ds,$$

(12)

where $H(t,\gamma )$ is given. Let $\mathcal{P}={\left\{{\gamma }_k\right\}}_{k=1}^p$ be a set of distinct points in the closed interval $[0,{\tau }_n]$, where $0={\gamma }_0\le {\gamma }_1<\dots <{\gamma }_p\le {\tau }_n$.

Replacing the integral in (12) by a quadrature from Section 4.1 yields the collocation approximation

$$y(t+{\gamma }_j)=H(t,{\gamma }_j)+{\displaystyle \sum _{k=0}^j}{\omega }_{jk}f(t+{\gamma }_k,y(t+{\gamma }_k)).$$

This can be written as

$$Y\left({\gamma }_j\right)=H\left({\gamma }_j\right)+{\displaystyle \sum _{k=0}^j}{\omega }_{jk}F({\gamma }_k,Y\left({\gamma }_k\right)),$$

(13)

where $Y\left(\gamma \right)$$=y(t+\gamma )$, $F(\gamma ,y)=f(t+\gamma ,y)$ and $H\left(\gamma \right)$ substitutes for $H(t,\gamma )$. At $t=t_n$, the values $Y_{n,j}$ are defined explicitly by the formula

$$Y_{n,j}=H_n\left({\gamma }_j\right)+{\displaystyle \sum _{k=0}^{j-1}}{\omega }_{jk}f(t_{n,k},Y_{n,k})$$

or implicitly through the nonlinear system

$$Y_{n,j}=H_n\left({\gamma }_j\right)+{\displaystyle \sum _{k=0}^j}{\tilde{\omega }}_{jk}f(t_{n,k},Y_{n,k}).$$

To obtain higher-order schemes, we can apply a spectral deferred correction (SDC) scheme [23,24,25].

While the time is stepping from $t_n$ to $t_{n+1}=t_n+{\tau }_n$, an SDC method subdivides the interval $[t_n,t_{n+1}]$ into p substeps

$$t_n\le t_{n,1}<t_{n,2}<\dots <t_{n,p}\le t_{n+1}.$$

At each substep $t_{n,j}$, the method computes a provisional solution using a low-order method with initial condition $y_n$. We will refer to the provisional solution as

$$Y_{n,j}^{\left[1\right]}\approx y\left(t_{n,j}\right),j=1,2,\dots ,p.$$

To improve the accuracy of the provisional solution, we perform a process of iterated corrections of the form

$$Y_{n,j}^{[k+1]}\approx Y_{n,j}^{\left[k\right]}+E_{n,j}^{\left[k\right]},$$

(14)

where

$$E_{n,j}^{\left[k\right]}\approx y\left(t_{n,j}\right)-Y_{n,j}^{\left[k\right]}.$$

Consider again the integral formulation of the solution,

$$Y\left(\gamma \right)=H\left(\gamma \right)+{\displaystyle \underset{0}{\overset{\gamma }{\int }}}\frac{1}{\Gamma \left(\alpha \right)}{\left(\gamma -s\right)}^{\alpha -1}F(s,Y\left(s\right))ds.$$

(15)

Let $\tilde{Y}\left(\gamma \right)$ be an approximate solution to (15). Then, plugging $Y\left(\gamma \right)=E\left(\gamma \right)+\tilde{Y}\left(\gamma \right)$ into (15) yields the integral equation

$$E\left(\gamma \right)=-\tilde{Y}\left(\gamma \right)+H\left(\gamma \right)+{\displaystyle \underset{0}{\overset{\gamma }{\int }}}\frac{1}{\Gamma \left(\alpha \right)}{\left(\gamma -s\right)}^{\alpha -1}F(s,\tilde{Y}\left(s\right)+E\left(s\right))ds.$$

Consequently, the error function can be written as

$$\begin{array}{cc}\hfill E\left(\gamma \right)& =\underset{\mathrm{Error}\mathrm{Term}}{\underbrace{{\displaystyle \underset{0}{\overset{\gamma }{\int }}}\frac{1}{\Gamma \left(\alpha \right)}{\left(\gamma -s\right)}^{\alpha -1}G(s,E\left(s\right))ds}}\hfill \\ \hfill & +\underset{\mathrm{Residual}\mathrm{Term}R(\gamma ,H(\gamma \left)\right)}{\underbrace{\left[H\left(\gamma \right)+{\displaystyle \underset{0}{\overset{\gamma }{\int }}}\frac{1}{\Gamma \left(\alpha \right)}{\left(\gamma -s\right)}^{\alpha -1}F(s,\tilde{Y}\left(s\right))ds\right]-\tilde{Y}\left(\gamma \right)}},\hfill \end{array}$$

(16)

where

$$G(s,E\left(s\right))=F(s,\tilde{Y}\left(s\right)+E\left(s\right))-F(s,\tilde{Y}\left(s\right)).$$

To compute the error

$$E_{n,j}^{\left[k\right]}\approx y\left(t_{n,j}\right)-Y_{n,j}^{\left[k\right]},$$

we replace $\tilde{Y}\left(s\right)$ with the Lagrange interpolating polynomial $L^{\left[k\right]}\left(s\right)$,

$$L^{\left[k\right]}\left(t\right)={\displaystyle \sum _{i=1}^p}{Y}_{n,i}^{\left[k\right]}{l}_i\left(t\right),\mathrm{where}{l}_i\left(t\right)={\displaystyle \prod _{\underset{m\ne k}{m=1}}^p}\frac{t-t_{n,m}}{t_{n,i}-t_{n,m}}.$$

To approximate the error, we first use a quadrature formula with nodes at ${\left\{{\gamma }_j\right\}}_{j=1}^p$ to approximate the integral appearing in the residual term $R(\gamma ,H(\gamma \left)\right)$, and then we use a fractional implicit or explicit Euler approximation to approximate the remaining integral term in (16). This process leads to the implicit ($\eta =1$) or the explicit ($\eta =0$) time-stepping method

$$E_{n,j+1}^{\left[k\right]}=E_{n,j}^{\left[k\right]}+\frac{{\left(h_{n,j}^{}\right)}^{\alpha }}{\Gamma (1+\alpha )}{G}_{n,j+\eta }^{\left[k\right]}+R_{n,j}^{\left[k\right]},$$

(17)

where

$$\begin{array}{cc}\hfill h_{n,j}& =t_{n,j+1}-t_{n,j},\hfill \\ \hfill G_{n,j+\eta }^{\left[k\right]}& =f(t_{n,j+\eta },Y_{n,j+\eta }^{\left[k\right]}+E_{n,j+\eta }^{\left[k\right]})-f(t_{n,j+\eta },Y_{n,j+\eta }^{\left[k\right]}),\hfill \end{array}$$

$$R_{n,j}^{\left[k\right]}=Y_{n,j}^{\left[k\right]}-Y_{n,j+1}^{\left[k\right]}+H_{n,j+1}-H_{n,j}+{\displaystyle \sum _{i=1}^p}{\omega }_{ji}f(t_{n,i},Y_{n,i}^{\left[k\right]}),$$

and

$$H_{n,j}=H(t_n,{\gamma }_j).$$

During the iteration process, we always use the initial conditions,

$$t_{n,0}=t_n,Y_{n,0}^{\left[k\right]}=y_n,E_{n,0}^{\left[k\right]}=0.$$

From now on, we will refer to the implicit method as IM-SDC and the explicit method as EX-SDC.

Practically, we start by using the fractional Euler’s method of order $\alpha$ to compute the provisional solution and then repeat the iteration (17) followed by (14) a total of m times to obtain solutions with an increasingly high order of accuracy

$${\left\{Y_{n,j}^{\left[1\right]}\right\}}_{j=1}^p,{\left\{Y_{n,j}^{\left[2\right]}\right\}}_{j=1}^p,\dots ,{\left\{Y_{n,j}^{[m+1]}\right\}}_{j=1}^p.$$

The scheme can be rewritten under a more compact form

$$Y_{n,j+1}^{[k+1]}=Y_{n,j}^{[k+1]}+\frac{{\left(h_{n,j}^{}\right)}^{\alpha }}{\Gamma (1+\alpha )}\left(f_{n,j+\eta }^{[k+1]}-f_{n,j}^{\left[k\right]}\right)+\left(H_{n,j+1}-H_{n,j}\right)+{\displaystyle \sum _{i=1}^p}{\omega }_{ji}{f}_{n,i}^{\left[k\right]},$$

where $f_{n,i}^{\left[k\right]}=f(t_{n,i},Y_{n,i}^{\left[k\right]})$.

4.3. Order of Convergence

Suppose that the integrator used to construct the provisional solution has order $p_0$, the integrator used to compute the correction approximation at level k has order $p_k$, and m correction steps are taken. If we have p groups of N uniformly distributed nodes, then the SDC methods can attain order [26,27]

$$O\left(h^{min(P_m,p+1)}\right),$$

where $P_m={\displaystyle \sum _{k=0}^m}{p}_k$.

The convergence rate of our fully discrete time-stepping SDC method is given in the following theorem.

Theorem 1.

For any sufficiently smooth function $f:[0,T]\times \mathbb{R}\to \mathbb{R}$, the SDC method outlined in the previous section using $p+1$ distinct quadrature nodes and m correction steps converges to the exact solution $(Y\left({\gamma }_0\right),\dots ,Y\left({\gamma }_p\right))$ with order of accuracy

$$O\left(h^{min\left\{p+1,(m+1)\alpha \right\}}\right).$$

5. Adaptive Implementation

Adaptive marching and accuracy control play a significant role in practical applications involving numerical schemes for solving differential equations. The schemes proposed in this paper are essentially one-step methods, making adaptive implementation relatively simple. Additionally, the schemes have a high order of convergence, which makes accuracy control much easier.

In what follows in this section, we describe details of an adaptive implementation of the SDC schemes and illustrate a number of numerical applications we have performed.

5.1. Accuracy Control

Given a fractional IVP (1) on the interval $[0,T]$, an approximate solution $Y_n^{[m+1]}$,

$$Y_n^{[m+1]}=\left(Y_{n,1}^{[m+1]},Y_{n,2}^{[m+1]},\dots ,Y_{n,p}^{[m+1]}\right)\approx \left(y\left(t_{n,1}\right),y\left(t_{n,2}\right),\dots ,y\left(t_{n,p}\right)\right)$$

and a positive $ϵ$, we would like to determine whether

$$∥Y_n^{[m+1]}-y∥<ϵ.$$

(18)

The structure of the proposed SDC method provides us with a number of completely reliable conditions which can be used as stopping criteria:

• We verify that $∥E_n^{[m+1]}∥<ϵ$, where the vector $E_n^{[m+1]}$ (see (17)),

  $$E_n^{[m+1]}=\left(E_{n,1}^{[m+1]},E_{n,2}^{[m+1]},\dots ,E_{n,p}^{[m+1]}\right)$$

  (19)

  is obtained during the last correction.
• Once $Y_n^{[m+1]}$ has been computed; the value of the approximate solution at an arbitrary time $t\in$$[t_n,t_{n+1}]$ can be obtained by the Lagrange interpolant. We apply the interpolation process to both $Y_n^{[m+1]}$ and $Y_n^{\left[m\right]}$ and demand that the difference be less than $ϵ$. This indicates that both the correction process and the discretization have converged to precision.

We can also control the step-size h. We start with an arbitrary step-size h and attempt to apply the SDC scheme. The step-size is halved and the procedure is repeated if the obtained precision is insufficient (per criteria (18) and (19)). However, the step-size is unchanged if the precision is acceptable and doubled if the precision is acceptable two steps in a row.

5.2. Computational Details and Numerical Tests

In what follows, we present some numerical results of our numerical method applied to nonlinear fractional differential equations.

We solve all test problems using the proposed SDC methods: implicit SDC (IM-SDC) or explicit SDC (EX-SDC). While using the implicit scheme, we use Newton’s method to solve the resulting system of nonlinear equations. Moreover, if the exact solution to the given problem is not known, we use reference solutions computed by the IM-SDC with error tolerance $ϵ={10}^{-10}$ instead.

As a first example, we study the nonlinear fractional differential equation

$$\left\{\begin{array}{c}{D}_t^{\alpha }y\left(t\right)=\frac{40320}{\Gamma (9-\alpha )}{t}^{8-\alpha }-\frac{3\Gamma (5+\alpha /2)}{\Gamma (5-\alpha /2)}{t}^{4-\alpha /2}+{\left(\frac{3}{2}{t}^{\alpha /2}-t^4\right)}^3\\ -y{\left(t\right)}^{3/2}+\frac{9\Gamma (\alpha +1)}{4},\\ y\left(0\right)=0.\end{array}\right.$$

(20)

The exact solution is $y\left(t\right)=t^8-3t^{4+\alpha /2}+\frac{9}{4}{t}^{\alpha }$.

We apply the IM-SDC method to (20) with $\alpha =0.75$, $p=5$, $m=3$, $ϵ={10}^{-6}$, and initial step-size $h_0=2^{-4}$. When $m=0$, the local order of the method is $1+\alpha =3/2$. When we perform correction iterations, i.e., $m\ge 1$, we measure a high order of convergence. For example, with $m=3$, we obtain order three which is in agreement with the expected theoretical results. We present our numerical results in Figure 1. The results show that high accuracy can be obtained using few correction iterations after applying high-order SDC methods.

We also studied the fractional-order Chua’s circuit system [28]

$$\begin{array}{cc}\hfill D_t^{\alpha }x\left(t\right)& ={\beta }_1\left(y\left(t\right)-H\left(x\right(t\left)\right)\right),\hfill \\ \hfill D_t^{\alpha }y\left(t\right)& =x\left(t\right)-y\left(t\right)+z\left(t\right),\hfill \\ \hfill D_t^{\alpha }z\left(t\right)& =-{\beta }_{2}y\left(t\right),\hfill \end{array}$$

where ${\beta }_1$ and ${\beta }_2$ are constant parameters, and $H(·)$ is given by

$$H\left(v\right)=-\frac{1}{7}v+\frac{2}{7}{v}^3.$$

In this test, $\alpha =0.9$, ${\beta }_1=12.5$, ${\beta }_2=100/7$, $p=6$, $m=5$, $ϵ={10}^{-6}$, $\left(x_1\left(0\right),x_2\left(0\right),x_3\left(0\right)\right)=(0,0.1,0)$, and $T=50$. The approximations $Y_n^{[m+1]}$ are compared with a reference solution $Y_{ref}$ obtained with $ϵ={10}^{-10}$. Figure 2 shows the reference solutions $x_{ref},y_{ref},z_{ref}$ and the curve $\left(x_{ref},z_{ref}\right)$ in the $\left(x,z\right)$ plane.

In Figure 3, the accuracy is plotted as a function of the average time-step $h_{avg}$.

6. Conclusions

We proposed high-order adaptive time-stepping schemes for nonlinear FDEs by using suitable kernel compression and SDC methods. We discussed the accuracy and local convergence of the proposed schemes. We also discussed the extension of these schemes to construct other higher-order methods for nonlinear fractional initial-value problems.

We provided numerical examples to verify the efficiency of the proposed schemes, which shows that applying suitable correction iterations can increase the accuracy in a significant way.

In future work, we will focus on proposing high-order schemes for the numerical solutions of fractional boundary value problems and time-fractional diffusion equations.