# A shooting-Newton procedure for solving fractional terminal value problems

Luigi Brugnano 111L. Brugnano, Dipartimento di Matematica e Informatica “U.
Dini”, Università di Firenze, Italy<EMAIL_ADDRESS>https://orcid.org/0000-0002-6290-4107 Gianmarco Gurioli 222G. Gurioli,
Dipartimento di Matematica e Informatica “U. Dini”, Università di Firenze,
Italy<EMAIL_ADDRESS>https://orcid.org/0000-0003-0922-8119
Felice Iavernaro 333F. Iavernaro, Dipartimento di Matematica, Università di
Bari Aldo Moro, Italy<EMAIL_ADDRESS>https://orcid.org/0000-0002-9716-7370

###### Abstract

In this paper we consider the numerical solution of fractional terminal value
problems (FDE-TVPs). In particular, the proposed procedure uses a Newton-type
iteration which is particularly efficient when coupled with a recently-
introduced step-by-step procedure for solving fractional initial value
problems (FDE-IVPs), able to produce spectrally accurate solutions of FDE
problems. Some numerical tests are reported to make evidence of its
effectiveness.

Keywords: fractional differential equations, fractional integrals, terminal
value problems, Jacobi polynomials, Fractional Hamiltonian Boundary Value
Methods, FHBVMs.

MSC: 34A08, 65R20.

##  1 Introduction

Fractional differential equations have gained more and more importance in many
applications: we refer, e.g., to the classical references [13, 24] for an
introduction.

The present contribution is mainly motivated by the recent paper of Diethelm
and Uhlig [16], where an iterative procedure, based on the secant method, is
proposed for solving fractional terminal value problems (FDE-TVPs) that,
without loss of generality, we assume in the form:

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

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)

Usually, one solves the fractional initial value problem (FDE-IVP) (see, e.g.
[15, 20, 21, 22, 23, 25]):

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

whose solution is given by

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

However, under mild hypothesis on $f$, also the FDE-TVP is well-posed [17],
and its numerical solution has been considered by many authors (see, e.g.,
[16, 17]). In particular, the scalar case of (1) ($m=1$, with $f(t,y)$ in
place of $f(y)$) allows using a shooting procedure coupled with the bisection
method [14, 18] or, more recently, with the secant method [16]. In this paper,
we propose an alternative approach, based on a straight Newton procedure,
which takes advantage of a recently introduced method for solving FDE-IVPs,
able to obtain spectrally accurate approximations [6]. This latter approach
has been derived as an extension of HBVM methods, originally devised for
Hamiltonian problems (see, e.g., [8, 9]), and later extended along several
directions (see, e.g., [1, 3, 5, 6, 7, 10]), including the numerical solution
of FDEs. A main feature of HBVMs is the fact that they can gain spectrally
accuracy, when approximating ODE-IVPs [2, 11, 12], and such a feature has been
recently extended to the FDE case [6].

With this premise, the structure of the paper is as follows: in Section 2 we
sketch the shooting-Newton procedure for solving (1); in Section 3 we recall
the main facts about the (possibly spectrally accurate) numerical solution of
FDE-IVPs proposed in [6], along with its extension for the shooting-Newton
procedure; in Section 4 we report a few numerical tests; at last, a few
conclusions are given in Section 5.

##  2 The shooting-Newton procedure

To begin with, let us introduce a perturbation result concerning the solution
of the FDE-IVP (4). In particular, let us denote by $y(t,\rho_{0})$ the
solution of this problem, in order to emphasize its dependence from the
initial condition. The following result holds true.

###### Theorem 1

For $t\in[0,T]$, one has:

$\frac{\partial}{\partial\rho_{0}}y(t,\rho_{0})=\Phi(t,\rho_{0}),$ (6)

the solution of the fractional variational problem 444As is usual,
$f^{\prime}(y)$ denotes the Jacobian matrix of $f(y)$.

$\Phi^{(\alpha)}(t,\rho_{0})=f^{\prime}(y(t,\rho_{0}))\Phi(t,\rho_{0}),\qquad
t\in[0,T],\qquad\Phi(0,\rho_{0})=I,$ (7)

explicitly given by:

$\Phi(t,\rho_{0})=I+\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-x)^{\alpha-1}f^{\prime}(y(x,\rho_{0}))\Phi(x,\rho_{0})\mathrm{d}x.$
(8)

