# A spectrally accurate step-by-step method for the numerical solution of
fractional differential equations

Luigi Brugnano Dipartimento di Matematica e Informatica “U. Dini”, Università
di Firenze, Italy<EMAIL_ADDRESS>Kevin Burrage School of
Mathematical Sciences, Queensland University of Technology, Australia,
<EMAIL_ADDRESS>Pamela Burrage School of Mathematical Sciences,
Queensland University of Technology, Australia<EMAIL_ADDRESS>Felice Iavernaro Dipartimento di Matematica, Università di Bari, Italy,
<EMAIL_ADDRESS>
###### Abstract

In this paper we consider the numerical solution of fractional differential
equations. In particular, we study a step-by-step graded mesh procedure based
on an expansion of the vector field using orthonormal Jacobi polynomials.
Under mild hypotheses, the proposed procedure is capable of getting spectral
accuracy. A few numerical examples are reported to confirm the theoretical
findings.

Keywords: fractional differential equations, fractional integrals, Jacobi
polynomials, Hamiltonian Boundary Value Methods, HBVMs.

AMS-MSC: 65L05, 65L03, 65L99.

##  1 Introduction

In recent years fractional differential equation modelling has become more and
more frequent - see, for example, the two monographs [37] and [19]. In fact
the use of fractional models has quite a long history and early applications
considered their behaviour in describing anomalous diffusion in a randomly
hetergeneous porous medium [31] and in water resources modelling [6]. In
mathematical biology, a number of researchers have considered fractional
reaction diffusion equations - see, for example, [16] who have studied via a
spectral approach the interfacial properties of the Allen-Cahn equation with a
quartic double well potential, a model often used for alloy kinetics and the
dynamics of phase boundaries. Henry, Langlands and Wearne [27] analysed Turing
pattern formation in fractional activator-inhibitor systems. In a
computational medicine setting, Henry and Langlands [26] developed fractional
cable models for spiny neuronal dendrites. Orsingher and Behgin [36] have
considered the dynamics of fractional chemical kinetics for unimolecular
systems using time change. [17] have made a series of arguments based on Riesz
potential theory and the wide range of heterogeneous cardiac tissues for the
use of fractional models in cardiac electrophysiology, while [18] have
explored the effects of fractional diffusion on the cardiac action potential
shape and electrical propagation. In addition, [34], [28] and [15] have
considered the calibration of fractional cardiac models through a fractional
Bloch-Torrey equation.

Due to these application advances there is clearly a need to develop advanced
numerical methods for fractional differential equations. In a purely time
fractional setting the standard approach has been to use a first order
approximation to the Caputo derivative. Lubich [32] has extended this and
developed higher order approximations based on an underlying linear multistep
method that are convergent of the order of this multistep method.

There are two important issues when designing effective numerical methods for
Fractional Differential equations: memory and non-locality. In the case of
memory, in principle one needs to evaluate the numerical approximations all
the way back to the initial time point at each time step and this becomes very
computationally intensive over long time intervals. Techniques such as fixed
time windowing have been proposed [37] and more recently a number of
approaches have been considered to compute the discrete convolutions to
approximate the fractional operator [39], [43]. Zeng et al. [42] developed
fast memory-saving time stepping methods in which the fractional operator is
decoupled into a local part with fixed memory length and the history part is
computed by Laguerre-Gauss quadrature, whose error is independent of the
stepsize.

In the second situation the solutions to FDEs are often non-smooth and may
have a strong singularity at $t=0$. In these cases the vector field may
inherit this singularity behaviour and so a constant timestep would be very
inefficient. Thus graded meshes have been proposed by a number of authors [30,
44, 40]. Lubich [33] on the other hand deals with singularities by introducing
certain correction terms. Other approaches have been considered such as Adams
methods [21], trapezoidal methods [23] and Bernstein splines [38] while
Garrappa [24] has given a survey and software tutorial on numerical methods
for FDEs.

In this paper, we consider a major improvement of the recent solution approach
described in [2, 10], based on previous work on Hamiltonian Boundary value
Methods (HBVMs) [12, 8, 9, 7] (also used as spectral methods in time [14, 13,
3, 5]), for solving fractional initial value problems in the form:

$y^{(\alpha)}(t)=f(y(t)),\qquad t\in[0,T],\qquad y(0)=y_{0}\in\mathbb{R}^{m},$
(1)

where, without loss of generality, we omit $t$ as a formal argument, at the
r.h.s. Further, for sake of brevity, the dimension $m$ of the state vector
will be also omitted, when unnecessary.

Here, for $\alpha\in(0,1]$, $y^{(\alpha)}(t)\equiv D^{\alpha}y(t)$ is the
Caputo fractional derivative:

$D^{\alpha}g(t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-x)^{-\alpha}\left[\frac{\mathrm{d}}{\mathrm{d}x}g(x)\right]\mathrm{d}x.$
(2)

The Riemann-Liouville integral associated to (2) is given by:

$I^{\alpha}g(t)=\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-x)^{\alpha-1}g(x)\mathrm{d}x.$
(3)

It is known that [37]:

$D^{\alpha}I^{\alpha}g(t)=g(t),~{}\quad
I^{\alpha}D^{\alpha}g(t)=g(t)-g(0),~{}\quad
I^{\alpha}t^{j}=\frac{j!}{\Gamma(\alpha+j+1)}t^{j+\alpha},\quad
j=0,1,2,\dots.$ (4)

Consequently, the solution of (1) is formally given by:

$y(t)=y_{0}+I^{\alpha}f(y(t)),\qquad t\in[0,T].$ (5)

We shall generalize the approach given in [2], by defining a step-by-step
procedure, whereas the procedure given in [2] was of a global nature. This
strategy, when combined with a graded mesh selection, enables better handling
of singularities in the derivative of the solution to (1) at the origin, which
occur, for example, when $f$ is a suitably regular function. Additionally, our
approach partially mitigates issues stemming from the non-local nature of the
problem. Indeed, at a given time-step of the integration procedure, it
requires the solution of a differential problem within the current time
interval along with a finite sequence of pure quadrature problems that account
for the contribution brought by the history term.

With this premise, the structure of the paper is as follows: in Section 2 we
recall some basic facts about Jacobi polynomials, for later use; in Section 3
we describe a piecewise quasi-polynomial approximation to the solution of (1);
in Section 4 the analysis of the new method is carried out; in Section 5 we
report a few numerical tests for assessing the theoretical achievements; at
last, in Section 6 a few concluding remarks are given.

##  2 Orthonormal Jacobi polynomials

Jacobi polynomials form an orthogonal set of polynomials w.r.t. a given
weighting function. In more detail, for $r,\nu>-1$:

$\displaystyle\bar{P}^{(r,\nu)}_{i}(x)\in\Pi_{i},\qquad\left(P_{0}^{(r,\nu)}(x)\equiv
1\right)$
$\displaystyle\int_{-1}^{1}(1-x)^{r}(1+x)^{\nu}\bar{P}_{i}^{(r,\nu)}(x)\bar{P}_{j}^{(r,\nu)}(x)\mathrm{d}x~{}=~{}\frac{2^{r+\nu+1}}{2i+r+\nu+1}\frac{\Gamma(i+r+1)\Gamma(i+\nu+1)}{\Gamma(i+r+\nu+1)i!}\delta_{ij},$
$\displaystyle\qquad i,j,=0,1,\dots,$

where, as is usual, $\Pi_{i}$ is the vector space of polynomials of degree at
most $i$ and $\delta_{ij}$ is the Kronecker symbol. Consequently, the shifted
and scaled Jacobi polynomials

$P_{i}^{(r,\nu)}(c)~{}:=~{}\sqrt{(2i+r+\nu+1)\frac{\Gamma(i+r+\nu+1)i!}{\Gamma(i+r+1)\Gamma(i+\nu+1)}}\bar{P}_{i}^{(r,\nu)}(2c-1),\qquad
i=0,1,\dots,$

are orthonormal w.r.t.

$\int_{0}^{1}(1-c)^{r}c^{\nu}P_{i}^{(r,\nu)}(c)P_{j}^{(r,\nu)}(c)\mathrm{d}c~{}=~{}\delta_{ij},\qquad
i,j=0,1,\dots.$

In particular, hereafter we shall consider the polynomials 111Hereafter, for
sake of brevity, we shall omit the upper indices for such polynomials.

$P_{i}(c)~{}:=~{}\frac{P_{i}^{(\alpha-1,0)}(c)}{\sqrt{\alpha}}~{}\equiv~{}\frac{\sqrt{2i+\alpha}}{\sqrt{\alpha}}\bar{P}_{i}^{(\alpha-1,0)}(2c-1),\qquad
i=0,1,\dots,$ (6)

with $\alpha\in(0,1]$ the same parameter in (1), such that:

$P_{0}(c)~{}\equiv~{}1,\qquad\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{i}(c)P_{j}(c)\mathrm{d}c~{}=~{}\delta_{ij},\qquad
i,j=0,1,\dots,$ (7)

where we have slightly changed the weighting function, in order that it has a
unit integral:

$\alpha\int_{0}^{1}(1-c)^{\alpha-1}\mathrm{d}c=1.$ (8)

We refer to [25] and the accompanying software, for their computation, and for
computing the nodes and weights of the Gauss-Jacobi quadrature of order $2k$.

###### Remark 1

As is clear, when in (6)–(8) $\alpha=1$, we obtain the scaled and shifted
Legendre polynomials, orthonormal in $[0,1]$.

As is well known, the polynomials (6) form an orthonormal basis for the
Hilbert space $L_{2}[0,1]$, equipped with the scalar product

$(f,g)=\alpha\int_{0}^{1}(1-c)^{\alpha-1}f(c)g(c)\mathrm{d}c,$ (9)

and the associated norm

$\|f\|=\sqrt{(f,f)}.$ (10)

Consequently, from the Parseval identity, it follows that any square summable
function can be expanded along such a basis, and the corresponding expansion
is convergent to the function itself.

That said, with reference to the vector field in (1), let us consider the
interval $[0,h]$, and the expansion

$f(y(ch))~{}=~{}\sum_{j\geq 0}\gamma_{j}(y)P_{j}(c),\qquad c\in[0,1],$ (11)

with the Fourier coefficients given by:

$\gamma_{j}(y)~{}=~{}\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)f(y(ch))\mathrm{d}c,\qquad
j=0,1,\dots.$ (12)

Consequently, on the interval $[0,h]$, (1) can be rewritten as

$y^{(\alpha)}(ch)~{}=~{}\sum_{j\geq 0}\gamma_{j}(y)P_{j}(c),\qquad
c\in[0,1],\qquad y(0)=y_{0},$ (13)

and, by virtue of (3)-(4), one obtains:

$y(ch)~{}=~{}y_{0}+h^{\alpha}\sum_{j\geq
0}\gamma_{j}(y)I^{\alpha}P_{j}(c),\qquad c\in[0,1].$ (14)

In particular, due to (7) and (8), one has:

$y(h)~{}=~{}y_{0}+\frac{h^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}(y).$ (15)

###### Remark 2

By considering that $P_{0}(c)\equiv 1$, from (12) one derives that

$\gamma_{0}(y)=\alpha\int_{0}^{1}(1-c)^{\alpha-1}f(y(ch))\mathrm{d}c.$

Consequently, because of (3), (15) becomes

$y(h)=y_{0}+I^{\alpha}f(y(h)),$

i.e., (5) with $t=h$. This clearly explains the use of the Jacobi basis for
the expansion (11).

Further, we observe that, when $\alpha=1$, one retrieves the approach
described in [12] for ODE-IVPs.

We also need the following preliminary results.

###### Lemma 1

Let $G:[0,h]\rightarrow V$, with $V$ a vector space, admit a Taylor expansion
at $t=0$. Then,

$\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)G(ch)\mathrm{d}c~{}=~{}O(h^{j}),\qquad
j=0,1,\dots.$

Proof By the hypothesis and (7), one has:

$\displaystyle\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)G(ch)\mathrm{d}c~{}=~{}\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)\sum_{\ell\geq
0}\frac{G^{(\ell)}(0)}{\ell!}(ch)^{\ell}\mathrm{d}c$ $\displaystyle=$
$\displaystyle\alpha\sum_{\ell\geq
0}\frac{G^{(\ell)}(0)}{\ell!}h^{\ell}\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)c^{\ell}\mathrm{d}c~{}=~{}\alpha\sum_{\ell\geq
j}\frac{G^{(\ell)}(0)}{\ell!}h^{\ell}\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)c^{\ell}\mathrm{d}c$
$\displaystyle=$ $\displaystyle O(h^{j}).\,\mbox{~{}$\Box{~{}}$}$

###### Remark 3

Concerning the above result, we observe that, by considering (9)-(10), from
the Cauchy-Schwartz theorem one derives that

$\left|\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)c^{\ell}\mathrm{d}c\right|=|(P_{j},c^{\ell})|\leq\|P_{j}\|\cdot\|c^{\ell}\|\leq
1,\qquad\forall\ell\geq j\geq 0.$

From Lemma 1, one has the following:

###### Corollary 1

Assume that the r.h.s. of problem (1) admits a Taylor expansion at $t=0$, then
the coefficients defined at (12) satisfy:

$\gamma_{j}(y)=O(h^{j}).$

The result of the previous lemma can be weakened as follows (the proof is
similar and, therefore, omitted).

###### Lemma 2

Let $G:[0,h]\rightarrow V$, with $V$ a vector space, admit a Taylor polynomial
expansion of degree $k$ with remainder at $t=0$. Then,

$\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)G(ch)\mathrm{d}c~{}=~{}O(h^{\min(j,k)}),\qquad
j=0,1,\dots.$

In such a case, the result of Corollary 1 is weakened accordingly.

###### Corollary 2

Assume that the r.h.s. of problem (1) admits a Taylor polynomial expansion of
degree $k$ with remainder at $t=0$. Then, the coefficients defined at (12)
satisfy:

$\gamma_{j}(y)=O(h^{\min(j,k)}).$

However, in general, at $t=0$ the solution may be not regular, since the
derivative may be singular (this is, indeed, the case, when $f$ is a regular
function). In such a case, we shall resort to the following result.

###### Lemma 3

Assume that the r.h.s. in (1) is continuous for all $t\in[0,h]$. Then, for a
convenient $\xi_{t}\in(0,t)$, one has:

$y(t)=y_{0}+\frac{y^{(\alpha)}(\xi_{t})}{\Gamma(\alpha+1)}t^{\alpha},\qquad
t\in[0,h].$

Proof In fact, from (4), and by using the weighted mean-value theorem for
integrals, one has:

$\displaystyle y(t)$ $\displaystyle=$ $\displaystyle
y_{0}+I^{\alpha}y^{(\alpha)}(t)~{}=~{}y_{0}+\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-x)^{\alpha-1}y^{(\alpha)}(x)\mathrm{d}x$
$\displaystyle=$ $\displaystyle
y_{0}+\frac{y^{(\alpha)}(\xi_{t})}{\Gamma(\alpha)}\int_{0}^{t}(t-x)^{\alpha-1}\mathrm{d}x~{}=~{}y_{0}+\frac{y^{(\alpha)}(\xi_{t})}{\Gamma(\alpha+1)}t^{\alpha}.\,\mbox{~{}$\Box{~{}}$}$

As a consequence, we derive the following weaker results concerning the
Fourier coefficients (12).

###### Corollary 3

In the hypotheses of Lemma 3, and assuming $f$ is continuously differentiable
in a neighborhood of $y_{0}$, one has:

$\gamma_{0}(y)=f(y_{0})+O(h^{\alpha}),\qquad\gamma_{j}(y)=O(h^{\alpha}),\quad
j=1,2,\dots.$

Proof In fact, from Lemma 3, one has:

$f(y(ch))=f\left(y_{0}+\frac{y^{(\alpha)}(\xi_{ch})}{\Gamma(\alpha+1)}(ch)^{\alpha}\right)=f(y_{0})+f^{\prime}(y_{0})\frac{y^{(\alpha)}(\xi_{ch})}{\Gamma(\alpha+1)}(ch)^{\alpha}+o((ch)^{\alpha}).$

From this relation, and from (12), the statement then follows.  $\Box{~{}}$

For later use, we also need to discuss the quadrature error in approximating
the first $s$ Fourier coefficients (12) by means of the Gauss-Jacobi formula
of order $2k$, for a convenient $k\geq s$, whose abscissae are the zeros of
$P_{k}$, as defined in (7), and with weights:

$P_{k}(c_{i})=0,\qquad
b_{i}=\alpha\int_{0}^{1}(1-c)^{\alpha-1}\ell_{i}(c)\mathrm{d}c,\qquad\ell_{i}(c)=\prod_{j=1,\,j\neq
i}^{k}\frac{c-c_{j}}{c_{i}-c_{j}},\qquad i=1,\dots,k.$ (16)

That is,

$\gamma_{j}(y)~{}\approx~{}\sum_{i=1}^{k}b_{i}P_{j}(c_{j})f(y(c_{i}h))~{}=:~{}\hat{\gamma}_{j}(y).$
(17)

Concerning the quadrature errors

$\Sigma_{j}^{\alpha}(y,h)~{}:=~{}\hat{\gamma}_{j}(y)-\gamma_{j}(y),\qquad
j=0,\dots,s-1,$ (18)

the following result holds true.

###### Theorem 1

Assume the function $G_{j,h}(c):=P_{j}(c)f(y(ch))$ be of class
$C^{(2k)}[0,h]$. Then, with reference to (18), one has, for a suitable
$\xi=(0,1):$

$\Sigma_{j}^{\alpha}(y,h)~{}=~{}K_{j}\frac{G_{j,h}^{(2k)}(\xi)}{(2k)!}~{}=~{}O(h^{2k-j}),\qquad
j=0,\dots,s-1,$ (19)

with the constants $K_{j}$ independent of both $G_{j,h}$ and $h$.

Proof The statement follows by considering that the quadrature is exact for
polynomials in $\Pi_{2k-1}$.  $\Box{~{}}$

However, as recalled above, when $f$ is a smooth function, the derivative of
the solution of problem (1) may have a singularity at $t=0$. In such a case,
estimates can be also derived (see, e.g., [20, 41] and [35, Theorem 5.1.8]).
However, for our purposes, because of (7), for $j=0,\dots,s-1$, from Corollary
3 we can easily derive the following one, assuming that $k\geq s$ and (13)
hold true:

$\displaystyle\hat{\gamma}_{j}(y)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f(y(c_{i}h))~{}=~{}\sum_{i=1}^{k}b_{i}P_{j}(c_{i})\sum_{\ell\geq
0}P_{\ell}(c_{i})\gamma_{\ell}(y)$ $\displaystyle=$
$\displaystyle\sum_{\ell=0}^{s-1}\underbrace{\left(\sum_{i=1}^{k}b_{i}P_{j}(c_{i})P_{\ell}(c_{i})\right)}_{=\,\delta_{j\ell}}\gamma_{\ell}(y)~{}+~{}O(h^{\alpha})~{}\equiv~{}\gamma_{j}(y)+\Sigma_{j}(y,h),\qquad
j=0,\dots,s-1.$

Consequently,

$\Sigma_{j}(y,h)~{}=~{}O(h^{\alpha}),\qquad j=0,\dots,s-1.$ (20)

In order to derive a polynomial approximation to (13), it is enough to
truncate the infinite series in (11) to a finite sum, thus obtaining:

$\sigma^{(\alpha)}(ch)~{}=~{}\sum_{j=0}^{s-1}\gamma_{j}(\sigma)P_{j}(c),\qquad
c\in[0,1],\qquad\sigma(0)=y_{0},$ (21)

with $\gamma_{j}(\sigma)$ formally given by (12) upon replacing $y$ by
$\sigma$. In so doing, (14) and (15) respectively become:

$\sigma(ch)~{}=~{}y_{0}+h^{\alpha}\sum_{j=0}^{s-1}\gamma_{j}(\sigma)I^{\alpha}P_{j}(c),\qquad
c\in[0,1],$ (22)

and

$\sigma(h)~{}=~{}y_{0}+\frac{h^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}(\sigma).$
(23)

It can be shown that Corollary 1, Corollary 2, Lemma 3, and Corollary 3
continue formally to hold for $\sigma$. Further, by considering the
approximation of the Fourier coefficients obtained by using the Gauss-Jacobi
quadrature of order $2k$,

$\hat{\gamma}_{j}(\sigma)~{}=~{}\sum_{i=1}^{k}b_{i}P_{j}(c_{j})f(\sigma(c_{i}h))~{}\equiv~{}\gamma_{j}(\sigma)+\Sigma_{j}^{\alpha}(\sigma,h),$
(24)

the result of Theorem 1 continues to hold. Moreover, we shall assume that the
expansion

$f(\sigma(ch))=\sum_{j\geq 0}P_{j}(c)\gamma_{j}(\sigma),\qquad c\in[0,1],$
(25)

holds true, similarly as (11), from which also (20) follows for the quadrature
error $\Sigma_{j}(\sigma,h)$, when $\sigma$ has a singular derivative at $0$.
In the next sections, we shall better detail, generalize, and analyze this
approach.

##  3 Piecewise quasi-polynomial approximation

To begin with, in order to obtain a piecewise quasi-polynomial approximation
to the solution of (1), we consider a partition of the integration interval in
the form:

$t_{n}=t_{n-1}+h_{n},\qquad n=1,\dots,N,$ (26)

where

$t_{0}=0,\qquad\sum_{n=1}^{N}h_{n}=T.$ (27)

Then, according to (21)–(23), on the first subinterval $[0,h_{1}]$ we can
derive a polynomial approximation of degree $s-1$ to (11), thus getting the
fractional initial value problem

$\sigma_{1}^{(\alpha)}(ch_{1})=\sum_{j=0}^{s-1}\gamma_{j}(\sigma_{1})P_{j}(c),\qquad
c\in[0,1],\qquad\sigma_{1}(0)=y_{0},$ (28)

where $\gamma_{j}(\sigma_{1})$ is formally given by (12), upon replacing $y$
by $\sigma_{1}$ at the right-hand side:

$\gamma_{j}(\sigma_{1})=\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)f(\sigma_{1}(ch_{1}))\mathrm{d}c,\qquad
j=0,\dots,s-1.$ (29)

The solution of (28) is a quasi-polynomial of degree $s-1+\alpha$, formally
given by:

$\sigma_{1}(ch_{1})=y_{0}+h_{1}^{\alpha}\sum_{j=0}^{s-1}\gamma_{j}(\sigma_{1})I^{\alpha}P_{j}(c),\qquad
c\in[0,1].$ (30)

However, in order to actually compute the Fourier coefficients
$\\{\gamma_{j}(\sigma_{1})\\}$, one needs to approximate them by using the
Gauss-Jacobi quadrature (16):222Hereafter, as a notational convention,
$\gamma_{j}^{i}:=\hat{\gamma}_{j}(\sigma_{i})$, $i=1,\dots,N$.

$\gamma_{j}^{1}~{}:=~{}\hat{\gamma}_{j}(\sigma_{1})~{}=~{}\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f(\sigma_{1}(c_{i}h_{1}))~{}\equiv~{}\gamma_{j}(\sigma_{1})+\Sigma_{j}^{\alpha}(\sigma_{1},h_{1}).$

In so doing, (30) now formally becomes

$\sigma_{1}(ch_{1})=y_{0}+h_{1}^{\alpha}\sum_{j=0}^{s-1}\gamma_{j}^{1}I^{\alpha}P_{j}(c),\qquad
c\in[0,1].$ (31)

Moreover, one solves the system of equations, having (block) dimension $s$
independently of $k$:

$\gamma_{j}^{1}~{}=~{}\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f\left(y_{0}+h^{\alpha}\sum_{\nu=0}^{s-1}\gamma_{\nu}^{1}\,I^{\alpha}P_{\nu}(c_{i})\right),\qquad
j=0,\dots,s-1.$ (32)

This kind of procedure is typical of HBVM$(k,s)$ methods, in the case of ODE-
IVPs [8]: the main difference, here, derives from the non-locality of the
operator. As is clear (compare with (23)), the approximation at $t=h_{1}$ will
be now given by: 333Hereafter, we shall in general denote by $\bar{y}_{n}$ the
approximation to $y(t_{n})$, since $y_{n}(t)$ will be later used to denote the
restriction of $y(t)$ to the subinterval $[t_{n-1},t_{n}]$, $n=1,\dots,N$.

$\bar{y}_{1}~{}:=~{}\sigma_{1}(h_{1})~{}=~{}y_{0}+\frac{h_{1}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}^{1}.$
(33)

###### Remark 4

For an efficient and stable evaluation of the integrals

$I^{\alpha}P_{0}(c_{i}),~{}I^{\alpha}P_{1}(c_{i}),~{}\dots,~{}I^{\alpha}P_{s-1}(c_{i}),\qquad
i=1,\dots,k,$ (34)

we refer to the procedure described in [4].

###  3.1 The general procedure

We now generalize the previous procedure for the subsequent integration steps.
For later use, we shall denote:

$y_{n}(ch_{n})~{}:=~{}y(t_{n-1}+ch_{n}),\qquad c\in[0,1],\qquad n=1,\dots,N,$
(35)

the restriction of the solution on the interval $[t_{n-1},t_{n}]$. Similarly,
we shall denote by $\sigma(t)\approx y(t)$ the piecewise quasi-polynomial
approximation such that:

$\sigma|_{[t_{n-1},t_{n}]}(t_{n-1}+ch_{n})~{}=:~{}\sigma_{n}(ch_{n})~{}\approx~{}y_{n}(ch_{n}),\qquad
c\in[0,1],\qquad n=1,\dots,N.$ (36)

Then, by using an induction argument, let us suppose one knows the quasi-
polynomial approximations

$\sigma_{i}(ch_{i})=\phi^{\alpha}_{i-1}(c,\sigma)+h_{i}^{\alpha}\sum_{j=0}^{s-1}\gamma_{j}^{i}I^{\alpha}P_{j}(c)~{}\approx~{}y_{i}(ch_{i}),\qquad
c\in[0,1],\qquad i=1,\dots,n,$ (37)

where $\phi^{\alpha}_{i-1}(c,\sigma)$ denotes a history term, to be defined
later, such that:

  * •

in the first subinterval, according to (31),
$\sigma_{0}^{\alpha}(c,\sigma)\equiv y_{0}$, $c\in[0,1]$;

  * •

for $i>1$, $\phi_{i-1}^{\alpha}(c,\sigma)$ only depends on
$\sigma_{1},\dots,\sigma_{i-1}$.

The corresponding approximations at the grid-points are given by

$\bar{y}_{i}~{}:=~{}\sigma_{i}(h_{i})~{}=~{}\phi^{\alpha}_{i-1}(1,\sigma)~{}+~{}\frac{h_{i}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}^{i}~{}\approx~{}y_{i}(h_{i})~{}\equiv~{}y(t_{i}),\qquad
i=1,\dots,n,$ (38)

and assume we want to compute

$\sigma_{n+1}(ch_{n+1}):=\phi^{\alpha}_{n}(c,\sigma)+h_{n+1}^{\alpha}\sum_{j=0}^{s-1}\gamma_{j}^{n+1}I^{\alpha}P_{j}(c)~{}\approx~{}y_{n+1}(ch_{n+1}),\qquad
c\in[0,1].$ (39)

Hereafter, we shall assume that the time-steps in (26) define a graded mesh.
In more detail, for a suitable $r>1$:

$h_{n}~{}=~{}rh_{n-1}~{}\equiv~{}r^{n-1}h_{1},\qquad n=1,\dots,N.$ (40)

