∎

11institutetext: L. Brugnano 22institutetext: Dipartimento di Matematica e
Informatica “U. Dini”
Università di Firenze, Italy
22email<EMAIL_ADDRESS>
https://orcid.org/0000-0002-6290-4107 33institutetext: G. Gurioli
44institutetext: Dipartimento di Matematica e Informatica “U. Dini”
Università di Firenze, Italy
44email<EMAIL_ADDRESS>
https://orcid.org/0000-0003-0922-8119 55institutetext: F. Iavernaro
66institutetext: Dipartimento di Matematica
Università di Bari Aldo Moro, Italy
66email<EMAIL_ADDRESS>
https://orcid.org/0000-0002-9716-7370

# Numerical solution of FDE-IVPs by using Fractional HBVMs: the fhbvm code

Luigi Brugnano  Gianmarco Gurioli  Felice Iavernaro

(Received: date / Accepted: date)

###### Abstract

In this paper we describe the efficient numerical implementation of Fractional
HBVMs, a class of methods recently introduced for solving systems of
fractional differential equations. The reported arguments are implemented in
the Matlab${}^{\copyright}\,$code fhbvm, which is made available on the web.
An extensive experimentation of the code is reported, to give evidence of its
effectiveness.

###### Keywords:

fractional differential equations fractional integrals Caputo derivative
Jacobi polynomials Fractional Hamiltonian Boundary Value Methods FHBVMs

###### MSC:

34A08 65R20 65-04

##  1 Introduction

Fractional differential equations have become a common description tool across
a variety of applications (see, e.g., the classical references Di2010 ; Po1999
for an introduction). For this reason, their numerical solution has been the
subject of many researches (see, e.g. AMN2014 ; DFF2002 ; DiFoFr2005 ; LYC16 ;
Lu1985 ; SOG17 ; ZK2014 ), with the development of corresponding software
(see, e.g., CCP2023 ; Ga2015 ; Ga2018 ). In this context, the present
contribution is addressed for solving initial value problems for fractional
differential equations (FDE-IVPs) 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, for the sake of brevity, we have omitted the argument $t$ for $f$.
Here, for $\alpha\in(0,1)$, $y^{(\alpha)}(t)\equiv D^{\alpha}y(t)$ is the
Caputo fractional derivative:111As is usual, $\Gamma$ denotes the Euler gamma
function, such that, for $x>0$, $x\Gamma(x)=\Gamma(x+1)$.

$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)

Consequently, the solution of (1) can be formally written as:

$y(t)=y_{0}+I^{\alpha}f(y(t))\equiv
y_{0}+\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-x)^{\alpha-1}f(y(x))\mathrm{d}x,\qquad
t\in[0,T].$ (4)

The numerical method we shall consider, relies on Fractional HBVMs (FHBVMs), a
class of methods recently introduced in BBBI2023 , as an extension of
Hamiltonian Boundary Value Methods (HBVMs), special low-rank Runge-Kutta
methods originally devised for Hamiltonian problems (see, e.g., BI2016 ;
BI2018 ), and later extended along several directions (see, e.g., ABI2019 ;
ABI2022-1 ; ABI2023 ; BBBI2023 ; BFCIV2022 ; BI2022 ), including the numerical
solution of FDEs. A main feature of HBVMs is the fact that they can gain
spectrally accuracy, when approximating ODE-IVPs ABI2020 ; BMIR2019 ; BMR2019
, and such a feature has been recently extended to the FDE case BBBI2023 .

With this premise, the structure of the paper is as follows: in Section 2 we
recall the main facts about the numerical solution of FDE-IVPs proposed in
BBBI2023 ; in Section 3 we provide full implementation details for the
Matlab${}^{\copyright}\,$code fhbvm used in the numerical tests; in Section 4
we report an extensive experimentation of the code, providing some comparisons
with another existing one; at last, a few conclusions are given in Section 5.

##  2 Fractional HBVMs

To begin with, in order to obtain a piecewise approximation to the solution of
the problem, we consider a partition of the integration interval in the form:

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

where $h_{n}>0$, $n=1,\dots,N$, and such that

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

Further, by setting

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

the restriction of the solution of (1) on the interval $[t_{n-1},t_{n}]$, and
taking into account (4) and (5)–(6), one obtains, for
$t\,\equiv\,t_{n-1}+ch_{n}$, $c\in[0,1]$,

$\displaystyle
y(t)~{}\equiv~{}y_{n}(ch_{n})~{}=~{}y_{0}+\frac{1}{\Gamma(\alpha)}\int_{0}^{t_{n-1}+ch_{n}}(t_{n-1}+ch_{n}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
(8) $\displaystyle=$
$\displaystyle~{}y_{0}+\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}\int_{t_{\nu-1}}^{t_{\nu}}(t_{n-1}+ch_{n}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle\qquad+~{}\frac{1}{\Gamma(\alpha)}\int_{t_{n-1}}^{t_{n-1}+ch_{n}}(t_{n-1}+ch_{n}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle=$
$\displaystyle~{}y_{0}+\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}\int_{0}^{h_{\nu}}(t_{n-1}-t_{\nu-1}+ch_{n}-x)^{\alpha-1}f(y_{\nu}(x))\mathrm{d}x$
$\displaystyle\qquad+~{}\frac{1}{\Gamma(\alpha)}\int_{0}^{ch_{n}}(ch_{n}-x)^{\alpha-1}f(y_{n}(x))\mathrm{d}x$
$\displaystyle=$
$\displaystyle~{}y_{0}+\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{t_{n-1}-t_{\nu-1}}{h_{\nu}}+c\frac{h_{n}}{h_{\nu}}-\tau\right)^{\alpha-1}f(y_{\nu}(\tau
h_{\nu}))\mathrm{d}\tau$
$\displaystyle\qquad+~{}\frac{h_{n}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}f(y_{n}(\tau
h_{n}))\mathrm{d}\tau$ $\displaystyle=$
$\displaystyle~{}y_{0}+\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}h_{\nu}^{\alpha}\int_{0}^{1}\left(\sum_{i=\nu}^{n-1}\frac{h_{i}}{h_{\nu}}+c\frac{h_{n}}{h_{\nu}}-\tau\right)^{\alpha-1}f(y_{\nu}(\tau
h_{\nu}))\mathrm{d}\tau$
$\displaystyle\qquad+~{}\frac{h_{n}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}f(y_{n}(\tau
h_{n}))\mathrm{d}\tau,\qquad c\in[0,1].$

To make manageable the handling of the ratios $h_{i}/h_{\nu}$,
$i=\nu,\dots,n$, $\nu=1,\dots,n-1$, we shall hereafter consider the following
choices for the mesh (5)–(6):

\- Graded mesh.

In order to cope with possible singularities in the derivative of the vector
field at the origin, we consider the graded mesh

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

where $r>1$ and $h_{1}>0$ satisfy, by virtue of (5)–(6),

$h_{1}\frac{r^{N}-1}{r-1}=T.$ (10)

As a result, one obtains that (8) becomes:

$\displaystyle y_{n}(ch_{n})~{}=~{}y_{0}~{}+$ (11)
$\displaystyle\frac{h_{1}^{\alpha}}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}r^{\alpha(\nu-1)}\int_{0}^{1}\left(\frac{r^{n-\nu}-1}{r-1}+cr^{n-\nu}-\tau\right)^{\alpha-1}f(y_{\nu}(\tau
r^{\nu-1}h_{1}))\mathrm{d}\tau$
$\displaystyle\qquad+~{}\frac{h_{1}^{\alpha}r^{\alpha(n-1)}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}f(y_{n}(\tau
r^{n-1}h_{1}))\mathrm{d}\tau$ $\displaystyle=:$
$\displaystyle~{}\phi_{n-1}^{\alpha}(c;h_{1},r,y)\,+\,\frac{h_{1}^{\alpha}r^{\alpha(n-1)}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}f(y_{n}(\tau
r^{n-1}h_{1}))\mathrm{d}\tau,$ $\displaystyle\qquad\qquad c\in[0,1].$