Proof In fact, from (2) and (4), one has:

$D^{\alpha}\frac{\partial}{\partial\rho_{0}}y(t,\rho_{0})\,=\,\frac{\partial}{\partial\rho_{0}}D^{\alpha}y(t)\,=\,\frac{\partial}{\partial\rho_{0}}f(y(t,\rho_{0}))\,=\,f^{\prime}(y(t,\rho_{0}))\frac{\partial}{\partial\rho_{0}}y(t,\rho_{0}),$

and

$\left.\frac{\partial}{\partial\rho_{0}}y(t,\rho_{0})\right|_{t=0}\,=\,\frac{\partial}{\partial\rho_{0}}\rho_{0}\,=\,I.$

Consequently, (6)-(7) follows and, therefore, also (8) follows from (3) and
(5).  $\Box{~{}}$

###### Remark 1

Hereafter, in order to guarantee the well-posedness of problem (1), if
$y(0)=\rho^{*}$ is the initial value of (4) fulfilling the FDE-TVP, i.e.,

$y(T,\rho^{*})=\eta,$ (9)

we shall assume that (see (6))

$\det\left(\Phi(T,\rho^{*})\right)\neq 0.$ (10)

Under regularity assumptions on $f(y)$, in turn (10) implies that

$\exists\,\delta>0~{}s.t.~{}\|\rho^{*}-\rho\|\leq\delta~{}\Rightarrow~{}\det\left(\Phi(T,\rho)\right)\neq
0.$ (11)

The previous results allow us stating the shooting-Newton procedure for
solving (1) sketched in Table 1, where a suitable stopping criterion has to be
adopted. Moreover, the starting approximation $\rho_{0}$ for the shooting-
Newton iteration has to be chosen in some way, possibly exploiting any
additional information; conversely, the choice $\rho_{0}=\eta$ (i.e., the
final value in (1)) can be considered, as proposed in [16].

Table 1: Algorithm 1 – the shooting-Newton procedure. $\displaystyle{\rm
fix~{}}\rho_{0}$ $\displaystyle{\rm for~{}~{}}\ell$ $\displaystyle=$
$\displaystyle 0,1,\dots$ $\displaystyle{\rm
solve:}~{}\,y^{(\alpha)}(t,\rho_{\ell})=f(y(t,\rho_{\ell})),\qquad\qquad~{}\,\,t\in[0,T],\qquad
y(0,\rho_{\ell})=\rho_{\ell}$ $\displaystyle{\rm
and}~{}~{}\quad\Phi^{(\alpha)}(t,\rho_{\ell})=f^{\prime}(y(t,\rho_{\ell}))\Phi(t,\rho_{\ell}),\quad
t\in[0,T],\qquad\Phi(0,\rho_{\ell})=I$ $\displaystyle{\rm
set:~{}}\quad\rho_{\ell+1}=\rho_{\ell}-\Phi(T,\rho_{\ell})^{-1}\left[y(T,\rho_{\ell})-\eta\right]$
$\displaystyle{\rm end~{}~{}}$

The following straightforward convergence result holds true.

###### Theorem 2

Under regularity assumptions on $f(y)$, the shooting-Newton procedure given in
Table 1 converges in a suitable neighborhood of $\rho^{*}$.

Proof In fact, from (9) it follows that $\rho^{*}$ is a fixed-point of the
corresponding iteration function,

$\Psi(\xi):=\xi-\Phi(T,\xi)^{-1}\left[y(T,\xi)-\eta\right],$

whose Jacobian (recall (6)) vanishes at $\xi=\rho^{*}$. Consequently, from the
Perron Theorem [26, Corollary 4.7.2], exponential convergence is granted, in a
suitable neighborhood of $\rho^{*}$.  $\Box{~{}}$

##  3 Implementing the algorithm

Following the approach in [6], let us now explain the way we compute
$y(T,\rho_{\ell})$ and $\Phi(T,\rho_{\ell})$ in Algorithm 1. From (5) and (8),
we have to compute:

$y(T,\rho_{\ell})~{}=~{}\rho_{\ell}+\frac{1}{\Gamma(\alpha)}\int_{0}^{T}(T-x)^{\alpha-1}f(y(x,\rho_{\ell}))\mathrm{d}x,$
(12)

and

