Mediterr. J. Math.

(2018) 15:113
https://doi.org/10.1007/s00009-018-1149-1
c Springer International Publishing AG,
part of Springer Nature 2018

Stability of Numerical Solutions

for Abel–Volterra Integral Equations
of the Second Kind
G. Izzo , E. Messina and A. Vecchio

Abstract. We analyze the stability of convolution quadrature methods

for weakly singular Volterra integral equations with respect to a linear
test equation. We prove that the asymptotic behavior of the numerical
solution replicates the one of the continuous problem under some re-
striction on the stepsize. Numerical examples illustrate the theoretical
results.
Mathematics Subject Classification. Primary 34A08 , 65L06; Secondary
65R20.
Keywords. Weakly singular integral equations, Numerical stability,
Convolution quadrature.

1. Introduction
We consider weakly singular Volterra integral equations (VIEs) of the type
 t
1
y(t) = f (t) + (t − s)α−1 g(s, y(s))ds, t ≥ 0. (1.1)
Γ(α) 0
The parameter α is assumed to take values in (0, 1] and defines the weak
singularity. We suppose that the forcing term f (t) is at least continuous for
t ≥ 0, and that g is differentiable. Existence and uniqueness results can be
found, for example, in [2].
The factor 1/Γ(α), where Γ is the Euler Gamma function, emphasizes
that (1.1) represents the integral formulation of a fractional differential equa-
tion
 α
D y(t) = G(t, y(t))
, t ≥ 0,
y(0) = y0
1
 t y (s)
where α is the order, Dα y(t) = Γ(1−α) 0 (t−s)α
ds is the fractional derivative
Caputo operator, G(t, y) = Dα f (t) + g(t, y) and y0 = f (0) (the monographs

G. Izzo, E. Messina and A. Vecchio are members of the INdAM Research group GNCS.

0123456789().: V,-vol
113 Page 2 of 14 G. Izzo et al. MJOM

of Diethelm [7] and Podlubny [23] provide a detailed analysis on the theory
and applications of fractional differential equations).
Our interest for Eq. (1.1) is mainly motivated by this equivalence and its
impact in the applications (see, e.g., [4–7,9,14], and the bibliography therein).
Since the stability analysis is crucial in this context, we mainly focus on
a comparative study about the asymptotic behavior of the analytical and
numerical solutions. Hence, we consider the following test equation
 t
1
y(t) = f (t) + (t − s)α−1 k(s)y(s)ds, t ≥ 0, (1.2)
Γ(α) 0
whose relevance in the applications stands also in the fact that it may be
representative of the behavior of non-linear problems of the kind (1.1) with
k(t) > 0 when, for example, |g(t, x)| ≤ k(t)|x|. We assume that
 t
1
sup (t − s)α−1 |k(s)|ds = η < 1, (1.3)
Γ(α) t≥t̄ t̄
for some t̄ ≥ 0, and we prove the stable behavior of the numerical solution
under some conditions on the discretization stepsize.
Although in [20,21] and [22] a similar analysis has already been carried out
for VIEs with regular kernels, the extension to Abel type equations is not
straightforward and deserves a separate investigation. In [12] the stability of
Eq. (1.2) has been analyzed with respect to (1.3) and t̄ = 0, for trapezoidal
product integration. Here we want to extend this result to the class of con-
volution quadrature methods [19] and for any non-negative value of t̄. While
for t̄ = 0 the asymptotic behavior of the analytical solution to (1.2) is well-
known (see [2]), for t̄ > 0 the boundedness of y(t) is not obvious. Thus, in
Sect. 2, we study the stability of the model problem for general t̄ ≥ 0. Most
significant for the rest of the paper is a property, that we prove in Sect. 3,
which relates the convolution quadrature method to the model problem. In
Sect. 4 we report our analysis on the numerical solution and state the sta-
bility result. Some numerical experiments are described in Sect. 5 to validate
the theoretical findings and in Sect. 6 an example of application is described.
Finally Sect. 7 concludes the paper.

