Adv Syst Sci Appl 2020; 02:56–70
Published online at https://ijassa.ipu.ru.

Numerical Methods for Constructing Solutions of Functional
Differential Equations of Pointwise Type

Armen Beklaryan1*

1National Research University Higher School of Economics, Moscow, Russia

Abstract: The article discusses construction of traveling wave type solutions for the Frenkel-
Kontorova model on the propagation of longitudinal waves. For the first time, based on the
existence and uniqueness theorem of traveling wave type solutions, as well as the approximation
theorem, a complete family of traveling wave type solutions is constructed in the form of
subfamilies of bounded solutions (horizontals) and unbounded solutions (verticals).

Keywords: traveling waves, functional differential equations, equations of mathematical
physics, splines, forward-backward differential equation

1. INTRODUCTION

The theory of differential equations with delay, which has developed rapidly in recent
decades [1–3], has in many ways acquired a complete form and is now actively used in
modeling various objects. Numerical algorithms for solving equations with delay of various
types were also developed [4, 5]. At the same time, numerical methods for equations with
advanced arguments (and, moreover, with mixed type deviations) are practically not studied,
although references to them have been encountered for a long time, mainly in connection
with the classification of equations with deviating argument [6]. As a rule, there is a parsing
of equations of a particular kind with further obtaining the existence and uniqueness theorems
based on the use of the properties of the right-hand side [7] and application of methods, such
as the study of the roots of the characteristic quasi-polynomial [8], collocation methods and
finite element scheme (expansion of a solution in terms of basis functions of some finite-
dimensional space) [9–11] or through the construction of a Hilbert space of the reproducing
kernels on the basis of boundary conditions [12].

For equations of a delayed type, as a rule, a well-known initial problem is considered when
the initial moment of time coincides with the left end of the interval of equation’s domain.
Such a statement is characterized by the fact that the initial-boundary conditions are local in
nature, and the equation itself can be integrated by the step method. In all other cases, and
even more so in the case of functional differential equations of pointwise type (FDEPT) with
mixed-type deviations, the initial-boundary conditions are nonlocal in nature.

One of the approaches to the construction of numerical methods for FDEPT is the
variational method based on solving the induced optimal control problem (OCP) as an
unconditional optimization problem for the residual functional. As noted earlier, the initial-
boundary conditions for FDEPT are not local, therefore, the solution of the variational
problem for the residual functional cannot be obtained by simple local improvements and
requires global optimization. Note that for problems of finite-dimensional optimization there

∗Corresponding author: abeklaryan@hse.ru



NUMERICAL METHODS FOR FUNCTIONAL DIFFERENTIAL EQUATIONS 57

are many approaches to finding a global extremum, sufficiently efficient algorithms are built,
there are a large number of publications that present precedents for successfully solved
problems [13, 14, for example].

At the same time, the task of finding global extremum for OCP remains one of the
most acute tasks in the extreme problems theory. Among classical approaches to this task
one cannot fail to mention the works [15] and [16]. Some works that apply search and
genetic algorithms are also well-known [17, 18]. But unfortunately these efforts have not yet
brought about effective algorithms that would be able to solve a wide spectrum of practical
problems [19].

One of the ideas in this area is the idea of reducing the optimal control problem to a
low-dimensional problem on the reachability set of a controlled system [20]. Unfortunately,
the problem of approximating the reachable set in the numerical solution is not much
simpler than the problem of finding a global extremum. Nevertheless, the properties of the
reachable set known from theoretical studies can be used to construct specialized optimization
algorithms for OCP [21]. In particular, the connectivity property of the reachability set allows
one to construct a scheme of continuous control variation that leads to the solution of the
problem. The idea of using the connectivity property of the reachability set when constructing
numerical optimization algorithms belongs to A.G. Chentsov [22].

In the majority of well-known approaches to the construction of non-convex optimization
methods, the solution of the problem is divided into two stages: “global”, where a wide scan
of the variable space is performed, and “local”, aimed on local refinement of the obtained
solution [14]. The combination of different methods at each stage, as well as the order of the
alternation of stages, determines the specific computational algorithm. For mentioned class
of problems, we consider a idea of combining both stages in one algorithm package and
creating a method that allows one to select and change both the “global scan” procedure and
local improvement of the obtained approximate solution at each iteration [23].

2. STATEMENT OF THE INITIAL BOUNDARY-VALUE PROBLEM

The most important goal in FDEPT is the study of the basic initial-boundary value problem

ẋ(t) = f(t, x(q1(t)), . . . , x(qs(t))), t ∈ BR, (2.1)
ẋ(t) = ϕ(t), t ∈ R\BR, ϕ(·) ∈ L∞(R,Rn), (2.2)
x(t̄) = x̄, t̄ ∈ R, x̄ ∈ Rn. (2.3)

where f : R× Rns −→ Rn is a mapping of the C(0) class, qj(·), j = 1, . . . , s, represent
diffeomorphisms of the line preserving orientation, andBR is either the closed interval [t0, t1],
or closed half-line [t0,+∞), or the whole line R.

In the case of a delayed-type equation, for BR = [t0, t1], [t0,+∞), t̄ = t0, the initial-
boundary-value problem becomes a well-known initial problem with the initial function
ζ(t) = x̄−

∫ t̄
t
ϕ(τ)dτ, t ≤ t̄.

The statement of the initial-boundary value problem (2.1)-(2.3) is correct without any
restrictions on the type of argument deviations.

Definition 2.1:
The solution to the basic initial-boundary value problem is any absolutely continuous solution
of the equation (2.1) satisfying the boundary condition (2.2) and the initial condition (2.3).

The aim of the FDE research is to study the solution space of the boundary value problem
(2.1)-(2.2) for each given boundary function ϕ, as well as to describe the obstacles due to
which the solutions of the boundary value problem (2.1)-(2.2) do not inherit the properties of
solutions of ordinary differential equations.

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