$\Phi(T,\rho_{\ell})~{}=~{}I+\frac{1}{\Gamma(\alpha)}\int_{0}^{T}(T-x)^{\alpha-1}f^{\prime}(y(x,\rho_{\ell}))\Phi(x,\rho_{\ell})\mathrm{d}x.$
(13)

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

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

where

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

In general, for coping with possible singularities in the derivative of the
vector field at the origin, we shall consider the following graded mesh,

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

where $r>1$ and $h_{1}>0$ satisfy, by virtue of (15)-(16),

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

Clearly, when a uniform mesh is considered then, in (16), $r=1$ and
$h_{1}=h:=T/N$, so that $h_{n}=h$, $n=1,\dots,N$.

By setting

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

the restriction of the solution of (12) on the interval $[t_{n-1},t_{n}]$, and
taking into account (14)–(16), one then obtains:

$\displaystyle
y(T,\rho_{\ell})~{}\equiv~{}y_{N}(h_{N},\rho_{\ell})~{}=~{}\rho_{\ell}+\frac{1}{\Gamma(\alpha)}\int_{0}^{T}(T-x)^{\alpha-1}f(y(x,\rho_{\ell}))\mathrm{d}x$
(19) $\displaystyle=$
$\displaystyle~{}\rho_{\ell}+\frac{1}{\Gamma(\alpha)}\sum_{n=1}^{N}\int_{t_{n-1}}^{t_{n}}(t_{N}-x)^{\alpha-1}f(y(x,\rho_{\ell}))\mathrm{d}x$
$\displaystyle=$
$\displaystyle~{}\rho_{\ell}+\frac{1}{\Gamma(\alpha)}\sum_{n=1}^{N}\int_{0}^{h_{n}}(t_{N}-t_{n-1}-x)^{\alpha-1}f(y_{n}(x,\rho_{\ell}))\mathrm{d}x$
$\displaystyle=$
$\displaystyle~{}\rho_{\ell}+\frac{1}{\Gamma(\alpha)}\sum_{n=1}^{N}h_{n}^{\alpha}\int_{0}^{1}\left(\frac{r^{N-n+1}-1}{r-1}-\tau\right)^{\alpha-1}f(y_{n}(\tau
h_{n},\rho_{\ell}))\mathrm{d}\tau.\qquad$

In case of a constant stepsize, $h=T/N$ is used, the previous equation
becomes:

$\displaystyle y(T,\rho_{\ell})$ $\displaystyle\equiv$ $\displaystyle
y_{N}(h,\rho_{\ell})$ (20) $\displaystyle=$
$\displaystyle\rho_{\ell}+\frac{h^{\alpha}}{\Gamma(\alpha)}\sum_{n=1}^{N}\int_{0}^{1}\left(N-n+1-\tau\right)^{\alpha-1}f(y_{n}(\tau
h,\rho_{\ell}))\mathrm{d}\tau.\quad$

Similarly, for (13), by setting

$\Phi_{n}(ch_{n},\rho_{\ell}):=\Phi(t_{n-1}+ch_{n},\rho_{\ell}),\qquad
c\in[0,1],\qquad n=1,\dots,N,$ (21)

the restriction of the solution on the interval $[t_{n-1},t_{n}]$, again by
virtue of (14)–(16), one obtains:

$\displaystyle\Phi(T,\rho_{\ell})~{}\equiv~{}\Phi_{N}(h_{N},\rho_{\ell})$
$\displaystyle=$ $\displaystyle
I+\frac{1}{\Gamma(\alpha)}\int_{0}^{T}(T-x)^{\alpha-1}f^{\prime}(y(x,\rho_{\ell}))\Phi(x,\rho_{\ell})\mathrm{d}x$
$\displaystyle=$ $\displaystyle
I+\frac{1}{\Gamma(\alpha)}\sum_{n=1}^{N}\int_{t_{n-1}}^{t_{n}}(t_{N}-x)^{\alpha-1}f^{\prime}(y(x,\rho_{\ell}))\Phi(x,\rho_{\ell})\mathrm{d}x$
$\displaystyle=$ $\displaystyle
I+\frac{1}{\Gamma(\alpha)}\sum_{n=1}^{N}\int_{0}^{h_{n}}(t_{N}-t_{n-1}-x)^{\alpha-1}f^{\prime}(y_{n}(x,\rho_{\ell}))\Phi_{n}(x,\rho_{\ell})\mathrm{d}x$
$\displaystyle=$ $\displaystyle
I+\frac{1}{\Gamma(\alpha)}\sum_{n=1}^{N}h_{n}^{\alpha}\int_{0}^{1}\left(\frac{r^{N-n+1}-1}{r-1}-\tau\right)^{\alpha-1}f^{\prime}(y_{n}(\tau
h_{n},\rho_{\ell}))\Phi_{n}(\tau h_{n},\rho_{\ell})\mathrm{d}\tau,$