2. Stability Analysis of the Model Problem

Here and in the following we will refer to (1.2) with (1.3) as to the model
problem, since it will be the subject of our investigations for the stability
analysis. From now on we assume that
sup |k(t)| = K < +∞, (2.1)
t≥0

and
sup |f (t)| = F < +∞. (2.2)
t≥0

Theorem 2.1. For problem (1.2), (1.3) assume that (2.1) and (2.2) hold, then
sup |y(t)| = Y < +∞. (2.3)
t≥0
MJOM Stability of Numerical Solutions for Abel–Volterra Integral Page 3 of 14 113

Proof. Consider 0 ≤ t < t̄. From Eq. (1.2), using (2.1) and (2.2), we get
 t
K
|y(t)| ≤ F + (t − s)α−1 |y(s)|ds,
Γ(α) 0
and thus, using Lemma 1.3.13 in [3], we get

|y(t)| ≤ Eα (K t̄α )F, ∀t ∈ [0, t̄), (2.4)

where Eα denotes the Mittag–Leffler function defined, for α > 0, by

∞
 zn
Eα (z) = .
n=0
Γ(nβ + 1)

Now let t ≥ t̄ and split Eq. (1.2) as follows

 t̄  t
1 1
y(t) = f (t) + (t − s)α−1 k(s)y(s)ds + (t − s)α−1 k(s)y(s)ds.
Γ(α) 0 Γ(α) t̄
Thus, using (2.4), (2.1) and (2.2), we get
 t̄  t
K
|y(t)| ≤ F + Eα (K t̄α )F (t − s)α−1 ds + sup |y(t)| (t − s)α−1 |k(s)|ds.
Γ(α) 0 t≥t̄ t̄

Since (1.3) holds,

 
F K α α
sup |y(t)| ≤ 1+ Eα (K t̄ )t̄ , (2.5)
t≥t̄ 1−η Γ(α + 1)

and (2.3) comes true. 

Remark 1. The bound (2.3) with (2.5) has a theoretical meaning and serves
to give a qualitative description of the solution to the model problem (1.2),
(1.3). It is clear that, when t̄ = 0, (2.5) reduces to the bound proved in [12],
supt≥0 |y(t)| ≤ F/(1 − η).
The following result on the limit of y(t) as t → +∞ is similar to the
one obtained in [15] for discrete equations and in [20] for general Volterra
equations on time scales. However, the weak singularity in the kernel of (1.2)
requires some additional details, which motivate us to report the proof.

Theorem 2.2. Consider Eqs. (1.2), (1.3) and assume that (2.1) holds, fur-
thermore let

lim f (t) = f∞ ,
t→+∞

and
  t 
 1 
lim lim sup  (t − s)α−1 k(s)ds − Ik  = 0, for some Ik ,
τ →+∞ t→+∞ Γ(α) τ

then limt→+∞ y(t) = y∞ , with y∞ = f∞ /(1 − Ik ).

113 Page 4 of 14 G. Izzo et al. MJOM

Proof. Since (1.3) holds, it is Ik ≤ η < 1. Let y∞ = f∞ + Ik y∞ , then from

Eq. (1.2) we have
 τ
1
|y(t) − y∞ | ≤ |f (t) − f∞ | + (t − s)α−1 |k(s)||y(s)|ds
Γ(α) 0
 t
1
+ (t − s)α−1 |k(s)||y(s) − y∞ |ds
Γ(α) τ
  t 
 1 
+  (t − s)α−1
k(s)ds − Ik  |y∞ |.
Γ(α) τ
We are in the hypotheses of Theorem 2.1, hence y(t) is bounded. Furthermore,
τ
for t ≥ τ̄ > τ, and 0 ≤ s ≤ τ, (t − s)α−1 is bounded so limt→+∞ 0 (t −
s)α−1 |k(s)||y(s)|ds = 0 and thus, for any τ > t̄,
lim sup |y(t) − y∞ | ≤ lim |f (t) − f∞ | + η sup |y(t) − y∞ |
t→+∞ t→+∞ t≥τ
  t 
 1 
