Higher-order fractional linear multi-step methods

In the literature on fractional differential and integral equations, there are several analytical methods to solve these equations. We can mention here the Fourier, Laplace, or Mellin transform methods and the Green function method. But these analytical methods have been used to solve very special types, mostly linear, of FDEs. Therefore, developing reliable and efficient numerical methods for solving general fractional differential equations, such as nonlinear types, is an attractive topic for researchers [1–4]. One can find numerical methods based on polynomial interpolation [5, 6], Gauss interpolation [7], Grunwald-Letnikov approximations, and generalizations of linear multistep methods. Firstly, these methods have been used for solving fractional differential equations containing the Riemann-Liouville derivative, and then natural generalizations of the above methods are introduced to solve equations with Caputo's derivative, the Riesz derivative, and other types of fractional derivatives. Some recent works related to the several computational methods to solve various types of fractional-order differential equations can be seen in ref. [8–11]. In [12], the authors proposed some novel analyses of multistep methods for fractional differential equations. In [13], a study on the fractional linear multistep methods for Abel-Volterra integral equations of the second kind has been proposed. Fractional linear multistep methods are a generalization of classical linear multistep methods for fractional differential equations.

A linear k-step method for solving initial value problem

$$\begin{eqnarray*}&&\left\{\begin{array}{l}y^{\prime} (t)=f(t,y(t)),\\ y({t}_{0})={y}_{0},\end{array}\right.\end{eqnarray*}$$

is of the form

$$\begin{eqnarray}&&\displaystyle \sum _{j=0}^{k}{\alpha }_{j}{y}_{n+j}=h\displaystyle \sum _{j=0}^{k}{\beta }_{j}{f}_{n+j}.\end{eqnarray} \tag{ 1 }$$

For the fractional initial value problem [14–16]

$$\begin{eqnarray}&&\left\{\begin{array}{l}{D}^{\beta }y(t)=f(t,y(t)),\\ \quad y({t}_{0})={y}_{0},\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where 0 < β < 1, we have

$$\begin{eqnarray}&&y(t)={y}_{0}+\displaystyle \frac{1}{{\rm{\Gamma }}(\beta )}{\int }_{{t}_{0}}^{t}{\left(t-s\right)}^{\beta -1}f(s,y(s)){ds}.\end{eqnarray} \tag{ 3 }$$

Using the following discrete convolution quadrature

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\rm{\Gamma }}(\beta )}{\int }_{{t}_{0}}^{{t}_{n}}\displaystyle \frac{f(s,y(s))}{{\left(t-s\right)}^{1-\beta }}{ds}\approx {h}^{\beta }\displaystyle \sum _{j=0}^{n}{\omega }_{j}f({t}_{n-j},{y}_{n-j})+{h}^{\beta }\displaystyle \sum _{j=0}^{s}{w}_{n,j}f({t}_{j},{y}_{j}),\end{eqnarray} \tag{ 4 }$$

we can obtain fractional linear multi-step methods to solve (2). Therefore, the general fractional form of a linear multi-step method is the form

$$\begin{eqnarray}&&{y}_{n}={y}_{0}+{h}^{\beta }\displaystyle \sum _{j=0}^{n}{\omega }_{j}f({t}_{n-j},{y}_{n-j})+{h}^{\beta }\displaystyle \sum _{j=0}^{s}{w}_{n,j}f({t}_{j},{y}_{j}).\end{eqnarray} \tag{ 5 }$$

Lubich [16, 17] proved that if ρ(ξ) and σ(ξ) are the first and second characteristic polynomials of a p-th order linear multistep method (1) and if the zeros of ρ(ξ) have absolute value less then 1, then the associated fractional form of the method, (5), would be convergent of order p. In (5), ω_j are the coefficients of the Taylor expansion of the so-called generating function ${\omega }^{\beta }(\xi )={\left(\tfrac{\sigma (\tfrac{1}{\xi })}{\rho (\tfrac{1}{\xi })}\right)}^{\beta },$ and the starting weights ω_n,j are chosen such that the method is exact when f(t) = t^μ , μ < p. It has been proved in [18], that the exact solution of (2) has the following form

$$\begin{eqnarray*}&&y(t)={y}_{0}+\displaystyle \sum _{\nu \in {{ \mathcal N }}^{* }}{y}_{\nu }{\left(t-{t}_{0}\right)}^{\nu }+{ \mathcal O }({\left(t-{t}_{0}\right)}^{q}),\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{{ \mathcal N }}^{* }:= \{\nu ={l}_{1}+\beta {l}_{2}| {l}_{1},{l}_{2}\in {\mathbb{N}}\,\,\,\,{and}\,\,\,\,\nu \in (0,q)\},\end{eqnarray*}$$

which is an asymptotic expansion of mixed powers of (t − t₀) and ${\left(t-{t}_{0}\right)}^{\beta }.$ Therefore, the starting weights w_n,j are computed by the assumption that (5) is exact for ${\left(t-{t}_{0}\right)}^{\nu },\,\nu \in { \mathcal A },\,{ \mathcal A }=\{\nu =j\,+\,l\beta :l\in \{0,1,\ldots \},\,\nu \leqslant p-1\},\,p\gt 0$ is the order of consistency and for the starting weights, we have a system of s + 1 linear equations and this number is equal to the cardinality of A, $s={Card}{ \mathcal A }-1.$ On the other hand, these forms of starting weights cause the errors y(t_n ) − y_n behave like h^k−ε(β), where 0 ≤ ε(β) < 1 − β, and so, due to the form of the exact solution, adding the terms of second summation in (5) increases the accuracy of the method. As one can see, in fractional generalizations of linear multi-step methods, we have to compute β-power of some functions. For this purpose, using the following Miller formula has great importance.