and, in case of a constant stepsize $h=T/N$,

$\displaystyle\Phi(T,\rho_{\ell})~{}\equiv~{}\Phi_{N}(h,\rho_{\ell})$ (23)
$\displaystyle=$ $\displaystyle
I+\frac{h^{\alpha}}{\Gamma(\alpha)}\sum_{n=1}^{N}\int_{0}^{1}\left(N-n+1-\tau\right)^{\alpha-1}f^{\prime}(y_{n}(\tau
h,\rho_{\ell}))\Phi_{n}(\tau h,\rho_{\ell})\mathrm{d}\tau.\qquad$

###  3.1 Piecewise quasi-polynomial approximation

The previous functions are then approximated via a piecewise quasi-polynomial
approximation, as described in [6], which we here briefly recall, and
generalize to the approximation of the fundamental matrix solution, too. In
more detail, with reference to (18) and (21), we shall look for
approximations:

$\displaystyle\sigma_{n}(ch_{n},\rho_{\ell})$ $\displaystyle\simeq$
$\displaystyle y_{n}(ch_{n},\rho_{\ell}),$ (24)
$\displaystyle\Psi_{n}(ch_{n},\rho_{\ell})$ $\displaystyle\simeq$
$\displaystyle\Phi_{n}(ch_{n},\rho_{\ell}),\qquad c\in[0,1],\qquad
n=1,\dots,N,$

and, consequently,

$y(T,\rho_{\ell})\simeq\sigma_{N}(h_{N},\rho_{\ell}),\qquad\Phi(T,\rho_{\ell})\simeq\Psi_{N}(h_{N},\rho_{\ell}).$
(25)

Following steps similar to those in [6, Section 2], we 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,$

resulting into a scaled and shifted family of Jacobi polynomials:555Here,
$\bar{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}}\bar{P}^{(a-1,0)}(2x-1),\qquad
x\in[0,1],\qquad j=0,1,\dots.$

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

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

with

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

The approximations are derived by truncating the infinite series to a finite
sum with $s$ terms. Consequently, for $n=1,\dots,N$, one obtains:666We refer
to [4] for an efficient procedure for computing the fractional integrals
$I_{j}^{\alpha}P_{j}(c)$, $j=0,\dots,s-1.$

$\sigma_{n}(ch_{n},\rho_{\ell})=\phi_{n-1}^{\alpha}(c,\rho_{\ell})+h_{n}^{\alpha}\sum_{j=0}^{s-1}\gamma_{j}(\sigma_{n},\rho_{\ell})I^{\alpha}P_{j}(c),\qquad
c\in[0,1],$ (26)

with

$\gamma_{j}(\sigma_{n},\rho_{\ell})=\alpha\int_{0}^{1}(1-\tau)^{\alpha-1}P_{j}(\tau)f(\sigma_{n}(\tau
h_{n},\rho_{\ell}))\mathrm{d}\tau,\qquad j=0,\dots,s-1.\quad$ (27)

and

$\displaystyle\phi_{n-1}^{\alpha}(c,\rho_{\ell})~{}=$ (28)
$\displaystyle\rho_{\ell}+\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}h_{\nu}^{\alpha}\sum_{j=0}^{s-1}J_{j}^{\alpha}\left(\frac{r^{n-\nu}-1}{r-1}+cr^{n-\nu}\right)\gamma_{j}(\sigma_{\nu},\rho_{\ell}),$

having set, for $x>1$, 777We refer to [6, Section 3.3] for the efficient
computation of such integrals.

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

If a constant stepsize $h=T/N$ is used, then (28) reads:

$\phi_{n-1}^{\alpha}(c,\rho_{\ell})=\rho_{\ell}+\frac{h^{\alpha}}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}\sum_{j=0}^{s-1}J_{j}^{\alpha}(n-\nu+c)\gamma_{j}(\sigma_{\nu},\rho_{\ell}),$
(30)

and similarly one modifies (26) and (27).