+ |y∞ | lim sup  (t − s)α−1 k(s)ds − Ik  .
t→+∞ Γ(α) τ
Passing to the limit as τ → +∞ we get the result. 

3. Continuous vs. Discrete

Convolution quadrature methods for problem (1.1) read
n

yn = fn + hα ωn−j kj g(tj , yj ), n ≥ 0, (3.1)
j=0
−1
with fn = f (tn ) + hα j=−n0 wnj kj g(tj , yj ). Here, ωn are the convolution
weights, tn = t0 + nh, n = 0, 1 . . . , with t0 = n0 h, and h positive constant,
kj = k(tj ) and yn ≈ y(tn ). The values y−n0 , . . . , y−1 are given and wnj are
the starting weights.
The class (3.1) of convolution quadrature methods includes, for example,
product integration methods [10,13,16], fractional linear multistep methods
[1,11,17]. For convergence and for stability w.r.t. basic test equation we refer
to [17–19].
In Sect. 2 we have proved that, if supt≥0 |f (t)| < +∞, the solution y(t)
to the model problem (1.2), (1.3) satisfies supt≥0 |y(t)| < +∞. To investigate
how hypothesis (1.3), which detects the boundedness of the analytical solu-
tion, influences the numerical solution, we solve the linear equation (1.2) by
the convolution quadrature (3.1)
n

yn = fn + hα ωn−j kj yj , n ≥ 0, n = 0, 1 . . . , (3.2)
j=0
−1
with fn = f (tn )+hα j=−n0 wnj kj yj , and emphasize the following properties
of the weights (see for example [19]) which will be crucial in the subsequent
analysis:
• wnj = O(nα−1 ) as n → ∞;
MJOM Stability of Numerical Solutions for Abel–Volterra Integral Page 5 of 14 113

• ωn are the convolution weights for which

∞

nα−1
ωn = + un , n≥1 with |un | < ∞. (3.3)
Γ(α) n=1

From here it is clear that, when 0 < α < 1, fn is bounded when f (t) is
bounded, limn→+∞ fn = limt→+∞ f (t) and that

sup |ωn | < W < +∞. (3.4)

n≥0

We prove the following result which links the continuous and discrete prob-
lems.

Lemma 3.1. Assume that, for Eqs. (1.2), (2.1) holds, furthermore let K̄ =
supt≥t̄ |k(t)|, t̄ ≥ 0. Then, we have
n
  t
α 1
sup h |ωn−j ||kj | ≤ sup (t − s)α−1 |k(s)|ds + hα A, (3.5)
n≥r j=r
Γ(α) t≥t̄ t̄

for any r ≥ 0 and h > 0 such that tr = t0 + rh = t̄, where

 t 
∞ 
1 

A= K̄ + sup |k (x)|dx + |un | + |ω0 | + W K̄. (3.6)
Γ(α) t≥t̄ t̄ n=1

Proof. Consider
 tj  tj
1 h 1
(tn − s)α−1 k(s)ds = (tn − tj )α−1 kj − (tn − tj )α−1 kj ds
Γ(α) tj−1 Γ(α) Γ(α) tj−1
 tj
1
+ (tn − s)α−1 k(s)ds
Γ(α) tj−1
h 1
= (tn − tj )α−1 kj −
Γ(α) Γ(α)
 tj

× (tn − tj )α−1 kj − (tn − s)α−1 k(s) ds
tj−1

h 1
= (tn − tj )α−1 kj −
Γ(α) Γ(α)
 tj  tj
 
× (tn − x)α−1 k(x) dxds,
tj−1 s

then, by exploiting the form (3.3) for the convolution weights, we have
 tj
1 1
(tn − s)α−1 k(s)ds = hα ωn−j kj − hα un−j kj −
Γ(α) tj−1 Γ(α)
 tj  tj
 
