Asen Dontchevt

Frank Lempio+

DIFFERENCE METHODS FOR DIFFERENTIAL INCLUSIONS: A SURVEY*

Keywords:

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.

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 inte1Val, Yo E lRn, and F be a map from I x ~n into the set of all subsets of~n.

Find an absolutely continuous function y( •) on I such that [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF] 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 [START_REF] Fiuppov | Differential equations with discontinuous right-hand side[END_REF], 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)))

8>0 ~-t(N)=O
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 [START_REF] Reissig | Erzwungene Schwingungen mit ziiher und trockener Reibung[END_REF], [START_REF] Andronov | Theory of Oscillators[END_REF], [START_REF]Differenzenverfahren for Schwingungen mit trockener und ziiher Reibung und for Regelungssysteme[END_REF], [START_REF] Pickering | The oscillations of a simple castor[END_REF], especially in vehicle dynamics for the description of locking phases during brake maneuvers [START_REF] Speer | Kontinuierliche Runge-Kutta-Approximationen filr steife Differentialgleichungen mit Anwendung in der Fahrzeugdynamik[END_REF], in elasto-plasticity [START_REF] Buhite | An ordinary differential equation from the theory of plasticity[END_REF], [15], variable structure systems [START_REF] Utkin | Sliding Regimes in Optimization and Control Problems[END_REF], and electrical circuits [START_REF] Raupp And 0 | On the corona effect: a lumped circuit model[END_REF], [START_REF] Cook | Nonlinear Dynamical Systems[END_REF]. 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., [START_REF] Rogers | A free boundary problem as diffusion with nonlinear absorption[END_REF], [START_REF] Cabannes | Mecanique des milieux continus-propagation des discontinuites dans les cordes vibrantes soumises a unfrottement solide[END_REF], [START_REF] Aubin | Set-Valued Analysis[END_REF]. In § §3 and 4 we use Example 3.3, cf. [START_REF] Reissig | Erzwungene Schwingungen mit ziiher und trockener Reibung[END_REF], [START_REF]Differenzenverfahren for Schwingungen mit trockener und ziiher Reibung und for Regelungssysteme[END_REF], 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. [START_REF] Aubin | Contingent derivatives of set-valued maps and existence of solutions to nonlinear inclusions and differential inclusions[END_REF], [START_REF] Aubin | Differential Inclusions[END_REF], [START_REF] Polyak | Introduction to Optimization[END_REF].

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) = j(t, y(t), u(t))

with feasible controls u(t) E U(t, y(t))

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 F(t, x) = {j(t, x, z) : z E U(t, x)}.

In §5, Example 5.3, we describe a special differential inclusion of this type, which was originally used in [START_REF] Kastner | Difference methods with selection strategies for differential inclusions[END_REF] 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 [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF]. 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. [START_REF] Clarke | Optimal control and the trne hamiltonian[END_REF] and [START_REF]Optimization and NonsmoothAnalysis[END_REF], Chapter 3, where a nice example of a nonsmooth problem in resource economics is presented, and [START_REF] Schilling | Simpliziale Algorithmen zur Berechnung von Fixpunkten mengenwerliger Operatoren[END_REF], [START_REF]An algorithm to solve boundary value problems for differential inclusions and applications in optimal control[END_REF], 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 [START_REF] Krasovski | The Control of Dynamic Systems[END_REF] for control synthesis of uncertain systems.

A number of other motivations and applications of differential inclusions can be found in [START_REF] Blagodatskikh | Differential inclusions and optimal control[END_REF] and in the books [START_REF] Aubin | Differential Inclusions[END_REF], [16], [START_REF]Optimization and NonsmoothAnalysis[END_REF], [START_REF]Differential Equations with Discontinuous Righthand Side, Mathematics and Its Applications[END_REF], [START_REF] Formalsky | Controllability and Stability of Systems with Restricted Ressources[END_REF], [START_REF] Pshenichny | Convex Analysis and Optimization[END_REF], 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 [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF]. 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

( N N N) "lo '"11 ' • • • '"'N (N EN'),
say, on a uniform grid

N N tN T to < t1 < ... < N = with stepsize T-to N N h= =t. -t. 1 N J J - (j = 1, • • •, N) ,
where, as usual, N' denotes a subsequence of N converging to infinity. Naturally, the question arises whether at least a subsequence of the sequence

