Electronic Journal of Differential Equations, Vol. 2023 (2023), No. 58, pp. 1–17.

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

DOI: 10.58997/ejde.2023.58

NON-LOCAL FRACTIONAL BOUNDARY VALUE PROBLEMS

WITH APPLICATIONS TO PREDATOR-PREY MODELS

MICHAL FEČKAN, KATERYNA MARYNETS

Abstract. We study a nonlinear fractional boundary value problem (BVP)

subject to non-local multipoint boundary conditions. By introducing an appro-
priate parametrization technique we reduce the original problem to an equiv-

alent one with already two-point restrictions. Using a notion of Chebyshev

nodes and Lagrange polynomials we construct a successive iteration scheme,
that converges to the exact solution of the non-local problem for particular

values of the unknown parameters, which are calculated numerically.

1. Introduction

Fractional calculus and fractional differential equations denote a separate and
dynamically developing branch in the classical theory of differential equations and
dynamical systems. The first mention of fractional calculus is dated by late 1695 and
belongs to Leibniz. However a real breakthrough in this field belongs to Riemann
and Liouville, who introduced a notion of the Riemann-Liouville fractional integral
and differential operators [15, 24]. Later followed also the Caputo type, Hilfer and
Hilfer-Prabhakar derivatives, all of which belong to the group of fractional operators
with a singular kernel [1, 15, 20]. Due to this singularity differential equations are
characterized by a non-local behavior, what makes analysis of the corresponding
dynamical systems even more challenging.

In this article we study a differential system with Caputo fractional derivative
of order p ∈ (0, 1). This differential operator is often used in modeling of the real-
world problems, because of its property to preserve values of the unknown function
and its derivatives, which in turn coincide with the integer order derivatives (for
more information on the applications we refer to [6, 15, 21, 22, 23]). This gives a
particular advantage in analysis of the corresponding initial value or boundary value
problems (IVPs or BVPs). To answer questions about existence and uniqueness of
solutions, their boundedness or stability one can apply tools from the operator and
fixed point theory, stability and functional analysis. However, less is possible to
find an explicit form of the exact solution to the studied problem. The well-known
techniques, such as the Laplace, Mellin or Fourier transforms, power series method
etc., have already their limitations when studying the linear systems. Since most

2020 Mathematics Subject Classification. 34A08, 34K07, 34K28.
Key words and phrases. Caputo derivative; non-local boundary conditions; Chebyshev nodes;

approximation of solutions; Lagrange polynomial interpolation; predator-prey model.
©2023. This work is licensed under a CC BY 4.0 license.

Submitted February 23, 2023. Published September 11, 2023.

1



2 M. FEČKAN, K. MARYNETS EJDE-2023/58

dynamical models, describing the real-world processes, are of the nonlinear nature,
development of efficient approximation techniques is one of the promising directions
in this field.

Recently a series of new results in approximation of solutions to the nonlinear
fractional differential systems (FDS) of the Caputo type was published. Authors
suggested a so-called “successive approximations technique” [19] to construct ap-
proximate solutions of the studied systems by incorporating also different types
of linear [2, 8, 9, 11] and nonlinear boundary conditions [10, 12]. This approach
demonstrates its effectiveness also when dealing with highly nonlinear systems, sys-
tems of a general order p ∈ (m,m+ 1) or mixed order equations (see discussions in
[3, 4, 8, 9]).

Motivated by these results, in this paper we analyze a fractional BVP (FBVP)
with a non-local behavior not only in the differential system, but also in the bound-
ary conditions. Moreover, we modify the aforementioned successive approximations
technique by coupling it with the polynomial interpolation method [5, 13, 14, 18].
The outline of the paper is the following. In Section 2 we introduce main notations,
used throughout the paper, and formulate necessary definitions and statements,
needed to prove our results. Section 3 contains the problem formulation and justifi-
cation of the parametrization techniques for elimination of the non-local boundary
conditions. In Section 4 we construct a “classical” and “improved” (polynomial)
successive approximations schemes and prove the main results of this paper. And
finally, in Section 5 we demonstrate efficiency of the developed method on a par-
ticular example of a non-local FBVP describing a predator-prey model with prey
refuge.

2. Notation, auxiliary statements and definitions

For a fixed n ∈ N and a bounded set D ⊂ Rn, the following notation apply:

• For any vector-column x = col(x1, x2, . . . , xn) ∈ Rn and n×n real matrix A
operations | · |, =, ≤, ≥, max, inf, and sup are understood component-wise;
• In is a unit n-dimensional matrix;
• On is a zero n-dimensional matrix;
• N0 := N ∪ {0};
• r(A) is the maximal (in modulus) eigenvalue of matrix A.

2.1. Auxiliary lemmas.

Definition 2.1 ([12]). For a set D ⊂ Rn, closed interval [a, b] ⊂ R, Caratheodory
function f : [a, b]×D → Rn and the n-dimensional square matrix K with non-
negative entries, we write f ∈ Lip(K,D), if

|f(t, u)− f(t, v)| ≤ K|u− v| (2.1)

holds, for all {u, v} ⊂ D and a.e. t ∈ [a, b].

Lemma 2.2 ([2]). Let f(t) be a continuous function for t ∈ [0, T ]. Then for all
t ∈ [0, T ] it holds∣∣∣ 1