× (tn − x)α−1 k(x) dxds.
tj−1 s
(3.7)
113 Page 6 of 14 G. Izzo et al. MJOM

From this,
 tj
α 1
h |ωn−j ||kj | ≤ (tn − s)α−1 |k(s)|ds + hα |un−j ||kj |
Γ(α) tj−1
 tj
h
+ (tn − x)α−1 |k  (x)|dx
Γ(α) tj−1
 tj 
h   
+  (tn − x)α−1  |k(x)|dx.
Γ(α) tj−1

Hence,
n
  tn−1 n−1

1
hα |ωn−j ||kj | ≤ (tn − s)α−1 |k(s)|ds + hα |kj ||un−j |
j=r
Γ(α) tr j=r+1
 tn−1
h
+ (tn − x)α−1 |k  (x)|dx
Γ(α) tr
 tn−1
h
+ (1 − α)(tn − x)α−2 |k(x)|dx
Γ(α) tr
+ hα |ω0 ||kn | + hα |ωn−r ||kr |.
Since x ∈ [0, tn−1 ], tn − x ≥ h, so (tn − x)α−1 ≤ hα−1 . Thus,

n  tn−1 
n−1
1
hα |ωn−j ||kj | ≤ (tn − s)α−1 |k(s)|ds + hα |kj ||un−j |
j=r
Γ(α) tr j=r+1
 tn−1
hα hα
+ |k (x)|dx + K̄ + hα |ω0 ||kn | + hα |ωn−r ||kr |.
Γ(α) tr Γ(α)
The result comes straightforward taking into account the form of A in (3.6).


4. Numerical Stability
Through Lemma 3.1 the main assumption (1.3) on the model problem trans-
lates into an analogous hypothesis on the discrete equation (3.2) if we assume
that  t
sup |k  (s)|ds < +∞, (4.1)
t≥t̄ t̄
and choose the stepsize h ∈ (0, h∗ ], where
h∗ = sup {h ≥ 0 s.t. ηh < 1} , (4.2)
α
with ηh = η + Ah , η and A defined in (1.3) and (3.6), respectively. Now, it
is possible to prove the following result on the boundedness of the numerical
solution yn .
Theorem 4.1. Assume that, for problem (1.2), hypotheses (1.3), (2.1), (2.2)
and (4.1) hold, then for the solution yn of (3.2), we have
sup |yn | ≤ Yh (4.3)
n≥0
MJOM Stability of Numerical Solutions for Abel–Volterra Integral Page 7 of 14 113

for any sufficiently small h.

Proof. Let 0 < h < h∗∗ , where
   α1 
1
h∗∗ ≤ min h∗ , , (4.4)
ω0 K
and h∗ defined in (4.2). Let n̄ = n̄(h) such that t̄ = n̄h, and consider 0 ≤ n <
n̄. From Eq. (3.2), using (2.1), (2.2) and (3.4), we have
n−1

F̄ α WK
|yn | ≤ α
+ h α
|yj |,
1 − h ω0 K 1 − h ω0 K j=0

where F̄ = supn≥0 |fn |. The discrete Gronwall inequality, stated, for example,
in [16, Th.7.1], yields
nhα W K F̄
|yn | ≤ e 1−hα ω0 K . (4.5)
1− hα ω 0 K, 0 ≤ n < n̄
For n̄ ≤ n ≤ N we split (3.2) to have
n̄−1
 n

yn = fn + hα ωn−j kj yj + hα ωn−j kj yj . (4.6)
j=0 j=n̄

We use (4.5) for the first sum in (4.6), thus

n̄−1
 n̄hα W K F̄
|ωn−j ||kj ||yj | ≤ n̄hα W Ke 1−hα ω0 K .
j=0
1 − hα ω0 K
By the previous expression and Lemma 3.1, for any fixed N ∈ N the following
bound for yn holds,
n̄hα W K
|yn | ≤ F̄ + F̄ e 1−hα ω0 K + ηh sup |yn |.
n̄≤n≤N
∗∗ ∗∗
For h ∈ (0, h ), with h defined in (4.4) we are then allowed to write
 