\- Uniform mesh.

This case is equivalent to allowing $r=1$ in (9). Consequently, from (5)–(6)
we derive

$h_{n}\equiv h_{1}:=\frac{T}{N},\quad n=1,\dots,N,\quad\Rightarrow\quad
t_{n}=nh_{1},\quad n=0,\dots,N.$ (12)

As a result, (8) reads:

$\displaystyle y_{n}(ch_{n})~{}\equiv~{}y_{n}(ch_{1})$ (13) $\displaystyle=$
$\displaystyle
y_{0}+\frac{h_{1}^{\alpha}}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}\int_{0}^{1}\left(n-\nu+c-\tau\right)^{\alpha-1}f(y_{\nu}(\tau
h_{1}))\mathrm{d}\tau$
$\displaystyle\qquad+~{}\frac{h_{1}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}f(y_{n}(\tau
h_{1}))\mathrm{d}\tau,$ $\displaystyle\equiv$
$\displaystyle~{}\phi_{n-1}^{\alpha}(c;h_{1},1,y)\,+\,\frac{h_{1}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}f(y_{n}(\tau
h_{1}))\mathrm{d}\tau$ $\displaystyle\qquad\qquad c\in[0,1].$

###### Remark 1

As is clear, in order to obtain an accurate approximation of the solution, it
is important to establish which kind of mesh (graded or uniform) is
appropriate. Besides this, also a proper choice of the parameters $h_{1}$,
$r$, and $N$ in (9)–(10) is crucial. Both aspects will be studied in Section
3.1.

###  2.1 Quasi-polynomial approximation

We now discuss a piecewise quasi-polynomial approximation to the solution of
(1),

$\sigma(t)\approx y(t),\qquad t\in[0,T],$

such that

$\sigma_{n}(ch_{n})\,:=\,\sigma(t_{n-1}+ch_{n}),\qquad c\in[0,1],\qquad
n=1,\dots,N,$ (14)

is the approximation to $y_{n}(ch_{n})$, defined in (7). According to BBBI2023
, such approximation will be derived through the following steps:

  1. 1.

expansion of the vector field, in each sub-interval $[t_{n-1},t_{n}]$,
$n=1,\dots,N$, (recall (5)–(6)) along a suitable orthonormal polynomial basis;

  2. 2.

truncation of the infinite expansion, in order to obtain a local polynomial
approximation.

Let us first consider the expansion of the vector field along the orthonormal
polynomial basis, w.r.t. the weight function

$\omega(x)=\alpha(1-x)^{\alpha-1},\qquad
s.t.\qquad\int_{0}^{1}\omega(x)\,\mathrm{d}x=1,$ (15)

resulting into a scaled and shifted family of Jacobi polynomials:222Here,
$\overline{P}_{j}^{\,(a,b)}(x)$ denotes the $j$-th Jacobi polynomial with
parameters $a$ and $b$, in $[-1,1]$.

$P_{j}(x):=\sqrt{\frac{2j+\alpha}{\alpha}}\,\overline{P}_{j}^{\,(\alpha-1,0)}(2x-1),\qquad
x\in[0,1],\qquad j=0,1,\dots,$

such that

$\int_{0}^{1}\omega(x)P_{i}(x)P_{j}(x)\mathrm{d}x=\delta_{ij},\qquad
i,j=0,1,\dots.$ (16)

In so doing, for $n=1,\dots,N$, one obtains:

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

with (see (15))

$\gamma_{j}(y_{n})=\int_{0}^{1}\omega(\tau)P_{j}(\tau)f(y_{n}(\tau
h_{n}))\mathrm{d}\tau,\qquad j=0,1,\dots.$ (18)

Consequently, the FDE (1), can be rewritten as:

$y_{n}^{(\alpha)}(ch_{n})=\sum_{j\geq 0}P_{j}(c)\gamma_{j}(y_{n}),\quad
c\in[0,1],\quad n=1,\dots,N,\qquad y_{1}(0)=y_{0}.$ (19)

The approximation is derived by truncating the infinite series in (17) to a
finite sum with $s$ terms, thus replacing (19) with a series of local
problems, whose vector field is a polynomial of degree $s$:

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

with $\gamma_{j}(\sigma_{n})$ defined similarly as in (18):

$\gamma_{j}(\sigma_{n})=\int_{0}^{1}\omega(\tau)P_{j}(\tau)f(\sigma_{n}(\tau
h_{n}))\mathrm{d}\tau,\qquad j=0,\dots,s-1.$ (21)

As a consequence, (11) will be approximated as:

$\displaystyle\sigma_{n}(ch_{n})~{}=$ (22) $\displaystyle
y_{0}\,+\,\frac{h_{1}^{\alpha}}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}r^{\alpha(\nu-1)}\int_{0}^{1}\left(\frac{r^{n-\nu}-1}{r-1}+cr^{n-\nu}-\tau\right)^{\alpha-1}\,\sum_{j=0}^{s-1}P_{j}(\tau)\gamma_{j}(\sigma_{\nu})\mathrm{d}\tau$
$\displaystyle\qquad+~{}\frac{h_{1}^{\alpha}r^{\alpha(n-1)}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}\sum_{j=0}^{s-1}P_{j}(\tau)\gamma_{j}(\sigma_{n})\mathrm{d}\tau~{}=~{}y_{0}\,+$
$\displaystyle\,h_{1}^{\alpha}\sum_{\nu=1}^{n-1}r^{\alpha(\nu-1)}\sum_{j=0}^{s-1}\left[\frac{1}{\Gamma(\alpha)}\int_{0}^{1}\left(\frac{r^{n-\nu}-1}{r-1}+cr^{n-\nu}-\tau\right)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau\right]\gamma_{j}(\sigma_{\nu})$
$\displaystyle\qquad+~{}h_{1}^{\alpha}r^{\alpha(n-1)}\sum_{j=0}^{s-1}\left[\frac{1}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau\right]\gamma_{j}(\sigma_{n})$
$\displaystyle=:~{}y_{0}\,+\,h_{1}^{\alpha}\sum_{\nu=1}^{n-1}r^{\alpha(\nu-1)}\sum_{j=0}^{s-1}J_{j}^{\alpha}\left(\frac{r^{n-\nu}-1}{r-1}+cr^{n-\nu}\right)\gamma_{j}(\sigma_{\nu})$
$\displaystyle\qquad+~{}h_{1}^{\alpha}r^{\alpha(n-1)}\sum_{j=0}^{s-1}I^{\alpha}P_{j}(c)\,\gamma_{j}(\sigma_{n})$
$\displaystyle=:~{}\phi_{n-1}^{\alpha,s}(c;h_{1},r,\sigma)~{}+~{}h_{1}^{\alpha}r^{\alpha(n-1)}\sum_{j=0}^{s-1}I^{\alpha}P_{j}(c)\,\gamma_{j}(\sigma_{n}),\qquad
c\in[0,1],$

where (see (3))

$I^{\alpha}P_{j}(c)=\frac{1}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau,\qquad
j=0,\dots,s-1,$ (23)

is the Riemann-Liouville integral of $P_{j}(c)$, and, for $x\geq 1$:

$J_{j}^{\alpha}(x)=\frac{1}{\Gamma(\alpha)}\int_{0}^{1}(x-\tau)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau,\qquad
j=0,\dots,s-1.$ (24)

The efficient numerical evaluation of the integrals (23) and (24) will be
explained in detail in Section 3.2.

###### Remark 2

We observe that, for $c=x=1$, by virtue of (15)–(16) one has:

$I^{\alpha}P_{j}(1)=J_{j}^{\alpha}(1)=\frac{\delta_{j0}}{\Gamma(\alpha+1)},\qquad
j=0,\dots,s-1.$ (25)

Similarly, when using a uniform mesh, by means of similar steps as above, (13)
turns out to be approximated by:

