On a randomized backward Euler method

On a randomized backward Euler method for
nonlinear evolution equations with
time-irregular coefficients

Monika Eisenmann Monika Eisenmann
Technische Universität Berlin
Institut für Mathematik, Secr. MA 5-3
Straße des 17. Juni 136
DE-10623 Berlin
Mihály Kovács Mihály Kovács
Department of Mathematical Sciences
Chalmers University of Technology and University of Gothenburg
SE-412 96 Gothenburg
Raphael Kruse Raphael Kruse
Technische Universität Berlin
Institut für Mathematik, Secr. MA 5-3
Straße des 17. Juni 136
DE-10623 Berlin
 and  Stig Larsson Stig Larsson
Department of Mathematical Sciences
Chalmers University of Technology and University of Gothenburg
SE-412 96 Gothenburg

In this paper we introduce a randomized version of the backward Euler method, that is applicable to stiff ordinary differential equations and nonlinear evolution equations with time-irregular coefficients. In the finite-dimensional case, we consider Carathéodory type functions satisfying a one-sided Lipschitz condition. After investigating the well-posedness and the stability properties of the randomized scheme, we prove the convergence to the exact solution with a rate of in the root-mean-square norm assuming only that the coefficient function is square integrable with respect to the temporal parameter.

These results are then extended to the numerical solution of infinite-dimensional evolution equations under monotonicity and Lipschitz conditions. Here we consider a combination of the randomized backward Euler scheme with a Galerkin finite element method. We obtain error estimates that correspond to the regularity of the exact solution. The practicability of the randomized scheme is also illustrated through several numerical experiments.

Key words and phrases:
Monte Carlo method, stratified sampling, evolution equations, ordinary differential equations, backward Euler method, Galerkin finite element method
2010 Mathematics Subject Classification:
65C05, 65L05, 65L20, 65M12, 65M60

1. Introduction

The aim of this paper is to introduce a new numerical scheme to approximate the solution of an ordinary differential equation (ODE) of Carathéodory type


for , and of a non-autonomous evolution equation


where is a possibly nonlinear operator that is defined on a Gelfand triple and satisfies suitable monotonicity and Lipschitz conditions with respect to the second argument.

We focus on the particular difficulty that the mappings and are irregular with respect to the temporal parameter. More precisely, we do not impose any continuity conditions but only certain integrability requirements with respect to . For a concise description of the general settings we refer to Sections 3 and 6, respectively. In particular, a precise statement of all conditions is given in Assumption 3.1 for (1.1) and in Assumption 6.1 for (1.2). To develop the idea of our scheme we mostly focus on the ODE problem (1.1) in this introduction. The derivation of the numerical scheme for the evolution equation (1.2) follows analogously and will be introduced in detail in Section 6.

When considering a right hand side that is only integrable, every deterministic algorithm can be ”fooled” if it only uses information provided by point evaluations on prescribed (deterministic) points. A simple, however academic, example for this would be a right hand side like


Using a deterministic scheme that only evaluates the data at rational points, the scheme never approximates the exact solution in a reasonable sense no matter how small the step size is chosen. In the same way, one can easily construct suitable fooling functions for more general classes of deterministic algorithms, for instance, based on adaptive strategies. We refer to the vast literature on the information-based complexity theory (IBC), which applies similar techniques to derive lower bounds for the error of certain classes of deterministic and randomized numerical algorithms. For instance, see [30, 35] for a general introduction into IBC and [20, 23, 24] for applications to the numerical solution of initial value problems.

One way to construct numerical methods for the solution of initial value problems with time-irregular coefficients consists of allowing the algorithm to use additional information of the right hand side as, for example, integrals of the form


This approach is often found in the existence theory of ODEs and PDEs when a numerical method is used to construct analytical solutions to the initial value problems (1.1) and (1.2) under minimal regularity assumptions. The complexity of such methods has also been studied in [23] (and the references therein) for the numerical solution of ODEs. It is also the state-of-the-art method in many recent papers for the numerical solution of evolution equations of the form (1.2). For example, we refer to [4, 11, 21, 29].