It can be shown (see [6]) that $\phi_{n-1}^{\alpha}(c,\rho_{\ell})$ is nothing
but the approximation of the memory term

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

such that, for all $c\in[0,1]$, and $n=1,\dots,N$:

$y_{n}(ch_{n},\rho_{\ell})=G_{n-1}^{\alpha}(c,\rho_{\ell})+\frac{h_{n}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}\left(c-\tau\right)^{\alpha-1}f(y_{n}(\tau
h_{n},\rho_{\ell}))\mathrm{d}\tau.$ (31)

As matter of fact, (19) corresponds to set $n=N$ and $c=1$ in (31).

Similarly, when a constant stepsize $h=T/N$ is used, then

$G_{n-1}^{\alpha}(c,\rho_{\ell})=\rho_{\ell}+\frac{h^{\alpha}}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}\int_{0}^{1}\left(n-\nu+c-\tau\right)^{\alpha-1}f(y_{n}(\tau
h,\rho_{\ell}))\mathrm{d}\tau,$

and (31) still formally holds, upon replacing $h_{n}$ with $h$. Consequently,
(20) corresponds again to set $n=N$ and $c=1$ in (31).

The Fourier coefficients (27) can be approximated up to machine precision by
using the Gauss-Jacobi formula of order $2k$ based at the zeros of $P_{k}(x)$,
$c_{1},\dots,c_{k}$, with corresponding weights $b_{1},\dots,b_{k}$, by
choosing $k$ large enough. As is explained in [6, Section 3], this allows
formulating the discrete problem for computing them as:888As 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}={\cal P}_{s}^{\top}\Omega\otimes
I_{m}f\left({\bm{\phi}}_{n-1}^{\alpha}+h_{n}^{\alpha}{\cal
I}_{s}^{\alpha}\otimes I_{m}{\bm{\gamma}}^{n}\right),$ (32)

with, by setting $\gamma_{j}^{n}(\rho_{\ell})$, $j=0,\dots,s-1$, the
approximation to $\gamma_{j}(\sigma_{n},\rho_{\ell})$ obtained by using the
Gauss-Jacobi quadrature formula,

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

and

${\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),~{}{\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)\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}.$

###### Remark 2

It is worth noticing that the discrete problem (32) has (block) dimension $s$,
independently of $k$. This, in turn, allows using relatively large values of
$k$, in order to have an accurate approximation of the Fourier coefficients,
without increasing too much the computational cost.

Moreover, the vector ${\bm{\phi}}_{n-1}^{\alpha}$ in (32) only depends on
known quantities, computed at the previous timesteps.

Further, we observe that also the matrices ${\cal P}_{s}$, ${\cal
I}_{s}^{\alpha}$, as well as all the required integrals (29), can be computed
in advance, once for all, and they can be used for each new approximation
$\rho_{\ell}$ in Algorithm 1.

Considering that

${\cal
I}^{\alpha}P_{j}(1)=\frac{1}{\Gamma(\alpha)}\int_{0}^{1}(1-x)^{\alpha-1}P_{j}(x)\mathrm{d}x=\frac{\delta_{j0}}{\Gamma(\alpha+1)},\qquad
j=0,\dots,s-1,$

from (25) and (26) one then obtains:

$y(T,\rho_{\ell})\simeq\sigma_{N}(h_{N},\rho_{\ell})\equiv\phi_{N-1}^{\alpha}(1,\rho_{\ell})+\frac{h_{N}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}^{N}(\rho_{\ell}).$
(33)

###### Definition 1

We shall refer to the method defined by the previous steps, as a Fractional
HBVM with parameters $k$ and $s$, in short FHBVM$(k,s)$.

###### Remark 3

When $\alpha=1$, the polynomials $\\{P_{j}\\}_{\geq 0}$, become the usual
Legendre polynomials orthomormal in $[0,1]$. Consequently, a FHBVM$(k,s)$
method reduces to a standard HBVM$(k,s)$ method, when $\alpha=1$.

In a similar way, for $n=1,\dots,N$:

$\Psi_{n}(ch_{n},\rho_{\ell})=\Theta_{n-1}^{\alpha}(c,\rho_{\ell})+h_{n}^{\alpha}\sum_{j=0}^{s-1}\Gamma_{j}(\sigma_{n},\rho_{\ell})I^{\alpha}P_{j}(c),\qquad
c\in[0,1],$ (34)

