MLMC diffusion

Multilevel Monte Carlo for random degenerate scalar convection diffusion equation

U. Koley
Institut für Mathematik,
Julius-Maximilians-Universität Würzburg,
Campus Hubland Nord, Emil-Fischer-Strasse 30,
97074, Würzburg, Germany.
N. H. Risebro
Centre of Mathematics for Applications (CMA)
University of Oslo
P.O. Box 1053, Blindern
N–0316 Oslo, Norway
Ch. Schwab
Seminar for Applied Mathematics
HG G. 57.1,
Rämistrasse 101, Zürich, Switzerland.
 and  F. Weber
Centre of Mathematics for Applications (CMA)
University of Oslo
P.O. Box 1053, Blindern
N–0316 Oslo, Norway
July 5, 2019

We consider the numerical solution of scalar, nonlinear degenerate convection-diffusion problems with random diffusion coefficient and with random flux functions. Building on recent results on the existence, uniqueness and continuous dependence of weak solutions on data in the deterministic case, we develop a definition of random entropy solution. We establish existence, uniqueness, measurability and integrability results for these random entropy solutions, generalizing [28, 29] to possibly degenerate hyperbolic-parabolic problems with random data. We next address the numerical approximation of random entropy solutions, specifically the approximation of the deterministic first and second order statistics. To this end, we consider explicit and implicit time discretization and Finite Difference methods in space, and single as well as Multi-Level Monte-Carlo methods to sample the statistics. We establish convergence rate estimates with respect to the discretization parameters, as well as with respect to the overall work, indicating substantial gains in efficiency are afforded under realistic regularity assumptions by the use of the Multi-Level Monte-Carlo method. Numerical experiments are presented which confirm the theoretical convergence estimates.

Research supported in part by ERC AdG247277 STAHDPDE



1. Introduction

Many problems in physics and engineering are modeled by nonlinear, possibly strongly degenerate, convection diffusion equation. The Cauchy problem for such equations takes the form


where with fixed, is the unknown function, the flux function, and the nonlinear diffusion. Regarding this, the basic assumption is that for all . When (1.1) is nondegenerate, i.e., , it is well known that (1.1) admits a unique classical solution [34]. This contrasts with the degenerate case where may vanish for some values of . A simple example of a degenerate equation is the porous medium equation

which degenerates at . This equation has served as a simple model to describe processes involving fluid flow, heat transfer or diffusion. Examples of applications are in the description of the flow of an isentropic gas through a porous medium, modelled by Leibenzon [27] and Muskat [32] around 1930, in the study of groundwater flow by Boussisnesq in 1903 [3] or in heat radiation in plasmas, Zel’dovich and collaborators around 1950, [38]. In general, a manifestation of the degeneracy in (1.1) is the finite speed of propagation of disturbances. If , and if at some fixed time the solution has compact support, then it will continue to have compact support for all later times.

By the term “strongly degenerate” we mean that there is an open interval such that if is in this interval. Hence, the class of equations under consideration is very large and contains the heat equation, the porous medium equation and scalar conservation laws. Independently of the smoothness of the initial data, due to the degeneracy of the diffusion, singularities may form in the solution . Therefore we consider weak solutions which are defined as follows.

Definition 1.1.

Set . A function

is a weak solution of the initial value problem (1.1) if it satisfies:

  1. .

  2. For all test functions


In view of the existence theory, the condition 1 is natural, and thanks to this we can replace (1.2) by

If is constant on a whole interval, then weak solutions are not uniquely determined by their initial data, and one must impose an additional entropy condition to single out the physically relevant solution. A weak solution satisfies the entropy condition if


for all convex, twice differentiable functions , where and are defined by

Via a standard limiting argument this implies that (1.3) holds for the Kružkov entropies for all constants . We call a weak solution satisfying the entropy condition an entropy solution.

For scalar conservation laws, the entropy framework (usually called entropy conditions) was introduced by Kružkov [25] and Vol’pert [36], while for degenerate parabolic equations entropy solution were first considered by Vol’pert and Hudajev [37]. Uniqueness of entropy solutions to (1.1) was first proved by Carrillo [4].