Γ(p)

[ ∫ t

0

(t− s)p−1f(s)ds−
( t
T

)p ∫ T

0

(T − s)p−1f(s)ds
]∣∣∣

≤ α1(t) max
t∈[0,T ]

|f(t)|,
(2.2)



EJDE-2023/58 NON-LOCAL FRACTIONAL BVPS WITH APPLICATIONS 3

where

α1(t) =
2tp

Γ(p+ 1)

(
1− t

T

)p
. (2.3)

Lemma 2.3 ([2]). Let {αm(t)}m∈Z0
be a sequence of continuous functions at the

interval [0, T ] given by

αm(t) := (Iαm−1)(t), m ∈ N, (2.4)

with

(Iy)(t) :=
1

Γ(p)

[ ∫ t

0

(
(t− s)p−1 − t

T
(T − s)p−1

)
y(s)ds

+
t

T

∫ T

t

(T − s)p−1y(s)ds
]
,

(2.5)

where α0(t) = 1 and α1(t) is defined by formula (2.3). Then

αm+1(t) ≤ Tmp

2m(2p−1)Γm(p+ 1)
α1(t) ≤ T (m+1)p

2(m+1)(2p−1)Γm+1(p+ 1)
, (2.6)

holds for all m ∈ Z0.

2.2. Some results from polynomial interpolation theory.

Definition 2.4 ([18]). For a given continuous vector-function y : [0, T ]→ Rn and
a natural number q, we denote by Lqy the n-dimensional Lagrange polynomials of
degree q, such that

(Lqy)(tq
i ) = y(tq

i ), i = 1, q + 1, (2.7)

where

tqi =
T

2

[
cos

(2i− 1)π

2(q + 1)
+ 1
]
, i = 1, q + 1 (2.8)

are Chebyshev nodes, translated from the interval (−1, 1) to (0, T ), and

Lqy := col(Lqy1, L
qy2, . . . L

qyn). (2.9)

Definition 2.5 ([18]). Denote by Pq a set of all polynomials of degree not higher
than q (q ≥ N) on [0, T ]. For any continuous function y : [0, T ]→ R, there exists a
unique polynomial p∗q ∈ Pq, for which

max
t∈[0,T ]

|y(t)− p∗q | = Eq(y),

where

Eq(y) := inf
p∈Pq

max
t∈[0,T ]

|y(t)− p(t)|. (2.10)

Here p∗q is a polynomial of the best uniform approximation of y in Pq, and the
number Eq(y) is called the error of the best uniform approximation.

Proposition 2.6 ([17]). For any q ∈ N and a continuous function y : [0, T ] → R,
the corresponding interpolation polynomial (2.7) with the Chebyshev nodes (2.8)
admits the estimate

|y(t)− (Lq)y(t)| ≤
( 2

π
ln q + 1

)
Eq(t), t ∈ [0, T ]. (2.11)



4 M. FEČKAN, K. MARYNETS EJDE-2023/58

Proposition 2.7 ([16]). If y ∈ C([0, T ],R) and q ∈ N, then

Eq(t) ≤ 6ω
(
y;
T

2q

)
, (2.12)

where

ω (y; δ) := sup{|y(t)− y(s)| : {t, s} ⊂ [0, T ], |t− s| ≤ δ}, ∀δ > 0

is the modulus of continuity of a continuous function y (for more information we
refer to [15]).

Definition 2.8 ([18]). A function y : [0, T ]→ R is said to satisfy the Dini-Lipschitz
condition if its modulus of continuity has the property

lim
δ→∞

ω(y; δ) ln δ = 0.

Remark 2.9. Note, that from (2.12) it follows that

lim
q→∞

Eq(y) ln q = 0, (2.13)

for any y satisfying the Dini-Lipschitz condition. Taking into account estimate
(2.11), equality (2.13) ensures the uniform convergence of the Lagrange interpola-
tion polynomials at Chebyshev nodes for this class of functions. In particular, every
α-Hölder continuous function y : [0, T ]→ R with α > 0 satisfies the Dini-Lipschitz
condition.

3. Non-local FBVP and its parametrization

Consider a FDS of order p ∈ (0, 1):

C
0 D

p
t x(t) = f(t, x(t)), t ∈ [0, T ], x, f ∈ Rn, (3.1)

subject to the non-local multipoint boundary constraints

φ(x) = γ, (3.2)

where C
0 D

p
t (·) denotes the Caputo fractional differential operator with the lower

limit at zero [15], x : [0, T ] → Rn and f : [0, T ] × Rn → Rn are continuous vector-
functions, φ : Rn → Rn is a continuous on [0, T ] non-linear functional, evaluated at
an arbitrary number of points ti ∈ [0, T ] (i = 1, n), and γ is a given n-dimensional
vector.

Together with the BVP (3.1), (3.2) we study the parametrized two-point FBVP

C
0 D

p
t x(t) = f(t, x(t)), t ∈ [0, T ], x, f ∈ Rn,

x(0) = ξ, x(T ) = η,
(3.3)

and a perturbed IVP:

C
0 D

p
t x(t) = f(t, x(t)) + ∆, t ∈ [0, T ], x, f ∈ Rn, (3.4)

x(0) = ξ, (3.5)