Lemma 1.1. [19] Let $\varphi (\xi )=1+{\sum }_{j=1}^{\infty }{a}_{j}{\xi }^{j}$ be a formal power series. Then for any $\beta \in {\mathbb{C}},$

$$\begin{eqnarray*}&&{\left(\varphi (\xi )\right)}^{\beta }=\displaystyle \sum _{n=0}^{\infty }{\nu }_{n}^{(\beta )}{\xi }^{n},\end{eqnarray*}$$

where coefficients ${\nu }_{n}^{(\beta )}$ can be recursively calculated as

$$\begin{eqnarray*}&&{\nu }_{0}^{(\beta )}=1,\,\,{\nu }_{n}^{(\beta )}=\displaystyle \sum _{j=1}^{n}\left(\displaystyle \frac{(\beta +1)j}{n}-1\right){a}_{j}{\nu }_{n-j}^{(\beta )},\quad n=1,2,\cdots .\end{eqnarray*}$$

$$\begin{eqnarray}&&{y}_{n}-{y}_{n-1}=h\displaystyle \sum _{j=0}^{k}{\bar{\gamma }}_{j}{{\rm{\nabla }}}^{j}{f}_{n},\end{eqnarray} \tag{ 6 }$$

where ∇ denotes the backward difference operator, ∇f_n = f_n − f_n−1, the coefficients of a k-step fractional Adams-Moulton method can be obtained by the following generating function [15]

$$\begin{eqnarray*}&&{g}_{\beta }(\xi )=\displaystyle \frac{{\left(G(1-\xi )\right)}^{\beta }}{{\left(1-\xi \right)}^{\beta }},\end{eqnarray*}$$

where $G(t)=\tfrac{-t}{\mathrm{ln}(1-t)}$ is the generating function of the coefficients $\bar{{\gamma }_{j}},$ in (6). Similar to (5), we can obtain these methods as follows

$$\begin{eqnarray}&&\displaystyle \sum _{j=0}^{n}{\omega }_{j}({y}_{n-j}-{y}_{0})+\displaystyle \sum _{j=0}^{s}{w}_{n,j}({y}_{j}-{y}_{0})={h}^{\beta }\displaystyle \sum _{j=0}^{k}{\bar{\gamma }}_{j}^{(k)}{f}_{n-j},\end{eqnarray} \tag{ 7 }$$

in which w_j can be evaluated by using ${\left(1-\xi \right)}^{\beta }={\sum }_{j=0}^{\infty }{\omega }_{j}{\xi }^{j}$ and

$$\begin{eqnarray}&&{\left(G(1-\xi )\right)}^{\beta }=\displaystyle \sum _{j=0}^{\infty }{\gamma }_{j}{\left(1-\xi \right)}^{j}=\displaystyle \sum _{j=0}^{\infty }{\bar{\gamma }}_{j}^{(k)}{\xi }^{j}.\end{eqnarray} \tag{ 8 }$$

In a similar way, in the case of explicit Adams-Bashforth methods

$$\begin{eqnarray*}&&{y}_{n}-{y}_{n-1}=\displaystyle \sum _{j=0}^{k}{\gamma }_{j}^{* }{{\rm{\nabla }}}^{j}{f}_{n-j},\end{eqnarray*}$$

where the generating function of the coefficients ${\gamma }_{j}^{* }$ is ${G}^{* }(t)=\tfrac{-t}{(1-t)\mathrm{ln}(1-t)},$ instead of using β-power of this function which presents some technical difficulties, Garrappa [14] considered the generating function

$$\begin{eqnarray*}&&{\delta }^{* }(\xi )={\left(1-\xi \right)}^{-\beta }\left({\gamma }_{0}+{\gamma }_{1}(1-\xi )+\cdots +{\gamma }_{k-1}{\left(1-\xi \right)}^{k-1}+{c}_{k}{\left(1-\xi \right)}^{k}\right),\end{eqnarray*}$$

where c_k = − (γ₀ + γ₁ + ... + γ_k−1) to retain the explicitness. Then the fractional form of Adams-Bashforth methods can be written as the following

$$\begin{eqnarray}&&\displaystyle \sum _{j=0}^{n}{\omega }_{j}({y}_{n-j}-{y}_{0})+\displaystyle \sum _{j=0}^{s}{w}_{n,j}({y}_{j}-{y}_{0})={h}^{\beta }\displaystyle \sum _{j=1}^{k}{\gamma }_{j}^{* (k)}{f}_{n-j},\end{eqnarray} \tag{ 9 }$$

where ω_j , w_n,j are the same as in (7) and

$$\begin{eqnarray*}&&{\gamma }_{j}^{* (k)}={\left(-1\right)}^{j}\left(\displaystyle \sum _{l=j}^{k-1}\left({}_{j}^{l}\right){\gamma }_{l}-\left({}_{j}^{k}\right)\displaystyle \sum _{l=0}^{k-1}{\gamma }_{l}\right),\,\,\,j=1,2,\cdots ,k,\end{eqnarray*}$$

such that γ_j are defined in (8).

Now let us consider the following result.

Theorem 1.2. [15] A k-step FAM method (7) is consistent and convergent with order $p=k+1$.

Previous research works on fractional linear multi-step methods were limited to the fractionalization of linear multi-step methods in Adams-Moulten or Adams-Bashforth forms [14–16]. By introducing the arrays of coefficients of Adams-type methods in this paper, it is possible to construct the fractional form of general linear multi-step methods. Since these new fractional methods are written as a linear combination of Adams methods, by properly choosing the value of the parameters appearing in them, the order can be increased and the stability can be improved.