$\displaystyle\sigma_{n}(ch_{n})~{}=$ (26) $\displaystyle
y_{0}\,+\,h_{1}^{\alpha}\sum_{\nu=1}^{n-1}\sum_{j=0}^{s-1}J_{j}^{\alpha}\left(n-\nu+c\right)\gamma_{j}(\sigma_{\nu})\,+\,h_{1}^{\alpha}\sum_{j=0}^{s-1}I^{\alpha}P_{j}(c)\,\gamma_{j}(\sigma_{n})$
$\displaystyle\equiv$
$\displaystyle~{}\phi_{n-1}^{\alpha,s}(c;h_{1},1,\sigma)~{}+~{}h_{1}^{\alpha}\sum_{j=0}^{s-1}I^{\alpha}P_{j}(c)\,\gamma_{j}(\sigma_{n}),\qquad
c\in[0,1].\qquad$

In both cases, the approximation at $t_{n}$ is obtained, by setting $c=1$ and
taking into account (25), as:

$\displaystyle\sigma_{n}(h_{n})$ $\displaystyle=$
$\displaystyle\phi_{n-1}^{\alpha,s}(1;h_{1},r,\sigma)\,+\,\frac{h_{1}^{\alpha}r^{\alpha(n-1)}}{\Gamma(\alpha+1)}\gamma_{0}(\sigma_{n})$
(27) $\displaystyle\equiv$
$\displaystyle\phi_{n-1}^{\alpha,s}(1;h_{1},r,\sigma)\,+\,\frac{h_{n}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}(\sigma_{n}),$

which formally holds also in the case of a uniform mesh (see (12)) by setting
$r=1$.

###  2.2 The fully discrete method

Quoting Dahlquist and Björk DB2008 , “as is well known, even many relatively
simple integrals cannot be expressed in finite terms of elementary functions,
and thus must be evaluated by numerical methods”. In our framework, this
obvious statement means that, in order to obtain a numerical method, at step
$n$ the Fourier coefficients $\gamma_{j}(\sigma_{n})$ in (21) need to be
approximated by means of a suitable quadrature formula. Fortunately enough,
this can be done up to machine precision by using a Gauss-Jacobi formula of
order $2k$ based at the zeros of $P_{k}(c)$, $c_{1},\dots,c_{k}$, with
corresponding weights (see (15))

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

by choosing a value of $k$, $k\geq s$, large enough. In other words,

$\gamma_{j}^{n}\,:=\,\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f(\sigma_{n}(c_{i}h_{n}))\,\doteq\,\gamma_{j}(\sigma_{n}),\qquad
j=0,\dots,s-1,$ (28)

where $\,\doteq\,$ means equal within machine precision. Because of this, and
for sake of brevity, we shall continue using $\sigma_{n}$ to denote the fully
discrete approximation (compare with (22), or with (26), in case
$r=1$):333Hereafter, we shall use $h_{n}$ in place of $h_{1}r^{n-1}$.

$\sigma_{n}(ch_{n})=\phi_{n-1}^{\alpha,s}(c;h_{1},r,\sigma)+h_{n}^{\alpha}\sum_{j=0}^{s-1}I^{\alpha}P_{j}(c)\gamma_{j}^{n},\qquad
c\in[0,1].$ (29)

As is clear, the coefficientes $\gamma_{j}^{\nu}$, $j=0,\dots,s-1$,
$\nu=1,\dots,n-1$, needed to evaluate
$\phi_{n-1}^{\alpha,s}(c;h_{1},r,\sigma)$, have been already computed at the
previous time-steps.

We observe that, from (28), in order to compute the (discrete) Fourier
coefficients, it is enough evaluating (29) only at the quadrature abscissae
$c_{1},\dots,c_{k}$. In so doing, by combining (28) and (29), one obtains a
discrete problem in the form:

$\displaystyle\gamma_{j}^{n}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f\left(\phi_{n-1}^{\alpha,s}(c_{i};h_{1},r,\sigma)+h_{n}^{\alpha}\sum_{\ell=0}^{s-1}I^{\alpha}P_{\ell}(c_{i})\gamma_{\ell}^{n}\right),$
(30) $\displaystyle j~{}=~{}0,\dots,s-1.$

Once it has been solved, according to (27), the approximation of $y(t_{n})$ is
given by:

$\bar{y}_{n}\,:=\,\sigma_{n}(h_{n})\,=\,\phi_{n-1}^{\alpha,s}(1;h_{1},r,\sigma)\,+\,\frac{h_{n}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}^{n}.$
(31)

It is worth noticing that the discrete problem (30) can be cast in vector
form, by introducing the (block) vectors

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

and 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_{0}(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)\in\mathbb{R}^{k\times k},$

as:

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

with the obvious notation for the function $f$, evaluated in a vector of
(block) dimension $k$, of denoting the (block) vector of dimension $k$
containing the function $f$ evaluated in each (block) entry. As observed in
BBBI2023 , the (block) vector

$Y\,:=\,{\bm{\phi}}_{n-1}^{\alpha,s}+h_{n}^{\alpha}{\cal
I}_{s}^{\alpha}\otimes I_{m}{\bm{\gamma}}^{n}\in\mathbb{R}^{km},$ (33)

appearing at the r.h.s. in (32) as argument of $f$, satisfies the equation

$Y~{}=~{}{\bm{\phi}}_{n-1}^{\alpha,s}+h_{n}^{\alpha}{\cal I}_{s}^{\alpha}{\cal
P}_{s}^{\top}\Omega\otimes I_{m}\,f(Y),$ (34)

obtained combining (32) and (33). Consequently, it can be regarded as the
stage vector of a Runge-Kutta type method with Butcher tableau

$\begin{array}[]{c|c}{\bm{c}}&{\cal I}^{\alpha}_{s}{\cal
P}_{s}^{\top}\Omega\\\ \hline\cr\\\\[-8.53581pt]
&{\bm{b}}^{\top}\end{array},\qquad{\bm{b}}=\left(\begin{array}[]{ccc}b_{1},&\dots,&b_{k}\end{array}\right)^{\top},\qquad{\bm{c}}=\left(\begin{array}[]{ccc}c_{1},&\dots,&c_{k}\end{array}\right)^{\top}.$
(35)

###### Remark 3

Though the two formulations (32) and (34) are equivalent each other,
nevertheless, the former has (block) dimension $s$ independently of the
considered value of $k$ ($k\geq s$), which is the (block) dimension of the
latter. Consequently, in the practical implementation of the method, the
discrete problem (32) is the one to be preferred, since it is independent of
the number of stages.

When $\alpha=1$, the Runge-Kutta method (35) reduces to a HBVM$(k,s)$ method
ABI2020 ; ABI2022 ; BFCIV2022 ; BI2016 ; BI2018 ; BI2022 ; BMIR2019 ; BMR2019
. Consequently, we give the following definition BBBI2023 .

###### Definition 1

The method defined by (35) (i.e., (31)–(32)) is called fractional HBVM with
parameters $(k,s)$. In short, FHBVM$(k,s)$.

The efficient numerical solution of the discrete problem (32) will be
considered in Section 3.4.

##  3 Implementation issues

In this section we report all the implementation details used in the
Matlab${}^{\copyright}\,$code fhbvm, which will be used in the numerical tests
reported in Section 4.

###  3.1 Graded or uniform mesh?

The first relevant problem to face is the choice between a graded or a uniform
mesh and, in the former case, also choosing appropriate values for the
parameters $r$, $N$, and $h_{1}$ in (9)–(10). We start considering a proper
choice of this latter parameter, i.e., $h_{1}$ which, in turn, will allow us
to choose the type of mesh, too. For this purpose, the user is required to
provide, in input, a convenient integer value $M>1$, such that, if a uniform
mesh is appropriate, then the stepsize is given by:

$h=\frac{T}{M}.$ (36)

It is to be noted that, since the code is using a method with spectral
accuracy in time, the value of $M$ should be as small as possible. That said,
we set $h_{1}^{0}=h$ and apply the code on the interval $[0,h_{1}^{0}]$ both
in a single step, and by using a graded mesh of 2 sub-intervals defined by a
value of $r:=\hat{r}\equiv 3$. As a result, the two sub-intervals, obtained by
solving (10) for $h_{1}$, with