F̄ n̄hα W K 1−h n̄hα W K
sup |yn | ≤ 1+ e αω K
0 . (4.7)
n̄≤n≤N 1 − ηh 1 − hα ω0 K
From here and the arbitrariness of N ∈ N, we get (4.3). 
Since Eq. (3.2) is linear, Theorem 4.1 gives also information about the
behavior of the global error of the discretization. In this case, fn plays the
role of the local truncation error Tn (h). The bound for the truncation error
due the discretization of the integral term in (1.1) by convolution quadrature
is well-known when integrating over a finite interval [17,18]. However, since
we are interested in the long-time behavior of the numerical solution, we can
prove that for problem (1.2), (1.3) it is |Tn (h)| ≤ chα , under hypotheses
t
(2.1), (2.2), (4.1) and supt≥t̄ t̄ |k(x)y  (x)|dx < +∞, t̄ ≥ 0. Thus, for any
h ∈ (0, h∗∗ ] the global error en satisfies
sup |en | ≤ Chα
n≥0

and can be made arbitrarily small with h.

113 Page 8 of 14 G. Izzo et al. MJOM

In analogy to Theorem 2.2 we have the following result.

Theorem 4.2. Consider Eqs. (1.2), (1.3) and assume that (2.1) and (4.1)
hold, furthermore let
lim f (t) = f∞ ,
t→+∞

and
 
 n 
 α 

lim lim sup h h
ωn−j kj − Ik  = 0, for some Ikh ,
r→+∞ n→+∞ 
 j=r 
h
then for the solution yn of (3.2), we have that limn→+∞ yn = y∞ , with
h h
y∞ = f∞ /(1 − Ik ).
The proof comes straightforward also taking into account the properties
of the weights mentioned at the beginning of this section.
Since
  t 
 1 
|Ik − Ikh | ≤ Ik − (t − s)α−1 k(s)ds
Γ(α) τ
  t   
 n    n 
 1   
+ (t − s)α−1 k(s)ds − hα ωn−j kj  + hα ωn−j kj − Ikh  ,
 Γ(α) τ   
j=r j=r

using (3.7) in Lemma 3.1 it is possible to prove that, assuming that the
K
hypotheses of Theorems 2.2 and 4.2 hold, |Ik − Ikh | ≤ hα (A + Γ(α+1) ). Thus,
h
y∞ → y∞ as h → 0, allowing a sharp estimate of the limit value of the
analytical solution.

5. Numerical Experiments
For the following experiments, we use the fractional trapezoidal rule (see
[10]). Our first test consists in Eq. (1.2) with
k(t) = p exp (−λt), and
f (t) = D1−α sin 2πωt. (5.1)
Here we choose α = 0.8, p = 10, λ = 5, and ω = 1. The main assumption
(1.3) is accomplished for t̄ > 0.15 and the result stated in Theorems 2.1 and
4.1 hold. The nature of the kernel allows us to choose t̄ arbitrarily large to
get an arbitrarily small A in (3.6) and, in fact, no restriction on h due to
h∗ in (4.2). On the contrary, we cannot ignore the effect of h∗∗ in (4.4). h∗∗
 1/α
is at least bounded by ω01K . Since K = 10 and ω0 ≈ 0.92, we choose
h < 0.11. Figure 1 shows the expected bounded behavior of the numerical
solution for different values of the stepsize h. Next we consider Eq. (1.2) with
p
k(t) = , and
(t + 1)2
t
f (t) such that y(t) = . (5.2)
1 + 2t
MJOM Stability of Numerical Solutions for Abel–Volterra Integral Page 9 of 14 113

18
h=0.05
16 h=0.025
h=0.0125
14

12

10
y

0
0 2 4 6 8 10
t

Figure 1. Problem (5.1), numerical solution for different val-