The next portion of this paper is framed as follows: In the next section 2, we propose two arrays to produce the coefficients of all fractional Adams-Moulton and fractional Adams-Bashforth methods. Also, we provide recursive relations to obtain the elements of these arrays. In section 3, the construction of higher-order fractional linear multi-step methods is given. In section 4, the analyses of linear stability are performed. In section 5, an example is given to justify the validity and correctness of the proposed method. Section 6 concludes the findings of the paper.

2.1. Fractional Adams-Moulton array

Let us consider a k-step fractional Adams-Moulton method of the form

$$\begin{eqnarray}&&\displaystyle \sum _{j=0}^{n}{\omega }_{j}({y}_{n-j}-{y}_{0})+\displaystyle \sum _{j=0}^{s}{w}_{n,j}({y}_{j}-{y}_{0})={h}^{\beta }\displaystyle \sum _{j=0}^{k}{V}_{j}^{(k)}{f}_{n-j}.\end{eqnarray} \tag{ 10 }$$

Firstly, we prove a recursive relation for the coefficients ${V}_{j}^{(k)}.$ For this purpose, let us start with

$$\begin{eqnarray}\begin{array}{rcl}{\left(G(t)\right)}^{\beta } & = & {\left(\displaystyle \frac{-t}{\mathrm{ln}(1-t)}\right)}^{\beta }={V}_{0}^{(0)}+{V}_{1}^{(0)}t+{V}_{2}^{(0)}{t}^{2}+\cdots \\ & = & \displaystyle \sum _{n=0}^{\infty }{V}_{n}^{(0)}{t}^{n},\end{array}\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray*}&&{V}_{0}^{(0)}=1,\,\,{V}_{n}^{(0)}=\displaystyle \sum _{j=1}^{n}\left(\displaystyle \frac{(\beta +1)j}{n}-1\right){V}_{n-j}^{(0)},\hspace{0.6cm}n=1,2,\cdots ,\end{eqnarray*}$$

and we can easily obtain the first terms in the expansion of ${\left(G(\xi )\right)}^{\beta }$ as

$$\begin{eqnarray*}&&{V}_{0}^{(0)}=1,\hspace{0.1cm}{V}_{1}^{(0)}=-\displaystyle \frac{\beta }{2},\hspace{0.1cm}{V}_{2}^{(0)}=\displaystyle \frac{1}{8}{\beta }^{2}-\displaystyle \frac{5}{24}\beta ,\hspace{0.1cm}\cdots .\end{eqnarray*}$$

Lemma 2.1. For the coefficients ${V}_{j}^{(k)}$ appeared in (10), one has the following recursive relation

$$\begin{eqnarray}\left\{\begin{array}{ll}{V}_{l}^{(k)}={V}_{l}^{(k-1)}+{\left(-1\right)}^{l}\left({}_{l}^{k}\right){V}_{k}^{(0)}, & l\lt k,\\ {V}_{l}^{(k)}={\left(-1\right)}^{l}\left({}_{l}^{k}\right){V}_{k}^{(0)}, & l=k.\end{array}\right.\end{eqnarray} \tag{ 12 }$$

Proof. As we mentioned in previous section, the coefficients $\bar{{\gamma }_{j}}$ can be obtained by

$$\begin{eqnarray*}&&{\left(G(1-\xi )\right)}^{\beta }=\displaystyle \sum _{j=0}^{\infty }\bar{{\gamma }_{j}}{\left(1-\xi \right)}^{j}.\end{eqnarray*}$$

Now, according to this relation and by choosing k-terms to obtain the coefficients of k-step Adams-Moulton method, we have

$$\begin{eqnarray*}&&{\left({G}_{k}(1-\xi )\right)}^{\beta }={V}_{0}^{(0)}+{V}_{1}^{(0)}(1-\xi )+\cdots +{V}_{k}^{(0)}{\left(1-\xi \right)}^{k},\end{eqnarray*}$$

which, expanding ${\left(1-\xi \right)}^{l}$ terms, we conclude

$$\begin{eqnarray*}\begin{array}{rcl}{\left({G}_{k}(1-\xi )\right)}^{\beta } & = & \left({V}_{0}^{(0)}+{V}_{1}^{(0)}+\cdots +{V}_{k}^{(0)}\right)+\left(-{V}_{1}^{(0)}-2{V}_{2}^{(0)}-\cdots -{{kV}}_{k}^{(0)}\right)\xi \\ & & +\left({V}_{2}^{(0)}+3{V}_{3}^{(0)}+\cdots \right){\xi }^{2}+\cdots +{\left(-1\right)}^{k}{V}_{k}^{(0)}{\xi }^{k}\\ & = & {V}_{0}^{(k)}+{V}_{1}^{(k)}\xi +\cdots +{V}_{k}^{(k)}{\xi }^{k},\end{array}\end{eqnarray*}$$

where ${V}_{0}^{(k)},{V}_{1}^{(k)},\ldots ,{V}_{k}^{(k)}$ are the coefficients of the k-step fractional Adams-Moulton method.

On the other hand, since

$$\begin{eqnarray*}&&{\left({G}_{k-1}(1-\xi )\right)}^{\beta }={V}_{0}^{(0)}+{V}_{1}^{(0)}(1-\xi )+\cdots +{V}_{k-1}^{(0)}{\left(1-\xi \right)}^{k-1},\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\left({G}_{k}(1-\xi )\right)}^{\beta }={\left({G}_{k-1}(1-\xi )\right)}^{\beta }+{V}_{k}^{(0)}{\left(1-\xi \right)}^{k},\end{eqnarray*}$$