58 A. BEKLARYAN

If we introduce notation for a fixed function ϕ(·)

fϕ(t, z1, . . . , zs) =

{
f(t, z1, . . . , zs), t ∈ BR,
ϕ(t), t /∈ BR,

the boundary value problem (2.1)-(2.3) becomes a Cauchy problem

ẋ(t) = fϕ(t, x(q1(t)), . . . , x(qs(t))), t ∈ R (2.4)
x(t̄) = x̄, t̄ ∈ R, x̄ ∈ Rn. (2.5)

Obviously, in the case BR = R, the equality fϕ = f holds.
In the study of the Cauchy problem (2.4)-(2.5), the role of the group < q1, . . . , qs > is

very important. Often, instead of such a group, it is useful to consider some wider finitely
generated group Q of diffeomorphisms of the line, that is, < q1, . . . , qs >⊆ Q. Moreover, we
assume that the condition hq = supt∈R |q(t)− t| < +∞ is satisfied for all elements q of the
group Q.

Let’s formulate a system of restrictions on the right-hand side f : R× Rn·s 7−→ Rn of
FDEPT (2.1), and on diffeomorphisms qj(·), j = 1, . . . , s:

(a) f(·) ∈ C(0)(R× Rn·s,Rn);
(b) for any t, zj, z̄j, j = 1, s, there is a quasilinear growth

‖f (t, z1, . . . , zs) ‖Rn ≤M0(t) +M1

s∑
j=1

‖zj‖Rn , M0(·) ∈ C(0)(R,R),

and a Lipschitz condition

‖f (t, z1, . . . , zs)− f (t, z̄1, . . . , z̄s) ‖Rn ≤ Lf

s∑
j=1

‖zj − z̄j‖Rn ;

(c) there exists µ∗ ∈ R+ such that

M0(·) ∈ Lnµ∗C(0)(R);

(d) the values
hqj = sup

t∈R
|t− qj(t)|, j = 1, . . . , s,

are finite;
(e) with the constant µ∗ ∈ R+ from condition (c) the family of functions

f̃q,z1,...,zs(t) = f(q(t), z1, . . . , zs)(µ
?)hq , q ∈ Q, z1, . . . , zs ∈ Rn,

is equicontinuous on any finite interval.

The continuity condition, growth conditions with respect to the phase variables and time
variable, and the Lipschitz condition (conditions (a)− (b)) are standard in the theory of
ordinary differential equations. In fact, in item (b), the first inequality in the form of the
growth condition with respect to the phase variables and time variable is a consequence of
the second condition in the form of the Lipschitz condition. Under the Lipschitz constant
Lf , the minimum value among the possible values of these constants should be understood.
Accordingly, we can assume that M1 = Lf . Moreover, we wrote the first inequality out
separately in order to formulate condition (c) for functionM0(·). Condition (c) for function f
is related to the study of solutions on the half-line and line, which requires certain restrictions

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



NUMERICAL METHODS FOR FUNCTIONAL DIFFERENTIAL EQUATIONS 59

on the time growth of the right-hand side. It should be noted that we can always satisfy
condition (d) by making a time change, but we may break condition (c). The last condition (e)
is necessary for the right-hand side of the induced infinite-dimensional ordinary differential
equation, with the phase space in a suitable Banach space, to easily establish the fact of its
Bochner integrability. In fact, condition (e) can be removed, but this leads to further technical
complications, as the right-hand side of the induced infinite-dimensional ordinary differential
equation will only be measurable, and it will be necessary to establish the fact of its Bochner
integrability [24].

The right-hand side f(·) of FDEPT will be considered as an element of the Banach space
Vµ∗(R× Rns,Rn) with Lipschitz norm

Vµ∗(R× Rns,Rn) =
{
f(·) : f(·) satisfies the conditions (a)-(d)

}
,

‖f(·)‖Lip = sup
t∈R
‖f(t, 0, . . . , 0)(µ∗)|t|‖Rn +

+ sup
(t,z1,...,zs,z̄1,...,z̄s)∈R1+2ns

‖f(t, z1, . . . , zs)− f(t, z̄1, . . . , z̄s)‖Rn∑s
j=1 ‖zj − z̄j‖Rn

,

where the parameter µ∗ ∈ R+ coincides with the corresponding constant from the condition
(c).

Obviously, for the function f(·) ∈ Vµ∗(R× Rns,Rn) the smallest value of the constant
Lf from the Lipschitz condition (the condition (b)) coincides with the value of the second
summand in the definition of the norm f(·). In what follows, speaking of the Lipschitz
condition, by the constant Lf we mean exactly its smallest value.

The right-hand side of the FDEPT (2.1) uniquely defines a pair (f(·), h), where f(·)
is an element of the Banach space Vµ∗(R× Rns,Rn), and h = (hq1 , . . . , hqs), hqj ≥ 0, j =
1, . . . , s, are maximum deviation of the argument from the condition (d), which we will
consider as parameters. Therefore, such a right-hand side of the equation will uniquely
determine the pair (Lf ;h).

We will seek a solution to the FDEPT with a quasilinear right-hand side in a one-
parameter family of Banach spaces of functions with a given exponential growth. The
exponent is the parameter of the selected family of functions, which is defined as follows

LnµC(k)(R) =

{
x(·) : x(·) ∈ C(k) (R,Rn) , max

0≤r≤k
sup
t∈R
‖x(r)(t)µ|t|‖Rn < +∞

}
,

‖x(·)‖(k)
µ = max

0≤r≤k
sup
t∈R
‖x(r)(t)µ|t|‖Rn , k = 0, 1, . . . , µ ∈ (0,+∞).

2.1. Existence and Uniqueness Theorem for the Basic Initial-Boundary Value Problem

In terms of the parameter µ ∈ (0, 1) and the space LnµC(0)(R) we formulate a condition
guaranteeing the existence and uniqueness of a solution for the initially-boundary problem.
Theorem 2.1 ( [25]):
Let function f and diffeomorphisms qj(·), j = 1, . . . , s, satisfy the conditions (a)− (e) from
the Section 2. If for some µ ∈ (0, µ∗) ∩ (0, 1) the inequality

Lf

s∑
j=1

µ−hqj < lnµ−1, (2.6)

is satisfied then for any fixed initial-boundary conditions x̄ ∈ Rn, ϕ(·) ∈ L∞(R,Rn)
there exists a solution (absolutely continuous) x(·) ∈ LnµC(0)(R) of the basic initial-
boundary value problem (2.1)-(2.3). Such a solution is unique and, as an element of the

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



60 A. BEKLARYAN

space LnµC(0)(R), depends continuously on the initial-boundary conditions x̄ ∈ Rn, ϕ(·) ∈
L∞(R,Rn) and the right-hand side of the equation (function f(·)).

The condition (2.6) is exact and not improvable: there are differential equations for which
violation of this condition leads to a violation of either the existence of the solution or
uniqueness. Note that if the inequality (2.6) has a solution, then there are µ1(Lf ;h), µ2(Lf ;h)
such that for each solution µ the following condition holds

µ1(Lf ;h) < µ < µ2(Lf ;h).

3. VARIATIONAL APPROACH TO CONSTRUCTING SOLUTIONS TO THE
BASIC INITIAL-BOUNDARY VALUE PROBLEM

We study absolutely continuous solutions of the system

Fi(t, ẋ(t), x(t+ τ1), . . . , x(t+ τs)) = 0, i = 1, k, t ∈ [tl, tr], (3.7)

under boundary conditions outside the system definition interval

ẋ(t) = ϕl(t), t ∈ [tll, tl], (3.8)
ẋ(t) = ϕr(t), t ∈ [tr, trr], (3.9)

as well as initial-boundary conditions on the system definition interval

Km(ẋ(ξ), x(ξ1), . . . , x(ξp)) = 0, (3.10)

where Fi : R× Rn × Rns −→ Rn, i = 1, k, is a mapping of the C(0) class; τj ∈ Z, j = 1, s;
tl, tr ∈ R; tll = tl + min{0, τ1, . . . , τs}, trr = tr + max{0, τ1, . . . , τs}; ϕl, ϕr : Rn −→ Rn

are mappings of theC(0) class, defining fixed boundary vector functions;Km : Rn × Rnp −→
Rn,m = 1, q, is a mapping of the C(0) class; ξ, ξi ∈ [tl, tr], i = 1, p, is a fixed set of points on
the system definition interval.

We give a formal statement of the optimization problem induced by the initial-boundary-
value problem.

Problem A. Find the trajectory x̂(t), which provides a minimum of residual functional

I(x̂(t)) = v(N)

(
k∑
i=1

∫ tr

tl

F 2
i (t, ˙̂x(t), x̂(t+ τ1), . . . , x̂(t+ τs))dt+∫ tl

tll

[ ˙̂x(t)− ϕl(t)]2dt+

∫ trr

tr

[ ˙̂x(t)− ϕr(t)]2dt
)

+

v(K)

q∑
m=1

K2
m( ˙̂x(ξ), x̂(ξ1), . . . , x̂(ξp)), (3.11)

where v(N), v(K) ∈ R+ are weighting factors.
The following statement holds.

Proposition 3.1:
Each solution to the initial-boundary value problem (3.7)-(3.10) is a solution to the
optimization Problem A.

In particular, under the conditions of Theorem 2.1, a solution to the initial-boundary value
problem (2.1)-(2.3) exists and is a solution to the corresponding Problem A. In the general

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



NUMERICAL METHODS FOR FUNCTIONAL DIFFERENTIAL EQUATIONS 61

case, the converse statement to Proposition 3.1 is false, but, nevertheless, taking into account
the uniqueness theorem and also the stability theorem of the solution, in the case of successful
construction of a numerical solution, we can speak with high confidence about constructing
a solution close to analytical one.

The proposed variational approach to constructing solutions of the basic initial-boundary
value problem lies at the basis of the numerical implementation of such solutions. The
numerical methods themselves are based on the Ritz method and spline collocation
constructions and were implemented in [23, 26]. In order to solve the problem of the class
under consideration, the trajectories of the system are discretized on a grid with a constant
step, and a generalized residual functional is formulated that includes both the weighted
residual of the original differential equation and the residual of the boundary conditions. A
spline differentiation technique is used to evaluate the derivatives of the desired trajectories
of the system, based on two spline approximation designs: using cubic natural splines and
using a special type of spline whose second derivatives at the edges are also controlled using
optimized parameters.

3.1. Software Package OPTCON-F
A set of algorithms for local and global optimization was implemented for solving the stated
finite-dimensional problems. The used technology includes: an algorithm for sequentially
increasing the accuracy of approximation by multiplying the number of nodes in the grid of
the discretization; algorithms for the difference evaluation of the derivatives of the functional
from the first to the sixth degree of accuracy inclusive; method of successively increasing the
precision of spline differentiation. The corresponding software complex (SC) OPTCON-F
was implemented in the language C under the control of operating systems OS Windows, OS
Linux and Mac OS using compilers BCC 5.5 and GCC. SC was designed to obtain a numerical
solution of boundary value problems, parametric identification problems and optimal control
for dynamical systems described by FDEPT [27].

Among the local algorithms included in the SC OPTCON-F there are: 1) PARTAN
method; 2) Powell-Brent’s method; 3) gradient method of confidence intervals; 4) Barzilai-
Borwein method; 5) Newton’s method with the difference estimate of the Hessian matrix;
6) generalized quasi-Newtonian method; 7) direct-dual method of gradient descent; 8)
differential Euler optimization method of the 2nd order; 9) differential Adams optimization
method of the 4th order and others.

