The implicit Euler scheme for one-sided Lipschitz differential inclusions

We propose a set-valued version of the implicit Euler scheme for 
relaxed one-sided Lipschitz differential inclusions and prove that 
the defining implicit inclusions have a well-defined 
solution. Furthermore, we give a convergence analysis based on stability theorems, which shows that the 
set-valued implicit Euler method inherits all favourable stability properties from the 
single-valued scheme. The impact of spatial discretization is discussed, a fully discretized 
version of the scheme is analyzed, and a numerical example is given.

Introduction. We consider the differential inclusioṅ in the Euclidean space Ê m , where F is continuous, has convex and compact values, and satisfies the relaxed one-sided Lipschitz condition (ROSL, see Definition 1 in Section 1 and assumption (A2) in Section 3). According to Theorem 2 in [4], the set of solutions of (1) is a nonempty compact subset of the space of continuous functions C([0, T ], Ê m ) equipped with the supremum norm.
In most articles which deal with the set-valued explicit Euler scheme, the righthand side of the differential inclusion is required to be Lipschitz continuous. One of the first attempts to give error estimates for the set-valued Euler scheme Φ exp h (t, x) := {y ∈ Ê m : y ∈ x + hF (t, x)} (2) was presented in [8]. The Euler method proposed in [9] only uses the extremal points of the right-hand side and adjusts the time discretization to this predefined discretization in the velocity space.
In [13], the spatial discretization is considered explicitly, and a more detailed analysis of spatial discretization effects has been given in [3]. The performance of the explicit Euler scheme in the presence of state constraints has been thoroughly investigated in the remarkable paper [2]. Recently, the error analysis of the setvalued Euler scheme has been extended to nonconvex-valued differential inclusions in [17]. The existence of attractors of differential inclusions and the corresponding Euler discretizations has been established in [12].
The explicit Euler scheme has been examined in the context of relaxed one-sided Lipschitz right-hand sides as well, see [6] and [7]. In [6], a ROSL version of the Filippov Theorem is proved which serves as the main tool for the error estimates, while [7] gives a detailed analysis of the dynamics of (1) and the corresponding Euler solutions.
The aim of the present paper is to establish an implicit set-valued Euler scheme for continuous and ROSL right-hand sides. In Section 1, we present a solvability theorem for ROSL maps which extends results from [5] and is needed in Section 2 to show that implicit Euler steps are well-defined and Lipschitz continuous with nonempty and compact values. For the system (1) a single implicit Euler step with stepsize h is given by Φ imp h (t, x) := {y ∈ Ê m : y ∈ x + hF (t + h, y)}. (3) In Section 3, we present error estimates for the implicit Euler scheme on bounded and unbounded intervals. We prove a version of the Filippov Theorem which is very similar to the one presented in [6] and a discrete stability theorem for the implicit Euler scheme from which all further estimates follow in a natural way.
Finally, we introduce a fully discretized implicit Euler scheme in Section 4 and show how a solution of the implicit inclusion (3) can be approximated. We apply this method to the well-known Michaelis-Menten model in order to demonstrate that for stiff differential inclusions the implicit Euler scheme (3) respects the stability properties of (1), while the explicit Euler fails to recognize them. In a way this shows that the favourable properties of A-stable and B-stable methods for stiff ODEs (see [10]) transfer to the implicit Euler method for differential inclusions.
The notation used in this paper is standard: The Euclidean norm is denoted by | · |, while A := sup a∈A |a| denotes the maximal norm of the elements of a set For a set-valued mapping F : Ê m → C(Ê m ), we define F (A) := ∪ a∈A F (a). Thus, the composition F •G of two set-valued mappings F and G is given by