###### Remark 5

As is clear, in the limit case where $r=1$, (40) reduces to a uniform mesh
with a constant time-step $h=T/N$.

In order to derive the approximation (39), we start computing the solution of
the problem in the subinterval $[t_{n},t_{n+1}]$. Then, for $t\equiv
t_{n}+ch_{n+1}$, $c\in[0,1]$, one has:

$\displaystyle
y(t)~{}\equiv~{}y_{n+1}(ch_{n+1})~{}=~{}y_{0}~{}+~{}\frac{1}{\Gamma(\alpha)}\int_{0}^{t_{n}+ch_{n+1}}(t_{n}+ch_{n+1}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle=$ $\displaystyle
y_{0}~{}+\frac{1}{\Gamma(\alpha)}\int_{0}^{t_{n}}(t_{n}+ch_{n+1}-x)^{\alpha-1}f(y(x))\mathrm{d}x~{}+\frac{1}{\Gamma(\alpha)}\int_{t_{n}}^{t_{n}+ch_{n+1}}(t_{n}+ch_{n+1}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle=$ $\displaystyle
y_{0}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}\int_{t_{\nu-1}}^{t_{\nu}}(t_{n}+ch_{n+1}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle~{}+~{}\frac{1}{\Gamma(\alpha)}\int_{t_{n}}^{t_{n}+ch_{n+1}}(t_{n}+ch_{n+1}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle=$ $\displaystyle
y_{0}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}\int_{0}^{h_{\nu}}(t_{n}-t_{\nu-1}+ch_{n+1}-x)^{\alpha-1}f(y_{\nu}(x))\mathrm{d}x$
$\displaystyle~{}+~{}\frac{1}{\Gamma(\alpha)}\int_{0}^{ch_{n+1}}(ch_{n+1}-x)^{\alpha-1}f(y_{n+1}(x))\mathrm{d}x$
$\displaystyle=$ $\displaystyle
y_{0}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}+cr^{n-\nu+1}-\tau\right)^{\alpha-1}f(y_{\nu}(\tau
h_{\nu}))\mathrm{d}\tau$
$\displaystyle~{}+~{}\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}f(y_{n+1}(\tau
h_{n+1}))\mathrm{d}\tau~{}\equiv~{}G_{n}^{\alpha}(c,y)~{}+~{}I^{\alpha}f(y_{n+1}(ch_{n+1})),$

having set the history term

$G_{n}^{\alpha}(c,y)~{}:=~{}y_{0}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}+cr^{n-\nu+1}-\tau\right)^{\alpha-1}f(y_{\nu}(\tau
h_{\nu}))\mathrm{d}\tau,$ (41)

which is a known quantity, since it only depends on $y_{1},\dots,y_{n}$ (see
(35)). Consequently, we have obtained

$y_{n+1}(ch_{n+1})=G_{n}^{\alpha}(c,y)+I^{\alpha}f(y_{n+1}(ch_{n+1})),\qquad
c\in[0,1],$ (42)

which reduces to (5), when $n=0$ and $t\in[0,h]$. Further, by considering the
expansion

$f(y_{n+1}(ch_{n+1}))=\sum_{j\geq 0}\gamma_{j}(y_{n+1})P_{j}(c),\qquad
c\in[0,1],$ (43)

with the Fourier coefficients given by:

$\gamma_{j}(y_{n+1})=\alpha\int_{0}^{1}(1-c)^{\alpha-1}P_{j}(c)f(y_{n+1}(ch_{n+1}))\mathrm{d}c,\qquad
j\geq 0,$ (44)

one obtains that (42) can be rewritten as:

$y_{n+1}(ch_{n+1})=G_{n}^{\alpha}(c,y)+\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha)}\sum_{j\geq
0}\gamma_{j}(y_{n+1})I^{\alpha}P_{j}(c),\qquad c\in[0,1],$ (45)

with the value at $t=h$ given by:

$y(t_{n+1})~{}\equiv~{}y_{n+1}(h_{n+1})~{}=~{}G_{n}^{\alpha}(1,y)+\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}(y_{n+1}).$
(46)

The corresponding approximation is obtained by truncating the series in (42)
after $s$ terms, and approximating the corresponding Fourier coefficients via
the Gauss-Jacobi quadrature of order $2k$:

$\sigma_{n+1}(ch_{n+1})~{}=~{}\phi^{\alpha}_{n}(c,\sigma)~{}+~{}\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha)}\sum_{j=0}^{s-1}\gamma_{j}^{n+1}I^{\alpha}P_{j}(c),\qquad
c\in[0,1],$ (47)

and

$\bar{y}_{n+1}~{}:=~{}\sigma_{n+1}(h_{n+1})~{}=~{}\phi^{\alpha}_{n}(1,\sigma)~{}+~{}\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}^{n+1},$
(48)

with $\phi_{n}^{\alpha}(c,\sigma)$ an approximation of $G_{n}^{\alpha}(c,y)$
in (41), defined as follows:444As anticipated above, from (49) it follows that
$\phi_{0}^{\alpha}(c,\sigma)\equiv y_{0}$, $c\in[0,1]$, so that (47) reduces
to (31), whereas $\phi_{n}^{\alpha}(c,\sigma)$, $n>1$, only depends on
$\sigma_{1},\dots,\sigma_{n}$.

$\displaystyle\phi_{n}^{\alpha}(c,\sigma)$ $\displaystyle:=$ $\displaystyle
y_{0}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}+cr^{n-\nu+1}-\tau\right)^{\alpha-1}\sum_{j=0}^{s-1}P_{j}(\tau)\gamma_{j}^{\nu}\mathrm{d}\tau$
(49) $\displaystyle=$ $\displaystyle
y_{0}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\sum_{j=0}^{s-1}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}+cr^{n-\nu+1}-\tau\right)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau\gamma_{j}^{\nu}$
$\displaystyle\equiv$ $\displaystyle
y_{0}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\sum_{j=0}^{s-1}J_{j}^{\alpha}\left(\frac{r^{n-\nu+1}-1}{r-1}+cr^{n-\nu+1}\right)\gamma_{j}^{\nu},$

having set

$J_{j}^{\alpha}(x):=\int_{0}^{1}(x-\tau)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau,\qquad
j=0,\dots,s-1,\qquad x>1,$ (50)

which, as we shall see in Section 3.3, can be accurately and efficiently
computed.

As is clear, in (49) we have used the approximation

$f(\sigma_{\nu}(\tau
h))~{}=~{}\sum_{j=0}^{s-1}P_{j}(\tau)\gamma_{j}^{\nu}~{}+~{}R_{\nu}^{\alpha}(\tau),$
(51)

with the error term given by (see (17)-(18) and (25)):

$\displaystyle R_{\nu}^{\alpha}(\tau)$ $\displaystyle:=$
$\displaystyle\sum_{j=0}^{s-1}P_{j}(\tau)\left[\gamma_{j}(\sigma_{\nu})-\gamma_{j}^{\nu}\right]~{}+~{}\overbrace{\sum_{j\geq
s}P_{j}(\tau)\gamma_{j}(\sigma_{\nu})}^{=:\,E_{\nu}^{\alpha}(\tau)}$ (52)
$\displaystyle\equiv$
$\displaystyle\underbrace{\sum_{j=0}^{s-1}P_{j}(\tau)\Sigma_{j}^{\alpha}(\sigma_{\nu},h_{\nu})}_{=:\,S_{\nu}^{\alpha}(\tau)}~{}+~{}E_{\nu}^{\alpha}(\tau)~{}\equiv~{}S_{\nu}^{\alpha}(\tau)+E_{\nu}^{\alpha}(\tau).$

Because of the results of Section 2, we shall hereafter assume that