and so substituting yields

$$\begin{eqnarray*}&&{V}_{0}^{(k)}+{V}_{1}^{(k)}\xi +\cdots +{V}_{k-1}^{(k)}{\xi }^{k-1}+{V}_{k}^{(k)}{\xi }^{k}={V}_{0}^{(k-1)}+{V}_{1}^{(k-1)}\xi +\cdots +{V}_{k-1}^{(k-1)}{\xi }^{k-1}+{V}_{k}^{(0)}{\left(1-\xi \right)}^{k},\end{eqnarray*}$$

which equating coefficients of ${\xi }^{l}$ gives the recursive formula (12). □

Therefore, knowing the coefficients (11) enables us to easily write an array of the coefficients, ${V}_{j}^{k},$ of a k-step fractional Adams-Moulton method. We shall christen this array the fractional Adams-Moulton arrary. The array is shown in table 1. In order to clarify the uses of the given relation (12) between the coefficients of table 1, we consider for example k = 2, the 2-step fractional Adams-Moulton coefficients are defined as follows:

$$\begin{eqnarray*}&&\displaystyle \sum _{j=0}^{n}{\omega }_{j}({y}_{n-j}-{y}_{0})+\displaystyle \sum _{j=0}^{s}{w}_{n,j}({y}_{j}-{y}_{0})={h}^{\beta }\left({V}_{0}^{(2)}{f}_{n}+{V}_{1}^{(2)}{f}_{n-1}+{V}_{2}^{(2)}{f}_{n-2}\right),\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{V}_{0}^{(2)} & = & {V}_{0}^{(1)}+\left({}_{0}^{2}\right){V}_{2}^{(0)}\\ {V}_{1}^{(2)} & = & {V}_{1}^{(1)}-\left({}_{1}^{2}\right){V}_{2}^{(0)}\\ {V}_{2}^{(2)} & = & \left({}_{2}^{2}\right){V}_{2}^{(0)}.\end{array}\end{eqnarray*}$$

Table 1. Fractional Adams-Moulton array.

	${V}_{0}^{(0)}$	${V}_{1}^{(0)}$	${V}_{2}^{(0)}$	⋯	${V}_{l-1}^{(0)}$	${V}_{l}^{(0)}$
k = 0	${V}_{0}^{(0)}$
k = 1	${V}_{0}^{(1)}$	${V}_{1}^{(1)}$
k = 2	${V}_{0}^{(2)}$	${V}_{1}^{(2)}$	${V}_{2}^{(2)}$
⋮	⋮	⋮	⋮
k = k	${V}_{0}^{(k)}$	${V}_{1}^{(k)}$	${V}_{2}^{(k)}$	⋯	${V}_{k-1}^{k}$	${V}_{k}^{(k)}$

Looking for another recursive relations to produce these methods, let us define a charactristic polynomial for a fractional Adams-Moulton method.

Definition 2.2. The characteristic polynomial of the k-th order, (k-1)-step, fractional Adams-Moulton method (10) is defined by

$$\begin{eqnarray}&&{\sigma }_{k}(\xi )={V}_{0}^{(k-1)}{\xi }^{k-1}+{V}_{1}^{(k-1)}{\xi }^{k-2}+\cdots +{V}_{k-2}^{(k-1)}\xi +{V}_{k-1}^{(k-1)}{\xi }^{0}.\end{eqnarray} \tag{ 13 }$$

Note that, here the subscript k denotes the order not the stepnumber of the method.

Proposition 2.3. The following property holds.

$$\begin{eqnarray}&&{\sigma }_{k+1}(\xi )=\xi {\sigma }_{k}(\xi )+{V}_{k}^{(0)}{\left(\xi -1\right)}^{k}.\end{eqnarray} \tag{ 14 }$$

Proof. We have

$$\begin{eqnarray*}\begin{array}{rcl}\xi {\sigma }_{k}(\xi )+{V}_{k}^{(0)}{\left(\xi -1\right)}^{k} & = & \displaystyle \sum _{j=0}^{k-1}{V}_{j}^{(k-1)}{\xi }^{k-j}+{V}_{k}^{(0)}\left(\displaystyle \sum _{j=0}^{k}{\left(-1\right)}^{j}\left({}_{j}^{k}\right){\xi }^{k-j}\right)\\ & = & \displaystyle \sum _{j=0}^{k-1}\left({V}_{j}^{(k-1)}+{\left(-1\right)}^{j}\left({}_{j}^{k}\right){V}_{k}^{(0)}\right){\xi }^{k-j}+{\left(-1\right)}^{k}{V}_{k}^{(0)}\\ & = & \displaystyle \sum _{j=0}^{k-1}{V}_{j}^{(k)}{\xi }^{k-j}+{V}_{k}^{(k)}\\ & = & {\sigma }_{k+1}(\xi ).\end{array}\end{eqnarray*}$$

□

Therefore, using (13), the fractional Adams-Moulton methods are generated as shown in table 2.

Table 2. Recursive construction of fractional Adams-Moulton methods.

