<figure><img src="image_1.png" /></figure>

HAL
open science

# Difference methods for differential inclusions: a survey

Asen Dontchev, Frank Lempio

► To cite this version:

Asen Dontchev, Frank Lempio. Difference methods for differential inclusions: a survey. SIAM Review, 1992, 34 (2), pp.263-294. 10.1137/1034050. hal-01571241

HAL Id: hal-01571241

https://hal.science/hal-01571241

Submitted on 1 Aug 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# DIFFERENCE METHODS FOR DIFFERENTIAL INCLUSIONS: A SURVEY*

ASEN DONTCHEV† AND FRANK LEMPIO‡

**Abstract.** The main objective of this survey is to study convergence properties of difference methods applied to differential inclusions. It presents, in a unified way, a number of results scattered in the literature and provides also an introduction to the topic.

Convergence proofs for the classical Euler method and for a class of multistep methods are outlined. It is shown how numerical methods for stiff differential equations can be adapted to differential inclusions with additional monotonicity properties. Together with suitable localization procedures, this approach results in higher-order methods.

Convergence properties of difference methods with selection strategies are investigated, especially strategies forcing convergence to solutions with additional smoothness properties.

The error of the Euler method, represented by the Hausdorff distance between the set of approximate solutions and the set of exact solutions is estimated. First- and second-order approximations to the reachable sets are presented.

**Key words.** differential inclusions, difference methods

**1. Introduction.** In this survey we consider the following initial value problem for ordinary differential inclusions.

**INITIAL VALUE PROBLEM 1.1.** Let $I = [t_0, T]$ be a finite interval, $y_0 \in \mathbb{R}^n$, and $F$ be a map from $I \times \mathbb{R}^n$ into the set of all subsets of $\mathbb{R}^n$.

Find an absolutely continuous function $y(\cdot)$ on $I$ such that

$$ (1) \qquad y(t_0) = y_0 \quad \text{and} \quad \dot{y}(t) \in F(t, y(t)) \quad \text{for almost all } t \in I, $$

where $\dot{y}(\cdot)$ is the derivative of $y(\cdot)$.

Just to give some motivation for studying differential inclusions, we briefly mention several applications where differential inclusions naturally occur.

A first motivation originates from differential equations with single-valued, discontinuous right-hand sides

$$ \dot{y}(t) = f(t, y(t)). $$

To get a sound notion of a solution, following [24], this problem has to be restated in the form

$$ \dot{y}(t) \in \bigcap_{\delta>0} \bigcap_{\mu(N)=0} \text{cl}(\text{conv}(f(t, \{z \in \mathbb{R}^n : \|z - y(t)\| \le \delta\} \setminus N))) $$

for almost all $t \in I$, where $\mu$ denotes Lebesgue measure on $\mathbb{R}^n$. Hence we arrive at a differential inclusion. Differential inclusions of this type occur in a variety of applications, e.g., in oscillating systems with combined dry and viscous damping [55], [2], [67], [51], especially in vehicle dynamics for the description of locking phases during brake

<footnote>This research was supported by the DFG (Deutsche Forschungsgemeinschaft) and by the Institute of Mathematics at the Bulgarian Academy of Sciences.</footnote>

<footnote>†Mathematical Reviews, 416 Fourth Street, Ann Arbor, Michigan 48107, on leave from the Institute of Mathematics, Bulgarian Academy of Sciences, Sofia, Bulgaria.</footnote>

<footnote>‡Lehrstuhl für Angewandte Mathematik, Universität Bayreuth, Postfach 10 12 51, D-8580 Bayreuth, Federal Republic of Germany.</footnote>
maneuvers [62], in elasto-plasticity [8], [15], variable structure systems [70], and electrical circuits [54], [19]. Moreover, the method of lines for nonlinear evolution equations with discontinuities, respectively, for evolutionary variational inequalities results in whole families of differential equations with discontinuous right-hand sides, cf., e.g., [56], [14], [6]. In §§3 and 4 we use Example 3.3, cf. [55], [67], describing forced vibrations with viscous and dry damping, as model problem for numerical tests.

A second motivation is given by differential inclusions of the following type

$$ \dot{y}(t) \in -\partial\varphi(y(t)) $$

where $\varphi: \mathbb{R}^n \to \mathbb{R}$ is a convex potential function with subdifferential $\partial\varphi$. Such inclusions have an important property: their equilibrium solutions minimize $\varphi$. Moreover, if $\varphi$ achieves a minimum at all, then for every initial value $y_0 \in \mathbb{R}^n$, the corresponding solution $y(t)$ as $t \to \infty$ converges to a minimizer of $\varphi$. Hence, there is an interesting connection between differential inclusions and subgradient methods for convex optimization problems, cf. [3], [4], [52].

A third motivation, naturally, is given by optimal control problems. Disregarding for the moment any objective function and the special structure of controls, the differential inclusions (1) could be obtained from a control system

$$ \dot{y}(t) = f(t, y(t), u(t)) $$

with feasible controls

$$ u(t) \in U(t, y(t)) $$

for almost all $t \in I$, where $f$ is single-valued and $U(t, x) \subset \mathbb{R}^m$ for all $t \in I$ and $x \in \mathbb{R}^n$, just by taking

$$ F(t, x) = \{f(t, x, z) : z \in U(t, x)\}. $$

In §5, Example 5.3, we describe a special differential inclusion of this type, which was originally used in [34] for the numerical test of several selection strategies. Such a reduction of an optimal control problem to a differential inclusion is especially appropriate when we consider Mayer's problem, which can be regarded as minimizing a given objective function on the reachable set at time $T$ of initial value problem (1). Furthermore, taking the full structure of control functions and even more general objective functions into account, the necessary optimality conditions can be analysed in terms of boundary value problems for differential inclusions, cf. [17] and [18], Chapter 3, where a nice example of a nonsmooth problem in resource economics is presented, and [59], [60], where simplicial fixed point algorithms for set-valued operators are investigated and used for the computation of optimal fishing strategies. Control systems with unknown but bounded disturbances can be described by differential inclusions; this observation is used in [36] for control synthesis of uncertain systems.

A number of other motivations and applications of differential inclusions can be found in [7] and in the books [4], [16], [18], [26], [27], [53], and their references.

The main objective of this survey is the study of difference methods for differential inclusions, which are motivated by difference methods for differential equations with single-valued right-hand sides.

Let $X$ be the set of solutions to (1). As a rule, the set $X$ consists of more than one element, that is we have a bundle of trajectories. Consequently, there are various closely connected approaches of approximating solutions $y(\cdot) \in X$.
The first approach uses a finite difference scheme together with suitable selection procedures resulting in a sequence of grid functions