with

$\displaystyle\Gamma_{j}(\sigma_{n},\rho_{\ell})$ $\displaystyle=$
$\displaystyle\alpha\int_{0}^{1}(1-\tau)^{\alpha-1}P_{j}(\tau)f^{\prime}(\sigma_{n}(\tau
h_{n},\rho_{\ell}))\Psi_{n}(\tau h_{n},\rho_{\ell})\mathrm{d}\tau,$ (35)
$\displaystyle\qquad\qquad j=0,\dots,s-1,$

and (see (29))

$\displaystyle\Theta_{n-1}^{\alpha}(c,\rho_{\ell})~{}=$ (36) $\displaystyle
I+\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}h_{\nu}^{\alpha}\sum_{j=0}^{s-1}J_{j}^{\alpha}\left(\frac{r^{n-\nu}-1}{r-1}+cr^{n-\nu}\right)\Gamma_{j}(\sigma_{\nu},\rho_{\ell}).$

Similarly as in (30), when a constant stepsize $h=T/N$ is used, then (36)
becomes:

$\Theta_{n-1}^{\alpha}(c,\rho_{\ell})=I+\frac{h^{\alpha}}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}\sum_{j=0}^{s-1}J_{j}^{\alpha}(n-\nu+c)\Gamma_{j}(\sigma_{\nu},\rho_{\ell}).$
(37)

As done for (27), by approximating the integrals in (35) by using the same
Gauss-Jacobi formula as before, from (34) and (35) one derives a discrete
problem in the form

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

where ${\bm{\phi}}_{n-1}^{\alpha}$ and ${\bm{\gamma}}^{n}$ have been already
computed in (32),

$f^{\prime}\left({\bm{\phi}}_{n-1}^{\alpha}+h_{n}^{\alpha}{\cal
I}_{s}^{\alpha}\otimes I_{m}{\bm{\gamma}}^{n}\right)\in\mathbb{R}^{km\times
km}$

is the block diagonal matrix whose diagonal blocks are given by the
corresponding evaluations of the Jacobian of $f$,

${\bm{\Theta}}_{n-1}^{\alpha}=\left(\begin{array}[]{c}\Theta_{n-1}^{\alpha}(c_{1},\rho_{\ell})\\\
\vdots\\\
\Theta_{n-1}^{\alpha}(c_{k},\rho_{\ell})\end{array}\right)\in\mathbb{R}^{km\times
m},$

and, by setting $\Gamma_{j}^{n}(\rho_{\ell})$, $j=0,\dots,s-1$, the
approximation to $\Gamma_{j}(\sigma_{n},\rho_{\ell})$ obtained through the
Gauss-Jacobi formula,

${\bm{\Gamma}}^{n}=\left(\begin{array}[]{c}\Gamma_{0}^{n}(\rho_{\ell})\\\
\vdots\\\
\Gamma_{s-1}^{n}(\rho_{\ell})\end{array}\right)\in\mathbb{R}^{sm\times m}.$

We observe that considerations similar to those made in Remark 2 for (32) can
be now repeated for (38), with the additional fact that (38) amounts to just
solving a linear system of equations.

At last, from (25) and (34) one eventually obtains:

$\Phi(T,\rho_{\ell})\simeq\Psi_{N}(h_{N},\rho_{\ell})\equiv\Theta_{N-1}^{\alpha}(1,\rho_{\ell})+\frac{h_{N}^{\alpha}}{\Gamma(\alpha+1)}\Gamma_{0}^{N}(\rho_{\ell}).$
(39)

###### Remark 4

By choosing values of $s$, and $k\geq s$, large enough, it can be seen, as is
shown in [6], that the approximations (33) and (39) provided by a FHBVM$(k,s)$
method can be accurate up to machine precision. As matter of fact, this
amounts to use the method as a spectral methods in time. This kind of
approximations will be considered in the implementation of the algorithm
listed in Table 1, which we shall use for the numerical tests reported in
Section 4.

###  3.2 Error estimation

It is worth mentioning that the procedure explained in the previous section
allows to derive, as a by-product, an estimate for the error in the computed
solution, due to the fact that, in Algorithm 1, the iteration is stopped when,
for a suitably small tolerance $tol$,

$|\rho_{\ell+1}-\rho_{\ell}|\leq tol.$ (40)