As algorithms of “closers” there are: 1) adaptive modification of the Hooke-Jeeves
method; 2) stochastic search methods in random subspaces of the indicated (2, 3, 4, or 5)
dimension; 3) local version of the curvilinear search method.

Non-local algorithms that form the basis of the SC are: 1) “parabolic” method – a
combination of coordinate-wise descent with a periodic multistart and a non-local one-
dimensional search by the parabola algorithm; 2) non-local method of curvilinear search;
3) Luus-Jaakola method; 4) “forest” method – multivariant adaptive method of random
multistart and others.

As part of the work on the SC, taking into account the specifics related to the qualitative
properties of FDEPT, the following key results were obtained:

• Algorithm of construction of “controllable splines” has been developed. The
problem of obtaining a high-precision approximation to the derivative of a function
of one variable over a set of values of the function itself, defined on a fixed grid, was
investigated. The performed computational experiments showed that in order to achieve
good accuracy in estimating derived trajectories for FDEPT systems, it is not possible
to use known types of splines (“natural”, with additional boundary conditions, Akima
splines, etc.). The proposed algorithm is based on the use of derivatives of cubic spline
functions.

• A technology has been developed for approximation of general FDEPT systems
using a finite-dimensional unconditional minimization problem. By analyzing the

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



62 A. BEKLARYAN