$r=3,\qquad N=2,\qquad T=h_{1}^{0},$

are 444It is to be noticed that the division by 4 is done without introducing
round-off errors.

$[0,h_{1}^{0}/4],\qquad[h_{1}^{0}/4,h_{1}^{0}].$

In so doing, we obtain two approximations to $y(h_{1}^{0})$, say $y_{1}$ and
$y_{2}$. According to the analysis in BBBI2023 , if 555Here, $./$ means the
componentwise division between two vectors, and $|y_{2}|$ denotes the vector
with the absolute values of the entries of $y_{2}$.

$\|(y_{1}-y_{2})./(1+|y_{2}|)\|_{\infty}\leq\mbox{\em tol},$ (37)

with tol a suitably small tolerance,666In view of the spectral accuracy of the
method, this tolerance is only slightly larger than the machine epsilon. then
this means the stepsize $h_{1}:=h_{1}^{0}$ is appropriate. Moreover, in such a
case, a uniform mesh with $N\equiv M$ turns out to be appropriate, too.

Conversely, the procedure is repeated on the sub-interval
$[0,h_{1}^{1}]\equiv[0,h_{1}^{0}/4]$, and so forth, until a convenient initial
stepsize $h_{1}$ is obtained. This process can be repeated up to a suitable
maximum number of times that, considering that at each iteration the current
guess of $h_{1}$ is divided by 4, allows for a substantial reduction of the
initial value (36). Clearly, in such a case, a graded mesh is needed. At the
end of this procedure, we have then chosen the initial stepsize $h_{1}$ which,
if the procedure ends at the $\ell$-th iteration, for a convenient $\ell>1$,
is given by:

$h_{1}=4^{1-\ell}h\equiv 4^{1-\ell}\frac{T}{M}.$ (38)

We need now to appropriately choose the parameters $r$ and $N$ in (9)–(10). To
simplify this choice, we require that the last stepsize in the mesh be equal
to (36). Consequently, we would ideally satisfy the following requirements:

$h_{1}\frac{r^{N}-1}{r-1}=T,\qquad h_{1}r^{N-1}=\frac{T}{M},$

which, by virtue of (38), become

$4^{1-\ell}\frac{r^{N}-1}{r-1}=M,\qquad 4^{1-\ell}r^{N-1}=1.$

Combining the two equations then gives:

$r=\frac{M-4^{1-\ell}}{M-1}>1,\qquad N=1+\log_{r}(4^{\ell-1}).$ (39)

As is clear, this value of $N$ is not an integer, in general, so that we shall
choose, instead:

$N=\lceil 1+\log_{r}(4^{\ell-1})\rceil.$ (40)

###### Remark 4

From (40), and considering that $Nh_{1}<T$ (conversely, a uniform mesh could
have been used) and, by virtue of (38), one has:

$2\leq N<4^{\ell-1}M.$ (41)

As is clear, using the value of $N$ in (40) in place of that in (39) implies
that we need to recompute $r$, so that the requirement (now $h_{1}$, $N$, and
$T$ are given)

$h_{1}\frac{r^{N}-1}{r-1}=T$ (42)

is again fulfilled. Equation (42) can be rewritten as

$r=\left[1+(r-1)\beta\right]^{1/N}=:\psi(r),\qquad\beta:=\frac{T}{h_{1}}>1,$
(43)

thus inducing the iterative procedure

$r_{0}>1\mbox{\quad(given)},\qquad r_{i+1}=\psi(r_{i}),\qquad i=0,1,\dots.$
(44)

As a convenient choice for $r_{0}$, one can use the guess for $r$ defined in
(39). The following result can be proved.

###### Theorem 3.1

There exists a unique $\bar{r}>1$ satisfying $\bar{r}=\psi(\bar{r})$, and the
iteration (44) globally converges to this value over the interval
$(1,+\infty)$.

###### Proof

A direct computation shows that:

$\displaystyle\psi(1)=1,\qquad\psi^{\prime}(r)>0,~{}\textrm{for~{}}r>1,$ (45)
$\displaystyle\psi^{\prime}(1)=\frac{\beta}{N}>1,~{}\lim_{r\rightarrow+\infty}\psi^{\prime}(r)=0,\qquad\psi^{\prime\prime}(r)<0,~{}\textrm{for~{}}r>1.$
(46)

From (45) we deduce that the (positive) mapping $\psi(r)$ is strictly
increasing and admits $(1,+\infty)$ as invariant set, since $r>1$ implies
$\psi(r)>\psi(1)=1$.

From (46) we additionally deduce that $\psi(r)$ is concave and $\psi(r)>r$ for
$r\in(1,1+\varepsilon)$, with $\varepsilon>0$ sufficiently small. These
properties and the fact that $\psi^{\prime}(r)\rightarrow 0$, for
$r\rightarrow+\infty$, imply that:

  * -

the equation $r=\psi(r)$ admits a unique solution $\bar{r}>1$ and
$\psi^{\prime}(\bar{r})<1$;

  * -

$\psi\big{(}(1,\bar{r})\big{)}\subset(1,\bar{r})$ and
$\psi\big{(}(\bar{r},+\infty)\big{)}\subset(\bar{r},+\infty)$ (since $\psi$ is
increasing);

  * -

$\psi(r)>r$ for $1<r<\bar{r}$ and $\psi(r)<r$ for $r>\bar{r}$.

From the two latter properties we conclude that, for any $r_{0}>1$, the
sequence generated by (44) converges monotonically to $\bar{r}$. ∎

Consequently, we have derived the parameters $h_{1}$ in (38), $N$ in (40), and
$r$ satisfying (42) of the graded mesh (9)–(10), the latter one obtained by a
few iterations of (44).

###### Remark 5

In the actual implementation of the code, we allow the use of a uniform mesh
also when $\ell=2$ steps of the previous procedure are required, provided that
$M$ has a moderate value (say, $M\leq 5$). Consequently, for the final mesh,
$h_{1}=h/4$, $r=1$, and $N=4M$.

###  3.2 Approximating the fractional integrals

The practical implementation of the method requires the evaluation of the
following integrals (recall (30) and the definitions (23)–(24)):

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

and

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

in case a graded mesh is used, or

$J_{j}^{\alpha}(\nu+c_{i}),\qquad i=1,\dots,k,\quad
j=0,\dots,s-1,\quad\nu=1\dots,N-1,$ (49)

in case a uniform mesh is used.

It must be emphasized that all such integrals ((47) and (48), or (47) and
(49)), can be pre-computed once for all, for later use. For their computation,
in the first version of the software, we adapted an algorithm based on ABI2022
which, however, required a quadruple precision, for the considered values of
$k$ and $s$. In this respect, the use of the standard vpa of
Matlab${}^{\copyright}\,$, which is based on a symbolic computation, turned
out to be too slow. For this reason, hereafter we describe two new algorithms
for computing the above integrals, which result to be quite satisfactory, when
using the standard double precision IEEE. Needless to say that, since vpa is
no more required, they turn out to be much faster than the previous ones.

Let us describe, at first, the computation of (47). One has, by considering
that $k\geq s$ and $c_{i}\in(0,1)$:

$\displaystyle
I^{\alpha}P_{j}(c_{i})~{}=~{}\frac{1}{\Gamma(\alpha)}\int_{0}^{c_{i}}(c_{i}-\tau)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau$
$\displaystyle=$
$\displaystyle\,\frac{c_{i}^{\alpha-1}}{\Gamma(\alpha)}\int_{0}^{c_{i}}\left(1-\frac{\tau}{c_{i}}\right)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau\,=\,\frac{c_{i}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{1}\left(1-\xi\right)^{\alpha-1}P_{j}(\xi
c_{i})\mathrm{d}\xi$ $\displaystyle=$
$\displaystyle\,\frac{c_{i}^{\alpha}}{\Gamma(\alpha+1)}\int_{0}^{1}\alpha\left(1-\xi\right)^{\alpha-1}P_{j}(\xi
c_{i})\mathrm{d}\xi\,\equiv\,\frac{c_{i}^{\alpha}}{\Gamma(\alpha+1)}\sum_{\ell=1}^{k}b_{\ell}P_{j}(c_{i}c_{\ell}),$
$\displaystyle\qquad\qquad j=0,\dots,s-1,$

where the last equality holds because the Jacobi quadrature formula has order
$2k$. Clearly, for each $i=1,\dots,k$, all the above integrals can be computed
in vector form in “one shot”, by using the usual three-term recurrence to
compute the Jacobi polynomials.

Concerning the integrals (48)-(49), let us now consider the evaluation of a
generic $J_{j}^{\alpha}(x)$, $j=0,\dots,s-1$, for $x>1$. In this respect,
there is numerical evidence that, for $x\geq 1.1$, a high-order Gauss-Legendre
formula is able to approximate the required integral to full machine
precision. Since we will use a value of $s=20$, we consider, for this purpose,
a Gauss-Legendre formula of order $60$, which turns out to be fully accurate.
Instead, for $x\in(1,1.1)$, one has:

$\displaystyle
J_{j}^{\alpha}(x)~{}=~{}\frac{1}{\Gamma(\alpha)}\int_{0}^{1}(x-\tau)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau$
$\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha)}\left[\int_{0}^{x}(x-\tau)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau~{}-~{}\int_{1}^{x}(x-\tau)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau\right]$
$\displaystyle=$
$\displaystyle\frac{x^{\alpha-1}}{\Gamma(\alpha)}\int_{0}^{x}\left(1-\frac{\tau}{x}\right)^{\alpha-1}P_{j}(\tau)\mathrm{d}\tau$
$\displaystyle\qquad~{}-~{}\frac{(x-1)^{\alpha-1}}{\Gamma(\alpha)}\int_{0}^{x-1}\left(1-\frac{c}{x-1}\right)^{\alpha-1}P_{j}(1+c)\mathrm{d}c$
$\displaystyle=$
$\displaystyle\frac{x^{\alpha}}{\Gamma(\alpha)}\int_{0}^{1}(1-\xi)^{\alpha-1}P_{j}(\xi
x)\mathrm{d}\xi-\frac{(x-1)^{\alpha}}{\Gamma(\alpha)}\int_{0}^{1}(1-\xi)^{\alpha-1}P_{j}(1+\xi(x-1))\mathrm{d}\xi$
$\displaystyle=$
$\displaystyle\frac{x^{\alpha}}{\Gamma(\alpha+1)}\int_{0}^{1}\alpha(1-\xi)^{\alpha-1}P_{j}(\xi
x)\mathrm{d}\xi$
$\displaystyle\qquad~{}-~{}\frac{(x-1)^{\alpha}}{\Gamma(\alpha+1)}\int_{0}^{1}\alpha(1-\xi)^{\alpha-1}P_{j}(1+\xi(x-1))\mathrm{d}\xi$
$\displaystyle\equiv$
$\displaystyle\frac{x^{\alpha}}{\Gamma(\alpha+1)}\sum_{\ell=1}^{k}b_{\ell}P_{j}(c_{\ell}x)~{}-~{}\frac{(x-1)^{\alpha}}{\Gamma(\alpha+1)}\sum_{\ell=1}^{k}b_{\ell}P_{j}(1+c_{\ell}(x-1)),$
$\displaystyle\qquad\qquad j=0,\dots,s-1,$

again, due to the fact that the quadrature is exact for polynomials of degree
at most $s-1$. Also in this case, for each fixed $x>1$, all the integrals can
be computed in “one shot” by using the three-term recurrence of the Jacobi
polynomials.

We observe that the previous expression is exact also for $x=1$, since
$J_{j}^{\alpha}(1)=\delta_{j0}/\Gamma(\alpha+1)$.

###  3.3 Error estimation

An estimate of the global error can be derived by computing the solution on a
doubled mesh. In other words, if (see (31))

$\bar{y}_{n}\,\approx\,y(t_{n}),\qquad n=0,\dots,N,$

with $t_{n}$ as in (5)–(6), are the obtained approximations, then

$e_{n}\,:=\,y(t_{n})-\bar{y}_{n}\,\approx\,\hat{y}_{2n}-\bar{y}_{n},\qquad
n=0,\dots,N,$

where $\hat{y}_{n}\approx y(\hat{t}_{n})$, $n=0,\dots,2N$, is the solution
computed on a doubled mesh.

When a uniform mesh (12) is used, the doubled mesh is simply given by:

$\hat{t}_{n}=n\frac{h_{1}}{2}\equiv n\frac{T}{2N},\qquad n=0,\dots,2N.$

Conversely, when a graded mesh (9)–(10) is considered, the doubled mesh is
given by:

$\hat{t}_{n}=\hat{t}_{n-1}+\hat{h}_{n},\qquad\hat{h}_{n}=\hat{r}\hat{h}_{n-1}\equiv\hat{r}^{n-1}\hat{h}_{1},\qquad
n=1,\dots,2N,$

with $\hat{t}_{0}=0$, and

$\hat{r}=\sqrt{r},\qquad\hat{h}_{1}=h_{1}\frac{\hat{r}-1}{r-1}.$

The choice of $\hat{h}_{1}$ is done for having

$\hat{h}_{1}\frac{\hat{r}^{2N}-1}{\hat{r}-1}=\hat{h}_{1}\frac{r^{N}-1}{\hat{r}-1}\equiv
h_{1}\frac{r^{N}-1}{r-1}=T,$

according to (10).

###  3.4 The nonlinear iteration

At last, we describe the efficient numerical solution of the discrete problem
(32), which has to be solved at the $n$-th integration step. As is clear, the
very formulation of the problem induces a straightforward fixed-point
iteration:

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

which can be conveniently started from ${\bm{\gamma}}^{n,0}=\bf 0$. The
following straightforward result holds true.

###### Theorem 3.2

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

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

###### Proof

See (BBBI2023, , Theorem 2).∎

Nevertheless, as is easily seen, also the simple equation

$y^{(\alpha)}=-\lambda y,\qquad t\in[0,T],\qquad y(0)=1,\qquad T,\lambda\gg
1,$

whose solution is almost everywhere close to 0, after an initial transient,
suffers from stepsize limitations, if the fixed-point iteration (50) is used,
since it has to be everywhere proportional to $\lambda^{-1/\alpha}$.

In order to overcome this drawback, a Newton-type iteration is therefore
needed. Hereafter, we consider the so-called blended iteration which has been
at first studied in a series of papers B2000 ; BM2002 ; BM2007 ; BM2009 . It
has been implemented in the Fortran codes BIM BM2004 , for ODE-IVPs, and BIMD
BMM2006 , for ODE-IVPs and linearly implicit DAEs, and in the Matlab code hbvm
BI2016 ; BIT2011 , for solving Hamiltonian problems. We here consider its
adaption for solving (32). By neglecting, for sake of brevity, the time-step
index $n$, we then want to solve the equation:

$G({\bm{\gamma}}):={\bm{\gamma}}-{\cal P}_{s}^{\top}\Omega\otimes
I_{m}f\left({\bm{\phi}}^{\alpha,s}+h^{\alpha}{\cal I}_{s}^{\alpha}\otimes
I_{m}{\bm{\gamma}}\right)={\bm{0}}.$ (51)

By setting $f_{0}^{\prime}$ the Jacobian of $f$ evaluated at the first entry
of ${\bm{\phi}}^{\alpha,s}$, $I=I_{s}\otimes I_{m}$, and

$X_{s}^{\alpha}:={\cal P}_{s}^{\top}\Omega{\cal I}_{s}^{\alpha},$ (52)

the application of the simplified Newton method then reads:

$\displaystyle\mathrm{solve:}$
$\displaystyle\left(I-h^{\alpha}X_{s}^{\alpha}\otimes
f_{0}^{\prime}\right)\Delta{\bm{\gamma}}^{\ell}=-G({\bm{\gamma}}^{\ell})\equiv{\bm{\eta}}^{\ell},$
(53) $\displaystyle\mathrm{set:}$
$\displaystyle{\bm{\gamma}}^{\ell+1}={\bm{\gamma}}^{\ell}+\Delta{\bm{\gamma}}^{\ell},\qquad\ell=0,1,\dots.$