ues of h

0.498

0.496

0.494

0.492

0.49
y

0.488

0.486
h=0.005
h=0.0025
0.484 h=0.00125
y(t)
0.482
40 50 60 70 80 90 100
t

Figure 2. Problem (5.2), numerical solution for different val-

ues of h

Here we choose α = 0.5, p = 5 and h < 0.08 to gain a stable behavior.

The numerical solution, computed for different values of the stepsize h, and
the analytical one y(t) are plotted in Fig. 2, allowing a comparison which,
according to Theorem 4.1, confirms numerical stability and, as expected from
Theorems 2.2 and 4.2, the same asymptotic limit, as h → 0. For the sake
of completeness, we report in Table 1 the number of correct digits of the
113 Page 10 of 14 G. Izzo et al. MJOM

Table 1. Error, correct digits and experimental order of the

numerical solution to problem (1.2), (5.2) with h =
0.05, 0.025, 0.0125

h Error Corr. dig. Exp. ord.

0.05 4.692e−2 1.33 –
0.025 6.861e−3 2.16 2.77
0.0125 1.175e−3 2.93 2.55

numerical solution with respect to the analytical one. It is clear that the
method behaves with order 2.
Our third experiment consists in a non-linear problem (1.1) with
y3
(a) g(t, y) = pe−λt ,
y2
+1
(b) g(t, y) = arctan(pe−λt y), (5.3)
and f (t) = 1/(1+t) . Both in case (a) and (b), |g(t, y)| ≤ pe |y|. So, pe−λt
2t −λt

can be assumed as the regular part k(t) of the kernel in the model problem
(1.2), (1.3). In the experiments we choose λ = 2, α = 0.5, p = 2 (case (a))
and p = 5 (case (b)). The bound for the stepsize h is 0.5 (case (a)) and 0.08
(case (b)); however, even larger stepsizes produce stable numerical solutions,
confirming the fact that the condition proved in Sect. 4 is only sufficient.
t
Furthermore, according to Theorem 2.2, since limt→+∞ τ (t − s)α−1 k(s)ds =
0 and thus Ik = 0, it is limt→+∞ y(t) = limt→+∞ f (t) = 0. Figure 3 shows
stable behavior of numerical solutions in the cases (a) and (b), respectively,
and the correct limit value of the numerical solution, which approximates
y∞ = f∞ = 0, is clear in the figure.

6. Applications
While, as already pointed out before, equations of the type (1.1) and (1.2)
represent significative mathematical models in a variety of fields, ranging
from mechanics, physics, control theory (see once again [7,14] and the bib-
liography therein), hypothesis (1.3) is quite restrictive and its applicability
to significative fields is scarce. However, the numerical stability with respect
to non-parametric test equations is quite an unexplored subject of study in
the context of weakly singular VIEs and the research described here may
be considered as a first attempt in this direction. Furthermore, the tech-
nique developed in the previous sections may be exploited in cases that are
closer to real life problems, as for example, in certain applications of frac-
tional calculus to viscoelasticity (see [24]), where one has that k(t) is sin-
gular at t = 0 and 0 < k(t) ≤ Kt−α , for t > 0, in (1.2). Hypothesis (2.1)
is not satisfied anymore, furthermore the uniqueness of the solution is not
guaranteed. To specify a unique trajectory, we assume that the value y(t̄)
at a certain time point t̄ > 0 is given (see [8]). Assume that f (t) ≥ 0 is
a smooth function such that supt≥t̄ ft(t)α < +∞ and consider t ≥ t̄, since
MJOM Stability of Numerical Solutions for Abel–Volterra Integral Page 11 of 14 113

2
y

0
0 200 400 600 800 1000
t
2

1.5

1
y

0.5

0
0 200 400 600 800 1000
t

Figure 3. Problem (5.3), a upper and b lower, numerical

solution with h = 0.05

1
 t̄