However, it is rarely discussed how a quantity such as in (1.4) is obtained in practice. Strictly speaking, since the computation of often requires the application of further numerical methods such as quadrature rules, algorithms relying on integrals such as (1.4) are, in general, not fully discrete solvers yet. More importantly, classical quadrature rules for the approximation of are again based on deterministic point evaluations of and may therefore be “fooled” in the same way as discussed above.

Instead of using linear functionals such as (1.4) we propose the following randomized version of the backward Euler method. For , a step size , and a temporal grid with for the randomized scheme for the numerical solution (1.1) is then given by


where is a uniformly distributed random variable with values in the interval . Note that we evaluate the right hand side at random points between the grid points. Since the evaluation points vary every time the algorithm is called, it is not possible to construct a fooling function as described above.

In the example with the indicator function in (1.3) the advantage of the randomization becomes immediately evident. Since the random variable will be an irrational number with probability one, the randomized scheme will almost surely converge to the exact solution. This result will be extended to general Carathéodory ODEs in Section 4. In particular, Theorem 4.7 shows that the numerical solution from (1.5) converges with order to the exact solution of (1.1), even if is only square integrable with respect to time. Due to the results in [20] this convergence rate is optimal in the sense that there exists no deterministic or randomized algorithm based on finitely many point evaluations of with a higher convergence rate within the class of all initial value problems satisfying Assumption 3.1.

The error analysis is based on the observation that the randomized scheme (1.5) is a hybrid of an implicit Runge–Kutta method and a Monte Carlo quadrature rule. In fact, if the ODE (1.1) is actually autonomous, that is, does not depend on , then we recover the classical backward Euler method. On the other hand, if is independent of the state variable , then the ODE (1.1) reduces to an integration problem and the randomized scheme (1.5) is the randomized Riemann sum for the approximation of given by

Observe that a randomized Riemann sum is a particular case of stratified sampling from Monte Carlo integration. It was introduced in [15], [16] together with further, higher order, quadrature rules. Our error analysis of the randomized scheme (1.5) combines techniques for the analysis of both time-stepping schemes and Monte Carlo integration. In particular, since we are interested in the discretization of evolution equations in later sections, we apply techniques for the numerical analysis of stiff ODEs developed in [18] and for stochastic ODEs in [3].

Before we give a more detailed account of the remainder of this paper, let us emphasize a few practical advantages of the randomized scheme (1.5):

  1. The implementation of the randomized scheme (1.5) is as difficult as for the classical backward Euler method in terms of the requirements of solving a nonlinear system of equations. On the other hand, the scheme (1.5) does not require integrals such as if the right hand side is time-irregular.

  2. This also holds true for the computational effort. Compared to the classical backward Euler method, the randomized scheme (1.5) only requires in each step the additional simulation of a single scalar-valued random variable. In general, the resulting additional computational effort is negligible compared to the solution of a potentially high-dimensional nonlinear system of equations. More importantly, due to the randomization we avoid the potentially costly computation of the integrals .

  3. In contrast to every deterministic method based on point evaluations of , the randomized scheme (1.5) is independent of the particular representation of an integrable function. To be more precise, let and be two representations of the same equivalence class . Then, it follows that with probability one, since almost everywhere.

We remark that the last item is only valid as long as the random variable is indeed uniformly distributed in . In practice, however, one usually applies a pseudo-random number generator which only draws values from the set of floating point numbers. Since this is a null set with respect to the Lebesgue measure, the argument given above is no longer valid. Of course, this problem affects any algorithm that uses the floating point arithmetic. Nevertheless, a randomized algorithm is often more robust regarding the particular choice of the representation of an equivalence class in and, hence, more user friendly. For instance, the mapping causes problems for the classical backward Euler method as it will evaluate the mapping in the singularity at . This problem does not occur for the randomized backward Euler method with probability one.

