Difference methods for differential inclusions: a survey

. 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 strate gies 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.

1. Introduction. In this survey we consider the following initial value problem for ordinary differential inclusions. Find an absolutely continuous function y( ·) on I such that (1) y (to) = Yo and iJ ( t) E F ( t, y ( t)) for almost all t E I, where iJ ( ·) is the derivative of y ( · ).
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 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 y(t) E n n cl (conv (f(t, {z ERn: liz-y(t)ll < 8} \ N))) for almost all t E I, where J-l denotes Lebesgue measure on JRn. 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 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 iJ(t) E -otp(y(t) ), where tp : Rn ----+ R is a convex potential function with subdifferential /J<p. Such inclusions have an important property: their equilibrium solutions minimize <p. Moreover, if <p achieves a minimum at all, then for every initial value y 0 E Rn, the corresponding solution y( t) as t ----+ oo converges to a minimizer of tp. 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

y(t), u(t))
with feasible controls for almost all t E J, where f is single-valued and U(t, x) c JRm for all t E I and x ERn, just by taking 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 areduction 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 timeT 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.
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(·) EX.
The first approach uses a finite difference scheme together with suitable selection procedures resulting in a sequence of grid functions say, on a uniform grid where, as usual, N' denotes a subsequence of N converging to infinity. Naturally, the question arises whether at least a subsequence of the sequence of, say, piecewise linear continuous interpolants of the grid functions converges to a solution y( ·) E 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(·) E 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 r;esult, for each h we obtain a set Xh 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 Xh to X ash -t 0. The error can be measured by a suitably defined distance between the sets Xh and X.
As an introductory example, in §2 we consider the classical Euler method. Following the first approach we prove directly that ifF is compact convex valued upper semicontinuous and with linear growth, then where lim sup is in the usual (Kuratowski) sense, i.e., limsupXh = {! E C(It: liminfdist(f,Xh) = o}' h~o h~o 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 "'N ( ·) E Xh of the difference inclusion with sufficiently small stepsize contains in its €-neighbourhood in C ( J) n a solution y( ·) E 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 where haus( ·, ·) 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.
Lett be a point in I. The "reachable set at the timet" of (1) is defined as That is, R(t) is the set of all points x that are ends oftrajectories of (1) on [to, t]. Section 7 shows that a sequence of sets can be determined from the Euler scheme, which is Hausdorff convergent to R(t) ash---+ 0, uniformly in t E J. 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].

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.
Let 'flo = y 0 , and for j = 0, · · ·, N -1 compute "li+l from (3) 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 such that where ('fJo, · · ·, 'fJN) is any grid function satisfying (3). Let Xh be the set of all solutions of (3) for given h.
The following theorem is present explicitely or implicitely in many works, in various forms, see, e.g., J.-P. Aubin  (ii) There exist constants k and a, such that l l z11 < kllxll +a Then every sequence (17N(·))NEN' with 'fJN(·) E Xhfor N E N' has a subsequence which converges as N--+ oo, uniformly in I, to some solution of the problem (1). In other words, (2) is fulfilled.