$$ (\eta_0^N, \eta_1^N, \dots, \eta_N^N) \quad (N \in \mathbb{N}'), $$

say, on a uniform grid

$$ t_0^N < t_1^N < \dots < t_N^N = T $$

with stepsize

$$ h = \frac{T - t_0}{N} = t_j^N - t_{j-1}^N \quad (j = 1, \dots, N), $$

where, as usual, $\mathbb{N}'$ denotes a subsequence of $\mathbb{N}$ converging to infinity.

Naturally, the question arises whether at least a subsequence of the sequence

$$ (\eta^N(\cdot))_{N \in \mathbb{N}'} $$

of, say, piecewise linear continuous interpolants of the grid functions converges to a solution $y(\cdot) \in X$ of (1). For linear multistep methods this question is investigated in §3, following the results of K. Taubert [66], [68].

Next, the closely connected question arises of how fast this subsequence converges, i.e., which order of convergence could be attained by a special sequence of difference approximations. This question is addressed to in §4. Stimulated by a paper of C. M. Elliott [23] on first-order convergence for a special class of methods A. Kastner-Maresch succeeded in adapting convergence proofs for numerical methods for stiff differential equations to differential inclusions satisfying a uniform one-sided Lipschitz condition; compare [32], [33]. This results in higher-order convergence on suitable subintervals of $I$, together with appropriate localization procedures in higher-order methods. Since the one-sided Lipschitz condition implies uniqueness of the solution of (1), only special classes of problems can be treated in this way, e.g., differential inclusions with additional monotonicity properties. At least for differential equations with single-valued discontinuous right-hand sides transformation algorithms like those proposed by D. Stewart [64], [65] lead to higher-order convergence under special assumptions.

Moreover, the question is interesting, of whether the limit function $y(\cdot) \in X$ has additional desirable properties. This question is treated in §5 for general differential inclusions, where convergence properties of difference methods with selection strategies are investigated forcing convergence to solutions with additional smoothness properties. We present algorithms and discuss the convergence of selections with minimal norm, minimal variation, and with respect to a given reference trajectory.

The second approach consists in approximating the whole solution set $X$ of (1). It exploits the fact that a difference scheme virtually describes a difference inclusion for each stepsize $h$ if no special selection procedure is used. As a result, for each $h$ we obtain a set $X_h$ of, say, piecewise linear continuous functions, approximating $X$. Apparently, as a definition of convergence of the method we can use an appropriate concept for convergence of sets $X_h$ to $X$ as $h \to 0$. The error can be measured by a suitably defined distance between the sets $X_h$ and $X$.

As an introductory example, in §2 we consider the classical Euler method. Following the first approach we prove directly that if $F$ is compact convex valued upper semicontinuous and with linear growth, then

$$ (2) \qquad \limsup_{h \to 0} X_h \subset X \quad \text{in } C(I)^n, $$
where $\limsup$ is in the usual (Kuratowski) sense, i.e.,

$$ \limsup_{h \to 0} X_h = \left\{ f \in C(I)^n : \liminf_{h \to 0} \mathrm{dist}(f, X_h) = 0 \right\}, $$

and $C(I)^n$ is endowed with the supremum norm.

If, moreover, $F$ is continuous and integrably Lipschitz in $x$ on bounded sets, then a partial result in the spirit of the second approach can be proved: every solution $\eta^N(\cdot) \in X_h$ of the difference inclusion with sufficiently small stepsize contains in its $\epsilon$-neighbourhood in $C(I)^n$ a solution $y(\cdot) \in X$ of the differential inclusion (1).

We consider again Euler method in §6 and show that if $F$ is Lipschitz in both $t$ and $x$, then we have

$$ \mathrm{haus}(X_h, X) = \mathcal{O}(h) \quad \text{in } C(I)^n, $$

where $\mathrm{haus}(\cdot, \cdot)$ denotes Hausdorff distance.

The set of solutions $X$ may be approximated not only as a set of functions. Sometimes it may be important to describe the values of all solutions at certain points. Thus we arrive at the problem of approximating the reachable sets.

Let $t$ be a point in $I$. The “reachable set at the time $t$” of (1) is defined as

$$ R(t) = \{x \in \mathbb{R}^n : x = y(t) \text{ for some } y(\cdot) \in X\}. $$

That is, $R(t)$ is the set of all points $x$ that are ends of trajectories of (1) on $[t_0, t]$. Section 7 shows that a sequence of sets can be determined from the Euler scheme, which is Hausdorff convergent to $R(t)$ as $h \to 0$, uniformly in $t \in I$. This follows from an older result due to A. I. Panasyuk and V. I. Panasyuk [50], concerning the so-called funnel equation. Error estimates obtained by M. S. Nikol'skiy [46]–[48] are discussed as well.

Section 8 presents some recent results due to V. M. Veliov, who found second-order approximations related to Runge-Kutta schemes, for both the trajectory bundle and the reachable sets.

All computer tests were made on the VAX cluster of the Computer Center of the University of Bayreuth, consisting of a DEC VAX 8600 and 6310. In this survey the numerical results are visualized by computer plots. For a more detailed presentation of these results compare [34] and [32].

**2. Euler method.** By far the simplest difference method for solving Initial Value Problem 1.1 is the classical Euler method which we present in the following as an introductory example.

**EULER METHOD 2.1.** For $N \in \mathbb{N}' \subset \mathbb{N}$ choose a grid

$$ t_0 < t_1 < \dots < t_N = T $$

with stepsize

$$ h = \frac{T - t_0}{N} = t_j - t_{j-1} \quad (j = 1, \dots, N). $$

Let $\eta_0 = y_0$, and for $j = 0, \dots, N-1$ compute $\eta_{j+1}$ from

$$ (3) \qquad \eta_{j+1} \in \eta_j + h F(t_j, \eta_j). $$
Where necessary, we use the exponent $N$ to emphasize the dependence on $N$. Only to avoid technical difficulties, we work with an equidistant grid.

As a solution of the difference inclusion (3) for a given stepsize $h$, it is convenient to consider any continuous and piecewise linear function

$$ \eta^N(\cdot) : I \to \mathbb{R}^n $$

such that

$$ \eta^N(t) = \eta_j + \frac{1}{h}(t - t_j)(\eta_{j+1} - \eta_j) \quad (t_j \le t \le t_{j+1}, \quad j = 0, \dots, N-1), $$

where $(\eta_0, \dots, \eta_N)$ is any grid function satisfying (3). Let $X_h$ be the set of all solutions of (3) for given $h$.

The following theorem is present explicitly or implicitly in many works, in various forms, see, e.g., J.-P. Aubin and A. Cellina [4, Lemma 1, p. 99], F. H. Clarke [18, Thm. 3.1.7, p. 118], or A. F. Filippov [26, Thm. 1, p. 77]. Its proof uses the idea of the classical Peano theorem to prove existence of solutions to differential equations.

**THEOREM 2.2.** Suppose that the set-valued map $F$ satisfies the conditions:

(i) $F$ is nonempty compact and convex valued upper semicontinuous in $I \times \mathbb{R}^n$.

(ii) There exist constants $k$ and $a$, such that

$$ \|z\| \le k\|x\| + a $$

whenever $z \in F(t, x)$, $x \in \mathbb{R}^n$, $t \in I$.

Then every sequence $(\eta^N(\cdot))_{N \in \mathbb{N}'}$ with $\eta^N(\cdot) \in X_h$ for $N \in \mathbb{N}'$ has a subsequence which converges as $N \to \infty$, uniformly in $I$, to some solution of the problem (1). In other words, (2) is fulfilled.

*Proof.* By (ii) we have

(4) $$ \|\eta_{j+1}\| \le (1 + kh)\|\eta_j\| + ha. $$

This implies

$$ \|\eta_{j+1}\| \le (1 + kh)^{j+1} \|y_0\| + \sum_{i=0}^{j} (1 + kh)^i ah \quad (j = 0, 1, \dots, N-1). $$

Hence

(5) $$ \limsup_{N \to \infty} \sup_{t \in I} \|\eta^N(t)\| < \infty. $$

By definition, the derivative $\dot{\eta}^N$ exists a.e. in $I$, and, moreover,

$$ \dot{\eta}^N(t) = \frac{1}{h} (\eta^N(t_{j+1}) - \eta_j^N) \in F(t_j, \eta_j^N) $$

for $t_j < t < t_{j+1}$.

Conditions (ii) and (5) imply that

(6) $$ \limsup_{N \to \infty} \sup_{t \in I} \|\dot{\eta}^N(t)\| < \infty. $$
By Arzela's theorem, $\eta^N(\cdot)$ has a subsequence uniformly convergent to some function $y(\cdot)$, moreover,

$$ \dot{\eta}^N(\cdot) \to z(\cdot) \quad \text{weakly in } L_\infty(I)^n $$

for some further subsequence. It is easy to observe that for all $t \in I$

$$ y(t) = y_0 + \int_{t_0}^{t} z(s)ds, $$

i.e., $y(\cdot)$ is absolutely continuous and $\dot{y}(t) = z(t)$ for almost every $t \in I$. Fix $\epsilon > 0$. From (6)

$$ \sup_{0 \le i \le N-1} \max_{t_i \le t \le t_{i+1}} \|\eta^N(t) - \eta^N(t_i)\| \to 0 \quad \text{as } h \to 0, $$

then by (i) there exists $N_1$ such that for $N > N_1$

$$ (7) \qquad \dot{\eta}^N(t) \in F(t, \eta^N(t)) + \epsilon B \quad \text{for a.e. } t \in I, $$

where $B$ is the unit ball. Pick some $q \in \mathbb{R}^n$ and let $\Delta$ be a measurable subset of $I$. Then from (7)

$$ \int_{\Delta} q^* \dot{\eta}^N(\tau) d\tau \le \int_{\Delta} [\text{supp}(q, F(\tau, \eta^N(\tau))) + \epsilon \|q\|] d\tau, $$

where $\star$ denotes transposition and $\text{supp}$ is the support function. By (i) the map $t \mapsto \text{supp}(q, F(t, \eta^N(t)))$ is upper semi-continuous, and from (ii) and (5) it is bounded above by a constant, independent of $N$. Then, for $N \to \infty$ for the corresponding subsequence, the dominated convergence theorem yields

$$ \int_{\Delta} q^* \dot{y}(\tau) d\tau \le \int_{\Delta} [\text{supp}(q, F(\tau, y(\tau))) + \epsilon \|q\|] d\tau. $$

Since $\Delta$ and $\epsilon$ are arbitrary and $F$ is convex valued, this means that $y(\cdot)$ is a solution of (1). $\square$

Under some more conditions we are able to prove that every solution $\eta^N(\cdot)$ of the discretized inclusion (3) with sufficiently small stepsize contains in its $\epsilon$-neighbourhood in $C(I)^n$ a solution of (1). For that purpose we need the following basic result of the theory of differential inclusions, often referred to as Gronwall-Filippov-Ważewski theorem; see A. F. Filippov [25] or J.-P. Aubin and A. Cellina [4, Thm. 2.4.1, p. 120], for a more general formulation.

**THEOREM 2.3** (Gronwall-Filippov-Ważewski). Let $\bar{y} : I \to \mathbb{R}^n$ be an absolutely continuous function with $\bar{y}(t_0) = y_0$, $\beta$ be a positive constant, and

$$ Q = \{(t, x) \in I \times \mathbb{R}^n : \|x - y(t)\| \le \beta\}. $$

Let $F : Q \to \mathbb{R}^n$ be nonempty closed valued and continuous and satisfy

$$ \text{haus}(F(t, x), F(t, z)) \le k(t) \|x - z\| $$

for all $(t, x)$ and $(t, z)$ from $Q$ with $k(\cdot) \in L_1(I)$.

Assume, moreover, that

$$ \text{dist}(\dot{y}(t), F(t, \bar{y}(t))) \le p(t) \quad \text{for almost all } t \in I, $$
for some $p(\cdot) \in L_1(I)$ such that

$$ \xi(t) = \int_{t_0}^{t} \exp \left[ \int_{s}^{t} k(\tau) d\tau \right] p(s) ds \le \beta \quad (t \in I). $$

Then there exists a solution $y(\cdot)$ to the initial value problem (1) such that

$$ \|y(t) - \bar{y}(t)\| \le \xi(t) $$

for all $t \in I$.

**THEOREM 2.4.** Suppose that the conditions of Theorem 2.2 hold and, moreover, that $F$ is continuous on $I \times \mathbb{R}^n$ and Lipschitz in $x$ on bounded sets in $\mathbb{R}^n$ with integrable in $I$ Lipschitz constant. Then for every $\epsilon > 0$ there exists $N_1$ such that for every $N > N_1$ and for every solution $\eta^N(\cdot)$ of the discrete inclusion (3) there exists a solution $y(\cdot)$ of the problem (1) such that

$$ \max_{t \in I} \|\eta^N(t) - y(t)\| \le \epsilon. $$

*Proof.* Let $\epsilon > 0$. As in the proof of Theorem 2.2 we show that the set $X_h$ of the discrete trajectories is bounded in $C(I)^n$ and, moreover, that there exists $N_1$ such that

$$ F(t_j, \eta^N(t_j)) \subset F(t, \eta^N(t)) + \epsilon B $$

for all $N > N_1$, for all $\eta^N(\cdot) \in X_h$ and for all $t$ with $t_j \le t \le t_{j+1}$, $j = 0, 1, \dots, N-1$. Hence for all $N > N_1$ every solution $\eta^N(\cdot) \in X_h$ will satisfy (7). Applying Gronwall-Filippov-Ważewski theorem we obtain that there exists a solution $y(\cdot)$ of (1) in the $\epsilon$-tube around $\eta^N(\cdot)$. $\square$

In §6, slightly strengthening the assumptions, we obtain an $O(h)$ estimate for the Hausdorff distance between the sets $X_h$ and the set of solutions $X$ of (1).

**3. Convergent multistep methods.** Beginning in 1973 K. Taubert investigated convergence properties of multistep methods for differential equations with discontinuous right-hand sides, later on he carried over his results to initial value problems for differential inclusions; compare [66], [68]. These methods are typically of the following form, where we deliberately suppress possible generalizations to multistage multistep methods.

**LINEAR MULTISTEP METHOD 3.1.** Let

$$ a_i, b_i \in \mathbb{R} \quad (i = 0, \dots, r) $$

with $a_r \ne 0$, $|a_0| + |b_0| > 0$,

$$ h = \frac{T - t_0}{N} \quad (N \in \mathbb{N}' \subset \mathbb{N}). $$

Let there be given starting values

$$ \eta_j \in \mathbb{R}^n \quad (j = 0, \dots, r-1) $$

and corresponding starting selections

$$ \zeta_j \in F(t_j, \eta_j) \quad (j = 0, \dots, r-1), $$
computed, e.g., by a linear $\hat{r}$-step method with $\hat{r} < r$ or by a 1-step method.

For $j = r, \dots, N$ compute $\eta_j$ from

$$ \frac{1}{h} \sum_{i=0}^{r} a_i \eta_{j-r+i} = \sum_{i=0}^{r} b_i \zeta_{j-r+i}, \quad \zeta_j \in F(t_j, \eta_j). $$

These methods are implicit in case $b_r \neq 0$. Generally, the solution of the inclusion

$$ \zeta_j \in F(t_j, \eta_j) $$

is not unique. Then it would be selected randomly or by a suitable optimization criterion. It can be obtained for example as the result of running appropriate phases of a constrained optimization algorithm, see, e.g., [61]. Typically, convergence results of the following form can be proven.

**CONVERGENCE THEOREM 3.2.** Let the following assumptions be satisfied:

(i) $F$ is nonempty closed and convex valued.

(ii) $F$ is bounded and upper semicontinuous in $I \times \mathbb{R}^n$.

(iii) The strong root condition is satisfied, i.e., all zeros $\lambda$ of the polynomial

$$ \sum_{i=0}^{r} a_i \lambda^i $$

have absolute value $|\lambda| < 1$ except the simple zero $\lambda = 1$.

(iv) The method is consistent, i.e.,

$$ \sum_{i=0}^{r} a_i = 0, \quad \sum_{i=0}^{r} ia_i = \sum_{i=0}^{r} b_i. $$

(v) The coefficients $b_i$ are nonnegative ($i = 0, \dots, r$).

(vi) The starting values satisfy

$$ \|\eta_{j+1}^N - \eta_j^N\| \le h M \quad (j = 0, \dots, r-2) $$

for all $N \in \mathbb{N}'$ with a constant $M$, which is independent of the stepsize $h = (T - t_0)/N$.

(vii) The approximations of the initial value $y_0$ satisfy $\lim_{N \to \infty} \eta_0^N = y_0$.

Then the sequence

$$ (\eta^N(\cdot))_{N \in \mathbb{N}'} $$

of piecewise linear continuous interpolants of the grid functions

$$ (\eta_0^N, \eta_1^N, \dots, \eta_N^N) $$

contains a subsequence which converges uniformly to a solution of the initial value problem (1).

Assumptions (i) and (ii) are Peano-type conditions, therefore no better convergence properties could be expected. Assumption (iii) is equivalent to inverse stability of the difference method with respect to the Spijker norm, compare [63, p. 81 and pp. 203–210]. This implies, that the family $(\eta^N(\cdot))_{N \in \mathbb{N}'}$ is uniformly Lipschitz continuous, which is the basis of the convergence proof.
Condition (v) is needed for certain operations in the calculus of sets. It restricts the class of feasible methods decisively. The proof could be extended to discrete Sobolev spaces; compare the papers [45], [44].

For the special selection of coefficients

$$r = 1: \begin{array}{l} a_0 = -1, \quad b_0 = 1, \\ a_1 = 1, \quad b_1 = 0, \end{array}$$

we get again the Euler method.

The following selections of coefficients were tested in [31]:

$$r = 2: \begin{array}{l} a_0 = 0, \quad b_0 = 0.5, \\ a_1 = -1, \quad b_1 = 0.5, \\ a_2 = 1, \quad b_2 = 0, \end{array}$$

respectively,

$$r = 3: \begin{array}{l} a_0 = -0.81, \quad b_0 = \frac{725}{1200}, \\ a_1 = -0.99, \quad b_1 = \frac{488}{1200}, \\ a_2 = 0.8, \quad b_2 = \frac{3119}{1200}, \\ a_3 = 1, \quad b_3 = 0, \end{array}$$

respectively,

$$r = 4: \begin{array}{l} a_0 = 0, \quad b_0 = 0.1, \\ a_1 = 0, \quad b_1 = 0.4, \\ a_2 = 0, \quad b_2 = 0.1, \\ a_3 = -1, \quad b_3 = 0.4, \\ a_4 = 1, \quad b_4 = 0. \end{array}$$

All these methods are consistent and strongly stable. The 3-step method would be consistent of order 3 for single-valued, sufficiently smooth right-hand sides, nevertheless it behaves badly for differential inclusions due to "almost instability." The program code together with some extensions to multistage multistep methods is contained in [1]. We conclude this section with the following simple differential equation with discontinuous right-hand side describing forced vibrations with viscous and Coulomb damping, cf. [55], [67].

*Example 3.3.* Find a function $z(\cdot)$ on $I = [0, T]$ with absolutely continuous derivative such that

$$(8) \qquad \ddot{z}(t) + z(t) + 0.2\dot{z}(t) + 4\operatorname{sgn}(\dot{z}(t)) = 2\cos \pi t$$

for almost all $t \in I$ and

$$z(0) = 3, \quad \dot{z}(0) = 4.$$
Naturally, this initial value problem for a second-order ordinary differential equation has to be reformulated as the following initial value problem for a first-order differential inclusion.

Find an absolutely continuous function $y(\cdot) : I \to \mathbb{R}^2$ such that

$$ 
\begin{aligned}
\dot{y}_1(t) &= y_2(t), \\
\dot{y}_2(t) &\in -y_1(t) - 0.2y_2(t) - 4\operatorname{Sgn}(y_2(t)) + 2 \cos \pi t
\end{aligned}
 $$

for almost all $t \in I = [0, T]$ and

$$ y_1(0) = 3, \quad y_2(0) = 4, $$

where the Sgn-function

$$ \operatorname{Sgn}(x_2) = \begin{cases} -1 & (x_2 < 0), \\ [-1, 1] & (x_2 = 0), \\ 1 & (x_2 > 0), \end{cases} $$

is the set-valued analogue of the usual sign function used in (8) in order to guarantee the existence of a generalized solution in the sense of Filippov. This generalized solution could be approximated with one of the above methods. The Euler method with fine stepsize gives some insight into the structure of the solution; compare the plots of the solution $z(\cdot)$ and its derivative $\dot{z}(\cdot)$ in Fig. 1 and the phase portrait shown in Fig. 2.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Solution $z(\cdot)$ and its derivative $\dot{z}(\cdot)$.</figcaption></figure>

Subintervals, where the solution is in fact a generalized one, are most interesting. The plot in Fig. 3 of the approximation of the derivative $\dot{z}(\cdot)$ by Euler method with step-size $h = 0.005$ shows the typical oscillations to be expected on these nonclassical intervals.

The computer tests show that it is not worthwhile to look for highly consistent standard methods; compare the approximation of $\dot{z}(\cdot)$ in Fig. 4 by the classical Runge-Kutta method with the same stepsize $h = 0.005$.

This additional effort pays only on subintervals where the right-hand side is single-valued and smooth.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Phase portrait of $z(\cdot)$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. Approximation of $\dot{z}(\cdot)$ by Euler method with $h = 0.005$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4. Approximation of $\dot{z}(\cdot)$ by Runge-Kutta method with $h = 0.005$.</figcaption></figure>
**4. One-sided Lipschitz condition and monotonicity.** We could get the impression that it may be hard to prove order of convergence results for general differential inclusions, if they really hold at all. Hence, we have to restrict our analysis to differential inclusions with special right-hand sides. In [41], [42] order of convergence results were proved for differential equations with discontinuous right-hand side, but there it is only allowed that the exact solution crosses directly the discontinuity manifolds; compare in this connection the recent work of R. Model [43]. D. Stewart has defended his Ph.D. thesis on the same problem class in 1990; compare [64]. His main objective is a skilful, but complex transformation of discontinuous right-hand sides into smooth ones on suitably selected submanifolds. This transformation will really pay if the auxiliary complementarity problems occuring in this connection can be analysed easily and if the resulting transformed classical differential equations are not stiff. In the following we analyse the class of right-hand sides which satisfy a uniform one-sided Lipschitz condition. On the one hand, this class is more general, since set-valued right-hand sides are admitted. On the other hand, it is more special, since the one-sided Lipschitz condition implies uniqueness of the solution of problem (1). Stimulated by a paper of C. M. Elliott [23] on first-order convergence for a special class of methods, A. Kastner-Maresch succeeded to adapt convergence proofs for numerical methods for stiff differential equations to this type of differential inclusions; compare [32] and [33]. In principle, in the following we could investigate the class of general linear methods; compare [13]. But for simplicity of presentation we restrict ourselves to the class of implicit s-stage Runge-Kutta methods.

**IMPLICIT s-STAGE RUNGE-KUTTA METHOD 4.1.** *Let there be given the scheme of coefficients*

$$ \begin{array}{c|cccc} c_1 & a_{11} & \dots & a_{1s} \\ \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & \dots & a_{ss} \\ \hline & b_1 & \dots & b_s \end{array} $$

with nonnegative $b_i$ ($i = 1, \dots, s$), a stepsize $h = (T - t_0)/N$ with $N \in \mathbb{N}'$, and an approximation $\eta_0$ of the initial value $y_0$.

For $j = 0, \dots, N-1$ solve the implicit system of inclusions

$$ \bar{\eta}_{\mu} = \eta_j + h \sum_{\nu=1}^{s} a_{\mu\nu} \zeta_{j\nu}, $$

$$ \zeta_{j\nu} \in F(t_j + c_{\nu}h, \bar{\eta}_{\nu}) \quad (\mu, \nu = 1, \dots, s) $$

and compute the next approximation

$$ \eta_{j+1} = \eta_j + h \sum_{\nu=1}^{s} b_{\nu} \zeta_{j\nu}. $$

Without further assumptions, we cannot expect reasonable convergence properties even for single-valued right-hand sides. Moreover, according to the well-known proceeding for stiff differential equations, we need a one-sided Lipschitz condition for set-valued maps. Therefore, we restrict ourselves to the following class of maps.

**DEFINITION 4.2.** Let $\mathbb{R}^n$ be equipped with the scalar product $(\cdot|\cdot)$ and the corresponding induced norm $\|\cdot\|$. A set-valued mapping

$$ F: I \times \mathbb{R}^n \Rightarrow \mathbb{R}^n $$
satisfies a one-sided Lipschitz condition with (one-sided) Lipschitz constant $L_F$ if

$$ (\zeta_1 - \zeta_2|x_1 - x_2|) \le L_F \|x_1 - x_2\|^2 $$

holds for all $x_1, x_2 \in \mathbb{R}^n$ and all $\zeta_1 \in F(t, x_1)$, $\zeta_2 \in F(t, x_2)$ uniformly for all $t \in I$.

An important subclass of all right-hand sides satisfying a one-sided Lipschitz condition consists of set-valued mappings $F: I \times \mathbb{R}^n \to \mathbb{R}^n$ possessing a decomposition

$$ (9) \qquad F(t, x) = f(t, x) - \beta(x) $$

with a single-valued function $f: I \times \mathbb{R}^n \to \mathbb{R}^n$ with (one-sided) Lipschitz constant $L_f$ and a monotone set-valued mapping $\beta: \mathbb{R}^n \to \mathbb{R}^n$. For example, the subdifferential $\partial\varphi$ of a convex functional $\varphi: \mathbb{R}^n \to \mathbb{R}$ defines a monotone mapping on $\mathbb{R}^n$. In [23] special right-hand sides of the type (9) were treated by a special class of methods. In fact, not this decomposition, but the validity of a one-sided Lipschitz condition is the basis for convergence and order of convergence proofs following proof structures for stiff differential equations, compare the books [20], [13] and the journal articles [12], [29], [28], [9], [10], especially the survey of K. Burrage [11]. For consistency proofs only smoothness properties of the solution of the differential inclusion are needed, not smoothness properties of the right-hand side as in the classical approach. For proofs of stability properties like C-stability, BS-stability, BSI-stability the fact is exploited that the above one-sided Lipschitz condition holds uniformly with respect to the selections.

A typical convergence result due to A. Kastner-Maresch [32] reads as the following.

**THEOREM 4.3.** Let the following assumptions be satisfied:

(i) Let $F$ be upper semicontinuous on $I \times \mathbb{R}^n$ and satisfy a one-sided Lipschitz condition.

(ii) The (necessarily unique) solution of initial value problem (1) is piecewise $(\mu + 1)$-times continuously differentiable with $\mu \ge 1$.

(iii) The simplifying conditions $B(\mu)$, i.e.,

$$ \sum_{j=1}^{s} b_j c_j^{k-1} = \frac{1}{k} \quad (k = 1, \dots, \mu), $$

and $C(\mu)$, i.e.,

$$ \sum_{j=1}^{s} a_{ij} c_j^{k-1} = \frac{1}{k} c_i^k \quad (k = 1, \dots, \mu; i = 1, \dots, s), $$

hold.

(iv) The method is BS-stable and C-stable.

(v) The initial approximations $\eta_0^N$ of $y_0$ satisfy

$$ \|\eta_0^N - y_0\| = \mathcal{O}(h). $$

Then the order of convergence is equal to one.

Crucial for the proof is the fact that the first derivative of the exact solution has only finitely many jump discontinuities because each of these discontinuities contributes a term of order one to the global discretization error. Hence, on intervals where the solution is smooth, we get higher-order convergence if the initial error is of higher-order as well.

**COROLLARY 4.4.** Let in Theorem 4.3 the exact solution be $(\mu + 1)$-times continuously differentiable on the whole interval $I$, and the initial approximations satisfy

$$ \|\eta_0^N - y_0\| = \mathcal{O}(h^\mu). $$
Then the order of convergence is equal to $\mu$.

This corollary serves as a basis for more or less sophisticated algorithms with built-in localization procedures for detecting possible discontinuities of the derivative of the exact solution. In practice, we have to localize numerically manifolds where the right-hand side is discontinuous. In principle, localization procedures in [41], [42], or [64] could be used, but they must be adapted to the above situation.

Especially, the implicit midpoint rule was tested numerically, which could be written as the following Runge-Kutta method

$$ 
\begin{aligned}
\bar{\eta} &= \eta_j + h \frac{1}{2} \zeta_j, \\
\zeta_j &\in F \left( t_j + \frac{1}{2} h, \bar{\eta} \right), \\
\eta_{j+1} &= \eta_j + h \zeta_j,
\end{aligned}
$$

with Butcher-array

$$ \frac{\frac{1}{2} \mid \frac{1}{2}}{1}, $$

or in more condensed form,

$$ \eta_{j+1} \in \eta_j + h F \left( t_j + \frac{1}{2} h, \frac{1}{2} (\eta_j + \eta_{j+1}) \right). $$

For this method the simplifying conditions B(1), B(2), and C(1) (see Theorem 4.3) hold. Hence the order of this method is equal to one if the solution is piecewise 2-times continuously differentiable; compare [23]. The order 2 can be proved on subintervals where the solution is 3-times continuously differentiable and the initial error is of order 2; here we have to exploit results in [35]. Together with a suitable localization procedure, we get the order 2 on the whole interval.

Applying, e.g., the implicit midpoint rule with stepsize $h = 0.005$ and localization procedure to Example 3.3 yields the result plotted in Fig. 5.

Naturally, the numerical effort is somewhat higher than for classical Runge-Kutta method with the same stepsize. But the results are much better; compare Fig. 4. Even with stepsize $h = 0.001$, the classical Runge-Kutta method would yield worse results than implicit midpoint rule with stepsize $h = 0.005$; compare Fig. 6.

In fact, for such a fine stepsize $h = 0.001$ the classical Runge-Kutta method needs more CPU-time and function evaluations than the implicit midpoint rule for $h = 0.005$, i.e., the higher-order of convergence of the implicit method at last beats the explicit method.

**5. Selection strategies.** Contrary to the differential inclusions treated in §4, general differential inclusions normally have a whole bundle of solutions. Hopefully, the qualitative properties of the approximated solution could be improved by more elaborate selection strategies, as compared with the random selection of Linear Multistep Method 3.1.

ALGORITHM 5.1 (selection with minimal norm). Choose a fixed element $z \in \mathbb{R}^n$ and minimize for $j = r, \dots, N$

$$ \|z - \zeta_j\|, $$
<figure><img src="image_7.png" /><figcaption>FIG. 5. Approximation of $\dot{z}(\cdot)$ by implicit midpoint rule with $h = 0.005$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 6. Approximation of $\dot{z}(\cdot)$ by Runge-Kutta method with $h = 0.001$.</figcaption></figure>

subject to

$$ \frac{1}{h} \sum_{i=0}^{r} a_i \eta_{j-r+i} = \sum_{i=0}^{r} b_i \zeta_{j-r+i}, \\ \zeta_j \in F(t_j, \eta_j). $$

For $z = 0_{\mathbb{R}^n}$, in each step we get the selection with minimal norm. In any case, the explicit difference method ($b_r = 0$) for the differential inclusion can be regarded as linear multistep method for the differential equation

$$ \dot{y}(t) = f(t, y(t)) $$

with

$$ f(t, x) = \mathrm{pr}_{F(t,x)}(z) \quad ((t, x) \in I \times \mathbb{R}^n), $$
where $\text{pr}_{F(t,x)}(z)$ is the projection of $z$ onto $F(t, x)$. Hence, the properties of $f(t, x)$ determine the properties of the algorithm. The calculation of starting values $\eta_0, \dots, \eta_{r-1}$ and starting selections $\zeta_0, \dots, \zeta_{r-1}$ has to be adjusted accordingly. Convergence is proved by Convergence Theorem 3.2. Using arguments in [25] we get the following.

**THEOREM 5.2.** In addition to the assumptions in Convergence Theorem 3.2, let $F$ be Hausdorff continuous on $I \times \mathbb{R}^n$ and $b_r = 0$. Then the sequence

$$ (\eta^N(\cdot))_{N \in \mathbb{N}} $$

of piecewise linear continuous interpolants of the grid functions corresponding to minimal norm selections contains a subsequence which converges uniformly to a continuously differentiable solution of the differential inclusion.

The following test problem from [34] is treated numerically as an example for the various selection strategies of this section.

**Example 5.3.** Find an absolutely continuous function $y: I \to \mathbb{R}^2$ such that

$$ \dot{y}_1(t) \in -y_2(t) + 0.1y_1(t)(9 - y_1(t)^2 - y_2(t)^2) + [-1, 1], $$

$$ \dot{y}_2(t) \in y_1(t) + 0.1y_2(t)(9 - y_1(t)^2 - y_2(t)^2) + [-1, 1] $$

for almost all $t \in I = [0, T]$ and

$$ y_1(0) = 5, \quad y_2(0) = 0. $$

To avoid all technical difficulties with the calculation of starting values, we apply Euler method. Note that no order of convergence result is available; therefore, the use of real multistep methods is justified anyway only on subintervals where the right-hand side degenerates to a single-valued function which is sufficiently smooth. The approximate solution for the selection with minimal norm is plotted in Fig. 7 with stepsize $h = 0.005$. The phase portrait is given in Fig. 8; the chosen selections are plotted in Fig. 9.

<figure><img src="image_9.png" /><figcaption>FIG. 7. Approximate solution $\eta^N(\cdot)$ for selection with minimal norm.</figcaption></figure>

In the following selection strategy we choose selections which have minimal variation in a certain sense, cf. the proof of Theorem 2 on page 115 in [4].

**ALGORITHM 5.4 (selection with minimal variation).** For $j = r, \dots, N$ minimize

$$ \|v_{j-1} - v_j\|, $$
<figure><img src="image_10.png" /><figcaption>FIG. 8. Phase portrait of $\eta^N(\cdot)$ for selection with minimal norm.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9. Selection $\zeta^N(\cdot)$ for selection with minimal norm.</figcaption></figure>

subject to

$$ \frac{1}{h} \sum_{i=0}^{r} a_i \eta_{j-r+i} = v_j, \\ v_j \in \sum_{i=0}^{r} b_i F(t_{j-r+i}, \eta_{j-r+i}). $$

If, in addition to the assumptions in the Convergence Theorem 3.2, the set-valued mapping $F$ is Lipschitz continuous in $I \times \mathbb{R}^n$; then a similar argument as in the proof of Theorem 3.2 applies to the selections $v_j$, and we get for explicit methods ($b_r = 0$).

**THEOREM 5.5.** In addition to the assumptions in Convergence Theorem 3.2, let $F$ be Lipschitz continuous on $I \times \mathbb{R}^n$, $b_r = 0$, and the starting selections satisfy

$$ \|v_{j+1}^N - v_j^N\| \le hM \quad (j = 0, \dots, r-2) $$

for all $N \in \mathbb{N}'$.
$$ (\eta^N(\cdot), v^N(\cdot))_{N \in \mathbb{N}} $$

of piecewise linear continuous interpolants of the grid functions and selections contains a subsequence which converges uniformly to a pair $(y, v)$, where $y$ is a solution of the differential inclusion with Lipschitz continuous derivative $v$.

The approximate solution of Example 5.3 for the selection with minimal variation, computed by Euler method with stepsize $h = 0.005$, is plotted in Fig. 10, the phase portrait given in Fig. 11, the according selections with minimal variation in Fig. 12.

<figure><img src="image_12.png" /><figcaption>FIG. 10. Approximate solution $\eta^N(\cdot)$ for selection with minimal variation.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 11. Phase portait of $\eta^N(\cdot)$ for selection with minimal variation.</figcaption></figure>

We get a further selection strategy by comparing the approximations at each grid point with a reference trajectory.

**ALGORITHM 5.6 (selection with respect to reference trajectory).** Choose a Lipschitz continuous reference trajectory

$$ \bar{y}: I \to \mathbb{R}^n. $$
<figure><img src="image_14.png" /><figcaption>FIG. 12. Selection $v^N(\cdot)$ for selection with minimal variation.</figcaption></figure>

For $j = r, \dots, N$, minimize

$$ \| \bar{y}(t_j) - \eta_j \|, $$

subject to

$$ \frac{1}{h} \sum_{i=0}^{r} a_i \eta_{j-r+i} = \sum_{i=0}^{r} b_i \zeta_{j-r+i}, $$

$$ \zeta_{j-r+i} \in F(t_{j-r+i}, \eta_{j-r+i}) \quad (i = 0, \dots, r). $$

For this selection strategy with reference trajectory $\bar{y}(t) \equiv 0_{\mathbb{R}^2}$, the Euler method with stepsize $h = 0.005$ gives the approximate solution of Example 5.3 plotted in Fig. 13, with phase portrait in Fig. 14.

<figure><img src="image_15.png" /><figcaption>FIG. 13. Approximate solution $\eta^N(\cdot)$ for selection with respect to reference trajectory.</figcaption></figure>

Most interesting is the comparison of the chosen selections; compare Fig. 15 with the results of the other selection strategies.

Contrary to the selections with minimal norm, cf. Fig. 9, and the selections with minimal variation, cf. Fig. 12, the selections with respect to reference trajectory now
<figure><img src="image_16.png" /><figcaption>FIG. 14. Phase portrait of $\eta^N(\cdot)$ for selection with respect to reference trajectory.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 15. Selection $\zeta^N(\cdot)$ for selection with respect to reference trajectory.</figcaption></figure>

indicate jump discontinuities of the derivatives of the solution, related to the bang-bang principle in optimal control. Accordingly, no smoothness properties can be expected for the approximated solution, but an error estimate for this selection strategy will be obtained in the next section.

At any rate, it would be worthwhile to extend the above qualitative results by estimations of the order of convergence, eventually by exploiting concepts of numerical methods for optimal control problems with nonstandard objective functions.

**6. Error estimates.** We are interested in estimating the distance between the sets of solutions of Initial Value Problem 1.1 and of the discrete inclusions (3), respectively, in a reasonable sense, by the stepsize $h$.

The first result in this direction probably is due to B. N. Pshenichny [53], who proved the following.

**THEOREM 6.1 ([53]).** Let $F: \mathbb{R}^n \to \mathbb{R}^n$ be a compact convex valued mapping and $y$ be a solution of
$\dot{y}(t) \in F(y(t)) \quad \text{for a.e. } t \in I, \quad y(t_0) = y_0$.

Suppose that there exists $\epsilon > 0$ such that $F$ is Lipschitz in the set $\{x \in \mathbb{R}^n : \|x - y(t)\| \le \epsilon$ for some $t \in I\}$. Then there exist $c$ and $N_1$ such that for all $N > N_1$ there exists a solution $\eta^N(\cdot)$ of the discretized inclusion

$$ (10) \qquad \eta_{j+1} \in \eta_j + hF(\eta_j), \quad j = 0, 1, \dots, N-1, \qquad \eta_0 = y_0 $$

such that

$$ \max_{0 \le j \le N} \|\eta^N(t_j) - y(t_j)\| \le ch. $$

In the proof, B. N. Pshenichny used the following construction of the discrete trajectory $\eta^N(\cdot)$: Let $\eta_0^N = y_0$ and define $\eta_{j+1}^N$ as the projection of $y(t_{j+1})$ on the set $\eta_j^N + hF(\eta_j^N)$. Then from (10), from the inclusion

$$ y(t_{j+1}) \in y(t_j) + hF(y(t_j)) + \mathcal{O}(h^2)B, $$

and from the choice of $\eta_{j+1}^N$, using the Lipschitz continuity of $F$, we obtain for some $a_1$ and $a_2$ that

$$ \|\eta_{j+1}^N - y(t_{j+1})\| \le (1 + a_1 h)\|\eta_j^N - y(t_j)\| + a_2 h^2 $$

for $j = 0, 1, \dots, N-1$. This implies the desired estimate.

The construction of B. N. Pshenichny is extended in [22] in the following way: Let $\bar{y} : I \to \mathbb{R}^n$ be a Lipschitz function. We obtain $\eta^N(\cdot)$ by the following procedure:

$$ (11) \qquad \begin{array}{ll} \eta_0^N = y_0, & \eta_{j+1}^N \text{ is the projection of } \bar{y}(t_{j+1}) \\ & \text{on the set } \eta_j^N + hF(t_j, \eta_j^N), \quad j = 0, 1, \dots, N-1. \end{array} $$

Combining the result of B. N. Pshenichny with Gronwall-Filippov-Ważewski theorem, we obtain that

$$ (12) \quad \max_{0 \le j \le N} \|\eta_j^N - \bar{y}(t_j)\| \le c \left( \|\bar{y}(t_0) - y_0\| + \int_{t_0}^T \operatorname{dist}(\dot{y}(t), F(t, \bar{y}(t))) dt + h \right) $$

for some $c$ independent of $h$.

This estimate was obtained in [22] by a direct proof. More precisely, the following theorem is proved.

**THEOREM 6.2 ([22]).** Suppose that the map $F : I \times \mathbb{R}^n \Rightarrow \mathbb{R}^n$ satisfies the following conditions.

$F$ is nonempty compact and convex valued; there exist $k$ and $a$ such that $\|f\| \le k\|x\|+a$ whenever $f \in F(t, x)$, $t \in I$, $x \in \mathbb{R}^n$; $F$ is Lipschitz in $x$ on bounded sets uniformly in $t$; $F$ is of bounded variation in $t$ uniformly in $x$ on bounded sets, i.e., for any bounded set $U \subset \mathbb{R}^n$

$$ \sup \left\{ \sum_{i=1}^{\nu} \sup_{x \in U} \{ \text{haus}(F(t_{i+1}, x), F(t_i, x)) \} : t_0 \le t_1 \le \dots \le t_{\nu+1} \le T, \nu \in \mathbb{N} \right\} < \infty. $$
Then there exists a constant c such that if $\eta^N(\cdot)$ is obtained by (11), then for all $N$ (12) is fulfilled.

In particular, if $\bar{y}$ is a solution of Initial Value Problem 1.1 then $\eta^N(\cdot)$ provides a first-order approximation to $\bar{y}$. The procedure (11) can be regarded as a "reference trajectory" strategy which can be used for numerical calculations; compare the plots in Fig. 13, 14, and 15. Some additional numerical examples are contained in [22].

Let $\eta^N(\cdot)$ be a piecewise linear function that solves (3). If $F$ is Lipschitz in both $x$ and $t$ then

$$ \dot{\eta}^N(t) \in F(t, \eta^N(t)) + hkB \quad \text{for a.e. } t \in [t_0, T] $$

for some $k$ and for all $N \in \mathbb{N}$. Then one can apply Gronwall-Filippov-Ważewski theorem, obtaining that there exists a solution $\tilde{y}^N$ of Initial Value Problem 1.1 such that

$$ (13) \qquad \max_{t_0 \le t \le T} \|\tilde{y}^N(t) - \eta^N(t)\| = \mathcal{O}(h). $$

In fact, to obtain (13) it is sufficient to assume the conditions of Theorem 6.2. Namely, we have

**THEOREM 6.3 ([22]). On the conditions of Theorem 6.2**

$$ \text{haus}(X, X_h) = \mathcal{O}(h) \quad \text{in } C(I)^n. $$

In [22] an averaged modulus of smoothness for set-valued maps is introduced and its properties are investigated. This helps to weaken the requirements for the map $F$.

Sometimes it may be convenient to use an approximation of the set-valued map $F$ in (10), e.g., by polygons. Let $\eta^N(\cdot)$ satisfy

$$ \eta_{j+1} \in \eta_j + hP_j, \quad j = 0, 1, \dots, N-1, \quad \eta_0 = y_0, $$

where

$$ \text{haus}(P_j, F(\eta_j)) = \mathcal{O}(h), \quad j = 0, 1, \dots, N-1. $$

Then on the assumptions of Theorem 6.2, using the same argument as above, we obtain first-order convergence of this approximation.

In [34], Theorem 6.2 is extended for the Algorithm 5.6 (selection with respect to reference trajectory), using Linear Multistep Method 3.1.

**THEOREM 6.4.** Assume that in addition to the assumptions in Convergence Theorem 3.2, $b_r = 0$, and the following conditions hold:

(i) $F$ is Lipschitz in $x \in \mathbb{R}^n$ uniformly in $t \in I$ and of bounded variation in $t \in I$ uniformly in $x$ on bounded sets.

(ii) the coefficients $a_i$ are nonnegative except the single coefficient $a_0$; or alternatively, exactly two of the coefficients $a_i$ ($i = 0, 1, \dots, r$) are nonzero.

(iii) The approximations of the initial value $y_0$ satisfy

$$ \|\eta_0^N - y_0\| = \mathcal{O}(h). $$

Then there exists a constant $c$ such that for all sufficiently small $h$, if $\eta^N(\cdot)$ is obtained by Algorithm 5.6, then (12) holds.
**7. Convergence of reachable sets.** Analogous to the reachable set of the differential inclusion (1), we can introduce a reachable set of its discrete approximation. Consider the Euler scheme (3) and let $i$ be a fixed number from the set $\{0, 1, \dots, N\}$. Then the reachable set, associated with (3) at the time $t_i$, will be

$$R_i^N = \{x \in \mathbb{R}^n : x = \eta^N(t_i) \text{ for some } \eta^N(\cdot) \in X_h\}.$$

**THEOREM 7.1.** Let $F$ be a compact convex and nonempty valued map that is continuous in both $t$ and $x$ in $I \times \mathbb{R}^n$ and let the reachable set $R(t)$ of Initial Value Problem 1.1 be contained in a bounded set $K$ in $\mathbb{R}^n$ for all $t \in I$. Let $F$ be Lipschitz in a neighbourhood of $K$ with respect to $x$ uniformly in $t \in I$. Then

$$\max_{0 \le i \le N} \text{haus} (R_i^N, R(t_i)) \to 0$$

as $N \to \infty$.

The proof of this theorem can be easily obtained from [50], for a more recent presentation see [49]. First we need some notation.

**DEFINITION 7.2 ([50]).** The closed valued and continuous map $R: I \to \mathbb{R}^n$ with $R(t_0) = y_0$ is an $R$-solution of Initial Value Problem 1.1 if and only if

$$(14) \qquad \frac{1}{h} \text{haus} \left( R(t+h), \bigcup_{x \in R(t)} (x+hF(t,x)) \right) \to 0$$

as $h \to 0$ uniformly in $t \in I$.

**THEOREM 7.3 ([50]).** On the conditions of Theorem 7.1 there exists a unique $R$-solution of Initial Value Problem 1.1 whose value at any $t \in I$ is the reachable set $R(t)$ of Initial Value Problem 1.1.

The relation (14) is known as *funnel equation*. The integral funnel (solution cone) of Initial Value Problem 1.1 is the union of the graphs of all solutions of Initial Value Problem 1.1. The reachable set $R(t)$ is the value at $t$ of the map, the graph of which is the integral funnel.

In [69] extensions of Theorem 7.3 for inclusions with right-hand side satisfying Carathoory conditions are obtained. A funnel equation to a linear differential inclusion with phase constraints is found in [38].

Let us prove Theorem 7.1. From (14) we have

$$(15) \qquad \text{haus} (R_1^N, R(t_1)) \le h\varphi(h)$$

and

$$(16) \qquad \text{haus} \left( R(t_2), \bigcup_{x \in R(t_1)} (x+hF(t_1,x)) \right) \le h\varphi(h),$$

where $\varphi(h) \to 0$ as $h \to 0$ uniformly in $t \in I$. If $k$ is the Lipschitz constant of $F$, it is easy to see from (15) that

$$(17) \qquad \text{haus} \left( \bigcup_{x \in R_1^N} (x+hF(t_1,x)), \bigcup_{x \in R(t_1)} (x+hF(t_1,x)) \right) \le (1+kh)h\varphi(h).$$
Combining (16) and (17) we have

$$ \text{haus}(R_2^N, R(t_2)) \le (1+kh)h\varphi(h) + h\varphi(h). $$

Proceeding in the same manner by induction we get

$$ \begin{aligned} \text{haus}(R_i^N, R(t_i)) &\le \sum_{\mu=1}^i (1+kh)^{i-\mu}h\varphi(h) \\ &\le \varphi(h)(T-t_0)e^{k(T-t_0)} \end{aligned} $$

for $i = 1, \dots, N$, which completes the proof of Theorem 7.1.

Consider now the Initial Value Problem 1.1 with $F$ independent of $t$. Then Euler approximation to the attainable set can be rewritten as

$$ \begin{aligned} R_1^N &= (I + hF)y_0, \\ R_2^N &= \bigcup_{x \in R_1^N} (x + hF(x)) = (I + hF)^2 y_0, \\ \vdots \\ R_i^N &= (I + hF)^i y_0, \end{aligned} $$

where the power of $(I + hF)$ is that of composition of set-valued maps. Then $R_N^N$ will converge to $R(T)$ if and only if

$$ R(T) = \lim_{N \to \infty} \left( I + \frac{T - t_0}{N} F \right)^N y_0. $$

More precisely, the following result is obtained in [74].

**THEOREM 7.4 ([74]).** Let $G \subset \mathbb{R}^n$ be open and $F : \mathbb{R}^n \Rightarrow \mathbb{R}^n$ be a map with nonempty compact values on $G$ locally Lipschitz on $G$. Fix $y_0 \in G$. Let

$$ \hat{T} = \sup \left\{ T : \text{cl} \left( \bigcup_{t_0 \le t \le T} R(t) \right) \text{ is compact} \right\}. $$

$$ (i) \text{ For } t_0 \le T < \hat{T} \\ \limsup_{N \to \infty} \left( I + \frac{T - t_0}{N} F \right)^N y_0 \subset \text{cl}(R(T)). $$

(ii) If, in addition, $F$ is assumed to have convex values, then for all $T \ge 0$,

$$ R(T) \subset \liminf_{N \to \infty} \left( I + \frac{T - t_0}{N} F \right)^N y_0. $$

On the assumptions of Theorem 6.2 we obtain that the Euler scheme provides an $O(h)$ estimate for the reachable set and the integral funnel. We note that the one-sided estimate (13) was observed (on slightly different conditions) independently in [22] and by M. S. Nikol'skiy [47] whose primary purpose was to estimate the Euler approximation to the reachable set. However, to obtain an approximation of a solution of Initial Value
Problem 1.1 by a solution of (3), he uses a procedure different from the projection strategy, giving an estimate of order $O(\sqrt{h})$ only. The paper [46] contains a similar result for a control system. In [48] the following modification of Euler scheme is proposed:

$$ 
\begin{aligned}
E_0^N &= \{y_0\}, U = a_1 B, \\
E_{i+1}^N &= \bigcup_{x \in E_i^N \cap U} (x + h\gamma(x)\co F(x) + ah^2 B), \quad i = 0, 1, \dots, N-1,
\end{aligned}
 $$

where the right-hand side $F$, independent of $t$, satisfies a growth condition and is locally Lipschitz and the function $\gamma$ and the constants $a, a_1$ are chosen in such a way that $\text{cl}(R(T)) \subset E_N^N \cap U$. Then from (13) an estimate of order $O(h)$ follows for the distance between $E_N^N \cap U$ and $R(T)$. For an extension, see [57].

**8. Higher-order approximations to reachable sets.** Generally speaking, a differential inclusion may have many “bad” solutions, e.g., nonsmooth and fast changing ones, to which it is unlikely to obtain higher-order approximations. Nevertheless, V. M. Veliov recently proposed in a series of papers [71]–[73] certain second-order approximations resulting from Runge-Kutta schemes, both for the trajectory bundle and the reachable sets.

To be specific, consider the simplest linear differential inclusion

$$ (18) \qquad \dot{y}(t) \in Ay(t) + U \quad \text{for a.e. } t \in I, \quad y(t_0) = y_0, $$

where $A$ is a $n \times n$-matrix, $U$ is a compact and convex set in $\mathbb{R}^n$. An absolutely continuous function $y$ is a solution of (18) if and only if there exists an integrable $n$-vector function $u$ with $u(t) \in U$ for almost every $t \in I$ such that

$$ (19) \qquad \dot{y}(t) = Ay(t) + u(t) \quad \text{for a.e. } t \in I. $$

It is clear that by applying a higher-order scheme to the equation (19), we could hardly expect higher-order accuracy since the function $u$ may be discontinuous. However, a higher-order approximation $\eta^N$ to $y$ may exist, being a solution of a discrete time inclusion resulting from the scheme, that corresponds possibly to another selection of the right-hand side of (18).

Let $N$ and $k$ be integers, $h = (T - t_0)/N$, $A_N^k$ be $n \times n$-matrices and $U_N^k$ be sets in $\mathbb{R}^n$. By a discrete approximation (of higher-order) to (18), we mean the following:

$$ (20) \qquad \eta_{i+1}^N \in A_N^k \eta_i^N + U_N^k, \quad i = 0, 1, \dots, N-1, \quad \eta_0^N = y_0. $$

Suppose that $A_N^k, U_N^k$ satisfy for $N \in \mathbb{N}'$,

$$ (21) \qquad \exp(Ah) = A_N^k + \mathcal{O}(h^{k+1}), $$

$$ (22) \qquad \int_0^h \exp(As)Uds = U_N^k + \mathcal{O}(h^{k+1}). $$

(i) Let $\eta^N = (\eta_0^N, \dots, \eta_N^N)$ be a solution of (20). Then, from (21) and (22),

$$ \eta_{i+1}^N \in A_N^k \eta_i^N + U_N^k \subset e^{A(t_{i+1}-t_i)} \eta_i^N + \int_{t_i}^{t_{i+1}} e^{A(t_{i+1}-s)} Uds + \mathcal{O}(h^{k+1}) $$

for $i = 0, 1, \dots, N-1$. Hence there exists a measurable function $u_i^N$ with $u_i^N(t) \in U$ for almost every $t \in I$, such that if $\bar{y}_i^N$ solves

$$ 
\begin{aligned}
\dot{\bar{y}}_i^N(t) &= A\bar{y}_i^N(t) + u_i^N(t) \quad \text{for a.e. } t \in [t_i, t_{i+1}], \\
\bar{y}_i^N(t_i) &= \eta_i^N,
\end{aligned}
 $$
then

$$
\bar{y}_i^N(t_{i+1}) = \eta_{i+1}^N + \mathcal{O}(h^{k+1}).
$$

Let $u^N$ be defined on $I$ as $u^N(t) = u_i^N(t)$ for $t \in [t_i, t_{i+1})$, $i = 0, 1, \dots, N-1$. If $\tilde{y}^N$ is a solution of

$$
\begin{align*}
\dot{\tilde{y}}^N(t) &= A\tilde{y}^N(t) + u^N(t) \quad \text{for a.e. } t \in I \\
\tilde{y}^N(t_0) &= y_0,
\end{align*}
$$

then it is easy to see that

$$
(23) \qquad \max_{0 \le i \le N} \| \tilde{y}^N(t_i) - \eta_i^N \| = \mathcal{O}(h^k).
$$

(ii) Now let $y$ be a solution of (18). Then

$$
\begin{align*}
y(t_{i+1}) &\in e^{Ah}y(t_i) + \int_{t_i}^{t_{i+1}} e^{A(t_{i+1}-s)}Uds \\
&\subset A_N^k y(t_i) + U_N^k + \mathcal{O}(h^{k+1}).
\end{align*}
$$

Let $\eta_0^N = y_0$. We determine $\eta_{i+1}^N$ on the basis of $\eta_i^N$ as projection of $y(t_{i+1})$ on the set $A_N^k \eta_i^N + U_N^k$. Then from

$$
y(t_{i+1}) \in A_N^k \eta_i^N + U_N^k - A_N^k (\eta_i^N - y(t_i)) + \mathcal{O}(h^{k+1})
$$

it follows that

$$
\| \eta_{i+1}^N - y(t_{i+1}) \| \le \| A_N^k \| \| \eta_i^N - y(t_i) \| + \mathcal{O}(h^{k+1}),
$$

which yields

$$
(24) \qquad \max_{0 \le i \le N} \| \eta_i^N - y(t_i) \| = \mathcal{O}(h^k).
$$

Thus, under the conditions (21) and (22), the Hausdorff distance between the sets of solutions of (18) and (20), in the sense of (23) and (24), is of order $\mathcal{O}(h^k)$.

The condition (21) will hold if we take

$$
A_N^k = \sum_{i=0}^{k} \frac{h^i}{i!} A^i.
$$

It turns out that the condition (22) is much more restrictive. It obviously holds for $k = 1$
(Euler scheme); moreover, we have always

$$
\int_0^h \exp(As) Uds \subset U_N^k + \mathcal{O}(h^{k+1}),
$$

provided that

$$
(25) \qquad U_N^k = h \sum_{i=0}^{k-1} \left( \frac{h^i A^i}{(i+1)} U \right).
$$
This means that (24) is satisfied, i.e., every solution of (18) can be approximated of order $O(h^k)$ by a solution of (20). The converse case needs, however, more conditions.

It follows from [73] that for $k = 2$ (25) implies (22) on the condition that $U$ is a strongly convex set, i.e., there exists a constant $\mu > 0$ such that $u_1, u_2 \in U$ implies

$$0.5(u_1 + u_2) + 0.25\mu\|u_1 - u_2\|l \in U \quad \text{for every } \|l\| \le 1.$$

The proof of this is based on the observation that for strongly convex $U$ the extremal (boundary) controls are equi-Lipschitz. This means that if $U$ is strongly convex, then the scheme

$$(26) \quad \eta_{i+1}^N \in \left( I + hA + \frac{1}{2}h^2 A^2 \right) \eta_i^N + h \left( I + \frac{h}{2} A \right) U, \quad i = 0, 1, \dots, N-1, \quad \eta_0^N = y_0$$

gives $O(h^2)$ approximation to the set of trajectories, and hence to the reachable set as well.

In [73] V. M. Veliov considers the following more general problem:

$$(27) \qquad \dot{y}(t) \in F(t, y(t)) \quad \text{for a.e. } t \in I, \quad y(t_0) \in R_0$$

approximated by

$$(28) \quad \eta_{i+1}^N \in \eta_i^N + 0.5h \{z + F(t_{i+1}, \eta_i^N + hz) : z \in F(t_i, \eta_i^N)\}, \quad \eta_0^N \in R_0$$

for $i = 0, 1, \dots, N-1$. This approximation is a generalization (but not the only possible one) of Euler scheme; it corresponds to the classical Euler-Cauchy method.

The map $F$ is defined on the compact and convex set $\Delta \times S \subset \mathbb{R}^{n+1}$ and satisfies the following conditions:

(a) $F$ is compact valued and there exists $\mu > 0$ such that $F(t, x)$ is strongly convex for every $t \in \Delta$ and $x \in S$. The support function $\tau(l, t, x) = \text{supp}(l, F(t, x))$ is differentiable with respect to $t$ and $x$, $\partial\tau/\partial x$ is Lipschitz in the set $(l, t, x) \in \mathcal{L} \times \Delta \times S$, $\partial\tau/\partial t$ is Lipschitz in $(l, t) \in \mathcal{L} \times \Delta$ uniformly in $x \in S$, where $\mathcal{L} = \{l \in \mathbb{R}^n : 0.5 \le \|l\| \le 2\}$; the set $R_0$ is compact, $I \subset \text{int}(\Delta)$, and all trajectories of (27) stay in $\text{int}(S)$.

(b) The function $\hat{y}(t, x, l)$ defined as

$$(29) \qquad l^* \hat{y}(t, x, l) = \text{supp}(l, F(t, x))$$

is Lipschitz in $t \in \Delta$ uniformly in $(l, x) \in \mathcal{L} \times S$.

**THEOREM 8.1 ([73]).** On the condition (a) there exist constants $c$ and $N_1$ such that if $N > N_1$ then

(i) For every solution $y$ of (27) there is a discrete solution $\eta^N = (\eta_0^N, \dots, \eta_N^N)$ of (28) such that

$$(30) \qquad \max_{0 \le i \le N} \|y(t_i) - \eta_i^N\| \le ch^{3/2}.$$

(ii) For every solution $\eta^N$ of (28) there exists a solution $y$ of (27) such that (30) holds. If in addition (b) is satisfied, then $h^{3/2}$ in (30) can be replaced by $h^2$.

In [73] it is also shown that if the map $F$ is

$$F(t, x) = f(t, x) + g(t, x)U$$
then the conditions (a) and (b) will be fulfilled if $U$ is strongly convex, $f$ and $g$ are sufficiently smooth and $g(t, x)$ is invertible for each $t \in \Delta, x \in S$. In this case the discrete inclusion (28) can be equivalently written as

$$ \eta_{i+1}^N \in \eta_i^N + 0.5h \{ \varphi(t_i, \eta_i^N, u) + \varphi(t_{i+1}, \eta_i^N + h\varphi(t_i, \eta_i^N, u)) : u \in U \}, $$

where $\varphi(t, x, u) = f(t, x) + g(t, x)u$. For fixed $u \in U$ this is exactly the second-order Euler-Cauchy scheme for right-hand side $\varphi(t, x, u)$. Another second-order approximation for the trajectory bundle without strong convexity of the right-hand side is obtained in [72] for the inclusion

$$ \dot{y}(t) \in f(t, y(t)) + \sum_{i=1}^{r} g_i(t, y(t))[0, 1], \quad y(t_0) \in \mathbb{R}_0. $$

It turns out that for linear control systems it is possible to obtain second-order approximations to the reachable set having, in the same time, only first-order convergence of the solutions, in general. Consider the differential inclusion

$$ (31) \qquad \dot{y}(t) \in A(t)y(t) + B(t)U \quad \text{for a.e. } t \in I, \quad y(t_0) \in X_0, $$

where $A(t)$ and $B(t)$ are $n \times n$ and $n \times r$ matrices, $U \subset \mathbb{R}^r$, on the following assumption:

(c) $A$ and $B$ are differentiable and their derivatives are Lipschitz continuous, and $X_0$ and $U$ are convex and compact.

Let $R(T)$ be the reachable set of (31) at the time $T$. Introducing a uniform grid in $I$, for each $N$ we define the set $X_N^N$ recurrently by the equation

$$ (32) \qquad X_{i+1}^N = A_i^N X_i^N + B_i^N U, \quad i = 0, 1, \dots, N-1, \quad X_0^N = X_0, $$

where

$$ A_i^N = I + 0.5h(A(t_i) + A(t_{i+1})) + 0.5h^2 A(t_{i+1})^2, $$

$$ B_i^N = 0.5h(B(t_i) + B(t_{i+1})) + 0.5h^2 A(t_{i+1})B(t_{i+1}). $$

Observe that if $B = I$ and $A, B$ are constant, from (32) we obtain the scheme (26).

**THEOREM 8.2 ([71]).** Under the condition (c) there exists a constant $c$ such that for every $N \in \mathbb{N}'$,

$$ \text{haus } (X_N^N, R(T)) \le c h^2. $$

In contrast to the strongly convex case (Theorem 8.1), the proof uses essentially an effect of *nonaccumulation of errors* at each step. Namely, the error at each step may be $O(h^2)$ while the global error is $O(h^2)$ as well. Furthermore, it is shown by an example in [72] that approximation order for the solution set may be only $O(h)$.

The paper [71] contains also a negative result showing that approximations better than of second-order cannot be obtained by means of a scheme of the type (32). More precisely, consider the inclusion (31) with constant $A$ and $B$, $n = 2, r = 1$, $X_0$ a singleton, and $U$ a nondegenerate segment. Let $X_N^N$ be obtained by the following recurrent formula

$$ X_{i+1}^N = P(h; A, B)X_i^N + \sum_{j=1}^q Q_j(h; A, B)U, $$
where $q$ is a fixed integer, $P$ and $Q$ are arbitrary matrices, $h = (T - t_0)/N$.

**THEOREM 8.3 ([71]).** For any controllable pair $(A, B)$ (i.e., $\text{rank}[B, AB] = 2$), there is a constant $c > 0$ such that

$$ \text{haus} (X_N^N, R(T)) \geq ch^2 $$

for all $N \in \mathbb{N}'$.

The proof uses the observation that $R(T)$ is strongly convex and $X_N^N$ is a convex polygon. On the other hand, a strongly convex set cannot be approximated by a polygon with order higher than 2.

Concerning numerical approximations to reachable sets, various approaches can be found in the literature. A common idea in this field is to use sets of simple structure as polyhedrons [27], [30], or ellipsoids [16], [37]. A different approach, based on a time-scale decomposition of the reachable set, is presented in [21].

**9. Concluding remarks.** In this survey we concentrated on difference methods for initial value problems, especially on error estimates and order of convergence results. Concluding, we want to summarize typical difficulties, influencing the numerical performance of these algorithms decisively, and some directions of future research.

In connection with implicit methods for problems satisfying a one-sided Lipschitz condition, originally formulated for classical stiff differential equations, efficient localization procedures have to be developed for detecting possible discontinuities of derivatives of the solution. Moreover, implicit discrete inclusions have to be solved at any grid point. Therefore, to reduce the complexity of these inclusions, it would be very desirable to find diagonally implicit general linear methods with the necessary consistency and stability properties. For differential equations with discontinuous right-hand sides, D. Stewart's transformation method [64], [65] would avoid at least the solution of implicit inclusions as long as the resulting transformed differential equations are not stiff. Without piecewise transformation to smooth right-hand sides, higher-order convergence of explicit difference methods for discontinuous ordinary differential equations cannot be expected. At most first-order convergence for special problem classes and special explicit methods can be proved [40].

Concerning selection strategies for general differential inclusions, naturally the relation to numerical methods for optimal control problems should be investigated further. Especially, for nonstandard selection strategies and corresponding nonstandard objective functions, it seems to be hard to prove order of convergence results. Imposing additional state constraints immediately leads to the question of viability of solutions; compare in this connection, e.g., [58].

Approximating the whole solution set of a general differential inclusion by the whole solution set of a difference inclusion until now led only to order of convergence results for special explicit methods. These results should be incorporated into the framework of general discretization theories. This would require at least a calculus of higher-order local approximations of set-valued mappings, an actual and interesting field of research. Naturally, this will result in efficient numerical algorithms only if the set of all solutions of the difference inclusions can be approximated of any prescribed order by simpler sets, e.g., by ellipsoids or simplicial complexes. This is an actual field of research as well, stimulated by the results of A. B. Kurzhansky and his group for special control systems.

Numerical methods for boundary value problems for differential inclusions are much less developed than desirable. Naturally, good methods for initial value problems make it worthwhile to attack boundary value problems by shooting methods. The additional difficulty with differential inclusions is that the finite-dimensional systems of equations to
be solved have very poor qualitative properties, especially for differential inclusions resulting from differential equations with discontinuous right-hand sides. Hence, efficient algorithms are needed for the solution of nonsmooth systems of equations. There are some other approaches to the direct solution of boundary value problems for differential inclusions, e.g., by simplicial fixed point algorithms, cf. [59], [60], or by difference methods, cf. [44]. But until now only convergence results, no order of convergence results are available.

Above all, more computer tests are necessary to get more insight into the algorithms and into the underlying problems. We are convinced that additional numerical experiments together with a sound mathematical analysis of the results will give the right inspiration for new algorithms and new ideas.

**Acknowledgments.** This survey was strongly influenced by discussions with our collaborators V. M. Veliov, E. Farkhi, Sofia, and A. Kastner-Maresch, Bayreuth.

## REFERENCES

[1] W. ALT AND F. LEMPIO, *Programmbibliothek des Lehrstuhls für Angewandte Mathematik der Universität Bayreuth*, Bayreuth, 1988.

[2] A. A. ANDRONOV, A. A. VITT, AND S. E. KHAIKIN, *Theory of Oscillators*, Pergamon Press, Oxford, 1966.

[3] J.-P. AUBIN, *Contingent derivatives of set-valued maps and existence of solutions to nonlinear inclusions and differential inclusions*, in *Advances in Mathematics. Supplementary Studies*, L. Nachbin, ed., Academic Press, New York, 1981, pp. 160–232.

[4] J.-P. AUBIN AND A. CELLINA, *Differential Inclusions*, Springer-Verlag, Berlin, 1984.

[5] J.-P. AUBIN AND H. FRANKOWSKA, *Set-Valued Analysis*, Vol. 2 of Systems and Control: Foundations and Applications, Birkhäuser, Boston, Basel, Berlin, 1990.

[6] V. BARBU, *Optimal Control of Variational Inequalities*, Vol. 100 of Research Notes in Mathematics, Pitner Advanced Publishing Program, Boston, London, Melbourne, 1984.

[7] V. I. BLAGODATSKIKH AND A. F. FILIPPOV, *Differential inclusions and optimal control*, in *Proceedings of the Steklov Institute of Mathematics*, 4, North-Holland, Amsterdam, 1986, pp. 199–259.

[8] J. L. BUHITE AND D. R. OWEN, *An ordinary differential equation from the theory of plasticity*, Arch. Rational Mech. Anal., 71 (1979), pp. 357–383.

[9] K. BURRAGE, (*k*, *l*)-algebraic stability of Gauss methods, IMA J. Numer. Anal., 7 (1987), pp. 251–259.

[10] ——, (*k*, *l*)-algebraic stability of Runge-Kutta methods, IMA J. Numer. Anal., 8 (1988), pp. 385–400.

[11] ——, *The dichotomy of stiffness: Pragmatism versus theory*, Appl. Math. Comp., 31 (1989), pp. 92–111.

[12] K. BURRAGE AND J. C. BUTCHER, *Non-linear stability of a general class of differential equation methods*, BIT, 20 (1980), pp. 185–203.

[13] J. C. BUTCHER, *The Numerical Analysis of Ordinary Differential Equations—Runge-Kutta and General Linear Methods*, John Wiley, Chichester, New York, Brisbane, Toronto, Singapore, 1987.

[14] H. CABANNES, *Mécanique des milieux continus—propagation des discontinuités dans les cordes vibrantes soumises à un frottement solide*, C. R. Acad. Sci. Paris, Ser. B, 289 (1979), pp. 127–130.

[15] S. P. CARRINGTON-PORTER, *Numerical solution of ordinary differential equations with discontinuous right-hand sides*. M.Sc. thesis, Imperial College, London, 1980.

[16] F. L. CHERNOUSKO, *Estimation of the Phase State of Dynamical Systems*, Nauka, Moscow, 1988. (In Russian.)

[17] F. H. CLARKE, *Optimal control and the true hamiltonian*, SIAM Rev., 21 (1979), pp. 157–166.

[18] ——, *Optimization and Nonsmooth Analysis*, Wiley-Interscience, New York, 1983.

[19] P. A. COOK, *Nonlinear Dynamical Systems*, Prentice-Hall, Englewood-Cliffs, NJ, 1986.

[20] K. DEKKER AND J. G. VERWER, *Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations*, Vol. 2 of CWI Monographs, North-Holland, Amsterdam, New York, Oxford, 1984.

[21] A. DONTCHEV, *Time-scale decomposition of the reachable set of constrained linear systems*, to appear.

[22] A. L. DONTCHEV AND E. M. FARKHI, *Error estimates for discretized differential inclusions*, Computing, 41 (1989), pp. 349–358.

[23] C. M. ELLIOTT, *On the convergence of a one-step method for the numerical solution of an ordinary differential inclusion*, IMA J. Numer. Anal., 5 (1985), pp. 3–21.
[24] A. F. FILIPPOV, *Differential equations with discontinuous right-hand side*, AMS Transl., 42 (1964), pp. 199-231.

[25] ——, *Classical solutions of differential equations with multivalued right-hand side*, SIAM J. Control Optim., 5 (1967), pp. 609-621.

[26] ——, *Differential Equations with Discontinuous Righthand Side*, Mathematics and Its Applications, Kluwer Academic Publishers, Dordrecht, Boston, London, 1988.

[27] A. M. FORMALSKY, *Controllability and Stability of Systems with Restricted Resources*, Nauka, Moscow, 1974. (In Russian.)

[28] R. FRANK, J. SCHNEID, AND W. UEBERHUBER, *Order results for implicit Runge-Kutta methods applied to stiff systems*, SIAM J. Numer. Anal., 22 (1985), pp. 515-534.

[29] ——, *Stability properties of implicit Runge-Kutta methods*, SIAM J. Numer. Anal., 22 (1985), pp. 497-514.

[30] B. H. GRAETTINGER AND B. H. KROGH, *Hyperplane method for reachable state estimation for linear time-invariant systems*, J. Optim. Theory Appl., 68 (1991), pp. 555-588.

[31] H. KAHL, *Direkte Methoden für optimale Steuerungsprobleme*, Diplomarbeit, Universität Bayreuth, Bayreuth, 1986.

[32] A. KASTNER-MARESCH, *Differenzenverfahren höherer Ordnung für Differentialinklusionen*, Ph.D. dissertation, Universität Bayreuth, 1990.

[33] ——, *Implicit Runge-Kutta methods for differential inclusions*, Numer. Funct. Anal. Optim., 11 (1991), pp. 937-958.

[34] A. KASTNER-MARESCH AND F. LEMPIO, *Difference methods with selection strategies for differential inclusions*, preprint, 1991.

[35] J. F. B. M. KRAAIJEVANGER, *B-convergence of the implicit midpoint rule and the trapezoidal rule*, BIT, 25 (1985), pp. 652-666.

[36] N. N. KRASOVSKI, *The Control of Dynamic Systems*, Nauka, Moscow, 1986. (In Russian.)

[37] A. B. KURZHANSKI AND I. VALYI, *The problem of control synthesis for uncertain systems: Ellipsoidal techniques*, in Modelling, Estimation and Control of Systems with Uncertainty, G. B. Di Masi, A. Gombani, and A. B. Kurzhansky, eds., Birkhäuser, Boston, Basel, Berlin, 1991, pp. 260-282.

[38] A. B. KURZHANSKY AND T. F. FILIPPOVA, *Description of the bundle of control-system trajectories remaining in a given set for a given time*, Differentsial'nye Uravneniya, 23 (1987), pp. 1303-1315. (In Russian.)

[39] F. LEMPIO, *The numerical treatment of differential inclusions*, in Proceedings of the Fourth Conference on Differential Equations and Applications, Russe, Bulgaria, 1989.

[40] ——, *Modified Euler methods for differential inclusions*, in Proceedings of the Workshop on Set-Valued Analysis and Differential Inclusions, Pamporovo, Bulgaria, 1990.

[41] R. MANNSHARDT, *Simulation of discontinuous systems by use of Runge-Kutta-methods combined with Newton iteration*, in Proceedings of the Sympos. SIMULATION, Montreux, 1977, pp. 163-167.

[42] ——, *One-step methods of any order for ordinary differential equations with discontinuous right-hand sides*, Numer. Math., 31 (1978), pp. 131-152.

[43] R. MODEL, *Zur Integration über Unstetigkeiten in gewöhnlichen Differentialgleichungen*, Z. Angew. Math. Mech., 68 (1988), pp. 161-169.

[44] H.-D. NIEPAGE, *Inverse stability and convergence of difference approximations for boundary value problems for differential inclusions*, Numer. Funct. Anal. Optim., 9 (1987), pp. 761-778.

[45] H.-D. NIEPAGE AND W. WENDT, *On the discrete convergence of multistep methods for differential inclusions*, Numer. Funct. Anal. Optim., 9 (1987), pp. 591-617.

[46] M. S. NIKOLSKIY, *On approximation of the attainability domain of control process*, Mat. Zametki., 41 (1987), pp. 71-76. (In Russian.)

[47] ——, *On the attainability domain of differential inclusions*, Vestnik Mosk. University, Vych. Mat. Kibern., 4 (1987), pp. 31-34. (In Russian.)

[48] ——, *On a method for approximation of attainable set for a differential inclusion*, Journ. of Vych. Mat. i Math. Phys., 28 (1988), pp. 1252-1254. (In Russian.)

[49] A. I. PANASYUK, *Equations of attainable set dynamics, Part 1: Integral funnel equations*, J. Optim. Theory Appl., 64 (1990), pp. 349-366.

[50] A. I. PANASYUK AND V. I. PANASYUK, *On one equation resulting from a differential inclusion*, Mat. Zametki, 27 (1980), pp. 429-445. (In Russian.)

[51] W. M. PICKERING AND D. M. BURLEY, *The oscillations of a simple castor*, Bull. IMA, 13 (1977), pp. 47-50.

[52] B. T. POLYAK, *Introduction to Optimization*, Translations Series in Mathematics and Engineering, Optimization Software, Publications Division, NY, 1987.

[53] B. N. PSHENICHNY, *Convex Analysis and Optimization*, Nauka, Moscow, 1980. (In Russian.)
[54] M. A. RAUPP AND O. R. BAIOCCHI, *On the corona effect: a lumped circuit model*, Math. Appl. Comp., 1 (1982), pp. 239-265.

[55] R. REISSIG, *Erzwungene Schwingungen mit zäher und trockener Reibung*, Math. Nachr., 11 (1954), pp. 345-384.

[56] J. C. W. ROGERS, *A free boundary problem as diffusion with nonlinear absorption*, J. Inst. Math. Appl., 20 (1977), pp. 261-268.

[57] E. N. RURENKO, *On the problem of approximate computation of the attainable set for a differential inclusion*, Vesn. Mosk. University, 2 (1989), pp. 81-83. (In Russian.)

[58] P. SAINT-PIERRE, *Approximation of slow solutions to differential inclusions*, Appl. Math. Optim., 22 (1990), pp. 311-330.

[59] K. SCHILLING, *Simpliziale Algorithmen zur Berechnung von Fixpunkten mengenwertiger Operatoren*, Wis-senschaftlicher Verlag Trier, Trier, 1986.

[60] ——, *An algorithm to solve boundary value problems for differential inclusions and applications in optimal control*, Numer. Funct. Anal. Optim., 10 (1989), pp. 733-764.

[61] K. SCHITTKOWSKI, *NLPQL: A FORTRAN subroutine solving constrained nonlinear programming problems*, Ann. Oper. Res., 5 (1985), pp. 485-500.

[62] T. SPEER, *Kontinuierliche Runge-Kutta-Approximationen für steife Differentialgleichungen mit Anwendung in der Fahrzeugdynamik*. Schwerpunktprogramm der DFG: Anwendungsbezogene Optimierung und Steuerung, TU München, Report No. 147, 1989.

[63] H. J. STETTER, *Analysis of Discretization Methods for Ordinary Differential Equations*, Vol. 23 of Springer Tracts in Natural Philosophy, Springer-Verlag, New York, 1973.

[64] D. STEWART, *High Accuracy Numerical Methods for Ordinary Differential Equations with Discontinuous Right-Hand Side*, Ph.D. thesis, University of Queensland, St. Lucia, Australia, 1989.

[65] ——, *A high accuracy method for solving ODEs with discontinuous right-hand side*, Numer. Math., 58 (1990), pp. 299-328.

[66] K. TAUBERT, *Über die Approximation nicht-klassischer Lösungen von Anfangswertaufgaben*, Ph.D. dissertation, Universität Hamburg, 1973.

[67] ——, *Differenzenverfahren für Schwingungen mit trockener und zäher Reibung und für Regelungssysteme*, Numer. Math., 26 (1976), pp. 379-395.

[68] ——, *Converging multistep methods for initial value problems involving multivalued maps*, Computing, 27 (1981), pp. 123-136.

[69] A. A. TOLSTONOGOV, *On the integral funnel equation of differential inclusions*, Mat. Zametki, 32 (1982), pp. 847-852. (In Russian.)

[70] V. I. UTKIN, *Sliding Regimes in Optimization and Control Problems*, Nauka, Moscow, 1981. (In Russian.)

[71] V. M. VELIOV, *Second order discrete approximations to linear differential inclusions*. SIAM J. Numer. Anal., 29 (1992), to appear.

[72] ——, *Approximations to differential inclusions by discrete inclusions*. IIASA Working Paper WP-89-017, 1989.

[73] ——, *Second order discrete approximations to strongly convex differential inclusions*, Systems and Control Lett., 13 (1989), pp. 263-269.

[74] P. R. WOLENSKI, *The exponential formula for the reachable set of a Lipschitz differential inclusion*, SIAM J. Control Optim., 28 (1990), pp. 1148-1161.