where the parameters ξ, η ∈ Rn and the perturbation term ∆ are to be defined.
Using the Caputo integral operator [15], solution of the IVP (3.4), (3.5) can be

written in an integral form that reads

x(t) = ξ +
1

Γ(p)

∫ t

0

(t− s)p−1f(s, x(s))ds+
∆tp

Γ(p+ 1)
. (3.6)



EJDE-2023/58 NON-LOCAL FRACTIONAL BVPS WITH APPLICATIONS 5

To find the perturbation term ∆ we require (3.6) to also satisfy the two-point
parametrized boundary conditions (3.3). Direct calculations show that

∆(ξ, η) =
1

T p

[
Γ(p+ 1)(η − ξ)− p

∫ T

0

(T − s)p−1f(s, x(s))ds
]
. (3.7)

Substitution of (3.7) into (3.6) yields

x(t) =
[
1−

( t
T

)p
]ξ +

( t
T

)p
η +

1

Γ(p)

[ ∫ t

0

(t− s)p−1f(s, x(s))ds

−
( t
T

)p ∫ T

0

(T − s)p−1f(s, x(s))ds
]
.

(3.8)

An appropriate choice of the vector parameters ξ, η, substituted into (3.8), leads
to the exact solution of the original problem (3.1), (3.2). However, the nonlinearity
f(·, x(·)) in the right hand-side of FDS (3.1), which depends also on the unknown
function x(·), enables the integral calculation in (3.8). But if instead of the exact
values of x(·) we consider their approximation, then formula (3.8) leads to an ex-
plicit approximate solution of problem (3.1), (3.2) for particular values of ξ and η.
In the next section we give a detailed justification of this approximation process.

4. Successive approximations technique and its modification

We assume, that BVP (3.1), (3.2) satisfies the following conditions:

(A1) f is bounded: there exists M > 0 such that

|f(t, x)| ≤M, ∀(t, x) ∈ Gf := [0, T ]×Dρ, (4.1)

and f ∈ Lip(K,Dρ).
(A2) The set D is defined by closed and bounded sets Dξ and Dη ⊂ Rn as

D := {(1− θ)ξ + θη : ξ ∈ Dξ, η ∈ Dη, θ ∈ [0, 1]}. (4.2)

Additionally, for any ρ ∈ Rn+ we define a componentwise ρ-neighborhood of
the set D by

Dρ := Ξρ(D), (4.3)

where

Ξρ(D) := ∪ξ∈DΞρ(ξ),

Ξρ(ξ) := {ξ ∈ Rn : |ξ − z| ≤ ρ}for all ξ.

(A3) The spectral radius of the matrix

Q :=
KT p

22p−1Γ(p+ 1)
(4.4)

satisfies

r(Q) < 1. (4.5)

Under conditions (A1)–(A3), we construct a classical and a modified iteration
processes for approximation of solutions of the original non-local BVP (3.1), (3.2),
prove their convergence and derive a determining system for calculation of the
numerical values of the unknown parameters ξ, η, introduced in (3.3).



6 M. FEČKAN, K. MARYNETS EJDE-2023/58

4.1. Classical numerical-analytic method and its convergence. Let us con-
nect with the exact solution of the parametrized FBVP (3.1), (3.3) a sequence of
functions in an operator form

x0(t, ξ, η) :=
[
1−

( t
T

)p]
ξ +

( t
T

)p
η,

xm(t, ξ, η) := x0(t, ξ, η) + (ΛNfxm−1(·, ξ, η))(t), t ∈ [0, T ], m ∈ N,
(4.6)

where Λ stands for a linear operator, acting on the space C([0, T ]), defined by

(Λy)(t) :=
1

Γ(p)

∫ t

0

(
(t− s)p−1

[
y(s)− 1

T p

∫ T

0

(T − τ)p−1y(τ)dτ
])
ds, (4.7)

for all t ∈ [0, T ] and Nf is the Nemytskii operator generated by the nonlinearity f
in the right hand-side of (3.1),

(Nfy)(t) := f(t, y(t)), t ∈ [0, T ]. (4.8)

Similarly to the results in [9, 10, 12], where authors studied FDSs with the
two-point linear and nonlinear boundary conditions, one can prove a uniform con-
vergence of the sequence (4.6) to a parametrized limit function x∞(·, ξ, η) and its
relation to the exact solution x(t) of the original FBVP (3.1), (3.2).

Theorem 4.1. Assume, that the FBVP (3.1), (3.3) satisfies conditions (A1)–(A3).
Then for each fixed (ξ, η) ∈ Dξ ×Dη:

(1) The sequence of functions (4.6) are continuous and satisfy the two-point
parametrized boundary conditions

xm(0, ξ, η) = ξ, xm(T, ξ, η) = η.

(2) The sequence of functions (4.6) for t ∈ [0, T ] converges uniformly as
m→∞ to a parameter-dependent limit function

x∞(t, ξ, η) = lim
m→∞

xm(t, ξ, η). (4.9)

(3) The limit function (4.9) satisfies the parametrized boundary conditions

x∞(0, ξ, η) = ξ, x∞(T, ξ, η) = η.

(4) The limit function (4.9) is a unique solution to the integral equation