Even though this iteration has the advantage of using a coefficient matrix
which is constant at each time-step, nevertheless, its dimension may be large,
when either $s$ or $m$ are large. To study a different iteration, able to get
rid of this problem, let us decouple the linear system into the various
eigenspaces of $f_{0}^{\prime}$, thus studying the simpler problem

$\displaystyle\mathrm{solve:}$ $\displaystyle\left(I_{s}-h^{\alpha}\mu
X_{s}^{\alpha}\right)\Delta\gamma^{\ell}=-g(\gamma^{\ell})\equiv\eta^{\ell},$
$\displaystyle\mathrm{set:}$
$\displaystyle\gamma^{\ell+1}=\gamma^{\ell}+\Delta\gamma^{\ell},\qquad\ell=0,1,\dots,$

with all involved vectors of dimension $s$, and $\mu\in\sigma(f_{0}^{\prime})$
a generic eigenvalue of $f_{0}^{\prime}$, and an obvious meaning of
$g(\gamma^{\ell})$. By setting $q=h^{\alpha}\mu$, the iteration then reads:

$\displaystyle\mathrm{solve:}$
$\displaystyle\left(I_{s}-qX_{s}^{\alpha}\right)\Delta\gamma^{\ell}=\eta^{\ell},$
(54) $\displaystyle\mathrm{set:}$
$\displaystyle\gamma^{\ell+1}=\gamma^{\ell}+\Delta\gamma^{\ell},\qquad\ell=0,1,\dots.$

Hereafter, we consider the iterative solution of the linear system in (54). A
linear analysis of convergence (in the case the r.h.s. is constant) is then
made, as at first suggested in HS1997 ; HS1997-1 ; HSS1996 , and later refined
in BM2002 ; BM2009 . Consequently, skipping the iteration index $\ell$, let us
consider the linear system to be solved:

$\left(I_{s}-qX_{s}^{\alpha}\right)\Delta\gamma=\eta,$

and its equivalent formulation, derived considering that matrix (52) is
nonsingular, and with $\xi>0$ a parameter to be later specified,

$\xi\left((X_{s}^{\alpha})^{-1}-qI_{s}\right)\Delta\gamma=\xi(X_{s}^{\alpha})^{-1}\eta=:\eta_{1}.$

Further, we consider the blending of the previous two equivalent formulations
with weights $\theta(q)$ and $I_{s}-\theta(q)$, where, by setting
$O\in\mathbb{R}^{s\times s}$ the zero matrix,