( 'f]N (.)) N EN'
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 [START_REF] Aubeiit | Ober die Approximation nicht-klassischer Losungen von Anfangswerlaufgaben[END_REF], [START_REF]Converging multistep methods for initial value problems involving multivalued maps[END_REF].

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 [START_REF] Ellioti | On the convergence of a one-step method for the numerical solution of an ordinary differential inclusion[END_REF] 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 [START_REF] Kas | 1NER-MARESCH, Differenzenverfahren hOherer Ordnung fii.r Differentialinklusionen[END_REF], [START_REF]Implicit Runge-Kutta methods for differential inclusions[END_REF]. 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 [START_REF] Stewart | High Accuracy Numerical Methods for Ordinary Differential Equations with Discontinuous Right-Hand Side[END_REF], [START_REF]A high accuracy method for solving ODEs with discontinuous right-hand side[END_REF] 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

(2) limsupXh C X in C(It, h--+0
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 [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF].

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

R(t) = {x ERn: x = y(t) for some y(•) EX}.
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 [START_REF] Nikoi | SKIY, On approximation of the attainability domain of control process[END_REF]- [START_REF]On a method for approximation of attainable set for a differential inclusion[END_REF] 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 [START_REF] Kastner | Difference methods with selection strategies for differential inclusions[END_REF] and [START_REF] Kas | 1NER-MARESCH, Differenzenverfahren hOherer Ordnung fii.r Differentialinklusionen[END_REF].

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 EN' C N choose a grid to < fi < • • • < t N = T with stepsize T -to h = N = t1 -t 1 _1 (j = 1, • • •, N) .
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

N 1 1] (t) = 'f}j + h (t-tj )('fJj+l -'f}j) J . -0 • • • N -1) -'
' '

where ('fJo, then by (i) there exists N 1 such that for N > N 1 [START_REF] Blagodatskikh | Differential inclusions and optimal control[END_REF] i]N (t) E F(t, 1JN (t)) + tB for a.e. t E I, where B is the unit ball. Pick some q E R"" and let ~be a measurable subset of I. Then from [START_REF] Blagodatskikh | Differential inclusions and optimal control[END_REF] 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 [START_REF]Classical solutions of differential equations with multivalued right-hand side[END_REF] or J.-P. Aubin and A. Cellina [4, Thm. 2.4.1, p. 120], for a more general formulation. THEOREM 2.3 (Gronwall-Filippov-Wai:ewski). Let y : I --+ JRn be an absolutely continuous function with f) ( t 0 ) = y 0 , f3 be a positive constant, and

Q = {(t, x) E I X Rn : llx-y(t)ll < /3}.
Let F: Q:::} Rn be nonempty closed valued and continuous and satisfy

haus(F(t,x),F(t,z)) < k(t)llx-zll 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 Wl = 1.: exp [J.' k(r)dr] p(s)ds < (3 (t E I).
Then there exists a solution y ( •) to the initial value problem ( 1) such that lly(t)-y(t)ll < e(t) 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 every solution 1JN ( •) of the discrete inclusion (3) there exists a solution y( •) of the problem [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF]

such that mruc II11N (t) -y(t) II < t. tEl 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 [START_REF] Blagodatskikh | Differential inclusions and optimal control[END_REF]. 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 [START_REF] Aubeiit | Ober die Approximation nicht-klassischer Losungen von Anfangswerlaufgaben[END_REF], [START_REF]Converging multistep methods for initial value problems involving multivalued maps[END_REF]. These methods are typically of the following form, where we deliberately suppress possible generalizations to multistage multistep methods.

(i = 0, .. •, r) (N EN' c N).
(j = 0, .. •, r -1)

(j = 0, • • •, r -1),
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., [START_REF] Schritkowski | NLPQL: A FORTRAN subroutine solving constrained nonlinear programming problems[END_REF]. 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 x lRn.

(iii) The strong root condition is satisfied, i.e., all zeros>.. of the polynomial have absolute value I >..I < 1 except the simple zero >.. = 1.

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

( i = 0, • • • , r ). (vi)
The starting values satisfy

(j = 0, .. •, r -2)
for all N EN' with a constant M, which is independent of the stepsize h = (Tt 0 )jN.

(vii) The approximations of the initial value y 0 satisfy lim N _. oo TJC' = y 0 .

NEN'

Then the sequence of piecewise linear continuous interpolants of the grid functions

( N N N) TJo , T/1 , • • • , TIN
contains a subsequence which converges uniformly to a solution of the initial value problem [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF]. 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. 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 [START_REF] Niepage | On the discrete convergence of multistep methods for differential inclusions[END_REF], [START_REF] Niepage | Inverse stability and convergence of difference approximations for boundary value problems for differential inclusions[END_REF].

For the special selection of coefficients r=1: ao=-1,

a1 = 1,
we get again the Euler method.

a2 = 0, b2 = 0.1, a3 = -1, b3 = 0.4, a4 = 1, b4 = 0.
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 [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF]. We conclude this section with the following simple differential equation with discontinuous right-hand side describing forced vibrations with viscous and Coulomb damping, cf. [START_REF] Reissig | Erzwungene Schwingungen mit ziiher und trockener Reibung[END_REF], [START_REF]Differenzenverfahren for Schwingungen mit trockener und ziiher Reibung und for Regelungssysteme[END_REF].

z(t) + z(t) + 0.2z(t) + 4sgn(z(t)) = 2cos7rt (8)
for almost all t E I and

z(O) = 3, z(O) = 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( •) :

I ---+ JR 2 such that YI ( t) Y2 ( t) ' iJ2(t) E -y1(t)-0.2y2(t)-4Sgn(y2(t)) + 2cos7rt
for almost all t E I= [0, T] and

YI(O) = 3, Y2(0) = 4,
where the Sgn-function -1

(x2 < 0), [-1, 1] (x2 = 0), 1 (x2 > 0),
is the set-valued analogue of the usual sign function used in [START_REF] Buhite | An ordinary differential equation from the theory of plasticity[END_REF] 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.

~----------------, It) ~~----~----~----~--~ o.o
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 [START_REF] Mannshardt | Simulation of discontinuous systems by use of Runge-Kutta-methods combined with Newton iteration[END_REF], [START_REF]One-step methods of any order for ordinary differential equations with discontinuous right-hand sides[END_REF] 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 [START_REF] Model | Zur Integration uber Unstetigkeiten in gewohnlichen Differentialgleichungen[END_REF]. D. Stewart has defended his Ph.D. thesis on the same problem class in 1990; compare [START_REF] Stewart | High Accuracy Numerical Methods for Ordinary Differential Equations with Discontinuous Right-Hand Side[END_REF]. 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 [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF]. Stimulated by a paper of C. M. Elliott [START_REF] Ellioti | On the convergence of a one-step method for the numerical solution of an ordinary differential inclusion[END_REF] 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 [START_REF] Kas | 1NER-MARESCH, Differenzenverfahren hOherer Ordnung fii.r Differentialinklusionen[END_REF] and [START_REF]Implicit Runge-Kutta methods for differential inclusions[END_REF]. In principle, in the following we could investigate the class of general linear methods; compare [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations-Runge-Kutta and General Linear Methods[END_REF]. But for simplicity of presentation we restrict ourselves to the class of implicit s-stage Runge-Kutta methods.

0 0 0 0 0 8 0 I g Y~------~--------~------~ 00
IMPLICIT s-STAGE RUNGE-KUTTA METHOD 4.1. Let there be given the scheme of coefficients and compute the next approximation 8

(J-l, v = 1' ... ' s) TJj+l = 7}j + h L bv (jv• v=l
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.

F(t, x) = f(t, x)-(3(x)
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 [START_REF] Ellioti | On the convergence of a one-step method for the numerical solution of an ordinary differential inclusion[END_REF] 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 [START_REF] Verwer | Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations[END_REF], [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations-Runge-Kutta and General Linear Methods[END_REF] and the journal articles [START_REF] Butcher | Non-linear stability of a general class of differential equation methods[END_REF], [START_REF]Stability properties of implicit Runge-Kutta methods[END_REF], [START_REF] Frank | Order results for implicit Runge-Kutta methods applied to stiff systems[END_REF], [START_REF] Burrage | k, l)-algebraic stability of Gauss methods[END_REF], [START_REF]k, l)-algebraic stability of Runge-Kutta methods[END_REF], especially the survey of K. Burrage [START_REF]The dichotomy of stiffness: Pragmatism versus theory[END_REF]. 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 [START_REF] Kas | 1NER-MARESCH, Differenzenverfahren hOherer Ordnung fii.r Differentialinklusionen[END_REF] reads as the following. THEOREM 4.3. Let the following assumptions be satisfied: (i) Let F be upper semicontinuous on I x JR.n and satisfy a one-sided Lipschitz condition.

(ii) The (necessarily unique) solution of initial value problem [START_REF] Lempio+ | Programmbibliothek des Lehrstuhls fUr Angewandte Mathematik der Universitiit Bayreuth[END_REF] is piecewise (J.t + 1)times continuously differentiable with J.l > 1.

(iii) The simplifying conditions B(J..t), i.e., and C(J.L), i.e.,

(k = 1, ... 'J..t), ( k-1 • • • II.' i = 1 • • • s) - ' ',..,, ' ' ' hold. (iv)
The method is BS-stable and C-stable.

( v) The initial approximations TJf/ of Yo satisfy

!ITJ~-Yo !I = 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. 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 [START_REF] Mannshardt | Simulation of discontinuous systems by use of Runge-Kutta-methods combined with Newton iteration[END_REF], [START_REF]One-step methods of any order for ordinary differential equations with discontinuous right-hand sides[END_REF], or [START_REF] Stewart | High Accuracy Numerical Methods for Ordinary Differential Equations with Discontinuous Right-Hand Side[END_REF] 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 For this method the simplifying conditions B(l), B(2), and C(l) (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 [START_REF] Ellioti | On the convergence of a one-step method for the numerical solution of an ordinary differential inclusion[END_REF]. 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 [START_REF] Kraaijevanger | B-convergence of the implicit midpoint rule and the trapezoidal rule[END_REF]. 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. 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

y(t) = f (t, y(t)) with f(t,x) = prF(t,x)(z)
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 [START_REF]Classical solutions of differential equations with multivalued right-hand side[END_REF] we get the following. THEOREM 5.2. 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 [START_REF] Kastner | Difference methods with selection strategies for differential inclusions[END_REF] is treated numerically as an example for the various selection strategies of this section.

Example 5.3. Find an absolutely continuous function y :

I --+ JR 2 such that -y2(t) + O.lyl(t)(9-Y1(t) 2 -Y2(t) 2 ) + [-1, 1), YI(t) + 0.1y2(t)(9-YI(t) 2 -Y2(t?) + [-1, 1]
for almost all t E I = [0, T] and YI(0)=5, Y2(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. 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 [START_REF] Aubin | Differential Inclusions[END_REF]. -----------------------~ ~,------------------------- If, in addition to the assumptions in the Convergence Theorem 3.2, the set-valued mapping F is Lipschitz continuous in I x JRn; then a similar argument as in the proof of Theorem 3.2 applies to the selections Vj, and we get for explicit methods (br = 0).

THEOREM 5.5. In addition to the assumptions in Convergence Theorem 3.2, let F be Lipschitz continuous on I x JRn, br = 0, and the starting selections satisfy llvf-n -vJ' II < h M (j = 0, • .. , r -2)

for all N EN'.

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. 12.0

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 0 ;, ------------------------- For j = r, • • • , N, minimize subject to l r r h 2..:.:: ai T/j-r+i = 2..:.: Most interesting is the comparison of the chosen selections; compare Fig. 15 with the results of the other selection strategies.

: bi (j-r+i , i=O i=O (j-r+i E F (tj-r+i' 'f/j-r+i) (i=O,•••,r).
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.

The first result in this direction probably is due to B. N. Pshenichny [START_REF] Pshenichny | Convex Analysis and Optimization[END_REF], who proved the following. THEOREM 6.1 ( [START_REF] Pshenichny | Convex Analysis and Optimization[END_REF]). Let F : Rn ~ IRn be a compact convex valued mapping andy be a solution of y(t) E F(y(t)) for a.e. t E I, y(to) =YO• Suppose that there exists£ > 0 such that F is Lipschitz in the set {x E lRn : llxy(t)ll < E for some t E I}. Then there exist c and N 1 such that for all N > N 1 there exists a solution 'fJN ( •) of the discretized inclusion [START_REF]k, l)-algebraic stability of Runge-Kutta methods[END_REF])

Tli+1 E 'r/j + hF(Tij), j = 0, 1, • • •, N -1,

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 [START_REF]k, l)-algebraic stability of Runge-Kutta methods[END_REF], 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.

The construction of B. N. Pshenichny is extended in [START_REF] Dontchev | Errorestimatesfordiscretized differential inclusions[END_REF] in the following way: Let f) : I ---+ Rn be a Lipschitz function. We obtain TJN ( •) by the following procedure: [START_REF]The dichotomy of stiffness: Pragmatism versus theory[END_REF] TJf/ =Yo, TJf-n is the projection of y(ti+l) on the set TJ;' + hF(tj, TJJ'),

j = 0, 1, • • •, N -1.
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 [START_REF] Dontchev | Errorestimatesfordiscretized differential inclusions[END_REF] by a direct proof. More precisely, the following theorem is proved. THEOREM 6.2 ( [START_REF] Dontchev | Errorestimatesfordiscretized differential inclusions[END_REF]). Suppose that the map F : I x Rn ===? Rn satisfies the following conditions.

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

U C Jltn
Then there exists a constant c such that if "'N ( •) is obtained by [START_REF]The dichotomy of stiffness: Pragmatism versus theory[END_REF], 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 [START_REF]The dichotomy of stiffness: Pragmatism versus theory[END_REF] 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 [START_REF] Dontchev | Errorestimatesfordiscretized differential inclusions[END_REF].

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 [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations-Runge-Kutta and General Linear Methods[END_REF] In fact, to obtain [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations-Runge-Kutta and General Linear Methods[END_REF] 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 haus(X, Xh) = O(h) in C(I)n.

In [START_REF] Dontchev | Errorestimatesfordiscretized differential inclusions[END_REF] 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 [START_REF]k, l)-algebraic stability of Runge-Kutta methods[END_REF], e.g., by polygons. Let "'N (•)satisfy

J •-o 1 ... N-1 -' ' ' ' "' o = Yo,
where haus (Pj, F(ryj)

) = O(h), j = 0, 1, • • •, 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 [START_REF] Kastner | Difference methods with selection strategies for differential inclusions[END_REF], 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:

(i) F 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.

(iii) The approximations of the initial value y 0 satisfy ll"' <f -Yoll = O(h).

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 as N---+ oo. max haus (Rf, R(ti)) ~ 0

OSiSN

The proof of this theorem can be easily obtained from [START_REF] Anasyuk | Equations of attainable set dynamics, Part 1: Integral funnel equations[END_REF], for a more recent presentation see [START_REF] Anasyuk | Equations of attainable set dynamics, Part 1: Integral funnel equations[END_REF]. First we need some notation. DEFINITION In [START_REF] Tolstonogov | On the integral funnel equation of differential inclusions[END_REF] 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 [START_REF] Kurzhansky | Description of the bundle of control-system trajectories remaining in a given set for a given time[END_REF].

Let us prove Theorem 7.1. From [START_REF] Cabannes | Mecanique des milieux continus-propagation des discontinuites dans les cordes vibrantes soumises a unfrottement solide[END_REF] we have (15) haus ( R{", R( t1)) :::; hrp( h) and ( 16)

haus (R(t2), U (x + hF(t1,x))) < hcp(h), xER (tt)
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 [START_REF] Clarke | Optimal control and the trne hamiltonian[END_REF] haus ( U (x + hF(t1, x)), U (x + hF(t1, x))) < (1 + kh)hcp(h).

xERf" xER(ti)

Combining ( 16) and ( 17) we have haus (Rf, R(t2)) < (1 + kh)h<p(h) + h<p(h). 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

( T-to )N R(T) = J~ I+ N F Yo•
More precisely, the following result is obtained in [START_REF] Wolenski | The exponential formula for the reachable set of a Lipschitz differential inclusion[END_REF]. (ii) If, in addition, F is assumed to have convex values, then for all T > 0, R(T) C ~~~ (J + T; to F) N Yo• 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 [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations-Runge-Kutta and General Linear Methods[END_REF] was observed (on slightly different conditions) independently in [START_REF] Dontchev | Errorestimatesfordiscretized differential inclusions[END_REF] and by M. S. Nikol'skiy [START_REF]On the attainability domain of differential inclusions[END_REF] 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 [START_REF] Nikoi | SKIY, On approximation of the attainability domain of control process[END_REF] contains a similar result for a control system. In [START_REF]On a method for approximation of attainable set for a differential inclusion[END_REF] the following modification of Euler scheme is proposed:

E~ {yo}, U = a1B, Eft_1 = U (x + hr(x)coF(x) + ah 2 B), i=Ol•••N-l ' ' ' ' xEEfnU
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 [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations-Runge-Kutta and General Linear Methods[END_REF] an estimate of order O(h) follows for the distance between E}J n U and R(T). For an extension, see [START_REF] Rurenko | On the problem of approximate computation of the attainable set fora differential inclusion[END_REF].

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 [START_REF] Veliov | Second order discrete approximations to linear differential inclusions[END_REF]- [START_REF]Second order discrete approximations to strongly convex differential inclusions[END_REF] 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 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 [START_REF]Optimization and NonsmoothAnalysis[END_REF].

Let Nand k be integers, h = (T-t 0 )/N, A~ be nxn-matrices and UJv be sets in JR.n. By a discrete approximation (of higher-order) to [START_REF]Optimization and NonsmoothAnalysis[END_REF], we mean the following:

(20) 'f/~1 E A~11f + UJv, i = 0, 1, .. •, N -1, 11~ =Yo• Suppose that A~, Ut satisfy for N EN', (21) exp(Ah)
fo' cxp(As)Uds UfV + 0 (h•+I) .

(i) Let 1JN = ( 11~, • • • , 11~) be a solution of [START_REF] Verwer | Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations[END_REF]. Then, from [START_REF] Dontchev | Time-scale decomposition of the reachable set of constrained linear systems[END_REF]

Thus, under the conditions ( 21) and [START_REF] Dontchev | Errorestimatesfordiscretized differential inclusions[END_REF], the Hausdorff distance between the sets of solutions of ( 18) and [START_REF] Verwer | Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations[END_REF], in the sense of ( 23) and [START_REF] Fiuppov | Differential equations with discontinuous right-hand side[END_REF], is of order O(hk).

The condition [START_REF] Dontchev | Time-scale decomposition of the reachable set of constrained linear systems[END_REF] 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 [START_REF]Classical solutions of differential equations with multivalued right-hand side[END_REF] This means that ( 24) is satisfied, i.e., every solution of ( 18) can be approximated of order O(hk) by a solution of [START_REF] Verwer | Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations[END_REF]. The converse case needs, however, more conditions.

It follows from [START_REF]Second order discrete approximations to strongly convex differential inclusions[END_REF] 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 [START_REF]Second order discrete approximations to strongly convex differential inclusions[END_REF] V. M. 28)

77~1 E 77{" + 0.5h { z + F (ti+l, 17{" + hz) : z E F(ti, 17{")}, 1]~ E Ro 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 I ax is Lipschitz in the set ( l' t, X) E .c X ~ X s' 8-r I at is Lipschitz in (l, t) E .C x ~uniformly in x E S, where .C = {l E Rn : 0.5 < IIlii < 2};

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 [START_REF]Stability properties of implicit Runge-Kutta methods[END_REF] l*y(t, x, l) = supp(l, F(t, x))

is Lipschitz in t E ~ uniformly in ( l, x) E £ x S.

THEOREM 8.1 ([73]). On the condition (a) there exist constants c and N 1 such that if N > N1 then (i) For every solution y of [START_REF] Formalsky | Controllability and Stability of Systems with Restricted Ressources[END_REF] there is a discrete solution 1JN = ( 1J/i, • • •, 17~) of (28) such that [START_REF] Graettinger | Hyperplane method for reachable state estimation for linear timeinvariant systems[END_REF] (ii) For every solution 1JN of [START_REF] Frank | Order results for implicit Runge-Kutta methods applied to stiff systems[END_REF] there exists a solution y of [START_REF] Formalsky | Controllability and Stability of Systems with Restricted Ressources[END_REF] such that [START_REF] Graettinger | Hyperplane method for reachable state estimation for linear timeinvariant systems[END_REF] holds. If in addition (b) is satisfied, then h 3 1 2 in ( 30) can be replaced by h 2 • In [START_REF]Second order discrete approximations to strongly convex differential inclusions[END_REF] 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 E ~, x E S. In this case the discrete inclusion ( 28) can be equivalently written as ryf' !t-1 Eryf +0.5h{t.p(ti,TJf,u)+<p(ti+t,TJf +h<p(ti,TJf,u)) :uE U},

where <p(t,x,u) = f(t,x) + g(t,x)u. For fixed u E U this is exactly the second-order Euler-Cauchy scheme for right-hand side <p(t, x, u). Another second-order approximation for the trajectory bundle without strong convexity of the right-hand side is obtained in [START_REF]Approximations to differential inclusions by discrete inclusions[END_REF] (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 Observe that if B = I and A, B are constant, from [START_REF] Kas | 1NER-MARESCH, Differenzenverfahren hOherer Ordnung fii.r Differentialinklusionen[END_REF] we obtain the scheme [START_REF]Differential Equations with Discontinuous Righthand Side, Mathematics and Its Applications[END_REF].

THEOREM 8.2 ([71 ]). Under the condition (c) there exists a constant c such that for every N EN', hans (X%, R(T)) < ch 2 . 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 [START_REF]Approximations to differential inclusions by discrete inclusions[END_REF] that approximation order for the solution set may be only O(h).

The paper [START_REF] Veliov | Second order discrete approximations to linear differential inclusions[END_REF] 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 [START_REF] Kas | 1NER-MARESCH, Differenzenverfahren hOherer Ordnung fii.r Differentialinklusionen[END_REF]. More precisely, consider the inclusion [START_REF] Kahl | Direkte Methoden fii.r optimale Steuerungsprobleme[END_REF] 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 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. [START_REF] Schilling | Simpliziale Algorithmen zur Berechnung von Fixpunkten mengenwerliger Operatoren[END_REF], [START_REF]An algorithm to solve boundary value problems for differential inclusions and applications in optimal control[END_REF], or by difference methods, cf. [START_REF] Niepage | Inverse stability and convergence of difference approximations for boundary value problems for differential inclusions[END_REF]. 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.

LINEAR MULTISTEP METHOD 3 . 1 .

 31 Let with ar # 0, laol + lbol > 0, h= _T_t_o N Let there be given starting values 1]j E }Rn and corresponding starting selections

  The coefficients bi are nonnegative

  81 and pp. 203-210].

bo = 1 ,r = 4 :

 14 bt = 0, The following selections of coefficients were tested in [31]: r = 2: ao = 0, ao = 0, b 0 = 0.1, at= 0, bt = 0.4,

Example 3 . 3 .

 33 Find a function z( •) on I = [0, T] with absolutely continuous derivative such that

FIG. 1 .

 1 FIG. 1. Solution z( •) and its derivative i( • ).

8 FIG. 2 .

 82 FIG. 2. Phase porlrait of z( • ).

  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. T}j + h 2:: ap.v cjv , v=l (jv E F (tj + Cvh, fjv)

  DEFINITION 4.2. Let !Rn be equipped with the scalar product (•I•) and the corresponding induced norm II •II• A set-valued mapping F : J X Rn ==:} Rn satisfies a one-sided Lipschitz condition with (one-sided) Lipschitz constant Lp if {(I-(21xi-x2) < LF!Ixi-x2!1 2 holds for all XI, x2 E JRn and all (I E F(t, xi), (2 E F(t, x2) uniformly for all t E I.An important subclass of all right-hand sides satisfying a one-sided Lipschitz condition consists of set-valued mappings F : I x JR.n => Rn possessing a decomposition[START_REF] Burrage | k, l)-algebraic stability of Gauss methods[END_REF]

COROLLARY 4 . 4 .

 44 Let in Theorem 4.3 the exact solution be (J..t + 1 )-times continuously differentiable on the whole interval I, and the initial approximations satisfy 11"1~-Yo !I= 0 (h~-').

FIG. 5 .-FIG. 6 .

 56 FIG. 5. Approximation of i( •)by implicit midpoint rnle with h= 0.005.

FIG. 7 .

 7 FIG. 7. Approximate solution TJN (•)for selection with minimal norm.

FIG. 8 .

 8 FIG. 8. Phase portrait of TJN (•)for selection with minimal norm. 0 0

FIG. 9 .

 9 FIG. 9. Selection ( N (•)for selection with minimal norm.

FIG. 10 .FIG. 11 .

 1011 FIG.[START_REF]k, l)-algebraic stability of Runge-Kutta methods[END_REF]. Approximate solution rJN (•)for selection with minimal variation.

FIG. 12 .

 12 FIG. 12. Selection vN (•)for selection with minimal variation.

ForFlo. 13 .

 13 this selection strategy with reference trajectory y(t) _ OR2, 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. Approximate solution TIN (•)for selection with respect to reference trajectory.

FIG. 14 .FIG. 15 .

 1415 FIG. 14. Phase portrait of TIN (•)for selection with respect to reference trajectory.0 0

THEOREM 7 . 4 (

 74 [START_REF] Wolenski | The exponential formula for the reachable set of a Lipschitz differential inclusion[END_REF]). Let G c Rn be open and F : Rn =::::? Rn be a map with nonempty compact values on G locally Lipschitz on G. Fix y 0 E G. Let T =sup {r: cl ( U R(t)) is compact}. to9:::;T (i) For t 0 < T < T ( T-to )N li:~~ I+ N F Yo c cl (R(T)).

  5h(A(ti) + A(tH 1 )) + 0.5h 2 A(tH 1 ) 2 , Bf 0.5h(B(ti) + B(ti+t)) + 0.5h 2 A(ti+t)B(ti+t) •

2.25 2.50 2.75 4.3 FIG

. 3. Approximation of i:(•) by Euler method with h = 0.005. 2.00 2.25 2.50 2.75 FIG. 4. Approximation of z(•) by Runge-Kutta method with h = 0.005.

  7.2 ([50]). The closed valued and continuous map R : I ===> Rn with R(to) = y 0 is an R-solution of Initial Value Problem 1.1 if and only if On the conditions of Theorem 7.1 there exists a unique R-solution of Initial Value Problem 1.1 whose value at any t E 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.

	(14)	~ haus (R(t +h), U (x + hF(t, x))) ~ 0
		xER(t)
	ash ---+ 0 uniformly in t E /.
		THEOREM 7.3 ([50]).

  Proceeding in the same manner by induction we get

	i haus (Rf, R(ti)) < L(l + kh)i-Jlh<p(h)
		Jl=l
		< <p(h)(T-to)ek(T-to)
	fori = 1, • • •, 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
	R'(	(I+ hF)yo,
	Rf	U (x + hF(x)) =(I+ hF) 2 y 0 ,
		xERf
	Rf	(I+ hF)iYo,

  Let uN be defined on I as uN (t) = uf (t) fort E [ti, ti+ 1 ), i = 0, 1, .. •, N -1. If yN is a solution of AyN (t) + uN (t) for a.e. t E I Let ry{j = y 0 • We determine ryf' !t. 1 on the basis of ryf as projection of y(ti+ 1 ) on the set Atryf + U~. Then from

	Yo,
	then it is easy to see that
	(23)
	(ii) Now let y be a solution of (18). Then
	1 ti+l y(ti+t) E eAhy(ti) + t, eA(t-;+ 1 -s)Uds
	it follows that
	which yields
	and (22),
	l rJft-1 E A~r,[" + ut c eA(ti+l-tdrJf + ti eA(tH1-s)Uds + 0 (hk+l) ti+l

fori= 0, 1, • • •, N -1. Hence there exists a measurable function uf' with uf' (t) E U for almost every t E I, such that if yfl solves yf" (t) Ayfl (t) + uf (t) for a.e. t E [ti , ti+l],

yf (ti) = ryf,

then y[" (ti+t) = 'TJ~l + 0 (hk+l) .

c Aty(ti) + ut + O(hk+t ).

  Veliov considers the following more general problem:

	(27)
	approximated by
	(

y(t) E F(t, y(t)) for a.e. t E I, y(to) E Ro

  for the inclusion E f(t,y(t)) + l::Yi(t,y(t))[O, 1], y(to) E Ro .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[START_REF] Kahl | Direkte Methoden fii.r optimale Steuerungsprobleme[END_REF] y(t) E A(t)y(t) + B(t)U for a.e. t E I, y(to) E Xo,where A(t) and B(t) are nxn and nxr matrices, U c lRr, on the following assumption:

r y(t) i=I

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

This research was supported by the DFG (Deutsche Forschungsgemeinschaft) and by the Institute of Mathematics at the

where q is a fixed integer, P and Q are arbitrary matrices, h = (Tt 0 )jN.

THEOREM 8.3 ([71 ]). For any controllable pair (A, B) (i.e., rank[B, ABJ = 2), there is a constant c > 0 such that haus (X%, R(T)) > ch 2 for all N E N'.

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 [START_REF] Formalsky | Controllability and Stability of Systems with Restricted Ressources[END_REF], [START_REF] Graettinger | Hyperplane method for reachable state estimation for linear timeinvariant systems[END_REF], or ellipsoids [16], [START_REF] Kurzhanski | The problem of control synthesis for uncertain systems: Ellipsoidal techniques, in Modelling, Estimation and Control of Systems with Uncertainty[END_REF]. A different approach, based on a timescale decomposition of the reachable set, is presented in [START_REF] Dontchev | Time-scale decomposition of the reachable set of constrained linear systems[END_REF].

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 [START_REF] Stewart | High Accuracy Numerical Methods for Ordinary Differential Equations with Discontinuous Right-Hand Side[END_REF], [START_REF]A high accuracy method for solving ODEs with discontinuous right-hand side[END_REF] 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 [START_REF]Modified Euler methods for differential inclusions[END_REF].

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., [START_REF] Saint-Pierre | Approximation of slow solutions to differential inclusions[END_REF].

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