xm(t, ξ, η) := x0(t, ξ, η) + (ΛNfxm−1(·, ξ, η))(t), t ∈ [0, T ],

i.e., it is a unique solution on t ∈ [0, T ] of the Cauchy problem for a per-
turbed system of FDEs,

C
0 D

p
t x(t) = f(t, x(t)) + ∆(ξ, η), t ∈ [0, T ], x, f ∈ Rn, (4.10)

x(0) = ξ, (4.11)

where ∆ : Dξ ×Dη → Rn is the mapping defined by

∆(ξ, η) :=
1

T p

[
Γ(p+ 1)(η − ξ)− p

∫ T

0

(T − s)p−1Nfx(s)ds
]
. (4.12)

(5) The following error estimate holds:

|x∞(t, ξ, η)− xm(t, ξ, η)| ≤ T p

22p−1Γ(p+ 1)
Qm(In −Q)−1M, (4.13)

where M and Q are defined by (4.1) and (4.4) respectively.



EJDE-2023/58 NON-LOCAL FRACTIONAL BVPS WITH APPLICATIONS 7

Consider now a Cauchy problem

C
0 D

p
t x(t) = f(t, x(t)) + µ, t ∈ [0, T ], (4.14)

x(0) = ξ, (4.15)

where µ ∈ Rn, which we call a control parameter, and ξ ∈ Dξ.

Theorem 4.2. Let ξ ∈ Dξ, η ∈ Dη and µ ∈ Rn be given vectors. Assume that all
conditions of Theorem 4.1 are satisfied for the FDS (3.1).

Then the solution x = x(·, ξ, η, µ) of IVP (4.14), (4.15) also satisfies boundary
conditions (3.3) if and only if

µ = ∆(ξ, η), (4.16)

where ∆(ξ, η) is given by (4.12), and in this case

x(t, ξ, η, µ) = x∞(t, ξ, η), for all t ∈ [0, T ]. (4.17)

Theorem 4.3. Assume that BVP (3.1), (3.2) satisfies conditions (A1)–(A3). Then
x∞(·, ξ∗, η∗) is a solution to the FDS (3.1) with nonlinear boundary conditions (3.2)
if and only if the point (ξ∗, η∗) is a solution to the determining system

∆(ξ∗, η∗) =
1

T p

[
Γ(p+ 1)(η − ξ)− p

∫ T

0

(T − s)p−1Nfx∞(s, ξ, η)ds
]

= 0,

φ(x∞) = γ,

(4.18)

where ∆(ξ, η) is given by (4.12) and the second equation comes from the non-local
boundary condition (3.2).

Since the proofs of Theorems 4.1-4.3 overlap with their analogues in [9, 12] and do
not contain any new techniques, we leave them to the reader. In the next section
we modify iterations (4.6) by additionally interpolating them via the Lagrange
polynomials, constructed at the Chebyshev nodes. One of the main advantages of
this approach is in the polynomial form of the approximate solution that is easier
to analyze further.

4.2. Polynomial interpolation method and its convergence. Assuming that
the FBVP (3.1), (3.3) satisfies conditions (A1)–(A3), let us fix a natural number q
and introduce a modified approximation scheme:

uq0(t, ξ, η) := Lqx0(t, ξ, η),

uqm(t, ξ, η) := uq0(t, ξ, η) + (ΛqNfu
q
m−1(·, ξ, η))(t), t ∈ [0, T ], m ∈ N,

(4.19)

with

(Λqy)(t) :=
1

Γ(p)
Lq
(∫ t

0

(
(t−s)p−1

[
y(s)− 1

T p

∫ T

0

(T − τ)p−1y(τ)dτ
])
ds
)
, (4.20)

where t ∈ [0, T ], x0(t, ξ, η) is defined by (4.6), and Lq(·) are the Lagrange polyno-
mials (2.9) satisfying a component-wise property (2.7).

Similarly to (2.9), for any continuous vector-function y : [0, T ]→ Rn we put

Eqy := col(Eqy1, Eqy2, . . . , Eqyn).

Additionally, if D ⊂ Rn is a closed domain and f : [0, T ]×D → Rn, then

lq,D(f) :=
( 2

π
ln q + 1

)
sup

p∈Pq,D

Eq(Nfp), (4.21)



8 M. FEČKAN, K. MARYNETS EJDE-2023/58

where

Pq,D := {y : y ∈ Pnq , y([0, T ]) ⊂ D}, Pnq := Pq × · · · × Pq︸ ︷︷ ︸
n

.

Let us prove convergence of the iterative scheme (4.19) and establish connection
between the original approximations (4.6) and the modified sequences of functions
(4.19).

Theorem 4.4. Let there exist a non-negative vector

ρq :=
T pM

2p−1Γ(p+ 1)

[12

qp

( 2

π
ln q + 1

)
+ 1
]
, (4.22)

such that for all fixed y ∈ Dρq it holds that

f ∈ Lip(K,Dρq ).

Then, for all fixed (ξ, η) ∈ Dξ ×Dη and q ∈ N we have:

(1) For any m ∈ N0, function uqm(t, ξ, η) in (4.19) is a vector polynomial of
degree q having values in Dρq and satisfying the parametrized boundary
conditions

lim
q→∞

uqm(0, ξ, η) = ξ, lim
q→∞