In fact, in such a case, one expects that $|\rho_{\ell+1}-\rho^{*}|\approx
tol$ as well. Consequently, by considering that at the mesh points, for
$n=1,\dots,N$:

$y(t_{n},\rho_{\ell})\simeq\sigma_{n}(h_{n},\rho_{\ell})\equiv\phi_{n-1}^{\alpha}(1,\rho_{\ell})+\frac{h_{n}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}^{n}(\rho_{\ell}),$
(41)

and, similarly,

$\Phi(t_{n},\rho_{\ell})\simeq\Psi_{n}(h_{n},\rho_{\ell})\equiv\Theta_{n-1}^{\alpha}(1,\rho_{\ell})+\frac{h_{n}^{\alpha}}{\Gamma(\alpha+1)}\Gamma_{0}^{n}(\rho_{\ell}),$

by virtue of the perturbation result of Theorem 1, one derives the estimates

$\|y(t_{n},\rho^{*})-y(t_{n},\rho_{\ell})\|\approx 2\cdot
tol\cdot\|\Psi_{n}(h_{n},\rho_{\ell})\|,\qquad n=1,\dots,N.$ (42)

##  4 Numerical Tests

We here report a few numerical tests to illustrate the theoretical findings.
For all tests, we use $k=22$ and $s=20$, so that we are going to use a
FHBVM(22,20) method. In other words, we use a local polynomial approximation
of degree $s-1=19$ for the vector field, coupled with a Gauss-Jacobi
quadrature formula of order $2k=44$ for approximating the Fourier coefficients
(27) and (35). We have used straightforward fixed-point iterations, derived
from (32) and (38), respectively, to solve the corresponding discrete
problems.999We have used a fixed point iteration also for solving (38),
despite the fact that it is just a linear system of equations. Namely,

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

starting from ${\bm{\gamma}}^{n,0}={\bm{0}}$, for (32), and similarly for
(38). The iterations are carried out until full machine accuracy is gained, so
that we expect full machine accuracy for the computed approximation (39) to
$\Phi(T,\rho_{\ell})$, as well as a corresponding fully accurate discrete
solution (41). We consider the same test problems in [16, Section 5]. For all
problems, (see (1)) the initial guess $\rho_{0}=\eta$ has been considered. All
numerical tests have been performed in Matlab© (Rel. 2022b) on a 3GHz Intel
Xeon W10 computer with 64GB of memory. The iteration of Algorithm 1 is stopped
by using a tolerance $tol=10^{-14}$ in (40).

###  4.1 Example 1

The first problem is given by:

$\displaystyle y^{(0.3)}$ $\displaystyle=$
$\displaystyle-|y|^{1.5}+\frac{8!}{\Gamma(8.7)}t^{7.7}-3\frac{\Gamma(5.15)}{\Gamma(4.85)}t^{3.85}+\left(\frac{3}{2}t^{0.15}-t^{4}\right)^{3}+\frac{9}{4}\Gamma(1.3),$
(43) $\displaystyle\qquad t\in[0,1],\qquad y(1)=\frac{1}{4},$

whose solution is

$y(t)=t^{8}-3\,t^{4.15}+\frac{9}{4}\,t^{0.3}.$ (44)

In this case, we use a uniform mesh with stepsize $h=1/50$. The method
converges in 4 iterations producing the following approximations:

$\ell$ | $\rho_{\ell}$
---|---
0 | 2.500000000000000e-01
1 | -6.974187065949766e-03
2 | -6.268019571770665e-06
3 | -5.034525913700941e-12
4 | 6.626717195461731e-15

It is possible to appreciate the quadratic convergence of the iteration in the
first iterations (in the last one, roundoff errors clearly dominate). The
maximum error on the final solution is $\approx 6\cdot 10^{-15}$, whereas the
estimated one, by using (42), is $2\cdot 10^{-14}$.

###  4.2 Example 2

The second problem is given by:

$\displaystyle y^{(0.3)}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\,y,\qquad t\in[0,7],$ $\displaystyle y(7)$
$\displaystyle=$ $\displaystyle\frac{14}{5}E_{0.3}\left(-1.5\cdot
7^{0.3}\right)\simeq 0.6476128469955936,$ (45)

with $E_{0.3}$ the Mittag-Leffler function of order 0.3, with solution

$y(t)=\frac{14}{5}E_{0.3}\left(-1.5\cdot t^{0.3}\right).$ (46)