Let us also mention that randomized algorithms for the numerical solution of initial value problems have already been studied in the literature. In the ODE case, the complexity and optimality of such algorithms is considered in [7, 20, 24] under various degrees of smoothness of . The time-irregular case studied in the present paper was first investigated in [33, 34]. See also [22, 26] for a more recent exposition of explicit randomized schemes.

The present paper extends the earlier results in several directions. In order to deal with possibly stiff ODEs we consider a randomized version of the backward Euler method and prove its well-posedness and stability under a one-sided Lipschitz condition. In addition, we obtain estimates on the local truncation error only under local Lipschitz conditions with respect to the state variable extending results from [26]. We also avoid any (local) boundedness condition on as, for example, in [7, 22].

The stability properties also qualify the randomized backward Euler method as a suitable temporal integrator for non-autonomous evolution equations with time-irregular coefficients. To the best of our knowledge, there is no work found in the literature that applies a randomized algorithm to the numerical solution of evolution equations of the form (1.2). Instead, the standard approach in the time-irregular case relies on the availability of suitable integrals of the right hand side as in (1.4). In particular, we mention [11, 21]. Further results on optimal rates under minimal regularity assumptions for linear parabolic PDEs can be found, e.g., in [4, 5, 17]. For semilinear parabolic problems optimal error estimates are also found in [29], where a discontinuous Galerkin method in time and space is considered.

This paper is organized as follows. In Section 2, we shortly introduce the notation and recapture some important concepts of stochastic analysis that are relevant for this paper. In the following Section 3, we state the assumptions imposed on the ODE (1.1). We also discuss existence and regularity of the solution. In Section 4, we then prove the well-posedness and convergence of the randomized backward Euler method in the root-mean-square sense. The ODE part of this paper is completed in Section 5 by examining a numerical example.

In Section 6, we begin by introducing the setting for the irregular non-autonomous evolution equation (1.2) that we consider in the second part of this paper. Under some additional regularity assumptions on the exact solution, we prove the convergence of a fully discrete method that combines the randomized backward Euler scheme with a Galerkin finite element method. The additional regularity assumption is then discussed in more detail in Section 7. In particular, it is shown that the regularity condition is fulfilled for rather general classes of linear and semilinear parabolic PDEs. Finally in Section 8, we demonstrate that this new randomized method for evolution equations works through a numerical example based on the finite element software package FEniCS [27].

2. Preliminaries

In this section, we explain the necessary tools from probability theory and recall some important inequalities that are needed. First, we start by fixing the notation used in this paper.

We denote the set of all positive integers by and the set of all real numbers by . In , , we denote the Euclidean norm by which coincides with the absolute value of a real number for . The standard inner product in is denoted by . For a ball of radius around the center we write .

In the following, we will consider different spaces of functions with values in general Hilbert spaces. To this end, let be a real Hilbert space and . Then we will denote the space of continuous functions on with values in by where the norm is given by

It will also be important to consider functions which are a little more regular. For we denote the space of Hölder continuous functions by with the norm given by

For , we introduce the Bochner–Lebesgue space

where the norm is given by

In the case we write .

The space of linear bounded operators from to a Banach space is denoted by and in the case of we write . The norm of this space is the usual operator norm given by

Since we are interested in a randomized scheme, we will briefly repeat the most important probabilistic concepts needed in this paper. To this end, we consider a probability space which consists of a measurable space together with a finite measure such that for every and . A mapping is called a random variable if it is measurable with respect to the -algebra and the Borel -algebra in , i.e., for every

is an element of . The integral with respect to the measure is often denoted by

The space of -measurable random variables such is finite is denoted by .

For our purposes it is important to consider the space of square integrable -measurable random variables. This space is often abbreviated by if it is clear from the context which -algebra and measure is used. The space is endowed with the norm

Equipped with this norm and the inner product

the space is a Hilbert space.