Over the years, there has been a growing interest in numerical approximation of entropy solutions to degenerate parabolic equations. Finite difference and finite volume schemes for degenerate equations were analysed by Evje and Karlsen [12, 11, 10, 13] (using upwind difference schemes), Holden et al[19, 20] (using operator splitting methods), Kurganov and Tadmor [26] (central difference schemes), Bouchut et al[2] (kinetic BGK schemes), Afif and Amaziane [1] and Ohlberger, Gallouët et al[33, 15, 16] (finite volume methods), Cockburn and Shu [7] (discontinuous Galerkin methods) and Karlsen and Risebro [24, 23] (monotone difference schemes). Many of the above papers show that the approximate solutions converge to the unique entropy solution as the discretization parameter vanishes. Rigorous estimates of the convergence rate of finite volume schemes for degenerate parabolic equations were proved in [21] (1-d) and [22] (multi-d).

This classical paradigm for designing efficient numerical schemes assumes that data for (1.1), i.e., initial data , convective flux and diffusive flux are known exactly.

In many situations of practical interest, however, these data are not known exactly due to inherent uncertainty in modelling and measurements of physical parameters such as, for example, the specific heats in the equation of state for compressible gases, or the relative permeabilities in models of multi-phase flow in porous media. Often, the initial data are known only up to certain statistical quantities of interest like the mean, variance, higher moments, and in some cases, the law of the stochastic initial data. In such cases, a mathematical formulation of (1.1) is required which allows for random data. The problem of random initial data was considered in [29], and the existence and uniqueness of a random entropy solution was shown, and a convergence analysis for MLMC FV discretizations was given. The MLMC discretization of balance laws with random source terms was investigated in [31].

In [29] a mathematical framework was developed for scalar conservation laws with random initial data. This framework was extended to include random flux functions in [28].

The aim of this paper is to extend this mathematical framework to include degenerate convection diffusion equations with random convective and diffusive flux functions with possibly correlated random perturbations. Its outline is as follows. In Section 2 we review notions from probability and from random variables taking values in separable Banach spaces. Section 3.1 is devoted to a review of convergence rates from [21, 22] on convergence rates for scalar, degenerate deterministic convection-diffusion problems. Particular attention is paid to the definition of entropy solutions and to existence-, uniqueness- and continuous dependence results, and to the definition of the random entropy solutions, and to sufficient conditions ensuring their measurability and integrability. In Section 4, we then address the discretization. First, again reviewing convergence rates of FD schemes for the deterministic case from [21, 22], which we then extend to Monte-Carlo as well as Multi-Level Monte-Carlo versions for the degenerate convection-diffusion problem with random coefficients and flux functions. The final Section 5 is then devoted to numerical experiments which confirm the theoretical convergence estimates and, in fact, indicate that they probably are pessimistic, at least in the particular test problems considered.

2. Preliminaries from Probability

We use the concept of random variables taking values in function spaces. To this end, we recapitulate basic concepts from [8, Chapter 1].

Let be a measurable space, with denoting the set of all elementary events, and a -algebra of all possible events in our probability model. If denotes a second measurable space, then an -valued random variable (or random variable taking values in ) is any mapping such that the set : for any , i.e., such that is a -measurable mapping from into .

Assume now that is a metric space; with the Borel -field , is a measurable space and we shall always assume that -valued random variables will be measurable. If is a separable Banach space with norm and (topological) dual , then is the smallest -field of subsets of containing all sets

Hence if is a separable Banach space, is an -valued random variable iff for every , is an -valued random variable. Moreover, by [8, Lemma 1.5, p.19] the norm is a measurable mapping.

The random variable is called Bochner integrable if, for any probability measure on the measurable space ,

A probability measure on is any -additive set function from into such that , and the measure space is called probability space. We shall assume that is complete.

If is a random variable, denotes the law of under , i.e.,

The image measure on is called law or distribution of .

A random variable taking values in is called simple if it can take only finitely many values, i.e., if it has the explicit form (with the indicator function of )

We set, for simple random variables taking values in and for any ,

By density, for such , and all ,

For any random variable which is Bochner integrable, there exists a sequence of simple random variables such that, for all , as . Therefore, (2) and (2) extend in the usual fashion by continuity to any -valued random variable. We denote the integral

by (“expectation” of ). We shall require for Bochner spaces of -summable random variables taking values in the Banach space . By we denote the set of all (equivalence classes of) integrable, -valued random variables , equipped with the norm

More generally, for , we define as the set of -summable random variables taking values in and equip it with norm