$\theta(q):=I_{s}(1-\xi q)^{-1}\left\\{\begin{array}[]{cc}\approx I_{s},&\quad
q\approx 0,\\\\[2.84526pt] \rightarrow
O,&q\rightarrow\infty,\end{array}\right.$

In so doing, one obtains the linear system

$M(q)\Delta\gamma=\eta_{1}+\theta(q)(\eta-\eta_{1})=:\eta(q),$ (55)

with the coefficient matrix,

$M(q)=\xi\left((X_{s}^{\alpha})^{-1}-qI_{s}\right)+\theta(q)\left[\left(I_{s}-qX_{s}^{\alpha}\right)-\xi\left((X_{s}^{\alpha})^{-1}-qI_{s}\right)\right],$

such that BM2002 :

$M(q)\approx\left\\{\begin{array}[]{cc}I_{s},&q\approx 0,\\\\[5.69054pt] -\xi
qI_{s},&|q|\gg 1.\end{array}\right.$

This naturally induces the splitting matrix

$N(q):=I_{s}(1-\xi q)\equiv\theta(q)^{-1},$ (56)

defining the blended iteration

$\Delta\gamma_{i}=\left[I_{s}-\theta(q)M(q)\right]\Delta\gamma_{i-1}+\theta(q)\eta(q),\qquad
i=1,2,\dots.$ (57)

This latter iteration converges iff the spectral radius of the iteration
matrix,

$\rho\left(I_{s}-\theta(q)M(q)\right)=:\rho(q)<1.$ (58)

The iteration is said to be $A$-convergent if (58) holds true for all
$q\in\mathbb{C}^{-}$, the left-half complex plane, and $L$-convergent if, in
addition, $\rho(q)\rightarrow 0$, as $q\rightarrow\infty$. Since BM2002

$\theta(0)M(0)=I_{s},\qquad\theta(q)M(q)\rightarrow I_{s},\quad
q\rightarrow\infty,$

the blended iteration is $L$-convergent iff it is $A$-convergent. For this
purpose, we shall look for a suitable choice of the positive parameter
$\xi>0$. Considering that $\theta(q)$ is well defined for all
$q\in\mathbb{C}^{-}$, the following statement easily follows from the maximum
modulus theorem.

###### Theorem 3.3

The blended iteration is $L$-convergent iff the maximum amplification factor,

$\rho^{*}:=\max_{x>0}\rho(\mathrm{i}x)\leq 1.$

The following result also holds true.

###### Theorem 3.4

The eigenvalues of the iteration matrix $I_{s}-\theta(q)M(q)$ are given by:

$\frac{q(\lambda-\xi)^{2}}{\lambda(1-q\xi)^{2}},\qquad\lambda\in\sigma(X_{s}^{\alpha}).$

Consequently, the maximum amplification factor is given by:

$\rho^{*}=\max_{\lambda\in\sigma(X_{s}^{\alpha})}\frac{|\lambda-\xi|^{2}}{2\xi|\lambda|}.$
(59)

###### Proof

See (BM2002, , Theorem 2 and Equation (25)), by considering that

$\max_{q\in\mathbb{C}^{-}}\frac{q(\lambda-\xi)^{2}}{\lambda(1-q\xi)^{2}}~{}\equiv~{}\max_{x>0}\frac{x|\lambda-\xi|^{2}}{|\lambda|(1+x^{2}\xi^{2})}$

is obtained at $x=\xi^{-1}$, so that (59) follows. ∎

Slightly generalizing the arguments in BM2002 , we then consider the following
choice of the parameter $\xi$,

$\xi=\mathrm{argmin}_{\mu\in\sigma(X_{s}^{\alpha})}\max_{\lambda\in\sigma(X_{s}^{\alpha})}\frac{|\lambda-|\mu||^{2}}{2|\mu||\lambda|},$
(60)

which is computed once forall, and always provides, in our experiments, an
$L$-convergent iteration. In particular, the code fhbvm uses, at the moment,
$k=22$ and $s=20$: the corresponding maximum amplification factor (59) is
depicted in Figure 1, w.r.t. the order $\alpha$ of the fractional derivative,
thus confirming this.

Figure 1: Maximum amplification factor (59)–(60), $k=22$ and $s=20$.

Coming back to the original problem (53), starting from the initial guess
$\Delta{\bm{\gamma}}={\bm{0}}$, and updating the r.h.s. as soon as a new
approximation to ${\bm{\gamma}}$ is available, one has that the iteration
(55)–(57) simplifies to:

$\displaystyle{\bm{\eta}}^{\ell}$ $\displaystyle=$
$\displaystyle-G({\bm{\gamma}}^{\ell})$ $\displaystyle{\bm{\eta}}_{1}^{\ell}$
$\displaystyle=$ $\displaystyle\xi(X_{s}^{\alpha})^{-1}\otimes
I_{m}\,{\bm{\eta}}^{\ell}$ (61) $\displaystyle{\bm{\gamma}}^{\ell+1}$
$\displaystyle=$
$\displaystyle{\bm{\gamma}}^{\ell}+I_{s}\otimes\Theta\left[{\bm{\eta}}_{1}^{\ell}+I_{s}\otimes\Theta\left({\bm{\eta}}^{\ell}-{\bm{\eta}}_{1}^{\ell}\right)\right],\qquad\ell=0,1,\dots,$

with $\xi$ chosen according to (60), and

$\Theta=\left(I_{m}-h^{\alpha}\xi f_{0}^{\prime}\right)^{-1}.$

Consequently, only the factorization of one matrix, having the same dimension
$m$ of the problem, is needed. Moreover, the initial guess
${\bm{\gamma}}^{0}={\bm{0}}$ can be conveniently considered in (61).

###### Remark 6

It is worth mentioning that, due to the properties of the Kronecker product,
the iteration (61) can be compactly cast in matrix form, thus avoiding an
explicit use of the Kronecker product. This implementation has been considered
in the code fhbvm used in the numerical tests.

Actually, according to Theorem 3.2, in the code fhbvm we automatically switch
between the fixed-point iteration (50) or the blended iteration (61),
depending on the fact that

$h^{\alpha}\|f_{0}^{\prime}\|\|{\cal P}_{s}^{\top}\Omega\|\|{\cal
I}_{s}^{\alpha}\|\leq tol,$

with $tol<1$ a suitable tolerance.

##  4 Numerical Tests

In this section we report a few numerical tests using the
Matlab${}^{\copyright}\,$code fhbvm: the code implements a FHBVM(22,20) method
using all the strategies discussed in the previous section. The calling
sequence of the code is:

[t,y,stats,err] = fhbvm( fun, y0, T, M )

where:

•

In input:

  * –

fun is the identifier (or the function handling) of the function evaluating
the r.h.s. of the equation (also in vector mode), its Jacobian, and the order
$\alpha$ of the fractional derivative (see help fhbvm for more details);

  * –

y0 is the initial condition;

  * –

T is the final integration time;

  * –

M is the parameter in (36) (it should be as small as possible);

•

In output:

  * –

t,y contain the computed mesh and solution;

  * –

stats (optional) is a vector containing the following time statistics:

    1. 1.

the pre-processing time for computing the parameters $h_{1}$, $r$, and $N$
(see Section 3.1) and the fractional itegrals (47), and (48) or (49);

    2. 2.

the time for solving the problem;

    3. 3.

the pre-processing time for computing the fractional itegrals (47), and (48)
or (49) for the error estimation;

    4. 4.

the time for solving the problem on the doubled mesh, for the error
estimation;

  * –

err (optional), if specified, contains the estimate of the absolute error.
This estimate, obtained on a doubled mesh, is relatively costly: for this
reason, when the parameter is not specified, the solution on the doubled mesh
is not computed.

For the first two problems, we shall also make a comparison with the
Matlab${}^{\copyright}\,$code flmm2 Ga2015 ,777In particular, the BDF2 method
is selected (method=3), with the parameters tol=1e-15 and itmax=1000. in order
to emphasize the potentialities of the new code. All numerical tests have been
done on a M2-Silicon based computer with 16GB of shared memory, using
Matlab${}^{\copyright}\,$R2023b.

The comparisons will be done by using a so called Work Precision Diagram
(WPD), where the execution time (in sec) is plotted against accuracy. The
accuracy, in turn, is measured through the mixed error significant computed
digits (mescd) testset , defined, by using the same notation seen in (37), as
888This definition corresponds to set atol=rtol in the definition used in
testset .

${\tt
mescd}:=-\log_{10}\max_{i=0,\dots,N}\|(y(t_{i})-\bar{y}_{i})./(1+|y(t_{i})|)\|_{\infty},$

being $t_{i}$, $i=0,\dots,N$, the computational mesh of the considered solver,
and $y(t_{i})$ and $\bar{y}_{i}$ the corresponding values of the solution and
of its approximation.

###  4.1 Example 1

The first problem Ga2018  is given by:

$\displaystyle y^{(\alpha)}$ $\displaystyle=$
$\displaystyle-|y|^{1.5}+\frac{8!}{\Gamma(9-\alpha)}t^{8-\alpha}-3\frac{\Gamma(5+\alpha/2)}{\Gamma(5-\alpha/2)}t^{4-\alpha/2}+\left(\frac{3}{2}t^{\alpha/2}-t^{4}\right)^{3}$
(62) $\displaystyle~{}+\frac{9}{4}\Gamma(\alpha+1),\qquad t\in[0,1],\qquad
y(0)=0,$

whose solution is

$y(t)=t^{8}-3\,t^{4+\alpha/2}+\frac{9}{4}\,t^{\alpha}.$

We consider the value $\alpha=0.3$, and use the codes with the following
parameters, to derive the corresponding WPD:

  * •

flmm2 : $h=10^{-1}2^{-\nu}$, $\nu=1,\dots,20$;

  * •

fhbvm : $M=2,3,4,5$.

Figure 2 contains the obtained results: as one may see, flmm2 reaches less
than 12 mescd, since by continuing reducing the stepsize, at the 15-th mesh
doubling the error starts increasing. On the other hand, the execution time
essentially doubles at each new mesh doubling. Conversely, fhbvm can achieve
full machine accuracy by employing a uniform mesh with stepsize $h_{1}=1/M$,
with $M$ very small, thus using very few mesh points. As a result, fhbvm
requires very short execution times.

Figure 2: Work-precision diagram for problem (62), $\alpha=0.3$.

###  4.2 Example 2

We now consider the following linear problem:

$\displaystyle y^{(0.5)}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{rr}-50&~{}0\\\
-49&~{}-1\end{array}\right)y,\qquad t\in[0,20],$ (65) $\displaystyle y(0)$
$\displaystyle=$ $\displaystyle(2,\,3)^{\top},$ (66)

having solution

$y(t)=\left(\begin{array}[]{c}2\,E_{0.5}\left(-50\cdot
t^{0.5}\right)\\\\[2.84526pt] 2\,E_{0.5}\left(-50\cdot
t^{0.5}\right)+E_{0.5}\left(-t^{0.5}\right)\end{array}\right),$

with $E_{0.5}$ the Mittag-Leffler function.999We have used the
Matlab${}^{\copyright}\,$function ml Ga2015ml  for its evaluation. We use the
codes with the following parameters, to derive the corresponding WPD:

  * •

flmm2 : $h=10^{-1}2^{-\nu}$, $\nu=1,\dots,20$;

  * •

fhbvm : $M=5,\dots,10$.

Figure 4 contains the obtained results, from which one deduces that flmm2
achieves about 5 mescd (with a run time of about 85 sec), whereas fhbvm has
approximately 13 mescd, with an execution time of about 1 sec. Further, in
Figure 4, we plot the true and estimated (absolute) errors for fhbvm in the
case $M=10$ (corresponding to a computational mesh made up of 251 mesh-points,
with an initial stepsize $h_{1}\approx 7.3\cdot 10^{-12}$, and a final
stepsize $h_{250}\approx 2$): as one may see from the figure, there is a
substantial agreement between the two errors.

Figure 3: Work-precision diagram for problem (66).

Figure 4: True and estimated absolute errors for fhbvm solving problem (66),
$M=10$.

###  4.3 Example 3

We now consider the following nonlinear problem BBBI2023 :

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

having solution

$y(t)=\left(\begin{array}[]{c}t^{2/3}+1\\\\[2.84526pt]
t^{4/3}\end{array}\right).$

This problem is relatively simple and, in fact, both flmm2 and fhbvm solve it
accurately. We use it to show the estimated error by using fhbvm with
parameter $M=2$, which produces a graded mesh with 41 mesh-points, with
$h_{1}\approx 1.8\cdot 10^{-12}$ and $h_{40}\approx 0.49$. The absolute errors
(true and estimated) for each component are depicted in Figure 5, showing a
perfect agreement for both of them.

In this case, the evaluation of the solution requires $\approx 0.04$ sec, and
the error estimation requires $\approx 0.11$ sec.

Figure 5: True and estimated absolute errors for fhbvm solving problem (67),
$M=2$.

###  4.4 Example 4

At last, we consider the following fractional Brusselator model:

$\displaystyle y_{1}^{(0.7)}$ $\displaystyle=$ $\displaystyle
1-4y_{1}+y_{1}^{2}y_{2},$ $\displaystyle y_{2}^{(0.7)}$ $\displaystyle=$
$\displaystyle 3y_{1}-y_{1}^{2}y_{2},\qquad t\in[0,5],$ $\displaystyle y(0)$
$\displaystyle=$ $\displaystyle(1.2,\,2.8)^{\top},$ (68)

By solving this problem using fhbvm with parameter $M=5$, a graded mesh of 46
points is produced, with $h_{1}\approx 6.1\cdot 10^{-5}$ and $h_{45}\approx
0.98$. The maximum estimated error in the computed solution is less than
$3.5\cdot 10^{-13}$, whereas the phase-plot of the solution is depicted in
Figure 6.

In this case, the evaluation of the solution requires $\approx 0.04$ sec, and
the error estimation requires $\approx 0.14$ sec.

Figure 6: Phase-plot of the computed solution by using fhbvm solving problem
(68), $M=5$.

##  5 Conclusions

In this paper we have described in full details the implementation of the
Matlab${}^{\copyright}\,$ code fhbvm, able to solving systems of FDE-IVPs. The
code is based on a FHBVM(22,20) method, as described in BBBI2023 . We have
also provided comparisons with another existing Matlab${}^{\copyright}\,$code,
thus confirming its potentialities. In fact, due to the spectral accuracy in
time of the FHBVM(22,20) method, the generated computational mesh, which can
be either a uniform or a graded one, depending on the problem at hand,
requires relatively few mesh points. This, in turn, allows to reduce the
execution time due to the evaluation of the memory term required at each step.

We plan to further develop the code fhbvm, in order to provide approximations
at prescribed mesh points, as well as to allow selecting different
FHBVM$(k,s)$, $k\geq s$, methods BBBI2023 . At last, we plan to extend the
code to cope with values of the fractional derivative, $\alpha$, greater than
1.

#### Declarations.

The authors declare no conflict of interests, nor competing interests. No
grants were received for conducting this study.

#### Data availability.

All data reported in the manuscript have been obtained by the Matlab© code
fhbvm, Rel. 2024-03-06, available at the url fhbvm .

#### Acknowledgements.

The authors are members of the Gruppo Nazionale Calcolo Scientifico-Istituto
Nazionale di Alta Matematica (GNCS-INdAM).

## References

  * (1) L. Aceto, C. Magherini, P. Novati. Fractional convolution quadrature based on generalized Adams methods. Calcolo 51 (2014) 441–463. https://doi.org/10.1007/s10092-013-0094-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. Arbitrarily high-order energy-conserving methods for Poisson problems. Numer. Algoritms 91 (2022) 861–894. https://doi.org/10.1007/s11075-022-01285-z
  * (5) 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
  * (6) 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
  * (7) L. Brugnano. Blended block BVMs (B3VMs): A family of economical implicit methods for ODEs. J. Comput. Appl. Math. 116 (2000) 41–62. https://doi.org/10.1016/S0377-0427(99)00280-0
  * (8) L. Brugnano, K. Burrage, P. Burrage, F. Iavernaro. A spectrally accurate step-by-step method for the numerical solution of fractional differential equations. arXiv:2310.10526 [math.NA] https://doi.org/10.48550/arXiv.2310.10526 (submitted)
  * (9) 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
  * (10) L. Brugnano, F. Iavernaro. Line Integral Methods for Conservative Problems. Chapman et Hall/CRC, Boca Raton, FL, USA, 2016.
  * (11) L. Brugnano, F. Iavernaro. Line Integral Solution of Differential Problems. Axioms 7(2) (2018) 36. https://doi.org/10.3390/axioms7020036
  * (12) 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
  * (13) L. Brugnano, F. Iavernaro, D. Trigiante. A note on the efficient implementation of Hamiltonian BVMs. J. Comput. Appl. Math. 236 (2011) 375–383. https://doi.org/10.1016/j.cam.2011.07.022
  * (14) L. Brugnano, C. Magherini. Blended Implementation of Block Implicit Methods for ODEs. Appl. Numer. Math. 42 (2002) 29–45. https://doi.org/10.1016/S0168-9274(01)00140-4
  * (15) L. Brugnano, C. Magherini. The BiM Code for the Numerical Solution of ODEs. J. Comput. Appl. Math. 164–165 (2004) 145–158. https://doi.org/10.1016/j.cam.2003.09.004
  * (16) L. Brugnano, C. Magherini. Blended Implicit Methods for solving ODE and DAE problems, and their extension for second order problems. J. Comput. Appl. Math. 205 (2007) 777–790. https://doi.org/10.1016/j.cam.2006.02.057
  * (17) L. Brugnano, C. Magherini. Recent Advances in Linear Analysis of Convergence for Splittings for Solving ODE problems. Appl. Numer. Math. 59 (2009) 542–557. https://doi.org/10.1016/j.apnum.2008.03.008
  * (18) L. Brugnano, C. Magherini, F. Mugnai. Blended Implicit Methods for the Numerical Solution of DAE Problems. J. Comput. Appl. Math. 189 (2006) 34–50. https://doi.org/10.1016/j.cam.2005.05.005
  * (19) 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
  * (20) 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
  * (21) A. Cardone, D. Conte, B. Paternoster. A Matlab code for fractional differential equations based on two-step spline collocation methods. In Fractional Differential Equations, Modeling, Discretization, and Numerical Solvers, A. Cardone et al. Eds., Springer INDAM Series, Vol. 50, 2023, pp. 121–146. https://doi.org/10.1007/978-981-19-7716-9_8
  * (22) G. Dahlquist, Å. Björk. Numerical Methods in Scientific Computing. SIAM: Philadelphia, PA, USA, 2008.
  * (23) 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.
  * (24) K. Diethelm, N.J. Ford, A.D. Freed. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. 29 (2002) 3–22. https://doi.org/10.1023/A:1016592219341
  * (25) 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
  * (26) 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
  * (27) R. Garrappa. Trapezoidal methods for fractional differential equations: theoretical and computational aspects. Mathematics and Computers in Simulation 110 (2015) 96–112. https://doi.org/10.1016/j.matcom.2013.09.012
  * (28) 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
  * (29) Z. Gu. Spectral collocation method for nonlinear Riemann-Liouville fractional terminal value problems. J. Compt. Appl. math. 398 (2021) 113640. https://doi.org/10.1016/j.cam.2021.113640
  * (30) Z. Gu, Y. Kong. Spectral collocation method for Caputo fractional terminal value problems. Numer. Algorithms 88 (2021) 93–111. https://doi.org/10.1007/s11075-020-01031-3
  * (31) P.J. van der Houwen, J.J.B. de Swart. Triangularly implicit iteration methods for ODE-IVP solvers. SIAM J. Sci. Comput. 18 (1997) 41–55. https://doi.org/10.1137/S1064827595287456
  * (32) P.J. van der Houwen, J.J.B. de Swart. Parallel linear system solvers for Runge-Kutta methods. Adv. Comput. Math. 7 (1–2) (1997) 157–181. https://doi.org/10.1023/A:1018990601750
  * (33) P.J. van der Houwen, B.P. Sommeijer, J.J. de Swart. Parallel predictor–corrector methods. J. Comput. Appl. Math. 66 (1996) 53–71. https://doi.org/10.1016/0377-0427(95)00158-1
  * (34) 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. https://doi.org/10.1016/j.jcp.2016.04.039
  * (35) 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
  * (36) 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.
  * (37) 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. https://doi.org/10.1137/16M1082329
  * (38) M. Zayernouri, G.E. Karniadakis. Exponentially accurate spectral and spectral element methods for fractional ODEs. J. Comput. Phys. 257 (2014) 460–480. https://doi.org/10.1016/j.jcp.2013.09.039
  * (39) F. Mazzia, C. Magherini. Test Set for Initial Value Problem Solvers. Release 2.4, February 2008, Department of Mathematics, University of Bari and INdAM, Research Unit of Bari, Italy, available at: https://archimede.uniba.it/~testset/testsetivpsolvers/
  * (40) https://people.dimai.unifi.it/brugnano/fhbvm/