A further important concept is the independence of events . We call the events independent if for every finite subset

holds. This concept can be transfered to families of -algebras. Such a family is called independent if for every finite subset it follows that every choice of events with are independent. Similarly, a family of -valued random variables is called independent if the generated -algebras

are independent.

A family of -algebras is called a filtration if for every the -algebra is a subset of and holds for . Thus a random variable can be measurable with respect to but not with respect to for . In some of the arguments in this paper it will be important to project an -measurable random variable to a smaller -algebra . To this end, we introduce the conditional expectation of with respect to : For a random variable we introduce the -measurable random variable which fulfills

for every where is the characteristic function with respect to . The random variable is uniquely determined by these postulations. An important property of the conditional expectation of is the tower property which states that for two -algebras and of the filtration with we obtain that

In particular, if is already measurable with respect to then holds. If is independent of we obtain that .

In the course of this paper, we will often use random variables which are uniformly distributed on a given temporal interval . To denote such a random variable we write .

For a deeper insight of the probabilistic background, we refer the reader to [25].

The following inequalities will be helpful in order to give suitable a priori bounds for the solution of a differential equation and the solution of a numerical scheme.

Lemma 2.1 (Discrete Gronwall lemma).

Let and be two nonnegative sequences which satisfy, for given and , that

Then, it follows that

Lemma 2.2 (Gronwall lemma).

Let be nonnegative functions which satisfy, for given , that


For a proof of the discrete Gronwall lemma, we refer the reader to [6]. A proof for Lemma 2.2 can be found in [19].

3. A Carathéodory type ODE under a one-sided Lipschitz condition

In this section, we introduce an initial value problem involving an ordinary differential equation with a non-autonomous vector field of Carathéodory type, that satisfies a one-sided Lipschitz condition. We give a precise statement of all conditions on the coefficient function in Assumption 3.1, which are sufficient to ensure the existence of a unique global solution. The same conditions will also be used for the error analysis of the randomized backward Euler method in Section 4. Further, we briefly investigate the temporal regularity of the solution .

Let . We are interested in finding an absolutely continuous mapping that is a solution to the initial value problem


where denotes the initial value. The following conditions on the coefficient function will ensure the existence of a unique global solution:

Assumption 3.1.

The mapping is measurable. Moreover, there exists a null set such that:

  • There exists such that

    for all and .

  • There exists a mapping with such that

  • For every compact set there exists a mapping with such that

    for all and .

First, we note that from Assumption 3.1 (i) and (ii) we immediately get


for all and .

Moreover, it is well-known that Assumption 3.1 (ii) and (iii) are sufficient to ensure the existence of a unique local solution to the initial value problem (3.1) with a local existence time , see for instance [19, Chap. I, Thm 5.3]. Here, we recall that a mapping is a (local) solution in the sense of Carathéodory to (3.1) if is absolutely continuous and satisfies


for all . Moreover, for almost all with we have

due to (3.2). Hence, by canceling from both sides of the inequality we obtain

for almost all with . After integrating this inequality from to it follows

which holds for all . An application of the Gronwall lemma (Lemma 2.2) yields


for all . In particular, since we deduce from (3.4) that is in fact the unique global solution with .

Finally, let us investigate the regularity of the solution . To this end, we define


Clearly, is a compact set, that contains the origin and the complete curve due to (3.4). Then, an application of Assumption 3.1 (iii) with yields


for all .

For arbitrary with it follows from (3.3) that

Furthermore, after inserting (3.6), we have

Then, an application of the Cauchy–Schwarz inequality yields


for all . This proves that is Hölder continuous with exponent .

4. Error analysis of the randomized backward Euler method

This section is devoted to the error analysis of the randomized backward Euler method. Our error analysis partly relies on variational methods developed in [11], that have recently been adapted to stochastic problems in [3].