σ1(ξ) (0-step FAM)			1
ξ σ1(ξ)		ξ
${V}_{1}^{(0)}(\xi -1)$		$-\tfrac{\beta }{2}\xi$	$\tfrac{\beta }{2}$
σ2(ξ)		$\left(1-\tfrac{\beta }{2}\right)\xi$	$\tfrac{\beta }{2}$
ξ σ2(ξ)	$\left(1-\tfrac{\beta }{2}\right){\xi }^{2}$	$\tfrac{\beta }{2}\xi$
${V}_{2}^{(0)}{\left(\xi -1\right)}^{2}$	$\left(\tfrac{1}{8}{\beta }^{2}-\tfrac{5}{24}\beta \right){\xi }^{2}$	$\left(-\tfrac{1}{4}{\beta }^{2}+\tfrac{5}{24}\beta \right)\xi$	$\tfrac{1}{8}{\beta }^{2}-\tfrac{5}{24}\beta$
σ3(ξ)	$\left(\tfrac{1}{8}{\beta }^{2}-\tfrac{17}{24}\beta +1\right){\xi }^{2}$	$\left(-\tfrac{1}{4}{\beta }^{2}+\tfrac{11}{12}\beta \right)\xi$	$\tfrac{1}{8}{\beta }^{2}-\tfrac{5}{24}\beta$
⋮

2.2. Fractional Adams-Bashforth array

Let us consider a k-step fractional Adams-Bashforth method of the form

$$\begin{eqnarray}&&\displaystyle \sum _{j=0}^{n}{\omega }_{j}({y}_{n-j}-{y}_{0})+\displaystyle \sum _{j=0}^{s}{w}_{n,j}({y}_{j}-{y}_{0})={h}^{\beta }\displaystyle \sum _{j=1}^{k}{V}_{j}^{* (k)}{f}_{n-j}.\end{eqnarray} \tag{ 15 }$$

Firstly, we prove a recursive relation for the coefficients ${V}_{j}^{* (k)}.$ For this purpose, similar to the fractional Adams-Moulton methods starting from the β-power of the generating function $G(t)=\tfrac{-t}{\mathrm{ln}(1-t)},$ as the first row of the Adams-Bashforth array and using the relation

$$\begin{eqnarray}&&{\delta }^{* }(\xi )={\left(1-\xi \right)}^{-\beta }\left({\gamma }_{0}+{\gamma }_{1}(1-\xi )+\cdots +{\gamma }_{k-1}{\left(1-\xi \right)}^{k-1}+{c}_{k}{\left(1-\xi \right)}^{k}\right),\end{eqnarray} \tag{ 16 }$$

where δ^*(ξ) is the generating function of the coefficients ${V}_{j}^{* (k)}$, and γ_j are the coefficients of Adams-Moulton methods, we can prove the following lemma.

Lemma 2.4. For the coefficients ${V}_{j}^{* (k)}$ appeared in (15), one has the following recursive relation