$\left|\Sigma_{j}^{\alpha}(\sigma_{\nu},h_{\nu})\right|\leq\left\\{\begin{array}[]{cc}O(h_{1}^{\alpha}),&\nu=1,\\\\[5.69054pt]
O(h_{\nu}^{2k-j}),&\nu>1,\end{array}\right.\qquad\left|E_{\nu}^{\alpha}(\tau)\right|\leq\left\\{\begin{array}[]{cc}O(h_{1}^{\alpha}),&\nu=1,\\\\[5.69054pt]
O(h_{\nu}^{s}),&\nu>1.\end{array}\right.$ (53)

Consequently, considering that $k\geq s$, $\alpha\in(0,1)$, and
$j=0,\dots,s-1$, one deduces:

$\left|S_{\nu}^{\alpha}(\tau)\right|\leq\left\\{\begin{array}[]{cc}O(h_{1}^{\alpha}),&\nu=1,\\\\[5.69054pt]
O(h_{\nu}^{s+1}),&\nu>1,\end{array}\right.\qquad\left|R_{\nu}^{\alpha}(\tau)\right|\leq\left\\{\begin{array}[]{cc}O(h_{1}^{\alpha}),&\nu=1,\\\\[5.69054pt]
O(h_{\nu}^{s}),&\nu>1.\end{array}\right.$ (54)

In so doing, see (47), the Fourier coefficients satisfy the system of
equations:

$\gamma_{j}^{n+1}~{}=~{}\sum_{\ell=1}^{k}b_{\ell}P_{j}(c_{\ell})f_{n+1}\left(\phi_{n}^{\alpha}(c_{\ell},\sigma)+h_{n+1}^{\alpha}\sum_{i=0}^{s-1}\gamma_{i}^{n+1}I^{\alpha}P_{i}(c_{\ell}h)\right),\quad
j=0,\dots,s-1,$ (55)

having (block)-size $s$, independently of the considered value of $k$.

Clearly,

$\gamma_{j}^{n+1}~{}=~{}\gamma_{j}(\sigma_{n+1})~{}+~{}\Sigma_{j}^{\alpha}(\sigma_{n+1},h_{n+1}),\qquad
j=0,\dots,s-1$

with $\gamma_{j}(\sigma_{n+1})$ formally defined as in (44), upon replacing
$y_{n+1}$ with $\sigma_{n+1}$, and

$\Sigma_{j}^{\alpha}(\sigma_{n+1},h_{n+1})=O(h_{n+1}^{2k-j}),\qquad
j=0,\dots,s-1,$

the corresponding quadrature errors.

###  3.2 Solving the discrete problems

The discrete problem (55), to be solved at each integration step, can be
better cast in vector form by introducing the matrices

${\cal P}_{s}=\left(\begin{array}[]{ccc}P_{0}(c_{1})&\dots&P_{s-1}(c_{1})\\\
\vdots&&\vdots\\\
P_{0}(c_{k})&\dots&P_{s-1}(c_{k})\end{array}\right),\quad{\cal
I}_{s}^{\alpha}=\left(\begin{array}[]{ccc}I^{\alpha}P_{0}(c_{1})&\dots&I^{\alpha}P_{s-1}(c_{1})\\\
\vdots&&\vdots\\\
I^{\alpha}P_{c}(c_{k})&\dots&I^{\alpha}P_{s-1}(c_{k})\end{array}\right)\quad\in\mathbb{R}^{k\times
s},$ $\Omega=\left(\begin{array}[]{ccc}b_{1}\\\ &\ddots\\\
&&b_{k}\end{array}\right),$

and the (block) vectors:

${\bm{\gamma}}^{n+1}=\left(\begin{array}[]{c}\gamma_{0}^{n+1}\\\ \vdots\\\
\gamma_{s-1}^{n+1}\end{array}\right)\in\mathbb{R}^{sm},\qquad{\bm{\phi}}_{n}^{\alpha}=\left(\begin{array}[]{c}\phi_{n}^{\alpha}(c_{1},\sigma)\\\
\vdots\\\
\phi_{n}^{\alpha}(c_{k},\sigma)\end{array}\right)\in\mathbb{R}^{km}.$

In fact, in so doing (55) can be rewritten as:555As is usual, the function
$f$, here evaluated in a (block) vector of dimension $k$, denotes the (block)
vector made up by $f$ evaluated in each (block) entry of the input argument.

${\bm{\gamma}}^{n+1}={\cal P}_{s}^{\top}\Omega\otimes
I_{m}f\left({\bm{\phi}}_{n}^{\alpha}+h_{n+1}{\cal I}_{s}^{\alpha}\otimes
I_{m}{\bm{\gamma}}^{n+1}\right),$ (56)

This formulation is very similar to that used for HBVMs in the case $\alpha=1$
[11]. The formulation (56) has a twofold use:

  * •

it shows that, assuming, for example, $f$ Lipschitz continuous, then the
solution exists and is unique, for all sufficiently small timesteps $h_{n+1}$;

  * •

it induces a straightforward fixed-point iteration:

$\displaystyle{\bm{\gamma}}^{n+1,0}$ $\displaystyle=$ $\displaystyle{\bm{0}},$
(57) $\displaystyle{\bm{\gamma}}^{n+1,\ell}$ $\displaystyle=$
$\displaystyle{\cal P}_{s}^{\top}\Omega\otimes
I_{m}f\left({\bm{\phi}}_{n}^{\alpha}+h_{n+1}{\cal I}_{s}^{\alpha}\otimes
I_{m}{\bm{\gamma}}^{n+1,\ell-1}\right),\qquad\ell=1,2,\dots,$

which we shall use for the numerical tests.

In fact, the following result holds true.

###### Theorem 2

Assume $f$ be Lipchitz with constant $L$ in in the interval $[t_{n},t_{n+1}]$.
Then, the iteration (57) is convergent for all timesteps $h_{n+1}$ such that

$h_{n+1}L\|{\cal P}_{s}^{\top}\Omega\|\|{\cal I}_{s}^{\alpha}\|<1.$ (58)

Proof In fact, one has:

$\displaystyle\|{\bm{\gamma}}^{n+1,\ell+1}-{\bm{\gamma}}^{n+1,\ell}\|$
$\displaystyle=$ $\displaystyle\|{\cal P}_{s}^{\top}\Omega\otimes
I_{m}\left[f\left({\bm{\phi}}_{n}^{\alpha}+h_{n+1}{\cal I}_{s}^{\alpha}\otimes
I_{m}{\bm{\gamma}}^{n+1,\ell}\right)-f\left({\bm{\phi}}_{n}^{\alpha}+h_{n+1}{\cal
I}_{s}^{\alpha}\otimes I_{m}{\bm{\gamma}}^{n+1,\ell-1}\right)\right]\|$
$\displaystyle\leq$ $\displaystyle h_{n+1}L\|{\cal
P}_{s}^{\top}\Omega\|\|{\cal
I}_{s}^{\alpha}\|\cdot\|{\bm{\gamma}}^{n+1,\ell}-{\bm{\gamma}}^{n+1,\ell-1}\|,$

hence the iteration function defined at (57) is a contraction, when (58) holds
true.  $\Box{~{}}$

A simplified Newton-type iteration, akin to that defined in [11] for HBVMs,
will be the subject of future investigations.

###### Remark 6

We observe that the discrete problem (56) can be cast in a Runge-Kutta type
form. In fact, the vector

$Y^{n+1}~{}:=~{}{\bm{\phi}}_{n}^{\alpha}+h_{n+1}{\cal I}_{s}^{\alpha}\otimes
I_{m}{\bm{\gamma}}^{n+1},$ (59)

in argument to the function $f$ at the r.h.s. in (56), can be regarded as the
stage vector of a Runge-Kutta method, tailored for the problem at hand.
Substituting (56) into (59), and using $Y^{n+1}$ as argument of $f$, then
gives

$Y^{n+1}~{}=~{}\phi_{n}^{\alpha}+h_{n+1}{\cal I}_{s}^{\alpha}{\cal
P}_{s}^{\top}\Omega\otimes I_{m}f(Y^{n+1}).$ (60)

It is worth noticing that (60) has (block) dimension $k$, instead of $s$.
However, by considering that usually $k\gg s$ (see Section 5), it follows that
solving (60) is less efficient than solving (56).

This fact is akin to what happens for a HBVM$(k,s)$ method [12, 8, 9], which
is the $k$-stage Runge-Kutta method obtained when $\alpha=1$. In fact, for
such a method, the discrete problem can also be cast in the form (56), then
having (block) dimension $s$, independently of $k$ (which is usually much
larger than $s$).

###  3.3 Computing the integrals $J_{j}^{\alpha}(x)$

From (55) and (49), it follows that one needs to compute the integrals

$J_{j}^{\alpha}\left(\frac{r^{n-\nu+1}-1}{r-1}+c_{i}r^{n-\nu+1}\right),\qquad
i=1,\dots,k,\qquad\nu=1,\dots,n,$ (61)

with $\\{c_{1},\dots,c_{k}\\}$ the abscissae of the Gauss-Jacobi quadrature
(16). As an example, in Figure 1 we plot the Gauss-Jacobi abscissae for
$\alpha=0.5$ and $k=5,10,15,20,25,30$.

Figure 1: Gauss-Jacobi abscissae for $\alpha=0.5$.

Further, there is numerical evidence that a sufficiently high-order Gauss-
Legendre formula can compute the integrals (50) up to round-off, when $x\geq
1.5$. Namely,

$J_{j}^{\alpha}\left(\frac{r^{n-\nu+1}-1}{r-1}+c_{i}r^{n-\nu+1}\right),\qquad
i=1,\dots,k,\qquad\nu=1,\dots,n-1,$ (62)

and

$J_{j}^{\alpha}\left(1+c_{i}r\right),\qquad{\rm s.t.}\qquad c_{i}r\geq
0.5,\qquad j=0,\dots,s-1.$ (63)

However, this is no more the case, when computing:

$J_{j}^{\alpha}\left(1+c_{i}r\right),\qquad{\rm s.t.}\qquad c_{i}r<0.5,\qquad
j=0,\dots,s-1.$ (64)

Consequently, there is the problem of accurately and efficiently computing
these latter integrals.

###### Remark 7

We observe that the evaluation of the integrals (62)-(63), $j=0,\dots,s-1$,
via a $2p$-order Gauss-Legendre formula is inexpensive. As matter of fact,
only the values $P_{j}(\xi_{i})$ need to be computed, $j=0,\dots,s-1$, with

$\\{\xi_{1},\dots,\xi_{p}\\}$ (65)

the Gauss-Legendre abscissae of the considered formula. Further, we observe
that, with reference to (61), only the integrals

$J_{j}^{\alpha}\left(\frac{r^{n}-1}{r-1}+c_{i}r^{n}\right),\qquad
j=0,\dots,s-1,\qquad i=1,\dots,k,$

need to be actually computed: as matter of fact, the remaining integrals,

$J_{j}^{\alpha}\left(\frac{r^{\nu}-1}{r-1}+c_{i}r^{\nu}\right),\qquad
j=0,\dots,s-1,\qquad i=1,\dots,k,\qquad\nu=1,\dots,n-1,$

are inherited from the previous steps.

The starting point to derive an efficient algorithm for computing the
integrals (64), is that the Jacobi polynomials (6) satisfy, as does any other
family of orthogonal polynomials, a three-term recurrence:

$P_{j+1}(c)=(a_{j}c-b_{j})P_{j}(c)-d_{j}P_{j-1}(c),\quad j=0,1,\dots,\qquad
P_{0}(c)\equiv 1,\quad P_{-1}(c)\equiv 0,$ (66)

with prescribed $a_{j},b_{j},d_{j}$, $j\geq 0$. In fact, by setting
$\phi(c)=c$, and using the scalar product (9), to enforce (7), for
$P_{0},\dots,P_{s-1}$, one obtains:

$a_{j}=\frac{1}{(P_{j+1},\phi\cdot P_{j})},\qquad b_{j}=\frac{(P_{j},\phi\cdot
P_{j})}{(P_{j+1},\phi\cdot P_{j})},\qquad d_{j}=\frac{(P_{j},\phi\cdot
P_{j-1})}{(P_{j+1},\phi\cdot P_{j})},\qquad j=0,\dots,s-2.$ (67)

Consequently, by recalling the definition (50), from (66) one has:

$J_{0}^{\alpha}(x)=\frac{x^{\alpha}-(x-1)^{\alpha}}{\alpha},\qquad
J_{-1}^{\alpha}(x)=0,\qquad\alpha>0,\qquad x>1,$

and

$\displaystyle J_{j+1}^{\alpha}(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}(x-\tau)^{\alpha-1}P_{j+1}(\tau)\mathrm{d}\tau$
$\displaystyle=$ $\displaystyle a_{j}\int_{0}^{1}(x-\tau)^{\alpha-1}\tau
P_{j}(\tau)\mathrm{d}\tau-b_{j}J_{j}^{\alpha}(x)-d_{j}J_{j-1}^{\alpha}(x)$
$\displaystyle=$ $\displaystyle-
a_{j}\int_{0}^{1}(x-\tau)^{\alpha}P_{j}(\tau)\mathrm{d}\tau+[a_{j}x-b_{j}]J_{j}^{\alpha}(x)-d_{j}J_{j-1}^{\alpha}(x)$
$\displaystyle=$ $\displaystyle-
a_{j}J_{j}^{\alpha+1}(x)+[a_{j}x-b_{j}]J_{j}^{\alpha}(x)-d_{j}J_{j-1}^{\alpha}(x),\qquad
j=0,\dots,j-2.$

From the last two formulae, one derives that Jjalfa(a,b,d,alfa,1+ci*r)
computes all the integrals in (64) corresponding to the abscissa ci, with
Jjalfa the Matlab© function listed in Table 1, and the vectors a,b,d
containing the scalars (67). An implementation of the previous function, using
the standard variable precision arithmetic (i.e., using vpa in Matlab©),
allows handling values of $s$ up to 20, at least.

Table 1:  Matlab© function Jjalfa

    function Ialfa = Jjalfa( a, b, d, alfa, x )
    %
    % Matlab function for computing the integrals J_j^alfa(x), j=0...s-1.
    %
    s1     = length(a);  % s1 == s-1
    Ialfa  = zeros(s1+1,1);
    valfa  = alfa+(0:s1);
    I1     = ( x.^valfa -(x-1).^valfa )./valfa;
    Ialfa(1) = I1(1);
    if s1>=1
       I2 = ( a(1)*x -b(1) )*I1(1:s1) -a(1)*I1(2:s1+1);
       Ialfa(2) = I2(1);
       for j = 2:s1
           I0 = I1; I1 = I2;
           I2 = ( a(j)*x -b(j) )*I1(1:s1-j+1) -a(j)*I1(2:s1-j+2) -d(j)*I0(1:s1-j+1);
           Ialfa(j+1) = I2(1);
       end
    end
    return

##  4 Analysis of the method

From (45) and (47), and with reference to (51)–(54), one derives:

$\displaystyle y_{n+1}(ch_{n+1})-\sigma_{n+1}(ch_{n+1})$ $\displaystyle=$
$\displaystyle\left[G_{n}^{\alpha}(c,y)-\phi_{n}^{\alpha}(c,\sigma)\right]~{}+~{}\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha)}\sum_{j=0}^{s-1}\left[\gamma_{j}(y_{n+1})-\gamma_{j}^{n+1}\right]I^{\alpha}P_{j}(c)$
$\displaystyle~{}+~{}\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha)}\sum_{j\geq
s}\gamma_{j}(y_{n+1})I^{\alpha}P_{j}(c)$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}+cr^{n-\nu+1}-\tau\right)^{\alpha-1}\left[f(y_{\nu}(\tau
h_{\nu}))-f(\sigma_{\nu}(\tau
h_{\nu})+R_{\nu}^{\alpha}(\tau)\right]\mathrm{d}\tau$
$\displaystyle+\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha)}\sum_{j=0}^{s-1}\left[\gamma_{j}(y_{n+1})-\gamma_{j}(\sigma_{n+1})+\Sigma_{j}(\sigma_{n+1},h_{n+1})\right]I^{\alpha}P_{j}(c)~{}+~{}\frac{h_{n+1}^{\alpha}}{\Gamma(\alpha)}E_{n+1}^{\alpha}(c).$

Assuming $f$ Lipschitz with constant $L$ in a suitable neighborhood of the
solution, and setting

$\displaystyle R_{\nu}^{\alpha}$ $\displaystyle:=$
$\displaystyle\max_{\tau\in[0,1]}|R_{\nu}^{\alpha}(\tau)|,\qquad\nu=1,\dots,N-1,$
$\displaystyle E_{n+1}^{\alpha}$ $\displaystyle:=$
$\displaystyle\max_{c\in[0,1]}|E_{n+1}^{\alpha}(c)|,\qquad n=1,\dots,N-1,$
$\displaystyle e_{n}$ $\displaystyle:=$
$\displaystyle\max_{c\in[0,1]}|y_{n}(ch_{n})-\sigma_{n}(ch_{n})|,\qquad
n=1,2,\dots,N,$

one then obtains that, for $h_{n+1}^{\alpha}$ sufficiently small, and a
suitable constant $K_{1}>0$:

$\displaystyle e_{n+1}$ $\displaystyle\leq$
$\displaystyle\frac{K_{1}L}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}-\tau\right)^{\alpha-1}\mathrm{d}\tau
e_{\nu}$
$\displaystyle~{}+~{}\frac{K_{1}}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}-\tau\right)^{\alpha-1}\mathrm{d}\tau
R_{\nu}^{\alpha}~{}+~{}g_{n+1},\qquad n=1,\dots,N-1,$

with

$g_{n+1}^{\alpha}=O(h_{n+1}^{s+\alpha}).$ (68)

Considering that

$\displaystyle\frac{1}{\Gamma(\alpha)}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}-\tau\right)^{\alpha-1}\mathrm{d}\tau$
(69) $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha+1)}\left[\left(\frac{r^{n-\nu+1}-1}{r-1}\right)^{\alpha}-\left(\frac{r^{n-\nu+1}-r}{r-1}\right)^{\alpha}\right]~{}\leq~{}\frac{1}{\Gamma(\alpha+1)},\qquad\nu=1,\dots,n,~{}$

and (see (54))

$h_{1}^{\alpha}R_{1}^{\alpha}\leq O(h_{1}^{2\alpha}),\qquad
h_{\nu}^{\alpha}R_{\nu}^{\alpha}\leq O(h_{\nu}^{s+\alpha}),\quad\nu>1,$

so that (see (40)), for a suitable constant $K_{2}>0$,

$\displaystyle\frac{K_{1}}{\Gamma(\alpha+1)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}R_{\nu}^{\alpha}$
$\displaystyle\leq$ $\displaystyle
K_{2}\left(h_{1}^{2\alpha}+\sum_{\nu=2}^{n}h_{\nu}^{s+\alpha}\right)~{}=~{}K_{2}\left(h_{1}^{2\alpha}+h_{n}^{s+\alpha}\sum_{\nu=0}^{n-2}(r^{s+\alpha})^{-\nu}\right)$
$\displaystyle\leq$ $\displaystyle
K_{2}\left(h_{1}^{2\alpha}+\frac{h_{n}^{s+\alpha}}{1-r^{-(s+\alpha)}}\right),$

one eventually obtains:

$e_{n+1}~{}\leq~{}\frac{K_{1}L}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}-\tau\right)^{\alpha-1}\mathrm{d}\tau\,e_{\nu}~{}+~{}\psi_{n}^{\alpha},\qquad
n=1,\dots,N-1,$

with (see (40) and (68))

$\psi_{n}^{\alpha}~{}:=~{}K_{2}\left(h_{1}^{2\alpha}+\frac{h_{n}^{s+\alpha}}{1-r^{-(s+\alpha)}}\right)+g_{n+1}^{\alpha}~{}=~{}O(h_{1}^{2\alpha}+h_{n+1}^{s+\alpha}),\qquad
n=1,\dots,N-1.$