In this section, we consider the following randomized version of the backward Euler method: Let denote the number of temporal steps and set as the temporal step size. For given and we obtain an equidistant partition of the interval given by , . Further, let be a family of independent and -distributed random variables on a complete probability space and let be the family of random variables given by for . Then the numerical approximation of the solution is determined by the recursion


Note that (4.1) is an implicit Runge–Kutta method with one stage and a randomized node. More precisely, in each step we apply one member of the following family of implicit Runge–Kutta methods determined by the Butcher tableaux


where the value of the parameter is determined by the random variable in the -th step.

Further, the resulting sequence consists of random variables, since we artificially inserted randomness into the numerical method. From a probabilistic point of view, is in fact a discrete time stochastic process, that takes values in and is adapted to the complete filtration . Here, is the smallest complete -algebra such that the subfamily is measurable. Note that , whenever . More precisely,


In particular, each -null set (and each subset of a -null set) is contained in every -algebra , .

Next, let us introduce the following set of square-integrable and adapted grid functions. For each this set is defined by

Take note that is an arbitrary deterministic initial value and that the condition ensures that is square-integrable as well as measurable with respect to the -algebra . First, we will show that the randomized backward Euler method (4.1) with a sufficiently large number of steps uniquely determines an element in .

We begin by proving the existence of a solution to the implicit scheme. First, we state two technical lemmas to prove the existence and measurability of a solution.

Lemma 4.1.

For let be continuous and fulfill the condition

Then there exists at least one such that .

A proof of Lemma 4.1 is found, for instance, in [12, Sec. 9.1].

Remark 4.2.

For a symmetric, positive definite Lemma 4.1 can be extended as follows. If a function where is given by

is continuous and fulfills

then there exists such that . This extension of Lemma 4.1 can be proved by exploiting that

can be rewritten as

using the transformation .

The next result is needed in order to prove the measurability of the sequence generated by the implicit numerical method (4.1). For a closely related result we refer to [13, Lem. 3.8]. The proof presented here follows an approach from [9, Prop. 1], that can easily be extended to more general situations.

Lemma 4.3.

Let be a complete sub -algebra of the -algebra , with and such that the following conditions are fulfilled.

  1. [label=()]

  2. The mapping is continuous for every .

  3. The mapping is -measurable for every .

  4. For every there exists a unique root of the function .

Define the mapping

where is the unique root of for and is arbitrary for . Then is -measurable.


Define the (multivalued) mapping

for . We first show for an arbitrary open set that the set

is an element of . To this end, first note that since is -measurable. Then, it follows that

It remains to verify the equality


It is clear that is a subset of .

To prove we consider two cases. If is a subset of then it is a null set and lies in due to the completeness of the -algebra. Else, we can assume that there exist and with . In particular, we note that the function is continuous, since . Further, observe that is an open neighborhood of and is an open neighborhood of . Since is open, the continuity of implies that the set

is an open set in with . Thus, is nonempty and open. Therefore, there exists such that . This implies and completes the proof of (4.4). Consequently, for each open set .

Next, recall that for each the image of is defined as the unique element of . Thus, the set

consists of a single element which coincides with . Therefore we obtain

which also implies for every open set due to the completeness of . From this the measurability of the mapping follows. ∎

Lemma 4.4.

Let Assumption 3.1 be satisfied. Then for each with there exists a unique solution to the implicit scheme (4.1).


The assertion is proved using an inductive argument for . Since the case is evident. Next, assuming exists, we define the set

where is the null set from Assumption 3.1. Since and the set fulfills . We define the function by


In the following we consider a fixed . Then the mapping is continuous by Assumption 3.1 (iii). Further we write

Thus, for each with this implies

Hence, by Lemma 4.1, for every there exists such that holds. This is always unique: Assume there exists and such that

holds. Then we can write for the difference

which implies . Thus, the function is -measurable in the first entry, continuous in the second and has a unique root for every . Then, Lemma 4.3 implies that the function

where is the unique root of for and for is -measurable.

It remains to prove that is finite with respect to the -norm. Using (3.2), it follows