For , we denote by the set of all -valued random variables which are essentially bounded. This set is a Banach space equipped with the norm

If and , , we write . Note that for any separable Banach space , and for any ,

In the following, we will be interested in random variables , , mapping from some probability space into subsets of the Banach spaces , , equipped with the Borel -algebra , where , for a closed and bounded interval , , and , . On , we choose the norm

on , we will use

and on ,

Furthermore, we will need the following special case of the fact that a continuous mapping is measurable:

Lemma 2.1.

Let be a measurable space, be Banach spaces equipped with the Borel -algebras , and let , be a random variable and be a continuous mapping, that is for ,

where is a continuous function which satisfies and which increases monotonically in .

Then the mapping is - measurable, i.e., it is a -valued random variable.


We have to show that for any , . Since and by the assumption that is a random variable for every , this amounts to showing that for any . Since the Borel -algebra is generated by the open sets and the inverse image of a mapping has the two fundamental properties

for a countable index set and any countable collection of sets , it is enough to verify this for an arbitrary open, nonempty set . That is an open set if is open, then follows by the continuity of . ∎

3. Degenerate Convection Diffusion Equation with Random Diffusive Flux

We develop a theory of random entropy solutions for degenerate convection diffusion equation with a class of random flux flunctions, proving in particular the existence and uniqueness of a random entropy solution. To this end, we first review classical results on degenerate convection diffusion equation with deterministic data.

3.1. Deterministic Scalar Degenerate Convection Diffusion Equation

We consider the Cauchy problem for degenerate convection diffusion equation of the form


3.2. Entropy Solutions

It is well-known that if is Lipschitz continuous and , then the deterministic Cauchy problem (3.1) admits, for each , a unique entropy solution (see, e.g., [17, 35, 9]). Moreover, for every , and several properties of the (nonlinear) data-to-solution operator

will be crucial for our subsequent development. To state these properties of , following [12] we introduce the set of admissible initial data


We collect next fundamental results regarding the entropy solution of (3.1) in the following theorem, for a proof see [37, 5],

Theorem 3.1.

Let and be locally Lipschitz continuous functions. Then

  • For every , the initial value problem (3.1) admits a unique BV entropy weak solution .

  • For every , the (nonlinear) data-to-solution map given by


    • For fixed , is a (contractive) Lipschitz map, i.e.,

    • For every ,

    • Lipschitz continuity in time: For any , ,


Point 1) of Theorem 3.1 is proved in [37] or [5, Thm 1.1], (3.3), (3.5) also follow from [5, Thm 1.1], (3.4) was proved in [5, Thm 1.2], and (3.6), (3.7), (3.8) were proved in [37]. In our convergence analysis of MC-FD discretizations of degenerate convection diffusion equation with random fluxes, we will need the following result regarding continuous dependence of with respect to and ([6, Thm. 3])

Theorem 3.2.

Assume , , and , , , with .

Then the unique entropy solutions and of (3.1) with initial data , , convective flux functions and and with diffusive flux functions and satisfy the Kružkov entropy conditions, and the à priori continuity estimate


where and . The above estimate holds for every .

Remark 3.3.

Using that for nonnegative numbers , ,

it follows from (3.9) that under the assumptions of Theorem 3.2,

hence the mapping , is continuous as a mapping between Banach spaces if restricted to initial data in and satisfying . Moreover, since for with the above properties and bounded derivatives it holds


it follows that is a continuous mapping between the separable Banach spaces and if restricted to initial data in the set .

3.3. Random Entropy Solutions

We are interested in the case where the initial data , the convective flux function and the diffusive flux function in (1.1) are uncertain. Existence and uniqueness for random initial data and random flux for was proved in [29, 28]. Based on Theorem 3.1, we will now formulate (1.1) for random initial data , random convective flux and random diffusive flux . To this end, we denote a probability space and consider random variables


where and , . In order to establish the appropriate framework for the random degenerate convection diffusion equation (3.19), we will restrict ourselves to random data which satisfy -a.s. the following assumptions:


Since and are separable, (3.11) is well defined. Moreover, by Lemma [8, Lemma 1.5, p.19] each of the expressions on the left hand sides of (3.12) - (3.17) is a random variable and we may impose for the -th moment condition:


where the Bochner spaces with respect to the probability measure are defined in Section 2. Then we are interested in random solutions of the random degenerate convection diffusion equation