At last, by considering that, for $n=1,\dots,N-1$,

$\displaystyle\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu+1}-1}{r-1}-\tau\right)^{\alpha-1}\mathrm{d}\tau$
(70) $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha+1)}\sum_{\nu=1}^{n}\left[h_{\nu}^{\alpha}\left(\frac{r^{n-\nu+1}-1}{r-1}\right)^{\alpha}-h_{\nu+1}^{\alpha}\left(\frac{r^{n-\nu}-1}{r-1}\right)^{\alpha}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha+1)}\left[\left(h_{1}\frac{r^{n}-1}{r-1}\right)^{\alpha}-h_{n+1}^{\alpha}\right]~{}\leq~{}\frac{T^{\alpha}}{\Gamma(\alpha+1)},$

setting

$\rho=\exp\left(\frac{K_{1}LT^{\alpha}}{\Gamma(\alpha+1)}\right),$

and considering that $e_{1}\leq O(h_{1}^{2\alpha})$, from the Gronwall lemma
(see, e.g., [29]) one eventually obtains, for a suitable $K>0$:

$\displaystyle e_{n+1}$ $\displaystyle\leq$ $\displaystyle
K\rho\left(h_{1}^{2\alpha}+\sum_{\nu=1}^{n}\psi_{n}^{\alpha}\right)~{}=~{}K\rho\left((n+1)h_{1}^{2\alpha}+\frac{h_{n+1}^{s+\alpha}}{1-r^{-(s+\alpha)}}\right)$
(71) $\displaystyle\leq$ $\displaystyle
K\rho\left(Nh_{1}^{2\alpha}+\frac{h_{N}^{s+\alpha}}{1-r^{-(s+\alpha)}}\right),\qquad
n=1,\dots,N-1.$

Considering that (see (27) and (40)) $N=\log_{r}\left(1+T(r-1)/h_{1}\right)$,
the error is optimal when $h_{1}^{2\alpha}\sim h_{N}^{s+\alpha}$.

###### Remark 8

In the case where the vector field of problem (1) is suitably smooth at $t=0$,
so that a constant timestep $h=T/N$ can be used, the estimate

$e_{n}\leq O(nh^{s+\alpha}),\qquad n=1,\dots,N,$ (72)

can be derived for the error, by using similar arguments as above.

##  5 Numerical Tests

We here report a few numerical tests to illustrate the theoretical findings.
For all tests, when not otherwise specified, we use $k=30$ for the Gauss-
Jacobi quadrature (16), and $p=30$ for the Gauss-Legendre quadrature formula
(65). Also, the following values of $s$ will be considered:

$s\in\\{1,2,3,4,5,6,7,8,9,10,20\\}.$

All the numerical tests have been performed in Matlab©. As anticipated, we use
the fixed-point iteration (57) to solve the generated discrete problems (56).
The iteration is carried out until full machine accuracy is gained.

The first problem, taken from [24], is given by:

$y^{(0.6)}=-10y,\qquad t\in[0,5],\qquad y(0)=1,$ (73)

whose solution is given by the following Mittag-Leffler function:666We refer
to [22] and the accompanying software, for its efficient Matlab©
implementation.

$E_{0.6}(t)=\sum_{j\geq 0}\frac{(-10t^{0.6})^{j}}{\Gamma(0.6j+1)}.$

In Table 2, we list the obtained results, in terms of maximum error, when
using different values of $s$ and $h_{1}$. We use the value $r=1.01$ in (40),
and the number of timesteps $N$ is chosen in order that $t_{N}$ is the closest
point to the end of the integration interval. The ***, in the line
corresponding to $s=1$, means that the solution is not properly evaluated:
this is due to the failure of the fixed-point iteration (57) in the last
integration steps. As is expected from (71), the error decreases as $s$
increases and the initial timestep $h_{1}$ decreases. Further, having a large
value of $s$ is not effective, if $h_{1}$ is not suitably small, and vice
versa (again from the bound (71)). Remarkably, by choosing $s$ large enough
and $h_{1}$ suitably small, full machine accuracy is gained (cf. the last 4
entries in the last column of the table, having the same error).

In Table 3 we list the results obtained by using $k=s$, which is the minimum
value allowed. In such a case, the accuracy is generally slightly worse, due
to the fact that the quadrature error is of the same order as the truncation
error. It is, however, enough choosing $k$ only slightly larger than $s$, to
achieve a comparable accuracy, as is shown in Table 4 for $k=s+5$.
Nevertheless, it must be stressed that choosing larger values of $k$ is not an
issue, since the discrete problem (56), to be solved at each integration step,
has (block) size $s$, independently of $k$.

Table 2:  Maximum error for Problem (73), $r=1.01$ and $k=30$. $h_{1}$ | $10^{-4}$ | $10^{-5}$ | $10^{-6}$ | $10^{-7}$ | $10^{-8}$ | $10^{-9}$
---|---|---|---|---|---|---
$s\,\backslash\,N$ | 626 | 857 | 1088 | 1320 | 1551 | 1783
1 | *** | *** | *** | *** | *** | ***
2 | 3.73e-06 | 2.43e-07 | 5.40e-08 | 5.38e-08 | 5.37e-08 | 5.37e-08
3 | 5.81e-07 | 3.78e-08 | 2.40e-09 | 1.52e-10 | 4.64e-11 | 4.64e-11
4 | 1.47e-07 | 9.60e-09 | 6.11e-10 | 3.86e-11 | 2.44e-12 | 1.56e-13
5 | 5.29e-08 | 3.46e-09 | 2.20e-10 | 1.39e-11 | 8.79e-13 | 5.37e-14
6 | 2.04e-08 | 1.33e-09 | 8.49e-11 | 5.37e-12 | 3.41e-13 | 2.32e-14
7 | 1.19e-08 | 7.81e-10 | 4.98e-11 | 3.15e-12 | 1.97e-13 | 1.08e-14
8 | 4.26e-09 | 2.77e-10 | 1.76e-11 | 1.11e-12 | 7.18e-14 | 7.91e-15
9 | 4.56e-09 | 3.03e-10 | 1.94e-11 | 1.22e-12 | 7.57e-14 | 7.91e-15
10 | 1.89e-09 | 1.22e-10 | 7.79e-12 | 4.90e-13 | 2.96e-14 | 7.91e-15
20 | 1.84e-09 | 1.22e-10 | 7.77e-12 | 4.90e-13 | 2.96e-14 | 7.91e-15
Table 3:  Maximum error for Problem (73), $r=1.01$ and $k=s$. $h_{1}$ | $10^{-4}$ | $10^{-5}$ | $10^{-6}$ | $10^{-7}$ | $10^{-8}$ | $10^{-9}$
---|---|---|---|---|---|---
$s\,\backslash\,N$ | 626 | 857 | 1088 | 1320 | 1551 | 1783
1 | *** | *** | *** | *** | *** | ***
2 | 7.52e-06 | 4.99e-07 | 5.78e-08 | 5.76e-08 | 5.76e-08 | 5.76e-08
3 | 2.27e-06 | 1.51e-07 | 9.68e-09 | 6.13e-10 | 4.55e-11 | 4.55e-11
4 | 9.67e-07 | 6.41e-08 | 4.09e-09 | 2.59e-10 | 1.64e-11 | 1.03e-12
5 | 4.91e-07 | 3.26e-08 | 2.08e-09 | 1.32e-10 | 8.33e-12 | 5.23e-13
6 | 2.82e-07 | 1.87e-08 | 1.19e-09 | 7.56e-11 | 4.77e-12 | 2.99e-13
7 | 1.76e-07 | 1.16e-08 | 7.44e-10 | 4.71e-11 | 2.97e-12 | 1.85e-13
8 | 1.16e-07 | 7.71e-09 | 4.92e-10 | 3.12e-11 | 1.97e-12 | 1.22e-13
9 | 8.07e-08 | 5.35e-09 | 3.42e-10 | 2.16e-11 | 1.36e-12 | 8.42e-14
10 | 5.82e-08 | 3.86e-09 | 2.46e-10 | 1.56e-11 | 9.83e-13 | 6.01e-14
20 | 6.63e-09 | 4.39e-10 | 2.80e-11 | 1.77e-12 | 1.11e-13 | 7.80e-15
Table 4:  Maximum error for Problem (73), $r=1.01$ and $k=s+5$. $h_{1}$ | $10^{-4}$ | $10^{-5}$ | $10^{-6}$ | $10^{-7}$ | $10^{-8}$ | $10^{-9}$
---|---|---|---|---|---|---
$s\,\backslash\,N$ | 626 | 857 | 1088 | 1320 | 1551 | 1783
1 | *** | *** | *** | *** | *** | ***
2 | 3.58e-06 | 2.33e-07 | 5.40e-08 | 5.38e-08 | 5.38e-08 | 5.37e-08
3 | 6.56e-07 | 4.29e-08 | 2.73e-09 | 1.73e-10 | 4.63e-11 | 4.63e-11
4 | 8.57e-08 | 5.43e-09 | 3.43e-10 | 2.16e-11 | 1.36e-12 | 8.76e-14
5 | 9.20e-08 | 6.12e-09 | 3.92e-10 | 2.48e-11 | 1.56e-12 | 9.70e-14
6 | 4.38e-08 | 2.87e-09 | 1.83e-10 | 1.16e-11 | 7.30e-13 | 4.41e-14
7 | 3.39e-08 | 2.28e-09 | 1.46e-10 | 9.25e-12 | 5.83e-13 | 3.51e-14
8 | 2.58e-08 | 1.70e-09 | 1.09e-10 | 6.87e-12 | 4.32e-13 | 2.53e-14
9 | 2.03e-08 | 1.35e-09 | 8.60e-11 | 5.44e-12 | 3.43e-13 | 1.97e-14
10 | 1.64e-08 | 1.09e-09 | 6.93e-11 | 4.38e-12 | 2.76e-13 | 1.54e-14
20 | 3.28e-09 | 2.17e-10 | 1.38e-11 | 8.74e-13 | 5.40e-14 | 8.19e-15

The next problem that we consider is from [21] (see also [24]):

$y^{(0.5)}=-y^{1.5}+\frac{40320}{\Gamma(8.5)}t^{7.5}-3\frac{\Gamma(5.25)}{\Gamma(4.75)}t^{3.75}+\left(\frac{3}{2}t^{0.25}-t^{4}\right)^{3}+\frac{9}{4}\Gamma(1.5),\quad
t\in[0,1],\quad y(0)=0,$ (74)

whose solution is

$y(t)=t^{8}-3\,t^{4.25}+\frac{9}{4}\,t^{0.5}.$ (75)

According to [24], “this problem is surely of interest because, unlike several
other problems often proposed in the literature, it does not present an
artificial smooth solution, which is indeed not realistic in most of the
fractional-order applications.” Despite this, a constant timestep $h=1/N$
turns out to be appropriate since, unlike the solution (75), the vector field
(74) is very smooth at the origin, as one may see in Figure 2. In Table 5 we
list the maximum error by using different values of $N$: as is clear, by using
$s\geq 8$, only 32 steps are needed to gain full machine accuracy for the
computed solution.

Figure 2: solution (continuous line) and vector field (dashed line) for problem (74). Table 5:  Maximum error for Problem (74), constant timestep $h=1/N$. $s\,\backslash\,N$ | 2 | 4 | 8 | 16 | 32
---|---|---|---|---|---
1 | 9.22e-01 | 5.65e-02 | 1.28e-02 | 1.35e-02 | 9.12e-03
2 | 7.48e-03 | 2.68e-03 | 5.15e-04 | 8.02e-05 | 1.91e-05
3 | 2.02e-03 | 1.96e-04 | 1.23e-05 | 2.04e-06 | 5.07e-07
4 | 2.29e-04 | 8.42e-06 | 2.72e-07 | 3.55e-08 | 3.70e-09
5 | 1.63e-05 | 3.52e-07 | 4.43e-09 | 3.44e-10 | 1.62e-11
6 | 7.61e-07 | 9.80e-09 | 6.57e-11 | 2.26e-12 | 1.47e-13
7 | 4.11e-08 | 3.71e-10 | 9.02e-12 | 3.46e-13 | 2.18e-14
8 | 1.24e-09 | 6.02e-11 | 1.87e-12 | 6.54e-14 | 4.22e-15
9 | 4.56e-10 | 1.44e-11 | 4.27e-13 | 1.65e-14 | 1.11e-15
10 | 1.40e-10 | 4.40e-12 | 1.33e-13 | 4.77e-15 | 8.88e-16
20 | 4.93e-14 | 1.33e-15 | 6.66e-16 | 6.66e-16 | 8.88e-16

Next, we consider the following problems, taken from [38]:

$y^{(1/3)}=\frac{t}{10}\left[y^{3}-(t^{2/3}+1)^{3}\right]+\frac{\Gamma(5/3)}{\Gamma(4/3)}t^{1/3},\qquad
t\in[0,1],\qquad y(0)=1.$ (76)

whose solution is $y(t)=t^{2/3}+1$, and

$y^{(1/3)}=\frac{1}{3}\left(y^{3}-t^{4}\right)+\Gamma(7/3)t,\qquad
t\in[0,1],\qquad y(0)=0,$ (77)

whose solution is $y(t)=t^{4/3}$.

We solve Problem (76) by using $h_{1}=10^{-11}$ and $r=1.2$. In so doing 130
timesteps are needed to cover approximately the integration interval. The
obtained results are listed in Table 6. Also in this case, we obtain full
accuracy starting from $s=8$.

Table 6:  Maximum error for Problems (76) and (LABEL:prob34), $r=1.2$ and $h_{1}=10^{-11}$. $s$ |  error (76) |  error (LABEL:prob34)
---|---|---
1 | 3.25e-02 | ***
2 | 8.86e-05 | 5.13e-04
3 | 8.36e-07 | 4.21e-06
4 | 1.41e-08 | 7.55e-08
5 | 3.03e-10 | 1.63e-09
6 | 7.54e-12 | 3.95e-11
7 | 3.46e-13 | 1.06e-12
8 | 2.09e-13 | 2.09e-13
9 | 2.09e-13 | 2.09e-13
10 | 2.09e-13 | 2.09e-13
20 | 2.09e-13 | 2.09e-13