$$\begin{eqnarray}\left\{\begin{array}{ll}{V}_{l}^{* (k)}={V}_{l}^{* (k-1)}+{\left(-1\right)}^{l}\left({}_{l-1}^{k-1}\right){c}_{k}{V}_{0}^{(0)}, & l\lt k,\\ {V}_{l}^{* (k)}={\left(-1\right)}^{l}\left({}_{l}^{k}\right){c}_{k}{V}_{0}^{(0)}, & l=k,\end{array}\right.\end{eqnarray} \tag{ 17 }$$

where ${c}_{k}=-{\sum }_{j=0}^{k-1}{V}_{j}^{(0)}.$

Therefore, taking (17) into account enables us to easily write an array of the coefficients, ${V}_{j}^{* (k)},$ of a k-step fractional Adams-Bashforth method. We call this array as the fractional Adams-Bashforth array. The array is shown in table 3.

Table 3. Fractional Adams-Bashforth array.

	${V}_{0}^{(0)}$	${V}_{1}^{(0)}$	${V}_{2}^{(0)}$	⋯	${V}_{l-1}^{(0)}$	${V}_{l}^{(0)}$
k = 0	${V}_{1}^{* (1)}$
k = 1	${V}_{1}^{* (2)}$	${V}_{2}^{* (2)}$
k = 2	${V}_{1}^{* (3)}$	${V}_{2}^{* (3)}$	${V}_{3}^{* (3)}$
⋮	⋮	⋮	⋮
k = k	${V}_{1}^{* (k)}$	${V}_{2}^{* (k)}$	${V}_{3}^{(k)}$	⋯	${V}_{k-1}^{* (k)}$	${V}_{k}^{* (k)}$

In order to clarify how to use the given relation (17) between the coefficients of table 3, we consider for example k = 3 in (17), the 3-step fractional Adams-Bashforth method coefficients are as follows:

$$\begin{eqnarray*}&&\displaystyle \sum _{j=0}^{n}{\omega }_{j}({y}_{n-j}-{y}_{0})+\displaystyle \sum _{j=0}^{s}{w}_{n,j}({y}_{j}-{y}_{0})={h}^{\beta }\left({V}_{1}^{* (3)}{f}_{n-1}+{V}_{2}^{* (3)}{f}_{n-2}+{V}_{3}^{* (3)}{f}_{n-3}\right),\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&\begin{array}{l}{V}_{1}^{* (3)}={V}_{1}^{* (2)}-{(}_{0}^{2}){c}_{3}{V}_{0}^{(0)},\\ {V}_{2}^{* (3)}={V}_{2}^{* (2)}+{(}_{1}^{2}){c}_{3}{V}_{0}^{(0)},\\ {V}_{3}^{* (3)}={\left(-1\right)}^{3}{(}_{3}^{3}){c}_{3}{V}_{0}^{(0)}.\end{array}\end{eqnarray*}$$

Now, like the previous subsection, we define characteristic polynomial for a fractional Adams-Bashforth method.

Definition 2.5. For the characteristic polynomial of the k-the order, fractional Adams-Bashforth method (15) is defined by

$$\begin{eqnarray}&&{\sigma }_{k}^{* }(\xi )={V}_{1}^{* (k)}{\xi }^{k-1}+{V}_{2}^{* (k)}{\xi }^{k-2}+\cdots +{V}_{k-1}^{* (k)}\xi +{V}_{k}^{* (k)}{\xi }^{0}.\end{eqnarray} \tag{ 18 }$$

Similar to the fractional Adams-Moulton methods, we can prove the following property.

Proposition 2.6. The following property holds.

$$\begin{eqnarray}&&{\sigma }_{k+1}^{* }(\xi )=\xi {\sigma }_{k}^{* }(\xi )-{c}_{k+1}{\left(\xi -1\right)}^{k}.\end{eqnarray} \tag{ 19 }$$

Therefore, using (18) the fractional Adams-Bashforth methods are generated as shown in table 4.

Table 4. Recursive construction of fractional Adams-Bashforth methods.

${\sigma }_{1}^{* }(\xi )$ (1-step FAB)			1
$\xi {\sigma }_{1}^{* }(\xi )$		ξ
c2(ξ − 1)		$(-1+\tfrac{\beta }{2})\xi$	$1-\tfrac{\beta }{2}$
${\sigma }_{2}^{* }(\xi )$		$(2-\displaystyle \frac{\beta }{2})\xi$	$-1+\tfrac{\beta }{2}$
$\xi {\sigma }_{2}^{* }(\xi )$	$(2-\displaystyle \frac{\beta }{2}){\xi }^{2}$	$(-1+\tfrac{\beta }{2})\xi$
${c}_{3}{\left(\xi -1\right)}^{2}$	$(-1+\displaystyle \frac{17}{24}\beta -\tfrac{1}{8}\beta ){\xi }^{2}$	$(2-\displaystyle \frac{17}{12}\beta +\tfrac{1}{4}{\beta }^{2})\xi$	$-1+\displaystyle \frac{17}{24}\beta -\tfrac{1}{8}{\beta }^{2}$
${\sigma }_{3}^{* }(\xi )$	$(3-\displaystyle \frac{29}{24}\beta +\tfrac{1}{8}{\beta }^{2}){\xi }^{2}$	$(-3+\displaystyle \frac{23}{12}\beta -\tfrac{1}{4}{\beta }^{2})\xi$	$1-\displaystyle \frac{17}{24}\beta +\tfrac{1}{8}{\beta }^{2}$
⋮

Now, we try to illustrate the application of tables 1 and 2 by using these to construct higher-order fractional linear multi-step methods in general form with extended stability regions. For this purpose, for example, we use the provided fractional Adams arrays to construct the corresponding fractional form of a 3-step linear multistep method in general form

$$\begin{eqnarray*}&&\displaystyle \sum _{j=0}^{3}{\alpha }_{j}{y}_{n-j}=h\displaystyle \sum _{j=0}^{3}{\beta }_{j}{f}_{n-j}.\end{eqnarray*}$$

Due to consistency the first characteristic polynomial of the method has a root at +1 and therefore it has the following form

$$\begin{eqnarray*}&&\rho (\xi )=(\xi -1)({\xi }^{2}+{A}_{1}\xi +{A}_{2}),\end{eqnarray*}$$

where A₁ and A₂ are the parameters to be obtained. Therefore, the left-side of the above mentioned method becomes

$$\begin{eqnarray}&&({y}_{n}-{y}_{n-1})+{A}_{1}({y}_{n-1}-{y}_{n-2})+{A}_{2}({y}_{n-2}-{y}_{n-3}).\end{eqnarray} \tag{ 20 }$$

Then, using 3-step Adams-Bashforth and 2-step Adams-Moulton methods for parentheses of (20) yields the following two-parameter family of explicit linear 4-step method

$$\begin{eqnarray}&&\begin{array}{l}({y}_{n}-{y}_{n-1})+{A}_{1}({y}_{n-1}-{y}_{n-2})+{A}_{2}({y}_{n-2}-{y}_{n-3})\\ =\displaystyle \frac{h}{12}(23{f}_{n-1}-16{f}_{n-2}+5{f}_{n-3})\\ \,+\displaystyle \frac{{A}_{1}h}{12}(23{f}_{n-2}-16{f}_{n-3}+5{f}_{n-4})\\ \,+\displaystyle \frac{{A}_{2}h}{12}(5{f}_{n-2}+8{f}_{n-3}-{f}_{n-4}).\end{array}\end{eqnarray} \tag{ 21 }$$

Note that the order of (21) is 3 and the error constant is

$$\begin{eqnarray*}&&{c}_{4}=\displaystyle \frac{3}{8}+\displaystyle \frac{3}{8}{A}_{1}-\displaystyle \frac{1}{24}{A}_{2}.\end{eqnarray*}$$

By choosing A₂ = 9(1 + A₁), the order of the method increase to 4, and the error constant is ${C}_{5}=\tfrac{251}{720}+\tfrac{251}{720}{A}_{1}-\tfrac{19}{720}{A}_{2}.$ Of course, we must choose the parameter A₁ such that the method is zero-stable and therefore $-1\lt {A}_{1}\lt -\tfrac{8}{9}.$ We choose ${A}_{1}=\displaystyle \frac{-17}{18}.$ Now, applying the results of previous section and using table 1 and table 2, we constract the fractional form of (21) as the following

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{j=0}^{n}{\omega }_{j}({y}_{n-j}-{y}_{0})+\displaystyle \sum _{j=0}^{n-1}{\omega }_{j}({y}_{n-1-j}-{y}_{0})+\displaystyle \sum _{j=0}^{n-2}{\omega }_{j}({y}_{n-2-j}-{y}_{0})\\ \quad ={h}^{\beta }\left({V}_{1}^{* (3)}{f}_{n-1}+({V}_{2}^{* (3)}+{A}_{1}{V}_{1}^{* (3)}+{A}_{2}{V}_{0}^{(2)}){f}_{n-2}\right.\\ \quad \left.+({V}_{3}^{* (3)}+{A}_{1}{V}_{2}^{* (3)}+{A}_{2}{V}_{1}^{(2)}){f}_{n-3}+({A}_{1}{V}_{3}^{* (3)}+{A}_{2}{V}_{2}^{(2)}){f}_{n-4}\right).\end{array}\end{eqnarray} \tag{ 22 }$$

Note that we used the 3-step fractional Adams-Bashforth method and the 2-step fractional Adams-Moulton method, which are consistent and convergent with order 3. Due to some lemmas in [16], we have the following result for our new method. Now, taking theorem 1.2 into account, we can easily prove the following theorem.

Theorem 3.1. The method (20) is consistent and convergent with order 4.

Let us consider the following test problem

$$\begin{eqnarray*}&&{}_{{t}_{0}}{D}^{\beta }y(t)=\lambda y(t),\,\,\,y({t}_{0})={y}_{0},\,\,\,\lambda \in {\mathbb{C}},\,\,\,0\lt \beta \lt 1,\end{eqnarray*}$$

where the exact solution takes the form $y(t)={E}_{\beta }(\lambda {\left(t-{t}_{0}\right)}^{\beta }){y}_{0}$ such that ${E}_{\beta }(x)={\sum }_{k=0}^{\infty }\tfrac{{x}^{k}}{{\rm{\Gamma }}(\beta k+1)}$ denotes the Mittag-Leffler function. It follows form the properties of E_β (x) [20] that the solution of (20) satisfy y(t) ⟶ 0, as t ⟶ ∞ , If $\lambda \in {S}^{* }=\{z\in {\mathbb{C}}:| {\arg }(z)-\pi | \lt (1-\tfrac{\beta }{2})\pi \}.$ We now try to deduce conditions in order that the numerical solutions of (20), obtained by our new method, has the same property. In this regard, we define the stability domain of a fractional numerical method for the set of all z = h^β λ for which the numerical solution vanishes as n ⟶ ∞ . It is easy to prove [14, 15] that the stability domain of fractional Adams-Moulton and Adams-Bashforth method can be written by

$$\begin{eqnarray}&&S={\mathbb{C}}\setminus \left\{\displaystyle \frac{1}{\delta (\xi )};| \xi | \leqslant 1\right\},\end{eqnarray} \tag{ 23 }$$

and since our new method is a combination of these methods the same argument is true.

It is obvious that the stability domain does not depend on the starting weights and it can be obtained only by the generating function. The generating function of (20) is

$$\begin{eqnarray}&&\delta (\xi )={\left(1-\xi \right)}^{-\beta }{\left(1+{A}_{1}\xi +{A}_{2}{\xi }^{2}\right)}^{-1}({c}_{1}\xi +{c}_{2}{\xi }^{2}+{c}_{3}{\xi }^{3}+{c}_{4}{\xi }^{4}),\end{eqnarray} \tag{ 24 }$$

such that

$$\begin{eqnarray*}\begin{array}{rcl}{c}_{1} & = & {V}_{1}^{* (3)},\\ {c}_{2} & = & {V}_{2}^{* (3)}+{A}_{1}{V}_{1}^{* (3)}+{A}_{2}{V}_{0}^{(2)},\\ {c}_{3} & = & {V}_{3}^{* (3)}+{A}_{1}{V}_{2}^{* (3)}+{A}_{2}{V}_{1}^{(2)},\\ {c}_{4} & = & {A}_{1}{V}_{3}^{* (3)}+{A}_{2}{V}_{2}^{(2)},\end{array}\end{eqnarray*}$$

where A₁ and A₂ are already defined. With the aim of (23), we can plot stability regions of the method (22) in figure 1. At the end of this section, we compare stability regions of our new method, FAB and FAM methods.

Zoom In Zoom Out Reset image size

For more in-depth description of stability charactristics of the method, let us first note the following theorems from [21].

Theorem 4.1. A convolution quadrature is strongly stable if and only if there exists $c\gt 0$ such that

$$\begin{eqnarray*}&&{\mathbb{R}}e({g}_{\beta }(\xi ))\geqslant -c,\end{eqnarray*}$$

and the open disk contained in the stability domain has center $-(\tfrac{1}{2c})$ and radius $R=(\tfrac{1}{2c}).$

Theorem 4.2. If $| \xi | \leqslant 1,$ then

$$\begin{eqnarray*}&&{\mathbb{R}}e(\delta (\xi ))\geqslant \delta (-1),\end{eqnarray*}$$

where δ is the generating function of (24).

Proof. Using (24), we know that the generating function $\delta (\xi )$ is of the following form

$$\begin{eqnarray}\begin{array}{rcl}\delta (\xi ) & = & {\delta }_{1}(\xi )+{A}_{1}{\delta }_{1}(\xi )+{A}_{2}{\delta }_{2}(\xi )\\ & = & (1+{A}_{1}){\delta }_{1}(\xi )+{A}_{2}{\delta }_{2}(\xi ),\end{array}\end{eqnarray} \tag{ 25 }$$

where ${\delta }_{1}(\xi )$ and ${\delta }_{2}(\xi )$ are generating functions of 3-step FAB and 2-step FAM method, respectively. Taking theorem 4.1 in to account, we have

$$\begin{eqnarray*}&&\begin{array}{l}(1+{A}_{1}){\delta }_{1}(\xi )\geqslant (1+{A}_{1}){\delta }_{1}(-1),\\ {A}_{2}{\delta }_{2}(\xi )\geqslant {A}_{2}{\delta }_{2}(-1).\end{array}\end{eqnarray*}$$

In other words

$$\begin{eqnarray*}&&\begin{array}{l}\delta (\xi )\geqslant (1+{A}_{1}){\delta }_{1}(-1)+{A}_{2}{\delta }_{2}(-1)\\ \quad =\delta (-1).\end{array}\end{eqnarray*}$$

□

In order to compare the method (22) with classical fractional Adams methods, we consider interval of stability which is a subset of real axis such that, for $\lambda \in {\mathbb{R}},$ choosing h^β λ from this interval makes the numerical solution tends to zero as n ⟶ ∞ . Due to the above theorem, we have (B_k (β), 0) ⫅ I_k (α), where ${B}_{k}(\beta )=\tfrac{1}{{\delta }_{k}(-1)}$ and I_k (β) denotes the interval of stability, k is the order of the method. In the following, we show that the stability interval of the method (22), I(β), is at least 2.23 times and at most 2.79 times larger than the stability interval of the 4-step (of order 4) fractional Adams-Bashforth method, for 0 < β < 1. Indeed since

$$\begin{eqnarray*}\begin{array}{rcl}\delta (\xi ) & = & {\left(1-\xi \right)}^{-\beta }{\left(1+{A}_{1}\xi +{A}_{2}{\xi }^{2}\right)}^{-1}({c}_{1}\xi +{c}_{2}{\xi }^{2}+{c}_{3}{\xi }^{3}+{c}_{4}{\xi }^{4})\\ & = & {\left(1-\xi \right)}^{-\beta }{\left(1+{A}_{1}\xi +{A}_{2}{\xi }^{2}\right)}^{-1}{M}_{\beta }(\xi ),\end{array}\end{eqnarray*}$$

and

$$\begin{eqnarray*}\begin{array}{rcl}{\delta }_{4}^{* }(\xi ) & = & {\left(1-\xi \right)}^{-\beta }({\gamma }_{1}^{* (4)}\xi +{\gamma }_{2}^{* (4)}{\xi }^{2}+{\gamma }_{3}^{* (4)}{\xi }^{3}+{\gamma }_{4}^{* (4)}{\xi }^{4})\\ & = & {\left(1-\xi \right)}^{-\beta }{N}_{\beta }(\xi ).\end{array}\end{eqnarray*}$$

Therefore, it is easy to prove the following relation

$$\begin{eqnarray*}&&{{LB}}_{4}(\beta )\lt B(\beta )\lt {{MB}}_{4}(\beta ),\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}L & = & -{\gamma }_{1}^{* (4)}+{\gamma }_{2}^{* (4)}-{\gamma }_{3}^{* (4)}+{\gamma }_{4}^{* (4)}\\ & = & -15+\displaystyle \frac{21}{2}\beta -\displaystyle \frac{7}{3}{\beta }^{2}+\displaystyle \frac{1}{6}{\beta }^{3},\end{array}\end{eqnarray*}$$