Definition 3.4.

A random field , i.e., a measurable mapping from to , is called a random entropy solution of (3.1) with random initial data , flux function and diffusive flux satisfying (3.11) and (3.12) – (3.18) for some , if it satisfies:

  • Weak solution: for -a.e. , satisfies

    for all test functions .

  • Entropy condition: For any pair consisting of a (deterministic) entropy and (stochastic) entropy flux and i.e., and are functions such that is convex and such that , and for -a.s. , satisfies the following integral identity:

    for all test functions .

We state the following theorem regarding the random entropy solution of (3.19):

Theorem 3.5.

Consider the degenerate convection diffusion equation (3.1) with random initial data , flux function and random diffusion operator , as in (3.11), and satisfying (3.12) – (3.17) and the -th moment condition (3.18) for some integer . Then there exists a unique random entropy solution which is “pathwise”, i.e., for , described in terms of a nonlinear mapping , depending only on the random flux and diffusion,

such that for every and for every


and such that we have -a.s.


and, with for as in (3.12),


For , we define, motivated by Theorem 3.1, for -a.e.  a random function by


By the properties of the solution mapping , see Theorem 3.1, the random field defined in (3.26) is well defined; for -a.e. , is a weak entropy solution of the degenerate diffusion equation (3.1). Moreover, we obtain from Theorem 3.1 that -a.s. all bounds (3.21)–(3.24) hold, with assumption (3.12) also (3.25). The measurability of the mapping , follows from Lemma 2.1, (3.10) and the assumption that the mapping is a random variable. Finally, (3.20) follows from (3.18) together with (3.5) in Theorem 3.1. ∎

Theorem 3.5 generalizes the existence of random entropy solutions for random initial data from [29] and random convective flux function [28]. It ensures the existence of a unique random entropy solution with finite -th moments provided that for some .

Remark 3.6.

All existence and continuous dependence results stated so far are formulated for the deterministic Cauchy problem (3.1). By the ‘usual arguments’, verbatim the same results will also hold for solutions defined in a bounded, axiparallel domain , provided that periodic boundary conditions in each coordinate are enforced on the weak solutions. Weak solutions for these periodic problems cannot coincide with weak solutions of the Cauchy problem (3.1) since the -periodic extension of these solutions belongs to , but does not belong to .

4. Numerical approximation of random degenerate convection diffusion equation

We wish to compute various quantities of interest, such as the expectation and higher order moments, of the solution to the random degenerate diffusion equation (3.19). We choose to split the approximation into two steps: On one hand, we need to approximate in the stochastic domain and on the other hand, since in general exact solutions to (1.1) are not available, we need an approximation in the physical domain . In this paper, we will consider a Multilevel Monte Carlo Finite Difference Method (MLMC-FDM), that is, a combination of the multilevel Monte Carlo method with a deterministic finite difference discretization. We will briefly review the two methods and mention some relevant results in the following sections.

4.1. Monte Carlo method

We view the Monte Carlo method as a “discretization” of the random degenerate diffusion equation data , , with respect to . We assume that satisfying in addition (3.12)–(3.17). We also assume (3.18), i.e., the existence of -th moments of for some , to be specified later. We shall be interested in the statistical estimation of the first and higher moments of , i.e., . For , . The Monte Carlo (MC) approximation of is defined as follows: Given independent, identically distributed samples , , of initial data, flux function and diffusion, the MC estimate of at time is given by


where denote the unique entropy solutions of the Cauchy problems (1.1) with initial data , flux function and diffusion operator . Since

we have for every and for every , by (3.5),

Using the i.i.d. property of the samples and therefore of , and the linearity of the expectation , we obtain the bound

As the sample size , the sample averages (4.1) converge and the convergence result from [29, 28] holds as well:

Theorem 4.1.

Assume that in (3.19) the random variable as in (3.11) satisfies (3.12) and , a.s.  and

Then the MC estimates in (4.1) converge as , to and, for any , , and we have the bound


The proof of this result proceeds completely analogous to the proof of [29, Thm. 4.1], using the measurability and square integrability (3.20) (with ) of Theorem 3.5.

So far, we addressed the MC estimation of the mean field or first moment. A similar result holds for the MC sample averages of the -th moment .

Theorem 4.2.

Consider the random degenerate advection diffusion equation (3.19) with random data as in (3.11) and satisfying (3.12) and