In addition, we solve Problem (77) by using a constant timestep $h=1/N$: in
fact, the vector field can be seen to be a polynomial of degree 1. The
obtained results are listed in Table 7: an order 1 convergence can be observed
for $s=1$ (which is consistent with (72)), whereas full machine accuracy is
obtained for $s\geq 2$, due to the fact that, as anticipated above, the vector
field of problem (77) is a polynomial of degree 1 in $t$ and, consequently,
(13) and (21) coincide, for all $s\geq 2$.

Table 7:  Maximum error for Problem (77), constant timestep $h=1/N$. $s\,\backslash\,N$ | 2 | 4 | 8 | 16 | 32 | 64
---|---|---|---|---|---|---
1 | *** | 1.56e-01 | 7.01e-02 | 3.59e-02 | 1.87e-02 | 9.75e-03
2 | 8.88e-16 | 1.33e-15 | 8.88e-16 | 8.88e-16 | 8.88e-16 | 8.88e-16
3 | 4.44e-16 | 6.66e-16 | 4.44e-16 | 4.44e-16 | 3.33e-16 | 4.44e-16
4 | 6.66e-16 | 6.66e-16 | 5.55e-16 | 2.22e-16 | 3.33e-16 | 5.55e-16
5 | 9.99e-16 | 9.99e-16 | 6.66e-16 | 5.55e-16 | 2.22e-16 | 5.55e-16
6 | 1.33e-15 | 7.77e-16 | 8.88e-16 | 6.66e-16 | 3.33e-16 | 5.55e-16
7 | 1.33e-15 | 8.88e-16 | 7.77e-16 | 4.44e-16 | 4.44e-16 | 7.77e-16
8 | 1.78e-15 | 1.11e-15 | 7.77e-16 | 6.66e-16 | 4.44e-16 | 6.66e-16
9 | 2.00e-15 | 1.22e-15 | 8.88e-16 | 8.88e-16 | 4.44e-16 | 6.66e-16
10 | 2.00e-15 | 1.22e-15 | 8.88e-16 | 8.88e-16 | 4.44e-16 | 7.77e-16
20 | 2.78e-15 | 1.89e-15 | 1.44e-15 | 1.11e-15 | 6.66e-16 | 8.88e-16

Finally, we reformulate the two problems (76) and (77) as a system of two
equations, having the same solutions as above, as follows:

$\displaystyle y_{1}^{(1/3)}$ $\displaystyle=$
$\displaystyle\frac{t}{10}\left[y_{1}^{3}-(y_{2}^{1/2}+1)^{3}\right]+\frac{\Gamma(5/3)}{\Gamma(4/3)}t^{1/3},\qquad
y_{1}(0)=1,$ $\displaystyle y_{2}^{(1/3)}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left(y_{2}^{3}-(y_{1}-1)^{6}\right)+\Gamma(7/3)t,~{}~{}\qquad\qquad
y_{2}(0)=0,\qquad t\in[0,1].$

We solve Problem (LABEL:prob34) by using the same parameters used for Problem
(76): $h_{1}=10^{-11}$, $r=1.2$, and 130 timesteps. The obtained results are
again listed in Table 6: it turns out that they are similar to those obtained
for Problem (76) and, also in this case, we obtain full accuracy starting from
$s=8$.

##  6 Conclusions

In this paper we have devised a novel step-by-step procedure for solving
fractional differential equations. The procedure, which generalizes that given
in [2], relies on the expansion of the vector field along a suitable
orthonormal basis, here chosen as the shifted and orthonormal Jacobi
polynomial basis. The analysis of the method has been given, along with its
implementation details. These latter details show that the method can be very
efficiently implemented. A few numerical tests confirm the theoretical
findings.

It is worth mentioning that systems with FDEs of different orders can be also
solved by slightly adapting the previous framework: as matter of fact, it
suffices using different Jacobi polynomials, corresponding to the different
orders of the FDEs. This, in turn, allows easily managing fractional
differential equations of order $\alpha>1$, by casting them as a system of
$\lfloor\alpha\rfloor$ ODEs, coupled with a fractional equation of order
$\beta\,:=\,\alpha-\lfloor\alpha\rfloor\in(0,1)$.

As anticipated in Section 3.2, a future direction of investigation will
concern the efficient numerical solution of the generated discrete problems,
which is crucial when using larger timesteps. Also the optimal choice of the
parameters for the methods will be further investigated, as well as the
possibility of adaptively defining the computational mesh. Further, the
parallel implementation of the methods could be also considered, following an
approach similar to [1].

## References

  * [1] P. Amodio, L. Brugnano. Parallel implementation of block boundary value methods for ODEs. J. Comput. Appl. Math. 78 (1997) 197–211. https://doi.org/10.1016/S0377-0427(96)00112-4
  * [2] P. Amodio, L. Brugnano, F. Iavernaro. Spectrally accurate solutions of nonlinear fractional initial value problems. AIP Conf. Proc. 2116 (2019) 140005. https://doi.org/10.1063/1.5114132
  * [3] P. Amodio, L. Brugnano, F. Iavernaro. Analysis of Spectral Hamiltonian Boundary Value Methods (SHBVMs) for the numerical solution of ODE problems. Numer. Algorithms 83 (2020) 1489–1508. https://doi.org/10.1007/s11075-019-00733-7
  * [4] P. Amodio, L. Brugnano, F. Iavernaro. A note on a stable algorithm for computing the fractional integrals of orthogonal polynomials. Applied Mathematics Letters 134 (2022) 108338. https://doi.org/10.1016/j.aml.2022.108338
  * [5] P. Amodio, L. Brugnano, F. Iavernaro. (Spectral) Chebyshev collocation methods for solving differential equations. Numer. Algoritms 93 (2023) 1613–1638. https://doi.org/10.1007/s11075-022-01482-w
  * [6] D.A. Benson, S.W. Wheatcraft, M.M. Meerschaert. Applications of a fractional advection-dispersion equation. Water Resources Research 36(6) (2000) 1403–1412.
  * [7] L. Brugnano, G. Frasca-Caccia, F. Iavernaro, V. Vespri. A new framework for polynomial approximation to differential equations. Adv. Comput. Math. 48 (2022) 76 https://doi.org/10.1007/s10444-022-09992-w
  * [8] L. Brugnano, F. Iavernaro. Line Integral Methods for Conservative Problems. Chapman et Hall/CRC, Boca Raton, FL, USA, 2016.
  * [9] L. Brugnano, F. Iavernaro. Line Integral Solution of Differential Problems. Axioms 7(2) (2018) 36. https://doi.org/10.3390/axioms7020036
  * [10] L. Brugnano, F. Iavernaro. A general framework for solving differential equations. Ann. Univ. Ferrara Sez. VII Sci. Mat. 68 (2022) 243–258. https://doi.org/10.1007/s11565-022-00409-6
  * [11] L. Brugnano, F. Iavernaro, D. Trigiante. A note on the efficient implementation of Hamiltonian BVMs. Jour. Comput. Appl. Math. 236 (2011) 375–383. https://doi.org/10.1016/j.cam.2011.07.022
  * [12] L. Brugnano, F. Iavernaro, D. Trigiante. A simple framework for the derivation and analysis of effective one-step methods for ODEs. Appl. Math. Comput. 218 (2012) 8475–8485. https://doi.org/10.1016/j.amc.2012.01.074
  * [13] L. Brugnano, J.I. Montijano, F. Iavernaro, L. Randéz. Spectrally accurate space-time solution of Hamiltonian PDEs. Numer. Algorithms 81 (2019) 1183–1202. https://doi.org/10.1007/s11075-018-0586-z
  * [14] L. Brugnano, J.I. Montijano, F. Iavernaro, L. Randéz. On the effectiveness of spectral methods for the numerical solution of multi-frequency highly-oscillatory Hamiltonian problems. Numer. Algorithms 81 (2019) 345–376. https://doi.org/10.1007/s11075-018-0552-9
  * [15] A. Bueno-Orovio, K. Burrage. Exact solutions to the fractional time-space Bloch-Torrey equation for magnetic resonance imaging. Commun. Nonlinear Sci. Numer. Simul. 52 (2017) 91–109.
  * [16] A. Bueno-Orovio, D. Kay, K. Burrage, Fourier-spectral methods for fractional in space reaction diffusion equations, BIT 54 (2014) 937–954.
  * [17] A. Bueno-Orovio, D. Kay, V. Grau, B. Rodriguez, K. Burrage. Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization. J. of Royal Soc. Interface 11(97) (2014) 20140352.
  * [18] N. Cusimano, A. Bueno-Orovio, I. Turner, K. Burrage. On the order of the fractional Laplacian in determining the spatio-temporal evolution of a space fractional model of cardiac electrophysiology. PLOS One 10(12) (2015) e0143938.
  * [19] K. Diethelm. The analysis of fractional differential equations. An application-oriented exposition using differential operators of Caputo type. Lecture Notes in Math., 2004. Springer-Verlag, Berlin, 2010.
  * [20] R. De Vore, L.R. Scott. Error bounds for Gaussian quadrature and weighted-$L^{1}$ polynomial approximation. SIAM J. Numer. Anal. 21, No. 2 (1984) 400–412. https://doi.org/10.1137/0721030
  * [21] K. Diethelm, N.J. Ford, A.D. Freed. Detailed error analysis for a fractional Adams method. Numer. Algorithms 36 (2004) 31–52. https://doi.org/10.1023/B:NUMA.0000027736.85078.be
  * [22] R. Garrappa, Numerical evaluation of two and three parameter Mittag-Leffler functions. SIAM J. Numer. Anal. 53, No. 3 (2015) 1350–1369. https://doi.org/10.1137/140971191
  * [23] R. Garrappa. Trapezoidal methods for fractional differential equations: Theoretical and computational aspects. Math. Comp. Simul. 110 (2015) 96–112. http://doi.org/10.1016/j.matcom.2013.09.012.
  * [24] R. Garrappa. Numerical solution of fractional differential equations: a survey and a software tutorial. Mathematics 6(2) (2018) 16. http://doi.org/10.3390/math6020016
  * [25] W. Gautschi. Orthogonal Polynomials Computation and Approximation. Oxford University Press, 2004.
  * [26] B.I. Henry, T.A.M. Langlands. Fractional cable models for spiny neuronal dendrites. Phys. Rev. Letts. 100(12) (2008) 128103.
  * [27] B.I. Henry, T. Langlands, S. Wearne. Turing pattern formation in fractional activator-inhibitor systems. Phys. Rev. E 72(2) (2005) 026101.
  * [28] M. Hori, I. Fukunaga, V. Masutani, T. Taoka, K. Kamagata, Y. Suzuki, et al. Visualising non Gaussian diffusion - clinical application of q-space imaging and diffusional kurtosis imaging of the brain, and spine. Magn. Reson. Med. Sc. 11 (2012) 221–233.
  * [29] V. Lakshmikantham, D. Trigiante. Theory of Difference Equations. Academic Press, Inc., Boston, MA, 1988.
  * [30] C. Li, Q. Yi, A. Chen. Finite difference methods with non-uniform meshes for nonlinear fractional differential equations. J. Comput. Phys. 316 (2016) 614–631.
  * [31] K. Lindenberg, S.B. Yuste. Properties of the reaction front in a reaction-subdiffusion process. Noise in Complex Systems and Stochastic Dynamics II 5471 (2004) 20–28.
  * [32] Ch. Lubich. Fractional Linear Multistep Methods for Abel-Volterra Integral Equations of the Second Kind. Math. Comp.  45, No. 172 (1985) 463–469. https://doi.org/10.1090/S0025-5718-1985-0804935-7
  * [33] Ch. Lubich. Discretized fractional calculus. SIAM J. Math. Anal. 17 (1986) 704–719.
  * [34] R. Magin, X. Feng, D. Baleanu. Solving the fractional order Bloch equation. Concepts in Magnetic Resonance, Part A, 34A (2009) 16–23.
  * [35] G. Mastroianni, G. Milovanovic. Interpolation processes. Basic Theory and Applications. Springer Monogr. Math., Springer-Verlag, Berlin 2008.
  * [36] E. Orsingher, L. Beghin. Fractional diffusion equations and processes with randomly varying time. Ann. Probab. 37(1) (2009) 206–249.
  * [37] I. Podlubny. Fractional differential equations. An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Academic Press, Inc., San Diego, CA, 1999.
  * [38] Z. Satmari. Iterative Bernstein splines technique applied to fractional order differential equations. Mathematical Foundations of Computing 6 (2023) 41–53. https://doi.org/10.3934/mfc.2021039
  * [39] A. Schädle, M. Lopez-Fernandez, Ch. Lubich. Fast and oblivious convolution quadrature. SIAM J. Sci. Comput. 28 (2006) 421–438.
  * [40] M. Stynes, E. O’Riordan, J.L. Gracia. Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM J. Numer. Anal. 55 (2017) 1057–1079.
  * [41] W. Themistoclakis. Some error bounds for Gauss-Jacobi quadrature rules. Appl. Numer. Math. 116 (2017) 286–293. https://doi.org/10.1016/j.apnum.2017.02.009
  * [42] F. Zeng, I. Turner, K. Burrage. A stable fast time-stepping method for fractional integral and derivative operators. J. Sci. Comput. 77 (2018) 283–307.
  * [43] F. Zeng, I. Turner, K. Burrage, G. Karniadiakis. A new class of semi-implicit methods with linear complexity for nonlinear fractional differential equations. SIAM J. Sci. Comput. 40, No. 5 (2018) A2986–A3011. https://doi.org/10.1137/18M1168169
  * [44] F. Zeng, Z. Zhang, G.E. Karniadiakis. Second order numerical methods for multi-term fractional differential equations. Comput. Methods Appl. Mech. Engrg. 327 (2017) 478–502.


