Electronic Journal of Differential Equations, Vol. 2025 (2025), No. 05, pp. 1–13.

ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu

DOI: 10.58997/ejde.2025.05

θ-SCHEME FOR SOLVING CAPUTO FRACTIONAL

DIFFERENTIAL EQUATIONS

THAI SON DOAN, PHAN THI HUONG, PETER E. KLOEDEN

Abstract. We formulate a θ-numerical scheme for solving Caputo fractional
differential equations (Caputo FDEs) of order α ∈ (0, 1), with vector fields

satisfying a standard Lipschitz continuity condition in the state variable and
a Hölder continuity condition in the time variable. The convergence rate is es-

tablished and a numerical example is given to illustrate the theoretical results.

The scheme obtained includes the explicit (θ = 0) and implicit (θ = 1) coun-
terparts of Euler-like schemes for Caputo FDEs known in the literature as the

Adams-Bashford and Adams-Moulton schemes, respectively, and essentially

linearly interpolates them.

1. Introduction

Fractional order differential equations are suitable for describing various complex
phenomena, in particular memory and hereditary properties of dynamical processes,
see the monographs [3, 11, 26] and the references therein. An extensive qualitative
theory of such equations is now available (see e.g., [3, 5, 11, 15, 24, 26]. Explicit
solutions are rarely known and then often have complicated expressions, which are
of little practical use. Numerical methods have been developed over many decades
for equations containing derivatives and integrals of non-integer order. More specif-
ically, numerical methods for fractional differential equations have been considered
by many authors [4, 6, 7, 8, 9, 10, 11, 14, 16, 17, 22, 23, 24, 25]. These results
and many more are discussed in the monographs [20, 21, 23]. Of particular interest
here are the fractional Adams-Bashforth/Moulton methods developed by Diethelm
et al [8, 9, 11], which are essentially the Caputo counterparts of the explicit/implicit
Euler schemes.

The Euler scheme is the simplest numerical scheme for a first-order ordinary
differential equation (ODEs) with the vector field f(t, x) and is easily derived due
to the geometric interpretation of the classical derivative. The θ-scheme essentially
linearly interpolates the explicit (θ = 0) and implicit (θ = 1) versions of the Euler
scheme, specifically

x(n)(tn+1) = x(n)(tn) + (1− θ)f(tn, x
(n)(tn))h+ θf(tn, x

(n)(tn+1))h (1.1)

2020 Mathematics Subject Classification. 34A05, 65L99, 65R20.
Key words and phrases. Caputo fractional differential equations; θ-scheme; Euler scheme;

Adams-Bashford scheme; Adams-Moulton scheme.
©2025. This work is licensed under a CC BY 4.0 license.

Submitted August 9, 2024. Published January 9, 2025.

1



2 T. S. DOAN, P. T. HUONG, P. E. KLOEDEN EJDE-2025/05

for a time step h > 0 and parameter 0 ≤ θ ≤ 1. Here, tn = nh. It is implicit except
when θ = 0. This scheme is now finding increasing use, especially when numerical
stability is an issue (see [19]).

The piecewise linear interpolation x(n)(t) of the iterates of the θ-scheme can be
expressed in integral form as

x(n)(t) = x0 + (1− θ)

∫ t

0

f(ζn(s), x
(n)(ζn(s))) ds+ θ

∫ t

0

f(λn(s), x
(n)(λn(s))) ds,

(1.2)

for t ∈ [0, T ], where ζn(s) = kT
n and λn(s) = (k+1)T

n for s ∈
(

kT
n , (k+1)T

n

]
, k =

0, . . . , n− 1.
Now consider a Caputo FDEs of order α ∈ (0, 1) on the interval [0, T ] of the

form

CDα
0+Y (t) = f(t, Y (t)), (1.3)

with the same vector field f : [0, T ] × Rd → Rd. This has the integral equation
representation

Y (t) = y0 +
1

Γ(α)

∫ t

0

f(s, Y (s))

(t− s)1−α
ds, (1.4)

where Γ is the Gamma function, i.e., Γ(α) :=
∫∞
0

sα−1e−s ds.
This and the integral expression (1.2) suggest the following definition of the

piecewise linear interpolation y(n)(t) of the θ-scheme for the Caputo FDEs (1.3):

Y (n)(t) = y0 +
(1− θ)

Γ(α)

∫ t

0

f(ζn(s), Y
(n)(ζn(s)))

(t− s)1−α
ds

+
θ

Γ(α)

∫ t

0

f(λn(s), Y
(n)(λn(s)))

(t− s)1−α
ds

(1.5)

for t ∈ [0, T ]. Indeed, this approach was used by the authors [14] to define the
Euler-Maruyama scheme for Caputo fractional stochastic differential equations.

The θ-scheme (1.5) can be written step by step on each interval
(
kT
n , (k+1)T

n

]
,

k = 0, . . . , n− 1. Write Y
(n)
k = Y (n)(kTn ) with Y

(n)
0 = y0. Then

Y
(n)
k+1

= Y
(n)
0 +

1

Γ(α)

k∑
j=0

[
(1− θ)f(tj , Y

(n)
j ) + θf(tj+1, Y

(n)
j+1)

] ∫ tj+1

tj

1

(tk+1 − s)1−α
ds

= Y
(n)
0 +

1

αΓ(α)

k∑
j=0

[
(1− θ)f(tj , Y

(n)
j ) + θf(tj+1, Y

(n)
j+1)

]
×
(
(tk+1 − tj)

α − (tk+1 − tj+1)
α
)

= Y
(n)
0 +

1

Γ(α+ 1)

k∑
j=0

[
(1− θ)f(tj , Y

(n)
j ) + θf(tj+1, Y

(n)
j+1)

]
×
(
(k + 1− j)α − (k − j)α

)(T
n

)α
.



EJDE-2025/05 θ-SCHEME 3

This gives the recursive formula

Y
(n)
k+1 = Y

(n)
0 +

k∑
j=0

γk−j

[
(1− θ)f(tj , Y

(n)
j ) + θf(tj+1, Y

(n)
j+1)

] (T
n

)α
(1.6)

for k = 0, . . . , n− 1, here

γj :=
1

Γ(α+ 1)
((j + 1)α − jα) .

Note that γj decreases from γ0 = 1
Γ(α+1) as j increases, i.e., the weighting of the

contribution of f(tj , Y
(n)
j ) and f(tj+1, Y

(n)
j+1) for a given j decreases as n increases.

This scheme (1.6) is known in the literature as the Adams-Bashford scheme
for θ = 0 and Adams-Moulton scheme for θ = 1, see [9, 10, 11]. These are the
counterparts of the explicit and implicit Euler schemes for Caputo FDEs. They
were also given without derivation or error bound in [16, p. 19].

When α = 1, in which case the weights γj ≡ 1, the expression (1.6) reduces to
θ-scheme (1.1) for ordinary differential equations (ODEs). For simplicity this will
be shown for the case θ = 0.

Y
(n)
k+1 = Y

(n)
0 +

k∑
j=0

f(tj , Y
(n)
j )

T

n
, k = 0, . . . , n− 1,

which gives

Y
(n)
1 = Y

(n)
0 + f(t0, Y

(n)
0 )

T

n
,

Y
(n)
2 = Y

(n)
0 + f(t0, Y

(n)
0 )

T

n
+ f(t1, Y

(n)
1 )

T

n
= Y

(n)
1 + f(t1, Y

(n)
1 )

T

n
,

and so on, culminating in

Y
(n)
k+1 = Y

(n)
k + f(tk, Y

(n)
k )

T

n
.

The case θ = 1 and the general case are similar.
This article is organized as follows: In Section 2, we give a setting of the problem

and state the main result of the paper (Theorem 2.1) on θ-scheme. Section 3 is
devoted to proving the main result. Finally, in Section 4, a simple example is
studied numerically.

2. Main result

Fix T > 0 arbitrarily and consider the Caputo FDEs (1.3) of order α ∈ (0, 1) on
the interval [0, T ] where f : [0, T ]× Rd → Rd satisfies the following conditions

(H1) Lipschitz continuity with respect to the second variable: There exists L1 >
0 such that for all y, y1 ∈ Rd, t ∈ [0, T ],

∥f(t, y)− f(t, y1)∥ ≤ L1∥y − y1∥.

(H2) Hölder continuity with respect to the first variable: There exists L2 > 0
and β ∈ [0, 1] such that for all y ∈ Rd, t, t1 ∈ [0, T ]

∥f(t, y)− f(t1, y)∥ ≤ L2|t− t1|β .



4 T. S. DOAN, P. T. HUONG, P. E. KLOEDEN EJDE-2025/05

As a consequence, the function f also satisfies a linear growth bound

∥f(t, y)∥ ≤ M + L1∥y∥, y ∈ Rd, t ∈ [0, T ], (2.1)

with M := maxt∈[0,T ] ∥f(t, 0)∥.
Then, (1.3) with the initial value Y (0) = y0 ∈ Rd has a unique solution on

[0, T ] denoted by Y (t), which satisfies the integral representation (1.4) (see e.g., [2,
Theorem 2] and [27, Theorem 6.4]). Notice that the kernel (t − s)α−1 of (1.4) is
singular at the point s = t, but integrable.

The first part of the following result shows that the θ-scheme (1.5) is uniquely
solvable when the step size is small enough. The second part establishes a global
discretization error between the numerical solution Y (n)(t) and the exact solution
Y (t) of the Caputo FDEs (1.3). The next theorem presents the solvability and
convergence of the θ-scheme (1.5) for Caputo FDEs.

Theorem 2.1. Suppose that Assumptions (H1) and (H2) hold and let N be the
smallest natural number satisfying that

N ≥ T α

√
2 θ L1

Γ(α+ 1)
, (2.2)

where L1 is the constant in Assumption (H1). Then, the following statements hold:

(i) For any n ≥ N , equation (1.5) has a unique solution Y (n)(tk+1) with
Y (n)(0) := y0 for each k = 0, . . . , n− 1.

(ii) There exists a constant C depending only on θ, y0, T, L1, L2,M, α, β such
that

sup
0≤t≤T

∥Y (n)(t)− Y (t)∥ ≤ C

nκ
for all n ≥ N, (2.3)

where κ := min
{
α, β

}
.

Remark 2.2. Similar results have been given by Diethelm [11, Appendix C] for the
Adams-Bashford scheme (θ = 0) and the Adams-Moulton scheme (θ = 1). When
α = 1, i.e., (1.3) becomes an ODEs, the convergence rate of the scheme in Theorem
2.1 coincides with the well-known convergence rate of the explicit and implicit Euler
schemes for ODEs (with θ = 0 and θ = 1, respectively) (see e.g., [1]).

3. Proof of the main result

First we prove solvability of the θ-scheme (1.5), which is implicit for 0 < θ ≤ 1.

Proof of Theorem 2.1 (i). Choose and fix an arbitrary natural number n such that

n ≥ T α

√
2θL1

Γ(α+ 1)
, (3.1)

where L1 is the constant given in Assumption (H1). The Y (n)(t0), . . . , Y
(n)(tn−1)

are obtained inductively by solving the equations

Y (n)(tk+1) = y0 + (1− θ)

k∑
j=0

γk−jf(tj , Y
(n)(tj))

(T
n

)α
+ θ

k∑
j=0

γk−jf(tj+1, Y
(n)(tj+1))

(T
n

)α
.

(3.2)



EJDE-2025/05 θ-SCHEME 5

Obviously, Y (n)(t0) = y0 , so assume that the Y (n)(t0), . . . , Y
(n)(tj) have been

determined for some j ∈ {0, . . . , n−2}. It needs to be shown that there is a unique
solution Y (n)(tk+1) of (3.2). For this purpose, define the map φ : Rd → Rd by

φ(x) = x− θγ0f(tk+1, x)
(T
n

)α − y0 − (1− θ)

k∑
j=0

γk−jf(tj , Y
(n)(tj))

(T
n

)α
− θ

k−1∑
j=0

γk−jf(tj+1, Y
(n)(tj+1))

(T
n

)α
.

For arbitrary x, y ∈ Rd, we have

∥φ(x)− φ(y)∥ = θγ0
(T
n

)α∥f(tk+1, x)− f(tk+1, y)∥

≤ θγ0
(T
n

)α
L1∥x− y∥,

where we used Assumption (H1) to obtain the last inequality. From (3.1), we derive
that

θγ0
(T
n

)α
L1 =

θ

Γ(α+ 1)

(T
n

)α
L1 < 1.

Therefore, the map φ is contractive. Using the Banach fixed point theorem, there
exists a unique fixed point of this map in Rd. □

In the rest of this section, we will prove the Theorem 2.1(ii). Our first task
is to deduce from (1.5) an upper bound for sup0≤t≤T ∥Y (n)(t)∥. Since Y (n)(t)
is constructed implicitly, the version of Gronwall’s inequality for Caputo FDEs
stated in [11, Lemma 6.19] (also in Dixon [12]) is not applicable. To overcome this
difficulty, we use the weight function Eα (µ tα), the increasing monotonicity of the
Mittag-Leffler function Eα (·) on [0,∞), and some skillful transformations. Here,
the Mittag-Leffler function Eα : R → R is defined by

Eα(t) :=
∞∑
k=0

tk

Γ(αk + 1)
, t ∈ R, α ∈ (0, 1).

An upper bound for the approximate solution is presented in the following lemma.

Lemma 3.1. Let n ∈ N satisfy n ≥ T α

√
2θ L1

Γ(α+1) and define

C1 :=
(
2∥y0∥+

2M(1 + θ)Tα

Γ(α+ 1)

)µEα (µTα)

µ− 2L1
. (3.3)

Then

sup
0≤t≤T

∥Y (n)(t)∥ ≤ C1. (3.4)



6 T. S. DOAN, P. T. HUONG, P. E. KLOEDEN EJDE-2025/05

Proof. Using (1.5), and the inequality ∥y1+y2+y3+y4∥ ≤ ∥y1∥+∥y2∥+∥y3∥+∥y4∥
for all y1, y2, y3, y4 ∈ Rd, we obtain∥∥∥Y (n)(tk+1)−

θ

Γ(α)

∫ tk+1

tk

f
(
tk+1, Y

(n)(tk+1)
)
− f(tk+1, 0)

(tk+1 − s)1−α
ds
∥∥∥

≤ ∥y0∥+
θ

Γ(α)

∫ tk+1

tk

∥f(tk+1, 0)∥
(tk+1 − s)1−α

ds

+
θ

Γ(α)

∫ tk

0

∥f(λn(s), Y
(n)(λn(s)))∥

(tk+1 − s)1−α
ds

+
(1− θ)

Γ(α)

∫ tk+1

0

∥f(ζn(s), Y (n)(ζn(s)))∥
(tk+1 − s)1−α

ds.

(3.5)

By using Assumption (H1) and condition (3.1), we arrive at∥∥∥Y (n)(tk+1)−
θ

Γ(α)

∫ tk+1

tk

f
(
tk+1, Y

(n)(tk+1)
)
− f(tk+1, 0)

(tk+1 − s)1−α
ds
∥∥∥

≥ ∥Y (n)(tk+1)∥ −
θ

Γ(α)

∫ tk+1

tk

L1∥Y (n)(tk+1)∥
(tk+1 − s)1−α

ds

≥ ∥Y (n)(tk+1)∥
(
1− θ L1

Γ(α+ 1)

(T
n

)α) ≥ 1

2
∥Y (n)(tk+1)∥.

Substituting this into (3.5), then apply estimate (2.1) leads to

∥Y (n)(tk+1)∥ ≤ 2∥y0∥+
2 θ M

Γ(α+ 1)

(T
n

)α
+

2θMTα

Γ(α+ 1)
+

2 (1− θ) M Tα

Γ(α+ 1)

+
2θL1

Γ(α)

∫ tk

0

∥Y (n)(λn(s))∥
(tk+1 − s)1−α

ds

+
2(1− θ)L1

Γ(α)

∫ tk+1

0

∥Y (n)(ζn(s))∥
(tk+1 − s)1−α

ds.

(3.6)

Next, to give an upper bound for sup0≤t≤T ∥Y (n)(t)∥, we need weight Eα(µ tαk+1),
where µ is a positive constant such that

µ > 2L1. (3.7)

From the increasing monotonicity of the functions Eα(·) on [0,∞) and f(s) =
(tk+1 − s)α−1, we have∫ tk

0

∥Y (n)(λn(s))∥
(tk+1 − s)1−αEα(µtαk+1)

ds

≤
k−1∑
i=0

∫ ti+1

ti

Eα(µt
α
i+1)

(tk+1 − s)1−αEα(µ tαk+1)

∥Y (n)(λn(s))∥
Eα(µ(λn(s))α)

ds

≤
k−1∑
i=0

∫ ti+2

ti+1

Eα(µs
α)

(tk+1 − (s− T
n ))

1−αEα(µtαk+1)
ds sup

0≤t≤T

∥Y (n)(t)∥
Eα(µtα)

≤
∫ tk+1

0

Eα(µs
α)

(tk+1 − s)1−αEα(µtαk+1)
ds sup

0≤t≤T

∥Y (n)(t)∥
Eα(µtα)

≤ Γ(α)

µ
sup

0≤t≤T

∥Y (n)(t)∥
Eα(µtα)

,

(3.8)



EJDE-2025/05 θ-SCHEME 7

here in the final step, we used the inequality

µ

Γ(α)

∫ t

0

(t− s)
α−1

Eα (µsα) ds ≤ Eα (µtα),

see [13, Lemma 5].
Observing that ζn(s) ≤ s, hence in a similar manner as above, we obtain∫ tk+1

0

∥Y (n)(ζn(s))∥
(tk+1 − s)1−αEα(µtαk+1)

ds ≤ Γ(α)

µ
sup

0≤t≤T

∥Y (n)(t)∥
Eα(µtα)

.

Inserting this and (3.8) into (3.6) leads to

sup
0≤t≤T

∥Y (n)(t)∥
Eα(µtα)

≤ 2θL1

µ
sup

0≤t≤T

∥Y (n)(t)∥
Eα(µtα)

+
2(1− θ)L1

µ
sup

0≤t≤T

∥Y (n)(t)∥
Eα(µtα)

+ 2∥y0∥+
2M(1 + θ) Tα

Γ(α+ 1)
.

Consequently,(
1− 2L1

µ

)
sup

0≤t≤T

∥Y (n)(t)∥
Eα(µ tα)

≤ 2∥y0∥+
2M(1 + θ) Tα

Γ(α+ 1)
.

The condition (3.7) gives

sup
0≤t≤T

∥Y (n)(t)∥ ≤
(
2∥y0∥+

2M(1 + θ) Tα

Γ(α+ 1)

)µEα (µTα)

µ− 2L1
,

which completes the proof. □

Secondly, we determine an upper bound on ∥Y (n)(t) − Y (n)(t1)∥ in terms of
|t− t1|α.

Lemma 3.2. Let

C2 :=
2(M + L1C1)

Γ(α+ 1)
, (3.9)

where C1 is given as in (3.3). Then, for all n ≥ N and t, t1 ∈ [0, T ],

∥Y (n)(t)− Y (n)(t1)∥ ≤ C2|t− t1|α.

Proof. Choose a fix t, t1 ∈ [0, T ] with t > t1. From (1.5), we obtain

Y (n)(t)− Y (n)(t1)

=
θ

Γ(α)

∫ t

t1

f(λn(s), Y
(n)(λn(s)))

(t− s)1−α
ds+

(1− θ)

Γ(α)

∫ t

t1

f(ζn(s), Y
(n)(ζn(s)))

(t− s)1−α
ds

+
θ

Γ(α)

∫ t1

0

( 1

(t− s)1−α
− 1

(t1 − s)1−α

)
f(λn(s), Y

(n)(λn(s))) ds

+
(1− θ)

Γ(α)

∫ t1

0

( 1

(t− s)1−α
− 1

(t1 − s)1−α

)
f(ζn(s), Y

(n)(ζn(s))) ds.

Using the inequality ∥y1 + y2 + y3 + y4∥ ≤ ∥y1∥ + ∥y2∥ + ∥y3∥ + ∥y4∥ for all
y1, y2, y3, y4 ∈ Rd, and by (2.1), we obtain∥∥∥Y (n)(t)− Y (n)(t1)

∥∥∥
≤ θ

Γ(α)

∫ t

t1

M + L1∥Y (n)(λn(s))∥
(t− s)1−α

ds+
(1− θ)

Γ(α)

∫ t

t1

M + L1∥Y (n)(ζn(s))∥
(t− s)1−α

ds



8 T. S. DOAN, P. T. HUONG, P. E. KLOEDEN EJDE-2025/05

+
θ

Γ(α)

∫ t1

0

( 1

(t1 − s)1−α
− 1

(t− s)1−α

)
(M + L1∥Y (n)(λn(s))∥) ds

+
(1− θ)

Γ(α)

∫ t1

0

( 1

(t1 − s)1−α
− 1

(t− s)1−α

)
(M + L1∥Y (n)(ζn(s))∥) ds.

Using Lemma 3.1, we obtain

∥Y (n)(t)− Y (n)(t1)∥

≤
∫ t

t1

M + L1C1

Γ(α) (t− s)1−α
ds+

1

Γ(α)

∫ t1

0

( 1

(t1 − s)1−α
− (t− s)1−α

)
(M + L1C1)ds.

A direct computation gives

∥Y (n)(t)− Y (n)(t1)∥ ≤ 2(M + L1C1)

Γ(α+ 1)
(t− t1)

α,

which completes proof. □

Proof of Theorem 2.1 (ii). Choose and fix y0 ∈ Rd. From (1.4) and (1.5) we have

Y (n)(t)− Y (t) =
θ

Γ(α)

∫ t

0

f(λn(s), Y
(n)(λn(s)))− f(s, Y (s))

(t− s)1−α
ds.

+
1− θ

Γ(α)

∫ t

0

f(ζn(s), Y
(n)(ζn(s)))− f(s, Y (s))

(t− s)1−α
ds.

Therefore,

∥Y (n)(t)− Y (t)∥ =
θ

Γ(α)

∫ t

0

∥f(λn(s), Y
(n)(λn(s)))− f(s, Y (s))∥
(t− s)1−α

ds.

+
1− θ

Γ(α)

∫ t

0

∥f(ζn(s), Y (n)(ζn(s)))− f(s, Y (s))∥
(t− s)1−α

ds.

(3.10)

In view of Assumptions (H1) and (H2) it is easily seen that

∥f(λn(s), Y
(n)(λn(s)))− f(s, Y (s))∥ ≤ L1∥Y (n)(λn(s))− Y (s)∥+ L2 |λn(s)− s|β .

(3.11)
By the definition of ρn, we have |λn(s)−s| ≤ T

n for s ∈ [0, T ]. This together Lemma
3.2 implies that

∥f(λn(s), Y
(n)(λn(s)))− f(s, Y (s))∥

≤ L1∥Y (n)(λn(s))− Y (n)(s)∥+ L1∥Y (n)(s)− Y (s)∥+ L2

(T
n

)β
≤ L1C2|λn(s)− s|α + L1∥Y (n)(s)− Y (s)∥+ L2

(T
n

)β
≤ L1C2T

α 1

nα
+ L2T

β 1

nβ
+ L1∥Y (n)(s)− Y (s)∥,

(3.12)

where C2 is given as in (3.9). Similarly, we have

∥f(ζn(s), Y (n)(ζn(s)))−f(s, Y (s))∥ ≤ L1C2T
α 1

nα
+L2T

β 1

nβ
+L1∥Y (n)(s)−Y (s)∥.

Inserting this, (3.11), and (3.12) into (3.10) results in

∥Y (n)(t)− Y (t)∥ ≤ L1

Γ(α)

∫ t

0

∥Y (n)(s)− Y (s)∥
(t− s)1−α

ds+
L1C2T

2α

Γ(α+ 1)

1

nα
+

L2T
α+β

Γ(α+ 1)

1

nβ
.



EJDE-2025/05 θ-SCHEME 9

Applying the Gronwall inequality for FDEs (see e.g., [11, Lemma 6.19]), we obtain

∥Y (n)(t)− Y (t)∥ ≤ (L1C2T
2α + L2T

α+β)Eα (L1T
α)

Γ(α+ 1)

1

nκ
.

With

C :=
(L1 C2T

2α + L2T
α+β)Eα (L1T

α)

Γ(α+ 1)

this completes the proof. □

4. A numerical example

We investigate numerically a simple Caputo FDEs for which we have an explicit
expression for the solution. In particular, we consider the simple scalar linear
Caputo FDEs

CDα
0+Y (t) = Y (t) (4.1)

on the interval [0, 1]. It has the exact solution for Y (0) = 1 is given by Y (1) = Eα(1)
(see [11, p. 135]). Moreover, Y (1) has the integral representation

Y (1) = 1 +
1

Γ(α)

∫ 1

0

Y (s)

(1− s)1−α
ds.

The convergence rates obtained are consistent with that stated in Theorem 2.1,
which are for the whole class of vector fields and can be better in individual cases.

4.1. Explicit case: θ = 0. The recursion expression (1.6) for the numerical solu-
tion simplifies to

Y
(n)
k+1 = 1 +

k∑
j=0

γk−j Y
(n)
j

1

nα
, k = 0, . . . , n− 1.

We denote the global discretization error between the numerical solution Y
(n)
n and

the exact solution Y (1) by

en = |Y (1)− Y (n)
n | = |Eα(1)− Y (n)

n |

and the experimental order of convergence (EOC) by

EOC := log2
en
e2n

= log2
|Eα(1)−Y

(n)
n |

|Eα(1)−Y
(2n)
2n |

,

see e.g. [17, p. 8].
The results for the global discretization error, EOC and computational cost of

the θ-scheme with θ = 0 are given in Table 1. All time values in the table are given
in seconds.

The slopes of the lines in the coordinate system with the horizontal axis being
log2 n and the vertical axis being log2(en) in Figure 1 give the convergence rates of
the θ-scheme with θ = 0: for α = 0.5 the slope is −0.9464465128045231, and for
α = 0.9 the slope is −0.9587708473867738.



10 T. S. DOAN, P. T. HUONG, P. E. KLOEDEN EJDE-2025/05

Table 1. Global error, EOC and computational cost for the θ-
scheme with θ = 0

α = 0.5 α = 0.9
en EOC Time en EOC Time

n = 8 1.4158 0.8262 0.0010 0.5202 0.8880 0.0003
n = 16 0.7985 0.9186 0.0010 0.2811 0.9419 0.0000
n = 32 0.4224 0.9680 0.0040 0.1463 0.9706 0.0016
n = 64 0.2160 0.9912 0.0050 0.0747 0.9853 0.0037
n = 128 0.1086 1.0008 0.0160 0.0377 0.9928 0.0156
n = 256 0.0543 1.0041 0.0588 0.0190 0.9964 0.0568
n = 512 0.0271 1.0047 0.2218 0.0095 0.9982 0.2702

α = 0.5, slope = −0.9464465128045231 α = 0.9, slope = −0.9587708473867738

Figure 1. Convergence rate for the θ-scheme with θ = 0.

4.2. Implicit case θ = 0.5. The recursion expression (1.6) simplifies to

Y
(n)
k+1 = Y

(n)
0 +

k∑
j=0

γk−j

[
(1− θ) Y

(n)
j + θ Y

(n)
j+1

] 1

nα
,

This is uniquely solvable for n ≥ α

√
2θ

Γ(α+1) and is linear, so can be solved alge-

braically to give

Y
(n)
k+1 = ΓnY

(n)
0 + γ0Γn(1− θ)Y

(n)
k

1

nα
+ Γn

k−1∑
j=0

γk−j

[
(1− θ) Y

(n)
j + θY

(n)
j+1

] 1

nα
,

for k = 0, . . . , n− 1, where Γn := 1
1−γ0

1
nα

.

Results on the global discretization error, EOC and computational costs of the
θ scheme with θ = 0.5 are given in Table 2

The slopes of the lines in Figure 2 give the convergence rates of the θ-scheme
with θ = 0.5: for α = 0.5 the slope is −1.0207618200185753, for α = 0.9 the slope
is −0.9937723963388096.



EJDE-2025/05 θ-SCHEME 11

Table 2. Global error, EOC and computational cost for θ-scheme
with θ = 0.5

α = 0.5 α = 0.9
en EOC Time en EOC Time

n = 32 0.1917 1.0234 0.0010 0.0928 0.9825 0.0010
n = 64 0.0943 1.0251 0.0066 0.0469 0.9913 0.0060
n = 128 0.0463 1.0221 0.0233 0.0236 0.9957 0.0215
n = 256 0.0228 1.0177 0.0884 0.0118 0.9979 0.0878
n = 512 0.0113 1.0136 0.3773 0.0059 0.9989 0.3826
n = 1024 0.0056 1.0102 1.4344 0.0030 0.9995 1.5066
n = 2048 0.0028 1.0075 5.9136 0.0015 0.9997 5.5054

α = 0.5, slope= −1.0207618200185753 α = 0.9, slope= −0.9937723963388096.

Figure 2. Convergence rate for θ-scheme with θ = 0.5.

4.3. Fully implicit case θ = 1. For n ≥ α

√
2 θ

Γ(α+1) the θ-scheme (1.6) with θ = 1

has a unique numerical solution Y (n) of

Y
(n)
k+1 = Y

(n)
0 +

k∑
j=0

γk−jY
(n)
j+1

1

nα
, k = 0, . . . , n− 1,

which is linear and thus it can be solved algebraically to give

Y
(n)
k+1 = ΓnY

(n)
0 + Γn

k−1∑
j=0

γk−jY
(n)
j+1

1

nα
, k = 0, . . . , n− 1,

where Γn := 1
1−γ0

1
nα

.

Results on the global discretization error, EOC and computational costs of the
θ scheme with θ = 1 are given in the Table 3.

The slopes of the lines in Figure 3 give the convergence rates of the θ-scheme
with θ = 1: for α = 0.5 the slope is −0.7737488258780804, and for α = 0.9 the
slope is −1.0136941007464315.

Acknowledgements. This research is funded by Vietnam National Foundation for
Science and Technology Development (NAFOSTED) under grant number 101.02-
2023.32. P. E. KLoeden was funded by NSFC Grant 12371242 the Ministerio de



12 T. S. DOAN, P. T. HUONG, P. E. KLOEDEN EJDE-2025/05

Table 3. Global error, EOC and computational cost for the θ-
scheme with θ = 1.

α = 0.5 α = 0.9
en EOC Time en EOC Time

n = 8 0.1182 0.53023 0.0010 0.1589 1.0327 0.0010
n = 16 0.0819 0.7122 0.0010 0.0777 1.0188 0.0000
n = 32 0.0500 0.8076 0.0030 0.0383 1.0109 0.0010
n = 64 0.0285 0.8666 0.0120 0.0190 1.0063 0.0070
n = 128 0.0157 0.9062 0.0389 0.0095 1.0036 0.0259
n = 256 0.0084 0.9336 0.1606 0.0047 1.0021 0.0898
n = 512 0.0044 0.9529 0.6361 0.0024 1.0012 0.3405

α = 0.5, slope= −0.7737488258780804 α = 0.9, slope= −1.0136941007464315.

Figure 3. Convergence rate for the θ-scheme with θ = 1.

Ciencia e Innovación (Spain) and FEDER (European Community) under grant
PID2021-122991NB-C21.

References

[1] K. E. Atkinson, W. Han, D. Stewart; Numerical Solution of Ordinary Differential Equations,
John Wiley & Sons, Inc., 2009.

[2] D. Bǎleanu, O. G. Mustafa; On the global existence of solutions to a class of fractional
differential equations, Computer and Mathematics with Applications, 17 (2010) (59), 1583-

1841.
[3] B. Bandyopadhyay, S. Kamal; Stabilization and Control of Fractional Order Systems: A

Sliding Mode Approach, Lecture Notes in Electrical Engineering, 317, Springer International
Publishing, Switzerland, 2015.

[4] L. Blank; Numerical treatment of differential equations of fractional order, Nonlinear World,
4 (1997), 473–490.

[5] N. D. Cong, H. T. Tuan, H. Trinh; On asymptotic properties of solutions to fractional differ-
ential equations, Journal of Mathematical Analysis and Applications, 484(2) (2020), 123759.

[6] K. Diethelm; An algorithm for the numerical solution of differential equations of fractional
order, Electronic Transactions on Numerical Analysis, 5 (1997), 1–6.

[7] K. Diethelm, Y. Luchko; Numerical solutions of linear multi-term initial value problems of
fractional order, J. Comput. Anal. Appl., 6 (2004), no. 3, 243-263.

[8] K. Diethelm, N. J. Ford, A. D. Freed; A predictor-corrector approach for the numerical

solution of fractional differential equations, Nonlinear Dynamics 29 (2002), 3–22.



EJDE-2025/05 θ-SCHEME 13

[9] K. Diethelm, N. J. Ford, A. D. Freed; Detailed error analysis for a fractional Adams method,

Numerical Algorithms, 36 (2004), 31–52.

[10] K. Diethelm, N. J. Ford, A. D. Freed, Y. Luchko; Algorithms for the fractional calculus: A
selection of numerical methods, Computer Methods in Applied Mechanics and Engineering,

194 (2005), 743–773.

[11] K. Diethelm; The Analysis of Fractional Differential Equations. An Application-oriented
Exposition Using Differential Operators of Caputo Type, Lecture Notes in Mathematics,

2004, Springer-Verlag, Berlin, 2010.

[12] J. Dixon; On the order of the error in discretization methods for weakly singular second kind
Volterra integral equations with non-smooth solutions, BIT, 25 (1986), 624–634.

[13] T. S. Doan, P. T. Huong, P. E. Kloeden, H. T. Tuan; Asymptotic separation between solutions

of Caputo fractional stochastic differential equations. Stoch. Anal. Appl., 36 (2018), no. 4,
654–664.

[14] T. S. Doan, P. T. Huong, P. E. Kloeden, V. A. My; Euler-Maruyama scheme for Caputo frac-
tional stochastic differential equations, Journal of Computational and Applied Mathematics,

380 (2020), 112989.

[15] T. S. Doan, P. E. Kloeden; Semi-dynamical systems generated by autonomous Caputo frac-
tional differential equations, Vietnam Journal of Mathematics, 49 (2021), 1305-1315.

[16] Y. Feng, L. Li, J. G. Liu, X. Xu; Continuous and discrete one dimensional autonomous

fractional ODEs, Discrete and Continuous Dynamical Systems Series B, 23(8) (2018), 3109-
3135.

[17] N. J. Ford, J. A. Connolly; Comparison of numerical methods for fractional differential equa-

tions, Communications on Pure and Applied Analysis, 5(2) (2006), 1-19.
[18] R. Gorenflo; Fractional calculus: Some numerical methods, CISM Courses and Lectures, no.

378 (1997), Springer, Wien, 277-290.

[19] D.J. Higham, P. E. Kloeden; An Introduction to Numerical Simulation of Stochastic Differ-
ential Equations, SIAM, Philadelphia, 2021.

[20] B. Jin, Z. Zhou; Numerical Treatment and Analysis of Time-Fractional Evolution Equations,
Springer-Verlag, Berlin, 2023.

[21] G. E. Karniadakis (Editor); Handbook of Fractional Calculus with Applications, Vol. 3, Nu-

merical Methods., De Gruyter, Berlin, 2019.
[22] C. Li, A. Chen, J. Ye; Numerical approaches to fractional calculus and fractional ordinary

differential equation, Journal of Computational Physics, 230 (2011), 3352-3368.

[23] Changpin L, Fanhai Zeng; Numerical Methods for Fractional Calculus Chapman and
Hall/CRC, New York, 2015.

[24] J. A. T. Machado; Handbook of Fractional Calculus with Applications, De Gruyter, Berlin,

2019.
[25] K. B. Oldham, J. Spanier; The Fractional Calculus, Academic Press, New York, 1974.

[26] S. G. Samko, A. A. Kilbas, O. I. Marichev; Fractional Integrals and Derivatives: Theory and

Applications, Gordon and Breach Science Publishers, 1993.
[27] C. C. Tisdell; On the application of sequential and fixed-point methods to fractional differen-

tial equations of arbitrary order, Journal of Integral Equations and Applications, 24 (2012),
no. 2, 283–319.

Thai Son Doan

Institute of Mathematics, Vietnam Academy of Science and Technology, 18 Hoang Quoc
Viet, Ha Noi, Vietnam

Email address: dtson@math.ac.vn

Phan Thi Huong
Department of Mathematics, Le Quy Don Technical University, 236 Hoang Quoc Viet,

Ha Noi, Vietnam

Email address: pthuong175@gmail.com

Peter E. Kloeden

Mathematiches Insitutes, Universität Tübingen, D-72076 Tübingen, Germany
Email address: kloeden@math.uni-frankfurt.de