uqm(T, ξ, η) = η.

(2) The sequences of functions {uqm(t, ξ, η)} and {uq∞(t, ξ, η)} converge uni-
formly to their corresponding limit functions, as m, q →∞:

uq∞(t, ξ, η) := lim
m→∞

uqm(t, ξ, η),

u∞(t, ξ, η) = x∞(t, ξ, η) := lim
q→∞

uq∞(t, ξ, η)
(4.23)

for each t ∈ [0, T ] and m ∈ N0, where x∞(t, ξ, η) is defined by (4.9).
(3) The limit function u∞(t, ξ, η) satisfies the parametrized two-point boundary

constraints

u∞(0, ξ, η) = ξ, u∞(T, ξ, η) = η.

(4) The estimate

|x∞(t, ξ, η)− uqm(t, ξ, η)| ≤ T p

22p−1Γ(p+ 1)
Qm(In −Q)−1M

+(In −Q)−1Big(lq(x0) + 24
( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M
)

+Qm+1lq(x0)

(4.24)

holds, for all t ∈ [0, T ] and m ∈ N0, where x∞(t, ξ, η) and uqm(t, ξ, η) are
given by (4.9) and (4.19) respectively.

Proof. Let us fix values of parameters ξ ∈ Dξ and η ∈ Dη and show that

{uqm(t, ξ, η) : t ∈ [0, T ]} ⊂ Dρq . (4.25)

From the first expression in (4.19) it follows that (4.25) holds for m = 0. Let us
now prove (4.25) for all m ≥ 1. For this purpose we first derive some intermediate
estimates.

By putting v := Nfu
q
m−1 we compute

|(Λv)(t)− (Λqv)(t)| = |(Λv)(t)− (Lq)(Λv)(t)| ≤ 6
( 2

π
ln q + 1

)
ω
(

Λv;
T

2q

)
(4.26)



EJDE-2023/58 NON-LOCAL FRACTIONAL BVPS WITH APPLICATIONS 9

Next we evaluate the difference |Λv(t2)− Λv(t1)|, for t1, t2 ∈ [0, T ], t1 < t2:

|Λv(t2)− Λv(t1)|

≤
∣∣ 1

Γ(p)

∫ t2

0

(t2 − s)p−1v(s)ds− 1

Γ(p)

∫ t1

0

(t1 − s)p−1v(s))ds
∣∣

+
∣∣ tp2 − tp1
T p

∫ T

0

(T − τ)p−1v(τ)dτ
∣∣ds

≤ 1

Γ(p)

∣∣ ∫ t2

t1

(t2 − s)p−1v(s)ds
∣∣+

1

Γ(p)

∣∣ ∫ t1

0

(
(t2 − s)p−1 − (t1 − s)p−1

)
v(s)ds

∣∣
+

tp2 − t
p
1

Γ(p+ 1)
‖v‖0

≤ 2
(t2 − t1)p + (tp2 − t

p
1)

Γ(p+ 1)
‖v‖0

≤ 4
(t2 − t1)p

Γ(p+ 1)
‖v‖0,

since

tp2 − t
p
1 ≤ (t2 − t1)p.

This leads to the following estimate of the term ω(Λv; T2q ) in (4.26):

ω
(

Λv;
T

2q

)
≤ 4

T p

2pqpΓ(p+ 1)
‖v‖0.

Consequently, (4.26) gives