LaTeX2ε SVMono Document Class Version 5.x

Reference Guide

for

Monographs

# *

$\copyright$ 2018, Springer Nature
All rights reserved.

## Contents

section.1.1 section.1.2 subsection.1.2.1 subsection.1.2.2 subsection.1.2.3
subsection.1.2.4 subsection.1.2.5 subsection.1.2.6 subsection.1.2.7
subsection.1.2.8 subsection.1.2.8 subsection.1.2.9 subsection.1.2.10
subsection.1.2.11 subsection.1.2.12 section*.4

##  1 Introduction

This reference guide gives a detailed description of the LaTeX2ε SVMono
document class Version 5.x and its special features designed to facilitate the
preparation of scientific books for Springer Nature. It always comes as part
of the SVMono tool package and should not be used on its own.

The components of the SVMono tool package are:

  * •

The Springer LaTeX class `SVMono.cls`, MakeIndex styles svind.ist, svindd.ist,
BibTeX styles spmpsci.bst, spphys.bst, spbasic.bst as well as the templates
with preset class options, packages and coding examples;

  * Tip: Copy all these files to your working directory, run LaTeX2ε, BibTeX and MakeIndex—as is applicable— and and produce your own example *.dvi file; rename the template files as you see fit and use them for your own input.

  * •

Author Instructions with style and coding instructions.

  * Tip: Follow these instructions to set up your files, to type in your text and to obtain a consistent formal style in line with the Springer Nature layout specifications; use these pages as checklists before you submit your manuscript data.

  * •

The Reference Guide describing SVMono features with regards to their
functionality.

  * Tip: Use it as a reference if you need to alter or enhance the default settings of the SVMono document class and/or the templates.

The documentation in the Springer SVMono tool package is not intended to be a
general introduction to LaTeX2ε or TeX. For this we refer you to [1–3].

Should we refer in this tool package to standard tools or packages that are
not installed on your system, please consult the Comprehensive TeX Archive
Network (CTAN) at [4–6].

SVMono was derived from the LaTeX2ε book.cls and article.cls.

The main differences from the standard document classes `article.cls` and
`book.cls` are the presence of

  * •

multiple class options,

  * •

a number of newly built-in environments for individual text structures like
theorems, exercises, lemmas, proofs, etc.,

  * •

enhanced environments for the layout of figures and captions, and

  * •

new declarations, commands and useful enhancements of standard environments to
facilitate your math and text input and to ensure their output is in line with
the Springer Nature layout standards.

Nevertheless, text, formulae, figures, and tables are typed using the standard
LaTeX2ε commands. The standard sectioning commands are also used.

Always give a `\label` where possible and use `\ref` for cross-referencing.
Such cross-references may then be converted to hyperlinks in any electronic
version of your book.

The `\cite` and `\bibitem` mechanism for bibliographic references is also
obligatory.

##  2 SVMono Class Features

###  2.1 Initializing the SVMono Class

To use the document class, enter

`\documentclass [`$\langle$options$\rangle$`] {svmono}`

at the beginning of your input.

###  2.2 SVMono Class Options

Choose from the following list of SVMono class options if you need to alter
the default layout settings of the SVMono document class. Please note that the
optional features should only be chosen if instructed so by the editor of your
book.

Page Style

[norunningheads]

default

twoside, single-spaced output, contributions starting always on a recto page

referee

produces double-spaced output for proofreading

footinfo

generates a footline with name, date, $\ldots$ at the bottom of each page

norunningheads

suppresses any headers and footers

N.B. If you want to use both options, you must type referee before footinfo.

Body Font Size

[11pt, 12pt]

default

10 pt

11pt, 12pt

are ignored

Language for Fixed LaTeX Texts

In the SVMono class we have changed a few standard LaTeX texts (e.g. Figure to
Fig. in figure captions) and assigned names to newly defined theorem-like
environments so that they conform with Springer Nature style requirements.

[francais]

default

English

deutsch

translates fixed LaTeX texts into their German equivalent

francais

same as above for French

Text Style

[graybox]

default

plain text

graybox

automatically activates the packages `color` and `framed`
and places a box with 15 percent gray shade in the background
of the text when you use the SVMono environment
`\begin{svgraybox}...\end{svgraybox}`, see Sects. 2.3, 2.4.

Equations Style

[vecarrow]

default

centered layout, vectors boldface (math style)

vecphys

produces boldface italic vectors (physics style)
when `\vec`-command is used

vecarrow

depicts vectors with an arrow above when `\vec`-commandis used

Numbering and Layout of Headings

[nochapnum]

default

all section headings down to subsubsection level are numbered, second and
subsequent lines in a multiline numbered heading are indented; Paragraph and
Subparagraph headings are displayed but not numbered; figures, tables and
equations are numbered chapterwise, individual theorem-like environments are
counted consecutively throughout the book.

nosecnum

suppresses any section numbering; figures, tables and equations are counted
chapterwise displaying the chapter counter, if applicable.

nochapnum

suppresses the chapter numbering only, subsequent section headings as well as
figures, tables and equations are numbered chapterwise but without chapter
counter.

nonum

suppresses any numbering of any headings; tables, figures, equations are
counted consecutively throughout the book.

$\backslash$chapter*

must not be used since all subsequent numbering will go bananas $\ldots$

††margin:  Warning !

Numbering of Figures, Tables and Equations

[numart]

default

chapter-wise numbering

numart

numbers figures, tables, equations consecutively (not chapterwise) throughout
the whole text, as in the standard article document class

Numbering and Counting of Built-in Theorem-Like Environments

[envcountresetchap]

default

each built-in theorem-like environment gets its own counter without any
chapter or section prefix and is counted consecutively throughout the book

envcountchap

Each built-in environment gets its own counter and is numbered chapterwise. To
be selected as default setting for a book with numbered chapters.

envcountsect

each built-in environment gets its own counter and is numbered sectionwise

envcountsame

all built-in environments follow a single counter without any chapter or
section prefix, and are counted consecutively throughout the book

envcountresetchap

each built-in environment gets its own counter without any chapter or section
prefix but with the counter reset for each chapter

envcountresetsect

each built-in environment gets its own counter without any chapter or section
prefix but with the counter reset for each section

N.B.1 When the option envcountsame is combined with the options envcount-
resetchap or envcountresetsect all predefined environments get the same
counter; but the counter is reset for each chapter or section.

N.B.2 When the option envcountsame is combined with the options envcountchap
or envcountsect all predefined environments get a common counter with a
chapter or section prefix; but the counter is reset for each chapter or
section.

N.B.3 We have designed a new easy-to-use mechanism to define your own
environments.

N.B.4 Be careful not to use layout options that contradict the parameter of
the selected environment option and vice versa. ††margin:  Warning !

Use the Springer class option

[nospthms]

nospthms

only if you want to suppress all defined theorem-like environments and use the
theorem environments of original LaTeX package or other theorem packages
instead. (Please check this with your editor.)

References

[chaprefs]

default

the list of references is set as an unnumbered chapter starting on a new recto
page, with automatically correct running heads and an entry in the table of
contents. The list itself is set in small print and numbered with ordinal
numbers.

sectrefs

sets the reference list as an unnumbered section, e.g. at the end of a chapter

natbib