1.
A solvability theorem for ROSL maps. The notion of relaxed one-sided Lipschitz (ROSL) set-valued mappings was coined by Tzanko Donchev, see e.g. [4]. It generalizes the concepts of Lipschitz continuity and the (strong) one-sided Lipschitz property. A detailed analysis of this property can be found in [5] and several other works of this author.
Definition 1. A mapping F : Ê m → CC(Ê m ) is called relaxed one-sided Lipschitz with constant l ∈ Ê if for every x, x ′ ∈ Ê m and y ∈ F (x) there exists some It was shown in [5], Theorem 3, that ROSL mappings with negative Lipschitz constant are onto for a quite general class of Banach spaces. The following results are more specific, because they provide information about the location of a point where a desired value is attained.
As F is usc, and there exists an α > 0 such that |z| 2 ≤ R 2 follows from (6). This means that for this fixed α, and H(·) is also usc. By the Kakutani Theorem, H and thus also G have a fixed point x R in B R (0), which implies that 0 ∈ F (x R ).
Corollary 3. Let F : Ê m → CC(Ê m ) be usc and ROSL with constant l < 0, and let x ∈ Ê m and y ∈ Ê m be given. Then there exists anx ∈ Ê m with y ∈ F (x) and 2.1. Existence and estimates of solutions. The following theorem is the central existence and stability statement for the implicit Euler method. It follows from Corollary 3, but the roles of x and y are interchanged because of the implicit character of the scheme.
Theorem 4. If F is usc and ROSL with constant l such that lh < 1, and if Proof. Consider the set-valued mapping Since F is ROSL, for any z, z ′ ∈ Ê m and ξ ∈ F (z) there exists ξ ′ ∈ F (z ′ ) such that In particular, which means that G h satisfies the ROSL condition with constant −(1 − lh) < 0.
According to Corollary 3, there exists aȳ ∈ Ê m such that −x ∈ G h (ȳ), i.e.ȳ ∈ Φ h (x), and Proof. Let y ∈ Φ h (x), i.e. y ∈ x + hF (y). For any x ′ ∈ Ê m , Theorem 4 guarantees the existence of a y ′ with y ′ ∈ x ′ + hF (y ′ ) such that The implicit Euler scheme is robust under perturbations of the map F .
Moreover, if G is usc and one-sided Lipschitz with the same constant and if it has convex and compact values, then holds for all x ∈ Ê m .
Lemma 7. Let F : Ê m → CC(Ê m ) be an usc and relaxed one-sided Lipschitz setvalued mapping with constant l. If hl < 1, the diameter of the values of the implicit Euler scheme is bounded by Proof. Let y, y ′ ∈ Φ h (x), i.e. y ∈ x + hF (y) and y ′ ∈ x + hF (y ′ ). Because of the ROSL property, there exists an η ∈ x + hF (y ′ ) such that