behavior of homeomorphisms on the initial (main) time interval, an extended time
interval is calculated, for which the discretized grid is constructed. To approximate
the initial continuous problem on a fixed grid over time, the Ritz method is used: the
trajectories are approximated using controlled spline functions, the coefficients of which
are selected by searching for the minimum of the residual functional.

• A specialized global optimization algorithm was developed, based on the idea of
curvilinear search. To search for the minimum of the non-convex residual functional,
a specialized global optimization algorithm has been developed, based on the idea
of curvilinear search on the reachable set of a control system produced using a
pair of randomly generated “support” controls. Using the property of connectedness,
variations in the space of control functions that do not violate the existing constraints
are constructed at iterations of the algorithm. This is achieved by direct projection on a
parallelepiped. To improve the current approximation, a modification of a non-local one-
dimensional search algorithm based on the “parabolic” method is applied [28]. Also,
as a globalizing mechanism, a non-local search in random directions is used, repeated
many times at each iteration of the algorithm. To solve auxiliary non-convex problems
of one-dimensional search, a modification of the stochastic P -algorithm, proposed by
A. Zhigljavsky and A. Žilinskas [14], is implemented [29]. To enhance reliability of
the proposed method, a periodic random multi-start was provided in the algorithm’s
construction.

• An algorithm for the numerical integration of FDEPT systems based on the
sequential discretization technique has been developed. To improve the efficiency of
calculations, tools have been implemented to build a sequence of approximative finite-
dimensional optimization problems with a growing number of variables. At the same
time, solutions obtained at the previous stages of calculations performed on the current
discretization grid are projected onto a new grid with an increased number of nodes
while preserving the qualitative and quantitative characteristics of the trajectories.

• The proposed numerical algorithms were tested on a wide range of tasks [27] with
using the principle of “the best of known solutions” [30]. In all considered problems,
the proposed algorithm allowed us to find the best known solution. The calculation
experiments that have been carried out demonstrate considerably high fidelity of the
proposed algorithms.

The above heuristic search algorithm for the solution x̂(t) can be justified on the
basis of the existence and uniqueness theorems for initial-boundary value problems for the
investigated FDEPT, as well as theorems on approximating solutions of such equations on
the whole line by solutions of the initial-boundary value problem on a sequence of expanding
intervals. A description of such equations and corresponding results are presented in the
following sections.

4. EXISTENCE AND UNIQUENESS THEOREM FOR SOLITON SOLUTIONS
IN THE PROBLEM OF LONGITUDINAL VIBRATIONS OF AN INFINITE
HOMOGENEOUS ROD

We consider a problem from the theory of plastic deformation, in which solutions of the
traveling wave type for the following system are studied

mÿi(t) = yi+1(t)− 2yi(t) + yi−1(t) + Φ(yi(t)), i ∈ Z, yi ∈ R, t ∈ R, (4.12)
yi(t+ τ) = yi+1(t), τ ≥ 0. (4.13)

We will study such a system under the most weak conditions on the potential Φ(·) in the form
of the Lipschitz condition. The Lipschitz constant will be denoted by LΦ. We formulate a
theorem on the existence and uniqueness of a traveling wave type solution (soliton solution),

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



NUMERICAL METHODS FOR FUNCTIONAL DIFFERENTIAL EQUATIONS 63

τ

µ

1

µ̂

τ̂

µ1(τ )

µ2(τ )

0

Fig. 4.1. Graphs of functions µ1(τ), µ2(τ)

based on the correspondence of traveling-wave type solutions of the infinite-dimensional
system (4.12)-(4.13) to the solutions of the induced FDEPT [24], reduced to a first-order
system

ẋ1(t) = x2(t), x ∈ R2, x = (x1, x2)′, t ∈ R,
ẋ2(t) = m−1

(
x1(t+ τ)− 2x1(t) + x1(t− τ) + Φ(x1(t))

)
, (4.14)

as well as the existence and uniqueness theorem for solutions to such equations
(Theorem 2.1).

To do this, we consider a transcendental equation with respect to two variables τ ∈
[0,+∞) and µ ∈ (0, 1)

Cτ
(
2µ−1 + 1

)
= lnµ−1, (4.15)