sorts reference entries in the author-year system (make sure that you have the
natbib package by Patrick W. Daly installed. Otherwise it can be found at the
Comprehensive TeX Archive Network (CTAN…tex-
archive/macros/latex/contrib/supported/natbib/), see [4–6]

Use the Springer class option

[chaprefs]

oribibl

only if you want to set reference numbers in square brackets without automatic
TOC entry etc., as is the case in the original LaTeX bibliography environment.
But please note that most page layout features are nevertheless adjusted to
Springer Nature requirements. (Please check usage of this option with your
editor.)

###  2.3 Required and Recommended Packages

SVMono document class has been tested with a number of Standard LaTeX tools.
Below we list and comment on a selection of recommended packages for preparing
fully formatted book manuscripts for Springer Nature. If not installed on your
system, the source of all standard LaTeX tools and packages is the
Comprehensive TeX Archive Network (CTAN) at [4–6].

Font Selection

default |  Times font family as default text body font together with Helvetica clone as sans serif and Courier as typewriter font.
---|---
newtxtext.sty and newtxmath.sty |  Supports roman text font provided by a Times clone, sans serif based on a Helvetica clone, typewriter faces, plus math symbol fonts whose math italic letters are from a Times Italic clone

If the packages ‘newtxtext.sty and newtxmath.sty’ are not already installed
with your LaTeX they can be found at https://ctan.org/tex.archive/ fonts/newtx
at the Comprehensive TeX Archive Network (CTAN), see [4–6].

If Times Roman is not available on your system you may revert to CM fonts.
However, the SVMono layout requires font sizes which are not part of the
default set of the computer modern fonts.

[type1cm.sty]

type1cm.sty

The type1cm package enhances this default by enabling scalable versions of the
(Type 1) CM fonts. If not already installed with your LaTeX it can be found at
../tex-archive/macros/latex/contrib/type1cm/ at the Comprehensive TeX Archive
Network (CTAN), see [4–6].

Body Text

When you select the SVMono class option [graybox] the packages framed and
color are required, see Sect. 2.2

[framed.sty]

framed.sty

makes it possible that framed or shaded regions can break across pages.

color.sty

is part of the graphics bundle and makes it possible to selct the color and
define the percentage for the background of the box.

Equations

A useful package for subnumbering each line of an equation array can be found
at ../tex-archive/macros/latex/contrib/supported/subeqnarray/ at the
Comprehensive TeX Archive Network(CTAN), see [4–6].

[subeqnarray.sty]

subeqnarray.sty

defines the subeqnarray and subeqnarray* environments, which behave like the
equivalent eqnarray and eqnarray* environments, except that the individual
lines are numbered as 1a, 1b, 1c, etc.

Footnotes

[footmisc.sty]

footmisc.sty

used with style option [bottom] places all footnotes at the bottom of the page

Figures

[graphicx.sty]

graphicx.sty

tool for including graphics files (preferrably eps files)

References

[natbib.sty]

default

Reference lists are numbered with the references being cited in the text by
their reference number

natbib.sty

sorts reference entries in the author–year system (among other features). N.B.
This style must be installed when the class option natbib is used, see Sect.
2.2

cite.sty

generates compressed, sorted lists of numerical citations: e.g. [8,11–16];
preferred style for books published in a print version only

Index

[multicol.sty]

makeidx.sty

provides and interprets the command `\printindex` which “prints” the
externally generated index file *.ind.

multicol.sty

balances out multiple columns on the last page of your subject index, glossary
or the like

N.B. Use the MakeIndex program together with one of the following styles

[svindd.ist]

svind.ist

for English texts

svindd.ist

for German texts

to generate a subject index automatically in accordance with Springer Nature
layout requirements. For a detailed documentation of the program and its usage
we refer you to [1].

###  2.4 SVMono Commands and Environments in Text Mode

Use the environment syntax

`\begin{dedication}`
---
$\langle text\rangle$
`\end{dedication}`

to typeset a dedication or quotation at the very beginning of the in preferred
Springer layout.

Use the new commands  `\foreword` `\preface`

to typeset a Foreword or Preface with automatically generated runnings heads.

Use the new commands  `\extrachap{`$\langle heading\rangle$`}`
`\Extrachap{`$\langle heading\rangle$`}`

to typeset — in the front or back matter of the book—an extra unnumbered
chapter with your preferred heading and automatically generated runnings
heads.

`\Extrachap` furthermore generates an automated TOC entry.

Use the new command

`\partbacktext{`$\langle text\rangle$`}`

to typeset a text on the back side of a part title page.

Use the new command

`\chapsubtitle[`$\langle subtitle\rangle$`]`

to typeset a possible subtitle to your chapter title. Beware that this
subtitle is not tranferred automatically to the table of contents.

The command must be placed before the `\chapter` command.

Alternatively use the `\chapter`-command to typeset your subtitle together
with the chapter title and separate the two titles by a period or an en-dash.
††margin:  Alternative !

The command must be placed before the `\chapter` command.

Use the new command

`\chapauthor[`$\langle name\rangle$`]`

to typeset the author name(s) beneath your chapter title. Beware that the
author name(s) are not tranferred automatically to the table of contents.

The command must be placed before the `\chapter` command.

Alternatively, if the book has rather the character of a contributed volume as
opposed to a monograph you may want to use the SVMono package with features
that better suit the specific requirements. ††margin:  Alternative !

Use the new commands

`\chaptermark{}`
---
`\sectionmark{}`

to alter the text of the running heads.

Use the new command `\motto{`$\langle{\it text}\rangle$`}`

to include special text, e.g. mottos, slogans, between the chapter heading and
the actual content of the chapter in the preferred Springer layout.

The argument `{`$\langle{\it text}\rangle$`}` contains the text of your
inclusion. It may not contain any empty lines. To introduce vertical spaces
use `\\[height]`.

If needed, the you may indicate an alternative widths in the optional
argument.

N.B. The command must be placed before the relevant `heading`-command.

Use the new commands

`\abstract{`$\langle{\it text}\rangle$`}`
---
`\abstract*{`$\langle{\it text}\rangle$`}`

to typeset an abstract at the beginning of a chapter.

The text of `\abstract*` will not be depicted in the printed version of the
book, but will be used for compiling `html` abstracts for the online
publication of the individual chapters `www.SpringerLink.com`.

Please do not use the standard LaTeX environment ††margin:  Warning !!!

`\begin{abstract}...\end{abstract}` – it will be ignored when used with the
SVMono document class!

Use the new commands

`\runinhead[`$\langle{\it title}\rangle$`]`
---
`\subruninhead[`$\langle{\it title}\rangle$`]`

when you want to use unnumbered run-in headings to structure your text.

Use the new environment command  `\begin{svgraybox}` $\langle text\rangle$
`\end{svgraybox}`

to typeset complete paragraphs within a box showing a 15 percent gray shade.

N.B. Make sure to select the SVMono class option `[graybox]` in order to have
all the required style packages available, see Sects. 2.2, 2.3. ††margin:
Warning !

Use the new environment command

`\begin{petit}`
---
$\langle text\rangle$
`\end{petit}`

to typeset complete paragraphs in small print.

Use the enhanced environment command

`\begin{description}[`$\langle{\it largelabel}\rangle$`]`
---
`\item[`$\langle{\it label1}\rangle$`]` $\langle\textit{text1}\rangle$
`\item[`$\langle{\it label2}\rangle$`]` $\langle\textit{text2}\rangle$
`\end{description}`

for your individual itemized lists.

The new optional parameter `[`$\langle{\it largelabel}\rangle$`]` lets you
specify the largest item label to two levels to appear within the list. The
texts of all items are indented by the width of $\langle largelabel\rangle$
and the item labels are typeset flush left within this space. Note, the
optional parameter will work only two levels deep.

Use the commands

`\setitemindent{`$\langle{\it largelabel}\rangle$`}`
---
`\setitemitemindent{`$\langle{\it largelabel}\rangle$`}`

if you need to customize the indention of your “itemized” or “enumerated”
environments.

###  2.5 SVMono Commands in Math Mode

Use the new or enhanced symbol commands provided by the SVMono document class:

`\D` | upright d for differential d
---|---
`\I` | upright i for imaginary unit
`\E` | upright e for exponential function
`\tens` | depicts tensors as sans serif upright
`\vec` | depicts vectors as boldface characters instead of the arrow accent

N.B. By default the SVMono document class depicts Greek letters as italics
because they are mostly used to symbolize variables. However, when used as
operators, abbreviations, physical units, etc. they should be set upright.

All upright upper-case Greek letters have been defined in the SVMono document
class and are taken from the TeX alphabet.

Use the command prefix

`\var...`

with the upper-case name of the Greek letter to set it upright, e.g.
`\varDelta`.

Many upright lower-case Greek letters have been defined in the SVMono document
class and are taken from the PostScript Symbol font.

Use the command prefix

`\u...`

with the lower-case name of the Greek letter to set it upright, e.g. `\umu`.

If you need to define further commands use the syntax below as an example:

`\newcommand{\ualpha}{\allmodesymb{\greeksym}{a}}`

###  2.6 SVMono Theorem-Like Environments

For individual text structures such as theorems, definitions, and examples,
the SVMono document class provides a number of pre-defined environments which
conform with the specific Springer Nature layout requirements.

Use the environment command

`\begin{`$\langle{\it name~{}of~{}environment}\rangle$`}[`$\langle{\it
optional~{}material}\rangle$`]`
---
$\langle{\it text~{}for~{}that~{}environment}\rangle$
`\end{`$\langle{\it name~{}of~{}environment}\rangle$`}`

for the newly defined environments.

Unnumbered environments will be produced by

`claim` and `proof`.

Numbered environments will be produced by

case, conjecture, corollary, definition, example, exercise, lemma, note,
problem, property, proposition, question, remark, solution, and theorem.

The optional argument `[`$\langle{\it optional~{}material}\rangle$`]` lets you
specify additional text which will follow the environment caption and counter.

N.B. We have designed a new easy-to-use mechanism to define your own
environments.

Use the new symbol command `\qed`

to produce an empty square at the end of your proof.

In addition, use the new declaration

`\smartqed`

to move the position of the predefined qed symbol to be flush right (in text
mode). If you want to use this feature throughout your book the declaration
must be set in the preamble, otherwise it should be used individually in the
relevant environment, i.e. proof.

Example

`\begin{proof}`
`\smartqed`
`Text`
`\qed`
`\end{proof}`

Furthermore the functions of the standard `\newtheorem` command have been
enhanced to allow a more flexible font selection. All standard functions
though remain intact (e.g. adding an optional argument specifying additional
text after the environment counter).

Use the mechanism

`\spdefaulttheorem{`$\langle{\it env~{}name}\rangle$`}``{`$\langle
caption\rangle$`}``{`$\langle{\it cap~{}font}\rangle$`}``{`$\langle{\it
body~{}font}\rangle$`}`

to define an environment compliant with the selected class options (see Sect.
2.2) and designed as the predefined theorem-like environments.

The argument `{`$\langle{\it env~{}name}\rangle$`}` specifies the environment
name; `{`$\langle{\it caption}\rangle$`}` specifies the environment’s heading;
`{`$\langle{\it cap~{}font}\rangle$`}` and `{`$\langle{\it
body~{}font}\rangle$`}` specify the font shape of the caption and the text
body.

N.B. If you want to use optional arguments in your definition of a theoremlike
environment as done in the standard `\newtheorem` command, see below.

Use the mechanism `\spnewtheorem{`$\langle{\it
env~{}name}\rangle$`}``[`$\langle{\it
numbered~{}like}\rangle$`]``{`$\langle{\it caption}\rangle$`}``{`$\langle{\it
cap~{}font}\rangle$`}``{`$\langle{\it body~{}font}\rangle$`}`

to define an environment that shares its counter with another predefined
environment `[`$\langle{\it numbered~{}like}\rangle$`]`.

The optional argument `[`$\langle{\it numbered~{}like}\rangle$`]` specifies
the environment with which to share the counter.

N.B. If you select the class option “envcountsame” the only valid “numbered
like” argument is `[theorem]`.

Use the defined mechanism

`\spnewtheorem{`$\langle{\it env~{}name}\rangle$`}``{`$\langle{\it
caption}\rangle$`}``[`$\langle\langle{\it
within}\rangle\rangle$`]``{`$\langle{\it cap~{}font}\rangle$`}``{`$\langle{\it
body~{}font}\rangle$`}`

to define an environment whose counter is prefixed by either the chapter or
section number (use `[chapter]` or `[section]` for `[`$\langle{\it
within}\rangle$`]`).

Use the defined mechanism

`\spnewtheorem*{`$\langle{\it env~{}name}\rangle$`}``{`$\langle{\it
caption}\rangle$`}``{`$\langle{\it cap~{}font}\rangle$`}``{`$\langle{\it
body~{}font}\rangle$`}`

to define an unnumbered environment such as the pre-defined unnumbered
environments claim and proof.

Use the defined declaration `\nocaption`

in the argument `{`$\langle{\it caption}\rangle$`}` if you want to skip the
environment caption and use an environment counter only.

Use the defined environment  `\begin{theopargself}` `...` `\end{theopargself}`

as a wrapper to any theorem-like environment defined with the mechanism. It
suppresses the brackets of the optional argument specifying additional text
after the environment counter.

###  2.7 SVMono Commands for the Figure and Table Environments

Use the new declaration

`\sidecaption[`$\langle pos\rangle$`]`

to move the figure caption from beneath the figure (default) to the lower
lefthand side of the figure.

The optional parameter `[t]` moves the figure caption to the upper left-hand
side of the figure

N.B.1 (1) Make sure the declaration `\sidecaption` follows the
`\begin{figure}` command, and (2) remember to use the standard `\caption{}`
command for your caption text.

N.B.2 This declaration works only if the figure width is less than 7.8 cm. The
caption text will be set raggedright if the width of the caption is less than
3.4 cm.

Use the new declaration

`\samenumber`

within the figure and table environment – directly after the `\begin{`$\langle
environment\rangle$`}` command – to give the caption concerned the same
counter as its predecessor (useful for long tables or figures spanning more
than one page, see also the declaration `\subfigures` below.

To arrange multiple figures in a single environment use the newly defined
commands

`\leftfigure[`$\langle pos\rangle$`]` and `\rightfigure[`$\langle
pos\rangle$`]`

within a `{minipage}{\textwidth}` environment. To allow enough space between
two horizontally arranged figures use `\hspace{\fill}` to separate the
corresponding `\includegraphics{}` commands. The required space between
vertically arranged figures can be controlled with `\\[12pt]`, for example.

The default position of the figures within their predefined space is flush
left. The optional parameter `[c]` centers the figure, whereas `[r]` positions
it flush right – use the optional parameter only if you need to specify a
position other than flush left.

Use the newly defined commands

`\leftcaption{}` and `\rightcaption{}`

outside the `minipage` environment to put two figure captions next to each
other.

Use the newly defined command

`\twocaptionwidth{`$\langle width\rangle$`}``{`$\langle width\rangle$`}`

to overrule the default horizontal space of 5.4 cm provided for each of the
abovedescribed caption commands. The first argument corresponds to
`\leftcaption` and the latter to `\rightcaption`.

Use the new declaration

`\subfigures`

within the figure environment – directly after the `\begin{figure}` command –
to subnumber multiple captions alphabetically within a single figure-
environment.

N.B.: When used in combination with `\samenumber` the main counter remains the
same and the alphabetical subnumbering is continued. It works properly only
when you stick to the sequence `\samenumber\subfigures`.

If you do not include your figures as electronic files use the defined command

`\mpicplace{`$\langle width\rangle$`}{`$\langle height\rangle$`}`

to leave the desired amount of space for each figure. This command draws a
vertical line of the height you specified.

Use the new command `\svhline`

for setting in tables the horizontal line that separates the table header from
the table content.

###  2.8 SVMono Environments for Exercises, Problems and Solutions

Use the environment command

`\begin{prob}`
---
`\label{`$\langle problem{:}key\rangle$`}`
$\langle problem~{}text\rangle$
`\end{prob}`

to typeset and number each problem individually.

To facilitate the correct numbering of the solutions we have also defined a
solution environment, which takes the problem’s key, i.e. $\langle
problem{:}key\rangle$ (see above) as argument.

Use the environment syntax  `\begin{sol}{`$\langle problem{:}key\rangle$`}`
$\langle solution~{}text\rangle$ `\end{sol}`

to get the correct (i.e. problem =) solution number automatically.

###  2.9 SVMono Special Elements

Use the commands

`\begin{trailer}{`$\langle$Trailer Head$\rangle$`}`
---
`...`
`\end{trailer}`

If you want to emphasize complete paragraphs of texts in an `Trailer Head`.

Use the commands

`\begin{question}{`$\langle$Questions$\rangle$`}`
---
`...`
`\end{question}`

If you want to emphasize complete paragraphs of texts in an `Questions`.

Use the commands

`\begin{important}{`$\langle$Important$\rangle$`}`
---
`...`
`\end{important}`

If you want to emphasize complete paragraphs of texts in an `Important`.

Use the commands

`\begin{warning}{`$\langle$Attention$\rangle$`}`
---
`...`
`\end{warning}`

If you want to emphasize complete paragraphs of texts in an `Attention`.

Use the commands

`\begin{programcode}{`$\langle$Program Code$\rangle$`}`
---
`...`
`\end{programcode}`

If you want to emphasize complete paragraphs of texts in an `Program Code`.

Use the commands

`\begin{tips}{`$\langle$Tips$\rangle$`}`
---
`...`
`\end{tips}`

If you want to emphasize complete paragraphs of texts in an `Tips`.

Use the commands

`\begin{overview}{`$\langle$Overview$\rangle$`}`
---
`...`
`\end{overview}`

If you want to emphasize complete paragraphs of texts in an `Overview`.

Use the commands

`\begin{backgroundinformation}{`$\langle$Background Information$\rangle$`}`
---
`...`
`\end{backgroundinformation}`

If you want to emphasize complete paragraphs of texts in an `Background`
`Information`.

Use the commands

`\begin{legaltext}{`$\langle$Legal Text$\rangle$`}`
---
`...`
`\end{legaltext}`

If you want to emphasize complete paragraphs of texts in an `Legal Text`.

###  2.10 SVMono Commands for Styling References

The command

`\biblstarthook{`$\langle text\rangle$`}`

allows the inclusion of explanatory text between the bibliography heading and
the actual list of references. The command must be placed before the
`thebibliography` environment.

###  2.11 SVMono Commands for Styling the Index

The declaration

`\threecolindex`

sets the next index following the `\threecolindex` declaration in three
columns.

The Springer declaration

`\indexstarthook{`$\langle text\rangle$`}`

allows the inclusion of explanatory text between the index heading and the
actual list of references. The command must be placed before the `theindex`
environment.

###  2.12 SVMono Commands for Styling the Table of Contents

Use the command

`\setcounter{tocdepth}{number}`

to alter the numerical depth of your table of contents.

Use the macro

`\calctocindent`

to recalculate the horizontal spacing for large section numbers in the table
of contents set with the following variables:

`\tocchpnum` for the | chapter number
---|---
`\tocsecnum` | section number
`\tocsubsecnum` | subsection number
`\tocsubsubsecnum` | subsubsection
`\tocparanum` | paragraph number

Set the sizes of the variables concerned at the maximum numbering appearing in
the current document.

In the preamble set e.g:

`\settowidth{\tocchpnum}{36.\enspace}`
---
`\settowidth{\tocsecnum}{36.10\enspace}`
`\settowidth{\tocsubsecnum}{99.88.77}`
`\calctocindent`

## References

  1. [1]

L. Lamport: LaTeX: A Document Preparation System 2nd ed. (Addison-Wesley,
Reading, Ma 1994)

  2. [2]

M. Goossens, F. Mittelbach, A. Samarin: The LaTeX Companion (Addison-Wesley,
Reading, Ma 1994)

  3. [3]

D. E. Knuth: The TeXbook (Addison-Wesley, Reading, Ma 1986) revised to cover
TeX3 (1991)

  4. [4]

TeX Users Group (TUG), http://www.tug.org

  5. [5]

Deutschsprachige Anwendervereinigung TeX e.V. (DANTE), Heidelberg, Germany,
http://www.dante.de

  6. [6]

UK TeX Users’ Group (UK-TuG), http://uk.tug.org