‖Λv − Λqv‖0 ≤ 24
( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
‖v‖0, (4.27)

that is

‖Λ− Λq‖ ≤ 24
( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
.

Now, from the iterative formula (4.19), estimates (4.27) and (2.2), we obtain

|uqm − u
q
0| = |ΛqNfu

q
m−1| ≤ |ΛqNfu

q
m−1 − ΛNfu

q
m−1|+ |ΛNfum−1|

= |(Λq − Λ)Nfum−1|+ |ΛNfum−1|

≤ 24
( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M +

T p

2p−1Γ(p+ 1)
M

=
T pM

2p−1Γ(p+ 1)

[12

qp

( 2

π
ln q + 1

)
+ 1
]

= ρq.

(4.28)

Taking into account definition (4.3) for the set Dρq , inequality (4.28) proves that
(4.25) holds, for all m ∈ N0.

To prove that iterations (4.19) converge to (4.9) as q tends to infinity let us
estimate a difference,

|xm+1 − uqm+1|
≤ |x0 − uq0|+ |ΛNfxm − ΛqNfu

q
m|

≤ |x0 − uq0|+ |Λ(Nfxm −Nfuqm)|+ |(Λ− Λq)Nfu
q
m|

≤ lq(x0) +KI(|xm − uqm|) + 24
( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M,

(4.29)



10 M. FEČKAN, K. MARYNETS EJDE-2023/58

where we used inequality (4.27), notation (4.21) and Lemma 2.2 for

|Λ(Nfxm −Nfuqm)| ≤ I(|Nfxm −Nfuqm|) ≤ KI(|xm − uqm|).
The method of mathematical induction leads to the following results:
• for m = 0:

|x1 − uq1| ≤ lq(x0) + 24
( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M +Klq(x0)α1(t)

≤ lq(x0) + 24
( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M +Qlq(x0);

• for m = 1:

|x2 − uq2| ≤ lq(x0) + 24
( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M

+K
{
lq(x0) + 24

( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M
}
α1(t) +K2lq(x0)α2(t)

≤ Q0
{
lq(x0) + 24

( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M
}

+Q
{
lq(x0) + 24

( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M
}

+Q2lq(x0);

• and thus, for a general m:

|xm+1−uqm+1| ≤ (In−Q)−1
(
lq(x0)+24

( 2

π
ln q+1

) T p

2pqpΓ(p+ 1)
M
)

+Qm+1lq(x0).

Coupling the estimates above with inequality (4.13), we obtain

|x∞(t, ξ, η)− uqm(t, ξ, η)|
= |x∞(t, ξ, η)− xm(t, ξ, η)|+ |xm(t, ξ, η)− uqm(t, ξ, η)|

≤ T p

22p−1Γ(p+ 1)
Qm(In −Q)−1M

+ (In −Q)−1
(
lq(x0) + 24

( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M
)

+Qm+1lq(x0).

For large q from the second term in the inequality above we obtain

lim
q→∞

( 2

π
ln q + 1

) T p

2pqpΓ(p+ 1)
M = 0,

assuming conditions (A1)–(A3) and Proposition 2.6 hold. Thus, (4.19) is convergent
to (3.8) as m→∞ and its limit for q →∞ is near the solution x(t) = x∞(t, ξ∗, η∗)
of the original equation. This completes the proof. �

Since, according to Theorem 4.4, the limit functions of the sequences xm(·, ξ, η)
and uqm(·, ξ, η) converge to the same limit, Theorem 4.3 holds also for the iteration
process (4.19) and reads as follows.

Theorem 4.5. Let the non-local BVP (3.1), (3.2) satisfy conditions (A1)–(A3).
Then u∞(·, ξ∗, η∗) is a solution to the FDS (3.1) with non-local boundary conditions
(3.2), if and only if the point (ξ∗, η∗) is a solution to the determining system

∆q(ξ∗, η∗) = 0, (4.30)

φ(u∞) = γ, (4.31)



EJDE-2023/58 NON-LOCAL FRACTIONAL BVPS WITH APPLICATIONS 11

where

∆q(ξ∗, η∗) =
1

T p

[
Γ(p+ 1)(η − ξ)− p

∫ T

0

(T − s)p−1Nfu∞(·, ξ, η))(s)ds
]
,

and the second equation comes from the boundary condition (3.2).

Remark 4.6. (1) For practical reasons, instead of the exact system (4.30), (4.31)
we introduce an approximate determining system of the form

∆q
m(ξ, η) = 0,

φ(uqm) = γ,
(4.32)

where the map ∆q
m : Dξ ×Dη → Rn is defined as

∆q
m(ξ, η) :=

1

T p

[
Γ(p+ 1)(η − ξ)− p

∫ T

0

(T − s)p−1Nq
fum(·, ξ, η)(s)ds

]
.

(2) On every iteration step m, solutions (ξ̄m, η̄m) of the system (4.32) stand for
the m-th approximation to their exact values (ξ∗, η∗). By plugging (ξ̄m, η̄m) into the
sequence of functions (4.19) we obtain the m-th approximation to the exact solution
of the parametrized problem (3.1), (3.3), or of the original non-local problem (3.1),
(3.2). This solution is then given by

Xq
m(t) = uqm(t, ξ̄m, η̄m).

5. Predator-prey model with prey refuge

Predator-prey models are often used in modeling of biological and ecological
dynamical systems. They are usually written in the form of systems of ordinary
differential equations and are driven by a set of parameters that characterize a
physical system under consideration.

Recently the predator-prey models were modified to incorporate also a mem-
ory effect. For this purpose the integer order derivatives were replaced by their
fractional analogues. In this Section we apply the iteration technique, developed
earlier, to find an approximate solution to a fractional order predator-prey model
with prey refuge. For details about this model we refer the reader to [7].

We consider a dynamical system

C
0 D

p
t x(t) = rx

(
1− x

k

)
− c(1−m)xy,

C
0 D

p
t y(t) = ec(1−m)xy − dy,

(5.1)

for t ∈ [0, T ] and p ∈ (0, 1), where t stands for time, x(t), y(t) correspond to
the predator and prey densities, and r, k,m, e, d are positive constant parameters,
driving the system.

Under the set of physically relevant parameter values T = 4/5, r = 1.2, k = 40,
c = 1, d = 0.4, e = 0.2, m = 0.1, and p = 0.98, the FDS (5.1) is written as

C
0 D

0.98
t x(t) = 1.2x

(
1− x

40

)
− 0.9xy,

C
0 D

0.98
t y(t) = 0.18xy − 0.4y,

(5.2)

for t ∈ [0, 4/5].



12 M. FEČKAN, K. MARYNETS EJDE-2023/58

We aim to construct an approximate solution to system (5.2) that satisfies the
non-local boundary conditions

x2(0)− x
(4

5

)
y(0) = 2.54,

y2
(4

5

)
x(0) = 2.38

(5.3)

in the domain
Dρ = {(x, y) : 2.2 ≤ x ≤ 3.1, 1 ≤ y ≤ 1.5}. (5.4)

The parametrized constraints (3.3) in this case have the form:[
x(0)
y(0)

]
=

[
ξ1
ξ2

]
:= ξ,[

x(1)
y(1)

]
=

[
η1

η2

]
:= η

(5.5)

with

ξ ∈ Dξ := {ξ ∈ R2 : 2.2 ≤ ξ ≤ 2.5},
η ∈ Dη := {η ∈ R2 : 0.5 ≤ η ≤ 1.5}.

Calculations show that conditions (A1)–(A3) hold with

M =

[
0.6417
0.237

]
, K =

[
1.736 2.79
0.27 0.958

]
,

ρq =

[
1.819303410
0.1649011626

]
, r

(
Q =

KT p

22p−1Γ(p+ 1)

)
< 0.96.

Let us now use a modified version of the numerical-analytic method (4.6), where
the Lagrange polynomial interpolation (LPI) scheme is applied to obtain approxi-
mate solutions to the problem in the form of the q-th order polynomials. For this
purpose we fix the order of interpolation at q = 4 and compute the Chebyshev
nodes (2.8),

t1 = 0.7804226065, t2 = 0.6351141009, t3 = 0.4,

t4 = 0.1648858991, t5 = 0.0195773935.

On every iteration step we find numerical values of the unknown parameters
as roots of the approximate determining system (4.32). These values are given in
Table 1 and they determine approximate solutions to the problem (5.2), (5.3).

Table 1. Numerical values of the unknown parameters for m = 0, 1, 2, 3

Approx. Numerical value of ξ Approx. Numerical value of η

ξ1,0 2.303032253 η1,0 2.928500653
ξ1,1 2.308072335 η1,1 2.971104555
ξ1,2 2.308081308 η1,2 2.970853059
ξ1,3 2.308127313 η1,3 2.971164378
ξ2,0 0.9445718795 η2,0 1.016462968
ξ2,1 0.9389864677 η2,1 1.015336111
ξ2,2 0.9390368098 η2,2 1.015339844
ξ2,3 0.9390134110 η2,3 1.015328673



EJDE-2023/58 NON-LOCAL FRACTIONAL BVPS WITH APPLICATIONS 13

Below we also give explicit forms of these solutions and plot the estimates of our
computations. On the zero-th iteration step we obtain the following approximation
to the exact solution of the FBVP (5.2), (5.3):

Xq
0 (t) = −0.01956900249t4 + 0.0549367687t3 − 0.0646567198t2

+ 0.654212958t+ 2.303545749,

Y q0 (t) = −0.00224924845t4 + 0.00631442252t3 − 0.00743160017t2

+ 0.0751949749t+ 0.9446309006.

A comparison of the left and right hand-sides of the system (5.2) is given on Figure 1.

Figure 1. Comparison of the left and right hand-sides of (5.2) in
the zero-th approximation (LPI method)

Computations show that the first approximation to the exact solution of the
FBVP (5.2), (5.3) is given by a system of polynomials:

Xq
1 (t) = −0.02112540606t4 + 0.04090489985t3 − 0.0802533673t2

+ 0.722904954t+ 2.308639591,

Y q1 (t) = 0.00073293971t4 + 0.0000320947t3 + 0.05954191818t2

+ 0.0160367660t+ 0.9389937391.

A comparison of the left and right hand-sides of the system (5.2) in the first ap-
proximation is given on Figure 2.

Note, that the higher order approximations show an even better accuracy of our
computations. These results can be seen on Figures 3 and 4.

The error functions calculated on the second and the third iterations for the LPI
method, are plotted on Figures 5 and 6.

We additionally depict all 4 approximations to the exact solution of FBVP (5.2),
(5.3) (see Figure 7), from which the convergence behavior of our approximations
follows straightaway.



14 M. FEČKAN, K. MARYNETS EJDE-2023/58

Figure 2. Comparison of the left and right hand-sides of (5.2) in
the first approximation (LPI method)

Figure 3. Comparison of the left and right hand-sides of (5.2) in
the second approximation (LPI method)

Figure 4. Comparison of the left and right hand-sides of (5.2) in
the third approximation (LPI method)



EJDE-2023/58 NON-LOCAL FRACTIONAL BVPS WITH APPLICATIONS 15

Figure 5. Error functions in the second approximation (using the
LPI method)

Figure 6. Error functions in the third approximation (LPI method)

Figure 7. Approximations to the exact solution of (5.2), (5.3) for
m = 0, 1, 2, 3, obtained by the LPI method



16 M. FEČKAN, K. MARYNETS EJDE-2023/58

6. Conclusions

To summarize, we want to mention the main advantages of using the LPI method
and possibilities for future developments:

(1) it is easier to analyse the behaviour of polynomial solutions;
(2) increasing the number of nodes q improves the speed of convergence;
(3) the method is easy to implement using mathematical software;
(4) one could think of interpolation of functions under the integrals as well.

However, this approach needs other than Lagrange fundamental polynomi-
als due to a necessity to control the growth of the polynomial for large q.

Acknowledgments. This work is partially supported by the Slovak Grant Agency
VEGA No. 2/0127/20 and No. 1/0084/23, and by the National Scholarship Program
of the Slovak Republic.

References

[1] K. Diethelm; The Analysis of Fractional Differential Equations: An Application-Oriented

Exposition Using Differential Operators of Caputo Type, Springer, Berlin, 2010.

[2] M. Fečkan, K. Marynets; Approximation approach to periodic BVP for fractional differential
systems, Eur. Phys. J.: Special Topics 226 (2017), 3681–3692.

[3] M. Fečkan, K. Marynets; Approximation approach to periodic BVP for mixed fractional
differential systems, J. Comput. Appl. Math., 339 (2018), 208–217.

[4] M. Fečkan, K. Marynets, J. R. Wang; Periodic boundary value problems for higher order

fractional differential systems, Math. Methods Appl. Sci., 42 (2019), 3616–3632.
[5] A. A. Kilbas, S. A. Marzan; Cauchy problem for a differential equation with Caputo deriva-

tive, Fract. Calc. Appl. Anal., 7(3) (2004), 297–321.

[6] A. Kilbas, H. Srivastava, J. Trujillo; Theory and applications of fractional differential equa-
tions, Elsevier, Amsterdam, The Netherlands, 2006.

[7] H. L. Li, L. L. Zhang, C. Hu et al.; Dynamical analysis of a fractional-order predator-prey

model incorporating a prey refuge, J. Appl. Math. Comput. 54 (2017), 435–449.
[8] K. Marynets; On one interpolation type fractional boundary–value problem, Axioms: SI

Fractional Calculus, Wavelets and Fractals, 9(1) (2020), 1–21.

[9] K. Marynets, D. Pantova; Approximation approach to the fractional BVP with the Dirichlet
type boundary conditions. Differ. Equ. Dyn. Syst. (2022).

[10] K. Marynets, D. Pantova; Successive approximations and interval halving for fractional BVPs
with integral boundary conditions, J. of Comp. and Appl. Math., (2023), 1–20.

[11] K. Marynets; Solvability analysis of a special type fractional differential system, Comput.
Appl. Math. (2019), 13p.

[12] K. Marynets; Successive approximations technique in study of a nonlinear fractional boundary

value problem, Mathematics 9(7) (2021), 1–19.

[13] I. P. Natanson; Constructive function theory. Vol. III. Interpolation and approximation
quadratures, New York: Frederick Ungar Publishing Co., 1965.

[14] G. M. Phillips; Interpolation and approximation by polynomials, CMS Books in Mathematics,
Vol. 14, Springer-Verlag, New York, 2003.

[15] I. Podlubny; Fractional differential equations, Academic Press, New York, 1999.

[16] T. J. Rivlin; An introduction to the approximation of functions, Blaisdell Publishing Co.

Ginn and Co., Waltham, Mass.-Toronto, Ont.-Lindon, 1969.
[17] T. J. Rivlin; The Chebyshev polynomials, Wiley-Interscience, NY-London-Syndey, 1974.

[18] A. Rontó, M. Rontó, N. Shchobak; Parametrization for boundary value problems with tran-
scendental non-linearities using polynomial interpolation, Electron. J. Qual. Theory Differ.

Equ., (2018) No. 59 (2018), 1–22.

[19] M. I. Ronto, A. M. Samoilenko; Numerical-analytic methods in the theory of boundary value
problem, World Scientific Publishing Co., Inc., River Edge, 2000.

[20] T. Sandev, Ž. Tomovski, J. L. A. Dubbeldam at al.; Generalized diffusion-wave equation with

different memory kernels, J. Phys. A: Math. Theor., 52(1) (2019), 1–23.

[21] T. Sandev, Ž. Tomovski; Fractional Equations and Models, Springer Nature, 2019.



EJDE-2023/58 NON-LOCAL FRACTIONAL BVPS WITH APPLICATIONS 17

[22] Ž. Tomovski, J. L. A. Dubbeldam, J. Korbel; Applications of Hilfer-Prabhakar operator to

option pricing financial model, Fract. Calc. Appl. Anal., 23(4) (2020), 996–1012.

[23] Ž. Tomovski; Generalized Cauchy type problems for nonlinear fractional differential equations
with composite fractional derivative operator, Nonlinear Anal., Theory Methods Appl., 75(7)

(2012), 3364–3384.
[24] Y. Zhou; Basic Theory of Fractional Differential Equations, World Scientifc, Singapore

(2014).

Michal Fečkan

Department of Mathematical Analysis and Numerical Analysis, Comenius University in
Bratislava, Mlynská dolina, 842 48 Bratislava, Slovakia.

Institute of Mathematics, Slovak Academy of Sciences, Štefánikova 49, 814 73 Bratislava,

Slovakia
Email address: michal.feckan@fmph.uniba.sk

Kateryna Marynets
Delft Institute of Applied Mathematics, Faculty of Electrical Engineering, Mathe-

matics and Computer Science, Delft University of Technology, Mekelweg 4 2628CD

Delft, Netherlands
Email address: K.Marynets@tudelft.nl


	1. Introduction
	2. Notation, auxiliary statements and definitions
	2.1. Auxiliary lemmas
	2.2. Some results from polynomial interpolation theory

	3. Non-local FBVP and its parametrization
	4. Successive approximations technique and its modification
	4.1. Classical numerical-analytic method and its convergence
	4.2. Polynomial interpolation method and its convergence

	5. Predator-prey model with prey refuge
	6. Conclusions
	Acknowledgments

	References