where
C = max{1; 2m−1

√
L2

Φ + 2}.

The set of solutions of the equation (4.15) is described by functions µ1(τ), µ2(τ) given in
Fig. 4.1.

If we denote the right-hand side of the system (4.14) by f , then relations µ1(Lf ;h) =
τ
√
µ1(τ) and µ2(Lf ;h) = τ

√
µ2(τ) will be valid, where Lf = C, h = (τ, τ).

For each µ ∈ (0, 1) we define the phase space of the infinite-dimensional equation
(4.12)as a Hilbert space of sequences with the corresponding norm

K2
Z2µ = {κ : κ = {zi}+∞

−∞, zi ∈ R2, i ∈ Z,
∑
i∈Z

‖zi‖2
R2µ2|i| < +∞}.

Theorem 4.1 ( [24]):
Let the potential Φ satisfy the Lipschitz condition with constant LΦ. Then, for any initial
values ī ∈ Z, a, b ∈ R, t̄ ∈ R, and characteristics τ > 0 satisfying the condition

0 < τ < τ̂ ,

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



64 A. BEKLARYAN

for the initial system of differential equations (4.12) there exists a unique solution of the
traveling wave type {yi(·)}+∞

−∞ with characteristic τ such that it satisfies the initial conditions
yī(t̄) = a, ẏī(t̄) = b. For any parameter µ ∈ (µ1(τ), µ2(τ)) the vector function

ω(t) = {(yi(t), ẏi(t))
′}+∞
−∞

belongs to the space K2
Z2µ for any t ∈ R, and the function

ρ(t) = ‖ω(t)‖2µ

belongs to the space L1
τ
√
µC

(1)(R). Such a solution depends continuously on the initial values
a, b ∈ R, as well as on the mass m.

Theorem 4.1 not only guarantees the existence of a solution but also determines the
limitation of its possible growth both in time t and in coordinates i ∈ Z (over space). It is
obvious that for each 0 < τ < τ̂ the space K2

Z2(µ2(τ)−ε), for small ε > 0, is much narrower
than the space K2

Z2(µ1(τ)+ε). The theorem guarantees the existence of a solution in narrower
spaces and uniqueness in wider spaces.

5. APPROXIMATION OF THE SOLUTION OF A FUNCTIONAL DIFFERENTIAL
EQUATION DEFINED ON THE LINE BY SOLUTIONS OF AN INITIAL-
BOUNDARY VALUE PROBLEM WITH EXPANDING INTERVALS OF THE
DEFINITION

Next, we consider soliton solutions of the system, which will be implemented as solutions
of the induced functional differential equations defined on the whole line. For the numerical
integration of such an equation, one should be able to approximate the solutions of such an
equation by the solutions of initial-boundary value problems defined on an expanding family
of finite intervals. This section is devoted to this problem.

Let’s consider the initial boundary value problem

ẋ(t) = f(t, x(t+ τ1), . . . , x(t+ τs)), t ∈ BR, (5.16)
ẋ(t) = ϕ(t), t ∈ R\BR, ϕ(·) ∈ L∞(R,Rn), (5.17)
x(t̄) = x̄, t̄ ∈ R, x̄ ∈ Rn, (5.18)

where τj ∈ R. We formulate again the theorem on the continuous dependence on the initial-
boundary conditions. We define a Banach space of functions with weight

LnµL∞(R) =

{
ϕ(·) : ϕ(·) ∈ L∞ (R,Rn) , sup

t∈R
vrai‖ϕ(t)µ|t|‖Rn < +∞

}
, µ ∈ (0, 1)

and norm
‖ϕ(·)‖µ = sup

t∈R
vrai‖ϕ(t)µ|t|‖Rn .

If in the boundary condition (5.17), starting with some sufficiently large k ∈ Z+ with
the condition BR ⊂ [−k, k], the boundary condition ϕ(·) ∈ L∞(R) replace outside the
interval BR ⊂ [−k, k] so that the new boundary function ϕ̃(·) satisfies the condition ϕ̃(·) ∈
LnµL∞(R), then the corresponding solution x̃(·) of the boundary value problem on BR ⊂
[−k, k] will match the original solution x(·).

We are going to formulate a proposition on the approximation of solutions of an initial-
boundary value problem defined on the whole line by solutions of the initial-boundary value

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



NUMERICAL METHODS FOR FUNCTIONAL DIFFERENTIAL EQUATIONS 65

problem defined on the interval [−k, k] as k → +∞. We consider the initial-boundary value
problem on the whole line BR = R

ẋ(t) = f(t, x(t+ τ1), . . . , x(t+ τs)), t ∈ R, (5.19)
x(t̄) = x̄, t̄ ∈ R, x̄ ∈ Rn, (5.20)

where τj ∈ R, j = 1, . . . , s, j = 1, . . . , s. For each k ∈ Z we consider the initial-boundary
value problem on a finite interval BR = [−k, k]

ẋ(t) = f(t, x(t+ τ1), . . . , x(t+ τs)), t ∈ [−k, k], (5.21)
ẋ(t) = ϕ(t), t ∈ R\[−k, k], ϕ(·) ∈ Ln1L∞(R), (5.22)
x(t̄) = x̄, t̄ ∈ R, x̄ ∈ Rn. (5.23)

Theorem 5.1 ( [31]):
Let the map f(·) satisfies the conditions (a)− (d). If for µ ∈ (0, µ?) ∩ (0, 1) the inequality

Lf

s∑
j=1

µ−|τj | < lnµ−1 (5.24)