Γ(α) 0
(t̄ − s)α−1 k(s)y(s)ds = y(t̄) − f (t̄), it is simple to prove that, if
K
γK = < 1, (6.1)
Γ(α + 1)
then
y(t) F̃
α
≤ sup , (6.2)
t≥t̄ t 1 − γK
  t̄ 
with F̃ = supt≥t̄ ft(t) 1
α + Γ(α)t̄α 0
(t − s)α−1 k(s)y(s)ds , assuring that the
growth of the analytical solution is dominated by tα . This result holds for
any particular solution of Eq. (1.2).
Now we prove that the same holds for any numerical approximation defined
by Eq. (3.1). We first notice that Lemma 3.1 can be adapted to our case to
give,
n

α 1
h ωn−j ≤ (tn − t̄)α + hα B,
j=n̄
Γ(α + 1)

for
+∞ any n̄ ≥ 1 and h > 0 such that tn̄ = t0 + n̄h = t̄ > 0, where B =
n=1 |un | + |ω0 | + W .
Thus, for n ≥ n̄,
 
yn K (tn − t̄)α α KB yn yn
≤ F̄ + +h α sup α ≤ F̄ + γK (h) sup α ,
tαn Γ(α + 1) t α
n t n t
n≥n̄ n n≥n̄ tn
 n̄−1 
K fn
with γK (h) = Γ(α+1) + hα KB t̄α and F̄ = sup n≥n̄ α
tn + h αW
t̄α j=0 k(t j )y j .
Since (6.1) holds, it is always possible to choose the stepsize h such that
113 Page 12 of 14 G. Izzo et al. MJOM

1
0 10 20 30 40 50

Figure 4. Numerical solution yn to problem (1.2) with

α = 0.5, k(t) = sin2 (t)/tα and f (t) = tα , divided by tα
n

γK (h) < 1 and thus the following bound for the numerical approximation
holds
yn F̄
sup α ≤ ,
n≥n̄ tn 1 − γK (h)
that is the discrete equivalent to (6.2).
To observe numerically the behavior of the solution described in this section,
we have integrated (1.2) with α = 0.5,
k(t) = sin2 (t)/tα and f (t) = tα , by the fractional trapezoidal rule. In
this case k(t) < t−α , so according to the analysis reported above yn < C · tαn,
asymptotically, with 0 < C < +∞. This behavior is clear in Fig. 4, where we
have plotted the numerical solution yn divided by tα n.

7. Conclusions
In the analysis of numerical methods for solving weakly singular Volterra
equations, stability is generally studied with respect to the test equation
λ t
y(t) = f (t) + Γ(α) 0
(t − s)α−1 y(s)ds, t ≥ 0, with | arg λ − π| < 1 − 12 α π
(see for example [19]).
In this paper, we propose a different approach which involves the model
problem (1.2), (1.3), representative of a variety of problems that cannot be
treated through the classical test equation. The theoretical results consist
in the boundedness and the convergence of the analytical as well as of the
numerical solution and are illustrated by means of numerical experiments
both on linear and non-linear problems.
Acknowledgements
The research was supported by INdAM-GNCS.

References
[1] Aceto, L., Magherini, C., Novati, P.: Fractional convolution quadrature based
on generalized Adams methods. Calcolo 51(3), 441–463 (2014)
MJOM Stability of Numerical Solutions for Abel–Volterra Integral Page 13 of 14 113