We refer to [19] and the accompanying software ml.m, for an efficient Matlab©
implementation of the Mittag-Leffler function.

In this case a uniform mesh is not appropriate, since the vector field is
proportional to the solution, which has a singularity in the first derivative
at the origin. Consequently, we use a graded mesh, according to (16), with
$h_{1}=10^{-14}$ and $N=500$. Taking into account (17), this implies $r\simeq
1.064914852480467$. The method converges in 2 iterations producing the
following approximations:

$\ell$ | $\rho_{\ell}$
---|---
0 | 0.6476128469955936
1 | 2.799999999999973
2 | 2.799999999999984

The maximum error on the final solution is $\approx 2\cdot 10^{-13}$, whereas
the estimated one, by using (42), is $2\cdot 10^{-14}$ (in this case, the
maximum error is essentially close to the origin, where there is the
singularity of the derivative).

###  4.3 Example 3

The last problem is given by:

$\displaystyle y^{(0.7)}$ $\displaystyle=$
$\displaystyle\frac{1}{t+1}\sin(t\cdot y),\qquad t\in[0,20],$ $\displaystyle
y(20)$ $\displaystyle=$ $\displaystyle 0.8360565285776644,$ (47)

which corresponds to the initial value $y(0)=1$. In such a case, the solution
is not known in closed form, and the final value has been taken from a
reference solution computed by using the FHBVM(22,20) method with a constant
stepsize $h=0.02$ (i.e., by using 1000 timesteps). This solution is depicted
in Figure 1, and the estimated error (by using a doubled mesh) is $\approx
1.8\cdot 10^{-14}$.

For solving problem (47), we use a uniform mesh with stepsize $h=20/400=1/20$.
The method converges in 6 iterations, producing the following approximations:

$\ell$ | $\rho_{\ell}$
---|---
0 | 0.8360565285776644
1 | 1.115178544783083
2 | 1.057854760373078
3 | 1.006528883050734
4 | 0.9999714859685486
5 | 0.9999999991678460
6 | 0.9999999999999847

Also in the case, it is possible to appreciate a quadratic-like convergence of
the iteration. The maximum error in the final solution is $\approx 2\cdot
10^{-14}$, whereas the estimated one, by using (42), is $6\cdot 10^{-14}$.

##  5 Conclusions

In this paper we have described a novel shooting procedure which, coupled with
the Newton method, proves very appealing for numerically solving fractional
terminal value problems. The implementation details of the given procedure
have been thoroughly given, when the underlying numerical methods are FHBVMs.
These latter methods, when used as spectral methods in time, allow deriving
very accurate solutions, along with a suitable estimate of the error in the
computed solution. Numerical tests on scalar problems confirm the
effectiveness of the approach that, however, could be extended to cope with
vector two-point boundary value problems, too: this will be indeed the subject
of future investigations. Also the efficient numerical solution of the
generated discrete problems needs further investigation.

Figure 1: Reference solution for problem (47).

## References

  * [1] 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
  * [2] 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
  * [3] 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
  * [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] 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)
  * [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, 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
  * [12] 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
  * [13] 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.
  * [14] K. Diethelm. Increasing the efficiency of shooting methods for terminal value problems of fractional order. J. Comput. Phys. 293 (2015) 135–141. https://doi.org/10.1016/j.jcp.2014.10.054
  * [15] 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
  * [16] K. Diethelm, F. Uhlig. A new approach to shooting methods for terminal value problems of fractional differential equations. J. Sci. Comput. 97 (2023) 38. https://doi.org/10.1007/s10915-023-02361-9
  * [17] N.J. Ford, M.L. Morgado. Fractional boundary value problems: analysis and numerical algorithms. Fract. Calc. Appl. Anal. 14 (2011) 554–567. https://doi.org/10.2478/s13540-011-0034-4
  * [18] N.J. Ford, M.L. Morgado, M. Rebelo. High order numerical methods for fractional terminal value problems. Comput. Methods Appl. Math. 14 (2014) 55–70. https://doi.org/10.1515/cmam-2013-0022
  * [19] 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
  * [20] 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.
  * [21] 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
  * [22] 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
  * [23] 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
  * [24] 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.
  * [25] 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
  * [26] V. Lakshmikantham, D. Trigiante. Theory of Difference Equations: Numerical Methods and Applications. Academic Press Inc., Boston, 1988.