N-+oo tEl
By Arzela's theorem, TJN ( ·) has a subsequence uniformly convergent to some function y(·), moreover, r,N (·)--+ z(·) weakly in L 00 (I)"" for some further subsequence. It is easy to observe that for all t E I y(t) =Yo+ lt z(s)ds, to i.e., y(-) is absolutely continuous and y(t) = z(t) for almost every t E I.
then by (i) there exists N 1 such that for N > N 1 where B is the unit ball. Pick some q E R"" and let ~be a measurable subset of I. Then from (7) where* denotes transposition and supp is the support function. By (i) the map t 1---4 supp(q, F(t, ryN (t))) is upper semi-continuous, and from (ii) and (5) it is bounded above by a constant, independent of N. Then, for N--+ oo for the corresponding subsequence, the dominated convergence theorem yields i q*y(r) dr < i [supp(q, F(r, y(r))) + tllqll] dr.
Since ~ and t are arbitrary and F is convex valued, this means that y( ·) is a solution of (1). D Under some more conditions we are able to prove that every solution ryN ( ·) of the discretized inclusion (3) with sufficiently small stepsize contains in its €-neighbourhood in C(I)n a solution of (1). For that purpose we need the following basic result ofthe theory of differential inclusions, often referred to as Gronwall-Filippov-Wai:ewski theorem; see A. F. Filippov [25] or J.-P. Aubin   haus(F(t,x),F(t,z)) < k(t)llxzll for all (t, x) and (t, z)from Q with k(·) E L1(J).
Assume, moreover, that dist(y(t), F(t, y(t))) < p(t) for almost all t E I, for some p( ·) E L1 (I) such that Then there exists a solution y ( ·) to the initial value problem ( 1) such that for all t E I. THEOREM 2.4. Suppose that the conditions of Theorem 2.2 hold and, moreover, that F is continuous on I x lRn and Lipschitz in x on bounded sets in IRn with integrable in I Lipschitz constant. Then for every ~ > 0 there exists N 1 such that for every N > N 1 and for Proof. Let € > 0. As in the proof of Theorem 2.2 we show that the set Xh of the discrete trajectories is bounded in C(J)n and, moreover, that there exists N 1 such that for all N > Nb for al117N ( ·) E Xh and for all t with tJ < t < tJ+b j = 0, 1, · · · , N -1. Hence for all N > N 1 every solution 1JN (-) E Xh will satisfy (7). Applying Gronwall-Filippov-Wazewski theorem we obtain that there exists a solution y( ·) of (1) in the €-tube around 1JN (·). 0 In §6, slightly strengthening the assumptions, we obtain an () (h) estimate for the Hausdorff distance between the sets Xh and the set of solutions X of (1).

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.
computed, e.g., by a linear f-step method with f < r or by a 1-step method~ For j = r, · · · , N compute T/j from These methods are implicit in case br =f. 0. Generally, the solution of the inclusion 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. This implies, that the family ('TIN ( ·)) N EN' 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: ao=-1, we get again the Euler method.
The following selections of coefficients were tested in [31]: 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].
for almost all t E I and 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( ·) : I ---+ JR 2 such that where the Sgn-function 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(·) and its derivative i(·) in Fig. 1 and the phase portrait shown in Fig. 2. Subintervals, where the solution is in fact a generalized one, are most interesting. The plot in Fig. 3 of the approximation of the derivative i( ·) by Euler method with stepsize 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 z ( ·) 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 singlevalued and smooth.

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 b1 b8 with nonnegative bi (i = 1, · · · ,s), a stepsize h = (Tt 0 )/N with N E N', and an approximation TJo of the initial value y 0 • For j = 0, · · · , N -1 solve the implicit system of inclusions 8 r;p.
and compute the next approximation 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 with a single-valued function f : I x Rn -+ Rn with (one-sided) Lipschitz constant L 1 and a monotone set-valued mapping (3 : Rn =;.. JR.n. For example, the subdifferential 81.p of a convex functional <p : JR.n -+ 1R. defines a monotone mapping on JR.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. hold.
(iv) The method is BS-stable and C-stable.
( v) The initial approximations TJf/ of Yo satisfy 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. Then the order of convergence is equal to J.L This corollary serves as a basis for more or less sophisticated algorithms with builtin localization procedures for detecting possible discontinuities of the derivative of the exact solution. In practice, we have to localize numerically manifolds where the righthand 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 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 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. For z = O!Rn, in each step we get the selection with minimal norm. In any case, the explicit difference method (br = 0) for the differential inclusion can be regarded as linear multistep method for the differential equation where pr F( t,x) ( z) is the projection of z onto F( t, x). Hence, the properties off ( t, x) determine the properties of the algorithm. The calculation of starting values 'f}o, · · ·, 'T/r-1 and starting selections ( 0 , · · · , (r-1 has to be adjusted accordingly. Convergence is proved by Convergence Theorem 3.2. Using arguments in [25] we get the following. THEOREM

In addition to the assumptions in Convergence Theorem 3.2, let F be
Hausdorff continuous on I x lRn and br = 0. Then the sequence 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. 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.  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

Then the sequence 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. 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 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 worthwile to extend the above qualitative results byestimations 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.

TJo =Yo such that
In the proof, B. N. Pshenichny used the following construction of the discrete trajectory TJN (·): Let TJf: = y 0 and define TJJ' + 1 as the projection of y (ti+d on the set rJf + hF ( 1Jf). Then from (10), from the inclusion and from the choice of TJf-t-1 , using the Lipschitz continuity ofF, we obtain for some a 1 and a 2 that for j = 0, 1, · · ·, N -1. This implies the desired estimate.
Combining the result of B. N. Pshenichny with Gronwall-Filippov-Wazewski theorem, we obtain that (12) m.ax IITJ.f -fj(tj)ll <c(llii(to)-Yoii+1T dist(y(t),F(t,fj(t))) dt+h) 05.J5.N to for some c independent of h. This estimate was obtained in [22] by a direct proof. More precisely, the following theorem is proved. F is nonempty compact and convex valued; there exist k and a such that II f II < k llx II+ a whenever f E F(t, x), t E I, x E Rn; F is Lipschitz in x on bounded sets unifonnly in t; F is of bounded variation in t unifonnly in x on bounded sets, i.e., for any bounded set Then there exists a constant c such that if "'N ( ·) is obtained by (11 ), then for all N (12) is fulfilled.
In particular, if y is a solution of Initial Value Problem 1.1 then "'N ( ·) provides a first-order approximation to fi. 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 "'N ( ·) be a piecewise linear function that solves (3). IfF is Lipschitz in both x and t then i]N (t) E F(t, "'N (t)) + hkB for a.e. t E (to, T] for some k and for all N EN'. Then one can apply Gronwall-Filippov-Wazewski theorem, obtaining that there exists a solution fiN of Initial Value Problem 1.1 such that (13) In fact, to obtain (13) it is sufficient to assume the conditions of Theorem 6.2. Namely, we have 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 .
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, br = 0, and the following conditions hold:

is Lipschitz in x E JRn uniformly in t E I and of bounded variation in t E I uniformly in x on bounded sets.
(ii) the coefficients ai are nonnegative except the single coefficient a 0 ; or alternatively, exactly two of the coefficients ai ( i = 0, 1, · · · , r) are nonzero. Then there exists a constant c such that for all sufficiently small h, if "'N ( ·) 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, · · ·, N}. Then the reachable set, associated with (3) at the time ti, will be THEOREM 7.1. Let F be a compact convex and nonempty valued map that is continuous in both t and x in I x Rn and let the reachable set R(t) of Initial Value Problem 1.1 be contained in a bounded set K in Rn for all t E I. Let F be Lipschitz in a neighbourhood of K with respect to x unifonnly in t E I. Then The proof of this theorem can be easily obtained from [50], for a more recent presentation see [49]. First we need some notation. DEFINITION   In [69] extensions of Theorem 7.3 for inclusions with right-hand side satisfying Caratheodory conditions are obtained. A funnel equation to a linear differential inclusion with phase constraints is found in [38]. Let  where cp(h) ---+ 0 ash---+ 0 uniformly in t E /.If k is the Lipschitz constant ofF, it is easy to see from (15) that (17) haus ( U (x + hF(t1, x)), U (x + hF(t1, x))) < (1 + kh)hcp(h).
xERf" xER(ti) Combining (16) and (17) where the power of (I+ hF) is that of composition of set-valued maps. Then RIJ will converge to R(T) if and only if (ii) If, in addition, F is assumed to have convex values, then for all T > 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 CJ( v'ii) only. The paper [46] contains a similar result for a control system. In [ 48] the following modification of Euler scheme is proposed: where the right-hand side F, independent of t, satisfies a growth condition and is locally Lipschitz and the function r and the constants a, a 1 are chosen in such a way that cl(R(T)) c E}$ n U. Then from (13) an estimate of order O(h) follows for the distance between E}J n 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) where A is a nxn-matrix, U is a compact and convex set in ntn. 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) E U for almost every t E I such that (19) 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 1JN toy 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).
(i) Let 1JN = ( 11~, · · · , 11~) be a solution of (20). Then, from (21) and (22) 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 O(hk). The condition (21) will hold if we take It turns out that the condition (22) is much more restrictive. It obviously holds for k = 1 (Euler scheme); moreover, we have always provided that (25) This means that (24) is satisfied, i.e., every solution of (18) can be approximated of order O(hk) 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 J.1, > 0 such that u 1 , u 2 E U implies 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 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) y(t) E F(t, y(t)) for a.e. t E I, y(to) E Ro approximated by fori= 0, 1, · · ·, 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~ x S c JRn+l and satisfies the following conditions: (a) F is compact valued and there exists J.1, > 0 such that F(t, x) is strongly convex foreveryt E ~andx E S. Thesupportfunctionr(l,t,x) = supp(l,F(t,x))isdifferentiable with respect to t and x, 8-r the set R 0 is compact, I c int(~), and all trajectories of (27) stay in int(S).
(b) The function y( t, x, l) defined as (29) l*y(t, x, l) = supp(l, F(t, x)) is  (27) there is a discrete solution 1JN = ( 1J/i, · · ·, 17~) of (28) such that (30) (ii) For every solution 1JN of (28) there exists a solution y of (27) such that (30) holds. If in addition (b) is satisfied, then h 3 1 2 in ( 30) can be replaced by h 2 • In [73] it is also shown that if the map F is 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 E ~, x E S. In this case the discrete  <p(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 r y(t) E f(t,y(t)) + l::Yi(t,y(t))[O, 1], y(to) E Ro .
i=I 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) where A(t) and B(t) are nxn and nxr matrices, U c lRr, 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 timeT. Introducing an uniform grid in I, for each N we define the set X% recurrently by the equation In contrast to the strongly convex case (Theorem 8.1), the proof uses essentially an effect of nonaccumulation of e"ors at each step. Namely, the error at each step may be tJ(h 2 ) while the global error is tJ(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 ser.:ond-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% be obtained by the following recurrent formula q xi~l = P(h;A,B)xr + LQj(h;A,B)U, j=l where q is a fixed integer, P and Q are arbitrary matrices, h = (Tt 0 )jN. The proof uses the observation that R(T) is strongly convex and X~ 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 timescale 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.