[2] Becker, L.C.: Resolvents and solutions of weakly singular linear Volterra inte-
gral equations. Nonlinear Anal. 74(5), 1892–1912 (2011)
[3] Brunner, H.: Volterra Integral Equations. An Introduction to Theory and Ap-
plications. Cambridge University Press, Cambridge (2017)
[4] Burns, J.A., Cliff, E.M., Herdman, T.L.: A state-space model for an aeroelastic
system. In: Proc. 22nd IEEE Conference on Decision and Control, pp. 1074–
1077 (1983)
[5] Carlone, R., Figari, R., Negulescu, C.: The quantum beating and its numerical
simulation. J. Math. Anal. Appl. 450, 1294–1316 (2017)
[6] Choi, U.Jin, MacCamy, R.C.: Fractional order Volterra equations with appli-
cations to elasticity. J. Math. Anal. Appl. 139(2), 448–464 (1989)
[7] Diethelm, K.: The Analysis of Fractional Differential Equations. Springer,
Berlin (2010)
[8] Diogo, T., Edwards, J.T., Ford, N.J., Thomas, S.M.: Numerical analysis of a
singular integral equation. Appl. Math. Comput. 167, 372–382 (2005)
[9] Fedotov, S., Iomin, A., Ryashko, L.: Non-markovian models for migration-
proliferation dichotomy of cancer cells: anomalous switching and spreading
rate. Phys. Rev. E 84(6), 061131 (2011)
[10] Garrappa, R.: Trapezoidal methods for fractional differential equations: theo-
retical and computational aspects. Math. Comput. Simul. 110, 96–112 (2015)
[11] Garrappa, R., Galeone, L.: On multistep methods for differential equations of
fractional order. Mediterr. J. Math. 3(3–4), 565–580 (2006)
[12] Garrappa, R., Messina, E., Vecchio, A.: Effect of perturbation in the numerical
solution of fractional differential equations. Discrete Contin. Dyn. Syst. Ser. B
(2017). https://doi.org/10.3934/dcdsb.2017188
[13] Garrappa, R., Popolizio, M.: On accurate product integration rules for linear
fractional differential equations. J. Comput. Appl. Math. 235(5), 1085–1097
(2011)
[14] Gorenflo, R., Vessella, S.: Abel Integral Equations. Analysis and Applications,
Volume 1461 of Lecture Notes in Mathematics. Springer, Berlin (1991)
[15] Győri, I., Reynolds, D.W.: On admissibility of the resolvent of discrete Volterra
equations. J. Differ. Equ. Appl. 16(12), 1393–1412 (2010)
[16] Linz, P.: Analytical and Numerical Methods for Volterra Equations. SIAM,
Philadelphia (1985)
[17] Lubich, Ch.: Fractional linear multistep methods for Abel–Volterra integral
equations of the second kind. Math. Comput. 45(172), 463–469 (1985)
[18] Lubich, Ch.: Discretized fractional calculus. SIAM J. Math. Anal. 17(3), 704–
719 (1986)
[19] Lubich, Ch.: A stability analysis of convolution quadratures for Abel–Volterra
integral equations. IMA J. Numer. Anal. 6(1), 87–101 (1986)
[20] Messina, E., Vecchio, A.: Stability and convergence of solutions to Volterra
integral equations on time scales. Discrete Dyn. Nat. Soc. 2015(ID612156), 6
(2015)
[21] Messina, E., Vecchio, A.: Stability and boundedness of numerical approxima-
tions to Volterra integral equations. Appl. Numer. Math. 116, 230–237 (2017)
113 Page 14 of 14 G. Izzo et al. MJOM

[22] Messina, E., Vecchio, A.: A sufficient condition for the stability of direct quad-
rature methods for Volterra integral equations. Numer. Algorithms 74(4), 1223–
1236 (2017)
[23] Podlubny, I.: Fractional Differential Equations. Academic Press, San Diego
(1999)
[24] Surguladze, T.A.: On certain applications of fractional calculus to viscoelastic-
ity. J. Math. Sci. 112(5), 4517–4557 (2002)

G. Izzo and E. Messina

Dipartimento di Matematica e Applicazioni
Università degli Studi di Napoli “Federico II”
Via Cintia
80126 Naples
Italy
e-mail: [email protected]

E. Messina
e-mail: [email protected]

A. Vecchio
C.N.R. National Research Council of Italy, Institute for Computational Application
“Mauro Picone”
Via P. Castellino 111
80131 Naples
Italy
e-mail: [email protected]

Received: July 31, 2017.

Revised: April 17, 2018.
Accepted: April 27, 2018.