and

$$\begin{eqnarray*}\begin{array}{rcl}M & = & -{c}_{1}+{c}_{2}-{c}_{3}+{c}_{4}\\ & = & -\displaystyle \frac{118}{9}+\displaystyle \frac{353}{54}\beta -\displaystyle \frac{13}{18}{\beta }^{2},\end{array}\end{eqnarray*}$$

and B(β), B₄(β) denote stability boundaries of the method (22) and the 4-step fractional Adams-Bashforth method, respectively. Considering the ratio $\tfrac{L}{M}$ for different values of β, it follows that 2.23B₄(β) < B(β) < 2.79B₄(β). A comparison of the figures 1 and 2 confirms this.

Zoom In Zoom Out Reset image size

In this section, to show the high performance of the presented method (22), we use this method to solve a fractional problem and compare the results with the results of previous methods, especially with fractional classical Adams methods.

Example 5.1. Let us consider the following problem which has been considered in several papers [22–24]

$$\begin{eqnarray}&&{D}^{\beta }y(t)=\displaystyle \frac{2{t}^{2-\beta }}{{\rm{\Gamma }}(3-\beta )}+{t}^{2}-y(t),\end{eqnarray} \tag{ 26 }$$

with $y(0)=0,$ where the exact solution is $y(t)={t}^{2}.$

In figure 3, we provide resulting errors, E(N) = y(t_N ) − y_N , for stipulated methods. These results emphasize the convergence advantages of our method.

Zoom In Zoom Out Reset image size

In this article, we have introduced two arrays, each containing the coefficients of the fractional Adams-Bashforth and Adams-Moulton methods, and then used recursive relations to produce the elements of these arrays. Also, formulas for obtaining the generating function of k + 1-step method from k-step method were proved for both fractional Adams-Moulton and Adams-Bashforth methods. We have derived the application of proposed arrays to construct higher-order methods with extended stability regions for solving fractional differential equations. Apparent parameters in newly constructed methods enabled us to increase the order and improve the stability. An in-depth comparison between the order and stability regions of fractional Adams-type methods and our proposed methods was performed. An illustrative example has been given to justify the effectiveness of the proposed method compared with some previously published results.

This research is supported by the research grant of the University of Tabriz(number 210).

All data that support the findings of this study are included within the article (and any supplementary files).

This work does not have any conflicts of interest.