holds, and (µ1(Lf ;h), µ2(Lf ;h)) is the maximum solution interval for the inequality (5.24),
then for any x̄ ∈ Rn, ϕ(·) ∈ Ln1L∞(R), the solution x̂(·) of the initial-boundary value
problem (5.19)-(5.20), as an element of the space LnµC(0)(R), is approximated by solutions
x̂k(·) of the initial-boundary value problem (5.21)–(5.23) as k → +∞. Moreover, for any
arbitrarily small ε, 0 < ε < µ2 − µ1 there is Cfϕε such that the following estimate is true

‖x̂(·)− x̂k(·)‖(0)
µ ≤ Cfϕε

(
µ1(Lf ;h)

µ2(Lf ;h)− ε

)k
.

From the estimate it follows that the convergence rate will be the higher, the more the
µ1(Lf ;h) and µ2(Lf ;h) values will differ from each other, the difference between them is
inversely proportional to the value of Lf max1≤j≤s |τj|.

6. CONSTRUCTION OF THE COMPLETE FAMILY OF TRAVELING WAVE TYPE
SOLUTIONS IN THE FRENKEL-KONTOROVA MODEL

Considering periodic functionals for the model from the section 4, we obtain the Frenkel-
Kontorova model from the theory of plastic deformation. In this case, taking the functional
Φ(y) = A sin(By), A,B ∈ R

mÿi(t) = yi+1(t)− 2yi(t) + yi−1(t) + Asin(Byi(t)), i ∈ Z, yi ∈ R, t ∈ R, (6.25)
yi(t+ τ) = yi+1(t), τ > 0, (6.26)

we construct solutions of the traveling wave type with characteristic τ . As noted earlier, the
space of traveling wave type solutions for such a system coincides with the space of solutions
of the FDEPT system

ẋ1(t) = τx2(t), t ∈ R,
ẋ2(t) = τm−1[x1(t+ 1)− 2x1(t) + x1(t− 1) + A sin(Bx1(t))]. (6.27)

Moreover, there is a correspondence of solutions according to the rule

yi(t) = x(τ−1t+ i).

We will construct numerical solutions of the system (6.27) using the approximation Theorem
5.1.

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



66 A. BEKLARYAN

6.1. Numerical experiments
Next, the results of the computational experiments on the study of initial-boundary value
problems for systems of FDEPT using OPTCON-F software will be presented. Before
demonstrating the examples, it is necessary to make a number of significant observations:

(a) in the SC OPTCON-F, the possibility of satisfying the condition of uniform
boundedness of the derivative is realized. The maximum deviation from zero is
determined by the Lipschitz constant of the right side of the equation, by the parameter
µ, and also by the deviations of the argument;

(b) numerical integration on an interval with initial-boundary conditions is realized as an
integration process with given boundary conditions at the left end and procedures for
minimizing deviations from the boundary conditions at the right end for the solution
constructed with observance of the restrictions from the preceding item;

(c) in the obtained theorem on approximating solutions of the original equation on the
whole line by solutions on expanding finite intervals [−k, k], the only restriction on
the boundary conditions themselves is the condition for their uniform boundedness. In
particular, we choose the zero boundary condition;

(d) we note that the obtained condition for the existence of a solution of the traveling wave
type is just a sufficient condition. Therefore, solutions of the traveling wave type can
also be numerically constructed for τ > τ̂ . Nevertheless, many of the central conditions
of the presented theory are “exact”, that is, there are examples of equations for which
violation of the indicated conditions leads to the lack of solutions;

(e) in the SC OPTCON-F there is the possibility of sequential application of various
algorithms within the framework of constructing a solution for one task. Thus, the
constructed intermediate solution in the previous step becomes the starting solution
(“baseline”) for the following algorithm. In this case, such an implementation does
not prevent the global search algorithms from “popping out” of the local solution.
Separately, we note the presence of a programming module that allows predetermining
the order of application of algorithms, as well as the construction of complex chain of
steps (conditional statements, loops, etc.) depending on the current or historical values
of a number parameters (for example, error estimation or number of iterations). For
the examples presented below, in addition to the basic global optimization algorithm
described in Section 1, the following scheme was used in cycle: the generalized quasi-
Newtonian and Powell-Brent’s methods (with bi-directional line search along each
dimension) were used alternately, and after the error changed by less than 10−p the
adaptive modification of the Hooke-Jeeves method was used l times (p and l are
computable functions on the basis of the loop iteration number, as well as the Lipschitz
constants of the equation itself and a number of other technical characteristics). The
stopping criterion depended on the number of iterations in the first part of the cycle, as
well as the current error estimate and its dynamics;

(f) in view of all the points listed above, as well as stochastic elements in the applied
algorithms, the presented value of the residual functional (RF, i.e. error of a numerical
solution) can not be used to estimate the theoretical rate of convergence.