2.2.
Structure of the implicit Euler set. Estimate (12) for the diameters and estimates (43) and (45) for the location of the images of the method inherit the implicit character of the scheme, which is inconvenient for practical implementation, because only crude upper bounds for the diameters or the norm of the images of F can be used.
be an usc and relaxed one-sided Lipschitz set-valued mapping with constant l. If hl < 1, the values of the implicit Euler scheme are nonempty and compact.
Proof. Theorem 4 guarantees that the images of the implicit Euler scheme are nonempty. Let because the mapping y → y − hF (y) is usc. Consequently, x ∈ y − hF (y), and the images of Φ h are closed.
Assume that Φ h (x) is unbounded for some x ∈ Ê m . Without loss of generality, By definition, y k ∈ hF (y k ). As hF (·) is ROSL with constant hl, there exists a sequence {y ′ k } k∈AE in F (0) such that for all k ∈ AE, Below we will show that implicit Euler sets of a mapping F : Ê m → CC(Ê m ) that satisfies the ROSL condition need not be convex in dimension m > 1. However, in dimension one the answer is affirmative. Proof. The set-valued mapping F can be represented as where α(·) : Ê → Ê and β(·) : Ê → Ê are single-valued functions.
The ROSL condition states that for all x, x ′ ∈ Ê and y ∈ F (x) there exists a In this particular setting, this is equivalent with Because of (13), the implicit Euler set is Then either x > y 1 − hα(y 1 ) or x < y 1 − hβ(y 1 ). In the first case, , which yields hl > 1. In the second case, and Hence the image of the implicit Euler scheme in dimension one is convex whenever the scheme is well-defined according to Theorem 8.
In Ê m with m > 1, the ROSL condition does not seem to be as powerful as in dimension one, because it controls the expansion of the images in only one direction.
The following example provides a continuous map G : Ê 2 → CC(Ê 2 ) which satisfies the ROSL property and which has nonconvex implicit Euler sets.
Example 10. Let a set-valued mapping F : Ê × {0} → CC(Ê 2 ) be given by Consider its extension F : Ê 2 → CC(Ê 2 ) defined by A short discussion shows that . From the negative slopes in the first component of z ± one finds that for all Thus for any two points In the general case, given Finally, the mapping satisfies the ROSL condition with constant l = 3 4 . (Note that G satisfies the conditions of Lemma 12.) The implicit Euler sets associated with G are For x = 0 and h = 1 we have hl < 1 as in Theorem 4, but The following geometric characterization of convexity still holds.
is (either empty or) convex. In particular, this is satisfied if graph F (·) is convex.
If graph F (·) is convex then graph{x + hF (·)} and its intersection with the diagonal are also convex.
The question of connectedness of the implicit Euler sets could not be fully resolved. The following lemma gives an affirmative answer if the right-hand side F is parametrized by a continuous family of ROSL selections.
Lemma 12. Let F : Ê m → CC(Ê m ) be a set-valued mapping, let U ⊂ Ê d be any path-connected set, and let f : Ê m × U → Ê m be a single-valued function such that • there exists some l ∈ Ê such that lh < 1 and f (·, u) is ROSL with constant l for any u ∈ U ; • f (x, ·) is continuous for any x ∈ Ê m .
Then the images of the implicit Euler scheme are path-connected.
Proof. Theorem 4 implies that the implicit Euler scheme applied to the function f (·, u) has a solution. It is well-known that this solution is unique, because y = x + hf (y, u) and y ′ = x + hf (y ′ , u) imply and thus |y − y ′ | 2 = h f (y, u) − f (y ′ , u), y − y ′ ≤ lh|y − y ′ | 2 .
As lh < 1, it follows that y = y ′ . If u, u ′ ∈ U and y satisfies y = x + hf (y, u), then Theorem 4 and Now let y, y ′ ∈ Φ h (x) be arbitrary solutions of the implicit Euler scheme, and let u, u ′ ∈ U be the corresponding parameters such that y = x + hf (y, u) and y ′ = x + hf (y ′ , u ′ ). As U is path-connected, there exists a continuous path σ : (A2) There is a continuous function l : [0, T ] → Ê such that for t ∈ [0, T ], l(t)h < 1 and the mapping F (t, ·) is ROSL with constant l(t), i.e. for every x, x ′ ∈ Ê m and y ∈ F (t, x) there exists some y ′ ∈ F (t, x ′ ) such that Consider the differential inclusioṅ where T ∈ Ê + . The set of solutions to (18) will be denoted S(T, 0, x 0 ), while Consider also the discrete-time (but spatially continuous) system with step-size h := T /N , where N ∈ AE, t k := kh. The sets S h (T, 0, x 0 ) and for k ∈ AE with kh ≤ T are the discrete-time solution set and reachable set at time kh, respectively.
3.1. Stability theorems. The following theorem is a continuous version of the Filippov Theorem for ROSL right-hand sides. The proof follows the argument in [6], but existence of a solution is guaranteed by continuity instead of a growth condition. Proof. Define the set-valued mapping The values of H(·, ·) are nonempty: For given t and x, there exists some w ∈ F (t, y(t)) such that |ẏ(t) − w| = g(t). The ROSL condition implies the existence of some v ∈ F (t, x) with Thus, which means that v ∈ G(t, x).
Precisely the same argument yields the desired result on the infinite time interval.
The following theorem is a discrete version of the Filippov Theorem for the implicit Euler scheme.
Proof. Let {x j } N j=0 be given, and assume that a solution {y j } of (19) has already been constructed for j ∈ {0, . . . , k}, k < N . Because of Theorem 4, there exists a solution y k+1 ∈ Φ h (y k ) such that By induction,