We consider dynamical system in the following form:{
ẋ1(t) = 0.15x2(t),
ẋ2(t) = 0.15× 100−1(x1(t+ 1)− 2x1(t) + x1(t− 1) + 500 sin(0.1x1(t))),

t ∈ R,
initial conditions{

x1(0) = c,
x2(0) = 0.

(6.28)

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



NUMERICAL METHODS FOR FUNCTIONAL DIFFERENTIAL EQUATIONS 67

(a) Trajectories of x1(t) for the system (6.29) (b) Trajectories of x1(t) for the system (6.28)

Fig. 6.2. Trajectories of x1(t) at different c

Here, with respect to the system (6.27), we have A = 500, B = 0.1, τ = 0.15,m = 100, and
the equality (4.15) takes the form

0.003
√

2502(2µ−1 + 1) = lnµ−1

and has on the interval (0, 1) two solutions with approximate values µ1(0.15) = 0, 22 and
µ2(0.15) = 0, 424191 (the exact values are expressed in terms of the Lambert W -function
and can’t be written out in quadratures).

Taking into account the impossibility of considering the numerical solution of the system
on an infinite interval, we introduce the parameter k and the corresponding family of
expanding initial-boundary value problems{

ẋ1(t) = 0.15x2(t),
ẋ2(t) = 0.15× 100−1(x1(t+ 1)− 2x1(t) + x1(t− 1) + 500 sin(0.1x1(t))),

t ∈ [−k, k],
boundary conditions{

ẋ1(t) = 0,
ẋ2(t) = 0,

t ∈ (−∞,−k] ∪ [k,+∞),

initial conditions{
x1(0) = c,
x2(0) = 0.

(6.29)

According to the Theorem 5.1, convergence is guaranteed under any given essentially
bounded boundary conditions. In the system (6.29), the simplest boundary function is chosen,
the identity zero. Thus, the solution of the system (6.29) converges (according to the metric
of the space LnµC(0)(R) for µ ∈ (µ1(0.15), µ2(0.15))) to the solution of the system (6.28) as
k →∞.

Since the equation (6.27) is autonomous, the solution space of such equation is invariant
with respect to time-variable shifts. On the other hand, from the periodicity of the right-hand
side with respect to the x1(t) it follows that the solution space of such equation is invariant

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



68 A. BEKLARYAN

(a) Trajectories of x1(t) at different k

(b) Trajectories of x1(t) at different c

Fig. 6.3. Unbounded trajectories of x1(t)

with respect to a shift in x1(t) for a period equal to 2π
B

. Therefore, it suffices to consider a
family of solutions of the initial system (6.28) with a value of initial conditions x1(0) from
zero to the value of the period. Nevertheless, the stationary solutions x1(t) are repeated each
half-period π

B
. Since the right-hand side of the equation is an odd function with respect to the

x1(t), the solution space of such equation can withstand the reflection transformation with
respect to the axis t. Hence, it is sufficient to construct trajectories x1(t) in the strip from zero
to the half-period. Figure 6.2 shows the graphs x1(t) for different values of the parameter
c = x1(0) for both the system (6.29) and the initial system (6.28) (the values of c are reduced
to half-period). Note that Fig. 6.2 shows the complete family of bounded graphs x1(t) (up to
the above transformations) from the space LnµC(0)(R) for µ ∈ (µ1(0.15), µ2(0.15)).

At the same time, for some initial conditions, there are only unbounded x1(t). Figure 6.3a
demonstrates the evolution of such an unbounded graph x1(t) with increasing k for the initial
conditions x1(0) = π

B
, x2(0) = 50. Figure 6.3b, in turn, shows a one-parameter family of

unbounded graphs x1(t) for x1(0) = π
B

and different values of x2(0) = c.

7. CONCLUSION

The construction of numerical solutions of the traveling wave type for the Frenkel-Kontorova
model on the propagation of longitudinal waves using the developed software package was
demonstrated. For this model, a complete family of traveling wave solutions was constructed.

ACKNOWLEDGEMENTS
The reported study was funded by RFBR according to the research project 19-01-00147.

REFERENCES

1. Bellman, R. & Cooke, K. L. (1963) Differential-Difference Equations. New York:
Academic Press.

2. Hale, J. K. (1977) Theory of Functional Differential Equations. New York: Springer.

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)



NUMERICAL METHODS FOR FUNCTIONAL DIFFERENTIAL EQUATIONS 69

3. Myshkis, A. D. (1972) Linear Differential Equations with Retarded Arguments.
Moscow: Nauka, (in Russian).

4. Bellen, A. & Zennaro, M. (2003) Numerical Methods for Delay Differential Equations.
USA: Oxford University Press.

5. Maset, S. (2003) Numerical solution of retarded functional differential equations as
abstract cauchy problems. Journal of Computational and Applied Mathematics, 161,
259–282, https://dx.doi.org/10.1016/j.cam.2003.03.001.

6. El’sgol’ts, L. E. & Norkin, S. B. (1973) Introduction to the Theory and Application of
Differential Equations with Deviating Arguments. New York: Academic Press.

7. Baotong, C. (1995) Functional differential equations mixed type in banach spaces.
Rendiconti del Seminario Matematico della Università di Padova, 94, 47–54.

8. Ford, N. J. & Lumb, P. M. (2009) Mixed-type functional differential equations: A
numerical approach. Journal of Computational and Applied Mathematics, 229, 471–
479, https://dx.doi.org/10.1016/j.cam.2008.04.016.

9. Abell, K. A., Elmer, C. E., Humphries, A. R., & Van Vleck, E. S. (2005) Computation of
mixed type functional differential boundary value problems. SIAM Journal on Applied
Dynamical Systems, 4, 755–781, https://dx.doi.org/10.1137/040603425.

10. Lima, P. M., Teodoro, M. F., Ford, N. J., & Lumb, P. M. (2010) Ana-
lytical and numerical investigation of mixed-type functional differential equa-
tions. Journal of Computational and Applied Mathematics, 234, 2826–2837,
https://dx.doi.org/10.1016/j.cam.2010.01.028.

11. Lima, P. M., Teodoro, M. F., Ford, N. J., & Lumb, P. M. (2010) Finite element solution
of a linear mixed-type functional differential equation. Numerical Algorithms, 55, 301–
320, https://dx.doi.org/10.1007/s11075-010-9412-y.

12. Li, X. Y. & Wu, B. Y. (2014) A continuous method for nonlocal functional differential
equations with delayed or advanced arguments. Journal of Mathematical Analysis and
Applications, 409, 485–493, https://dx.doi.org/10.1016/j.jmaa.2013.07.039.

13. Evtushenko, Y., Posypkin, M., & Sigal, I. (2009) A framework for parallel large-scale
global optimization. Computer Science – Research and Development, 23, 211–215,
https://dx.doi.org/10.1007/s00450-009-0083-7.

14. Zhigljavsky, A. & Žilinskas, A. (2008) Stochastic Global Optimization. US: Springer,
https://dx.doi.org/10.1007/978-0-387-74740-8.

15. Bellman, R. (1957) Dynamic Programming. Princeton, New Jersey: Princeton
University Press.

16. Krotov, V. F. (1996) Global Methods in Optimal Control Theory. New York: Marcel
Dekker Inc.

17. Chachuat, B., Singer, A. B., & Barton, P. I. (2006) Global methods for dynamic
optimization and mixed-integer dynamic optimization. Industrial & Engineering
Chemistry Research, 45, 8373–8392, https://dx.doi.org/10.1021/ie0601605.

18. Lopez Cruz, I. L., Van Willigenburg, L. G., & Van Straten, G. (2003) Efficient
differential evolution algorithms for multimodal optimal control problem. Applied Soft
Computing, 3, 97–122, https://dx.doi.org/10.1016/S1568-4946(03)00007-3.

19. Floudas, C. A. & Gounaris, C. E. (2009) A review of recent advances
in global optimization. Journal of Global Optimization, 45, 3–38,
https://dx.doi.org/10.1007/s10898-008-9332-8.

20. Khrustalev, M. M. (1988) An accurate description of feasibility sets and conditions
for the global optimality of dynamic systems. ii. conditions for global optimality.
Automation and Remote Control, 49, 874–881.

21. Tolstonogov, A. (2000) Differential inclusions in a Banach space. Netherlands:
Springer, https://dx.doi.org/10.1007/978-94-015-9490-5.

22. Subbotin, A. I. & Chentsov, A. G. (1981) Guarantee Optimization in Control Problems.
Moscow: Nauka, (in Russian).

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)

https://dx.doi.org/10.1016/j.cam.2003.03.001
https://dx.doi.org/10.1016/j.cam.2008.04.016
https://dx.doi.org/10.1137/040603425
https://dx.doi.org/10.1016/j.cam.2010.01.028
https://dx.doi.org/10.1007/s11075-010-9412-y
https://dx.doi.org/10.1016/j.jmaa.2013.07.039
https://dx.doi.org/10.1007/s00450-009-0083-7
https://dx.doi.org/10.1007/978-0-387-74740-8
https://dx.doi.org/10.1021/ie0601605
https://dx.doi.org/10.1016/S1568-4946(03)00007-3
https://dx.doi.org/10.1007/s10898-008-9332-8
https://dx.doi.org/10.1007/978-94-015-9490-5


70 A. BEKLARYAN

23. Zarodnyuk, T. S., Anikin, A. S., Finkelshtein, E. A., Beklaryan, A. L., & Belousov,
F. A. (2016) The technology for solving the boundary value problems for nonlinear
systems of functional differential equations of pointwise type. Modern technologies.
System analysis. Modeling, 49, 19–26, (in Russian).

24. Beklaryan, L. A. (2007) Introduction to the theory of functional differential equations.
Group approach. Moscow: Factorial Press, (in Russian).

25. Beklaryan, L. A. (1991) A method for the regularization of boundary value problems
for differential equations with deviating argument. Soviet Math. Dokl., 43, 567–571.

26. Zarodnyuk, T. S., Gornov, A. Y., Anikin, A. S., & Finkelstein, E. A. (2017)
Computational technique for investigating boundary value problems for functional-
differential equations of pointwise type. In Proc. of the VIII International Conference
on Optimization and Applications (OPTIMA-2017), Petrovac, Montenegro, October 2-
7, CEUR-WS.org.

27. Gornov, A. Y., Zarodnyuk, T. S., Madzhara, T. I., Daneeva, A. V., & Veyalko,
I. A. (2013) A Collection of Test Multiextremal Optimal Control Problems, 257–274.
Springer.

28. Gornov, A. Y., Zarodnyuk, T. S., Finkelstein, E. A., & Anikin, A. S. (2016)
The method of uniform monotonous approximation of the reachable set bor-
der for a controllable system. Journal of Global Optimization, 66, 53–64,
https://dx.doi.org/10.1007/s10898-015-0346-8.

29. Gornov, A. Y. & Finkelstein, E. A. (2015) Algorithm for piecewise-linear approximation
of the reachable set boundary. Automation and Remote Control, 76, 385–393,
https://dx.doi.org/10.1134/S0005117915030030.

30. Floudas, C. A. & Pardalos, P. M. (1990) A Collection of Test Problems for
Constrained Global Optimization Algorithms. Berlin, Heidelberg: Springer-Verlag,
https://dx.doi.org/10.1007/3-540-53032-0.

31. Beklaryan, L. A. & Beklaryan, A. L. (2017) Traveling waves and functional differential
equations of pointwise type. what is common? In Proc. of the VIII International
Conference on Optimization and Applications (OPTIMA-2017), Petrovac, Montenegro,
October 2-7, CEUR-WS.org.

Copyright © 2020 ASSA. Adv Syst Sci Appl (2020)

https://dx.doi.org/10.1007/s10898-015-0346-8
https://dx.doi.org/10.1134/S0005117915030030
https://dx.doi.org/10.1007/3-540-53032-0

	Introduction
	Statement of the initial boundary-value problem
	Existence and Uniqueness Theorem for the Basic Initial-Boundary Value Problem

	Variational Approach to Constructing Solutions to the Basic Initial-Boundary Value Problem
	Software Package OPTCON-F

	Existence and Uniqueness Theorem for Soliton Solutions in the Problem of Longitudinal Vibrations of an Infinite Homogeneous Rod
	Approximation of the Solution of a Functional Differential Equation Defined on the Line by Solutions of an Initial-Boundary Value Problem With Expanding Intervals of the Definition
	Construction of the complete family of traveling wave type solutions in the Frenkel-Kontorova model
	Numerical experiments

	Conclusion