3.2.
The finite time interval. The boundedness of the continuous as well as of the numerical solutions follows immediately from the continuous and the discrete Filippov Theorems.

Corollary 15.
Any solution x(·) of (18) is bounded by Proof. The statement follows from Theorem 13 with y(t) = 0 for all t ∈ [0, T ]. Because of Corollaries 15 and 16, all solutions of (18) and (19) are contained in a ball B M (0) with sufficiently large radius M > 0. In particular, we may assume because of continuity that F is globally bounded, i.e. that there exists some P > 0 It is well-known that a continuous set-valued mapping F : Ê × Ê m → CC(Ê m ) is uniformly continuous on every compact set K ⊂ Ê× Ê m , i.e. for any ε > 0 there exists a δ > 0 such that dist H (F (x), F (x ′ )) ≤ ε whenever x, x ′ ∈ K and |x − x ′ | ≤ δ. Hence the moduli of continuity and τ (δ) := sup{dist H (F (t, x), F (t ′ , x)) : are well-defined.
Proof. By definition, y k+1 ∈ y k + hF (t k+1 , y k+1 ). Define the absolutely continuous Theorem 13, the continuous version of the Filippov Theorem, guarantees (piecewise, and thus globally) the existence of a solution x(·) of (18) such that for all t ∈ [0, T ] and in particular for all t ∈ {kh : 0 ≤ k < N }.

1−l(t)h
Proof. For k ≥ 1, let Then and the desired result follows by applying the discrete Filippov Theorem 14 with x k := x(t k ).
We have to assume that the relaxed one-sided Lipschitz constant is given by a Note that if l max := sup t∈[0,∞) l(t) < 0, the integrals and Riemann sums mentioned above are bounded indeed. This behaviour matches the shadowing results presented in [15], where it is shown that the errors of reasonable numerical approximations of differential inclusions with strictly contractive flow are bounded on [0, ∞). However, the above results are stronger than the shadowing theorem, because the boundedness condition on l(·) is by no means necessary.
3.4. The influence of spatial discretization. Corollary 5 implies the following estimate for errors induced by spatial discretization. Recall that compositions of set-valued maps are inductively defined as explained in the Introduction.
Theorem 19. Let F be a set-valued mapping which satisfies (A1) and (A2), and for some ρ ≥ 0. Then the error caused by spatial discretization is whenever dist H (X, X ρ ) ≤ ρ and kh ≤ T .
Proof. Assume that (36) already holds true for some k ∈ AE. Then, for any point By Corollary 5, Φ h (t k , ·) is Lipschitz continuous with Lipschitz constant Remark 20. The detailed analysis of the explicit Euler method given in [3] shows that it is necessary to choose ρ := h 2 in order to obtain convergence of order one for the fully discretized scheme. As the right-hand side of (36) behaves like the integral it is reasonable to adapt the same setting for ρ in the present context. 4. Implementation. Let ∆ ρ := ρ m be an equidistant grid in Ê m . For x ρ ∈ ∆ ρ , we define the fully discretized implicit Euler scheme Φ h,ρ : As in (30), let χ(·) denote the modulus of continuity of F w.r.t. the space coordinate. The reachable set of the fully discretized implicit Euler scheme is the set ). The following theorem shows that the tolerance ρ + hχ(ρ) in the right-hand side of (37) leads to a reasonable error estimate.
Theorem 21. If F : Ê m → CC(Ê m ) is continuous and ROSL with constant l such that lh < 1, then for any x ρ ∈ ∆ ρ .
The following lemmas provide information about the location of the image of the implicit Euler scheme. and hold.
Lemma 24. Let F : Ê m → CC(Ê m ) be continuous and ROSL with constant l such that lh < 1. Then and Proof. Because of Theorem 4, If y ∈ x + hF (y), then |y − x| ≤ h F (y) , which proves (45).  We have implemented a very simple, but very robust version of the implicit Euler scheme, which does not require any assumptions apart from (A1) and (A2). In every step, we estimate the region in which the reachable set R h,ρ ((k + 1)h, 0, x 0 ) must be located according to Lemma 23 or 24. Unfortunately, estimates (43) and (45) are of implicit nature and allow only a rough guess for the location based on crude upper bounds for the norm of the images of F . For every grid point y in this region, we compute an approximation of −y + hF (y) with a suitable precision and check for any grid point x ∈ R h,ρ (kh, 0, x 0 ), whether dist(y, x + hF (y)) = dist(−x, −y + hF (y)) ≤ ρ + hχ(ρ).
If the distance is below the tolerance, the point y is accepted as an element of R h,ρ ((k + 1)h, 0, x 0 ).
As a model example, we choose the Michaelis-Menten model which is of considerable importance in theoretical Biology. It describes a simple biochemical process, where an organic substrate molecule is changed into a product by means of an enzyme. At first, the substrate molecule forms a complex with an enzyme. In a second step, the substrate molecule reacts and becomes the product, which leaves the complex in a third step. The chemical equation can easily be expressed as a four dimensional system of differential equations which can be reduced to the system where the k i are rate constants, e 0 is the concentration of the enzyme, and s and c denote the concentrations of the substrate and substrate-enzyme complex, respectively. For a more detailed discussion of this model see e.g. [14].
If we choose the parameters e 0 = 10, k 1 = 1, k −1 = 0.6, and k 2 = 0.3, the eigenvalues of the Jacobian Df (0, 0) are λ 1 = −0.28 and λ 2 = −10.62, which means that the ODE is moderately stiff in a neighbourhood of zero. In such a setting, explicit numerical methods are known to yield uncontrollable errors if hλ j is outside the region of absolute stability (see [10]).
The solutions of the single-valued explicit and implicit Euler schemes applied to (46) with the above parameters are displayed in Figure 1. Whereas the solution of the explicit Euler scheme with step size h = 0.1 still converges to the stable equilibrium zero, the method produces an utterly erratic trajectory when performed with larger step sizes. The implicit Euler scheme, however, retains the stability properties of the underlying ODE as expected.
This differential inclusion is Lipschitz continuous and thus ROSL with constant l = L ≈ 10 in a neighbourhood of the origin. Figure 2 shows its reachable sets obtained using the explicit Euler scheme which we presented in [3] and those computed with the implicit Euler scheme described above. In all four computations, we choose ρ = 0.005 in order to provide images of equal optical quality. As in the single-valued case, the explicit Euler with step size h = 0.1 yields a reasonable result, but for larger step-sizes its images grow rapidly and diverge, whereas the implicit Euler scheme converges to the stable equilibrium. Its performance is Figure 2. Performance of the set-valued explicit and implicit Euler schemes applied to the Michaelis-Menten system better than predicted e.g. by Theorem 4, because it works well even for some values lh > 1.
Please note that neither the implicit nor the explicit Euler scheme (for h = 0.1) converge exactly to zero but attain a fixed set near the origin in finite time. This effect is due to spatial discretization: If the norm of the right-hand side becomes too small, the numerical scheme cannot reach the neighbouring grid points any more and eventually becomes constant.
In our computations, searching the spatial grid for elements of Φ h,ρ according to Lemma 24 lead to a better performance than according to Lemma 23, which is probably due to the fact that it is necessary to compute the explicit Euler image as a predictor in the latter case. However, the computational costs are still extremely high.
The main obstacle for the construction of more efficient algorithms is our limited knowledge about the properties of the images of the implicit Euler scheme. In dimension one, Lemma 9 guarantees that the images are convex, but for any higher dimensional space it is even unclear whether they are connected. Because of Lemma 12 there is hope that connectedness of the images might be true, which would be a crucial precondition for improvements.