Numerical analysis of equations of porous medium type

Robust numerical methods for local and nonlocal equations of porous medium type. Part I: Theory.


We develop the theory for monotone schemes of finite difference type for general nonlinear diffusion equations,

where is a general symmetric diffusion operator of Lévy type and is merely continuous and non-decreasing. This theory covers strongly degenerate Stefan problems and the full range of porous medium and fast diffusion equations. Examples of diffusion operators are the Laplacian , the fractional Laplacian , , discrete operators, and combinations of these. The observation that finite difference type operators are nonlocal Lévy operators, allows us to give a unified, compact, and efficient nonlocal theory for both local and nonlocal, linear and nonlinear diffusion equations. The theory includes well-posedness, stability, equicontinuity, compactness, and convergence of the methods under minimal assumptions including assumptions that lead to very irregular solutions. In other words, the schemes we introduce are robust in the sense that they converge under very unfavorable conditions. Such general results are possible because we are able to work in the largest class of weak solutions where uniqueness is known: distributional solutions in . As a byproduct of our numerical analysis, we obtain a new and general existence result for the original PDE, a result that was previously announced in [[28]]. We also present some numerical tests, but extensive numerical testing is deferred to the companion paper [[31]] along with a more detailed discussion of the numerical methods included in our theory.

Key words and phrases:
Numerical methods, finite differences, monotone methods, robust methods, convergence, stability, a priori estimates, nonlinear degenerate diffusion, porous medium equation, fast diffusion equation, Stefan problem, fractional Laplacian, Laplacian, nonlocal operators, distributional solutions, existence
2010 Mathematics Subject Classification:
65M06, 65M12, 35B30, 35K15, 35K65, 35D30, 35R09, 35R11, 76S05




1. Introduction

In this paper we develop the theory for monotone schemes of finite difference type for a large class of possibly degenerate nonlinear diffusion equations of porous medium type,


where is the solution, is a merely continuous and nondecreasing function, some right-hand side, and . The diffusion operator is given as


with local and nonlocal (anomalous) parts,


where , for and , and are the gradient and Hessian, is a characteristic function, and is a nonnegative symmetric Radon measure.

The assumptions we impose on and are so mild that many different problems can be written in the form (1.1). The assumptions on allow strongly degenerate Stefan type problems and the full range of porous medium and fast diffusion equations to be covered by (1.1). In the first case e.g. for and and in the second for any . Some physical phenomena that can be modelled by (1.1) are flow in a porous medium (oil, gas, groundwater), nonlinear heat transfer, phase transition in matter, and population dynamics. For more information and examples, we refer to Chapters 2 and 21 in [[59]] for local problems and to [[63], [51], [9], [60]] for nonlocal problems.

One important contribution of this paper is that we allow for a very large class of diffusion operators . This class coincides with the generators of the symmetric Lévy processes. Examples are Brownian motion, -stable, relativistic, CGMY, and compound Poisson processes [[5], [58], [3]], and the generators include the classical and fractional Laplacians and , (where ), relativistic Schrödinger operators , and surprisingly, also monotone numerical discretizations of . Since and may be degenerate or even identically zero, problem (1.1) can be purely local, purely nonlocal, or a combination.

A number of nonstandard and novel ideas on numerical methods for (1.1) and their analysis are presented in this paper. We will strongly use the fact that our (large) class of diffusion operators contain many of its own monotone approximations. This important observation from [[30]] is used to interpret discretizations of as nonlocal Lévy operators which again opens the door for powerful PDE techniques and a unified analysis of our schemes. Hence in this paper we consider discretizations of of the form

or equivalently with , where , the stencil points , the weights , and and . These discretizations are nonpositive in the sense that for any maximum point of , and as we will see, they include monotone finite difference quadrature approximations of . Our numerical approximations of (1.1) will then take the general form

where , , , and and are the discretization parameters in space and time respectively. By choosing in certain ways, we can recover explicit, implicit, -methods, and various explicit-implicit methods. In a simple one dimensional case,

an example of a discretization in our class is given by

Our class of schemes include both well-known discretizations and discretization that are new in context of (1.1). These new discretizations of (1.1) include higher order discretizations of the nonlocal operators, explicit schemes for fast diffusions, and various explicit-implicit schemes. See the discussion in Sections 2 and 3 and especially the companion paper [[31]] for more details.

The main contribution of this paper is to provide a uniform and rigorous analysis of such numerical schemes in this very general setting, a setting that covers local and nonlocal, linear and nonlinear, non-degenerate and degenerate, and smooth and nonsmooth problems. This novel analysis includes well-posedness, stability, equicontinuity, compactness, and -convergence results for the schemes, results which are completely new in some local and most nonlocal cases. Schemes that converge in such general circumstances are often said to be robust. Consistent numerical schemes are not robust in general, i.e. they need not always converge, or can even converge to false solutions. Such issues are seen especially in nonlinear, degenerate and/or low regularity problems. Our general results are therefore only possible because we have (i) identified a class of schemes with good properties (including monotonicity) and (ii) developed the necessary mathematical techniques for this general setting.

The main novelty of our analysis is that we are able to present the theory in a uniform, compact, and natural way. By interpreting discrete operators as nonlocal Lévy operators, and the schemes as holding in every point in space, we can use PDE type techniques for the analysis. This is possible because in recent papers [[30], [28]] we have developed a well-posedness theory for problem (1.1) which in particular allows for the general class of diffusion operators needed here. Moreover, the well-posedness holds for merely bounded distributional or very weak solutions. The fact that we can use such a weak notion of solution will simplify the analysis and make it possible to do a global theory for all the different problems (1.1) and schemes that we consider here. At this point the reader should note that if (1.1) has more regular (bounded) solutions (weak, strong, mild, or classical), then our results still apply because these solutions will coincide with the (unique) distributional solution.

The effect of the Lévy operator interpretation of the discrete operators, is that the main step of our analysis will be reduced to analyzing a semidiscrete in time approximation of (1.1) (cf. (2.5)). To do so we adapt some of the arguments of [[30]]. But since we now work with time discrete equations, we have to prove well-posedness and all a priori estimates anew. Unlike [[30]] where many results were inherited by approximation from other theories and papers, we now prove them from scratch because we have to and it is easier in most cases. The proofs of this paper are therefore more natural, elementary, and self-contained. A convergence result for the semi-discrete approximations is then obtained by adapting a compactness argument which is classical e.g. for scalar conservation laws, see Chapter 3 in [[44]]. We prove (i) uniform estimates in and and space/time translation estimates in /, (ii) compactness in via the Arzelà-Ascoli and Kolmogorov-Riesz theorems, (iii) limits of convergent subsequences are distributional solutions via stability results for (1.1), and finally (iv) full convergence of the numerical solutions by (ii), (iii), and uniqueness for (1.1).

To complete our proofs, we also need to connect the results for the semi-discrete scheme defined on the whole space with the fully discrete scheme defined on a spatial grid. We observe here that this part is easy for uniform grids where we prove an equivalence theorem under natural assumptions on discrete operators: Piecewise constant interpolants of solutions of the fully discrete scheme coincides with solutions of the corresponding semi-discrete scheme with piecewise constant initial data (see Proposition 2.10). Nonuniform grids is a very interesting case that we leave for future work.

The nonlocal approach presented in this paper gives a uniform way of representing local, nonlocal and discrete problems, different schemes and equations; compact, efficient, and easy to understand PDE type arguments that work for very different problems and schemes; new convergence results for local and nonlocal problems; and it is very natural since difference quadrature approximations are nonlocal by nature even when equation (1.1) is local.

We also mention that a consequence of our convergence and compactness theory is the existence of distributional solutions of the Cauchy problem (1.1). This result is new in this generality and was announced in [[28]]. Note that this approach to existence through semi-discrete in time approximation has superficial similarities with the Crandall-Ligget approximation and nonlinear semigroup theory [[20]]. But our approach uses much more elementary techniques and results to conclude.

Related work. In the local linear case, when and in (1.1), numerical methods and analysis can be found in undergraduate text books. In the nonlinear case there is a very large literature so we will focus only on some developments that are more relevant to this paper. For porous medium nonlinearities ( with ), there are early results on finite element and finite-difference interface tracking methods in [[56]] and [[33]] (see also [[53]]). There is extensive theory for finite volume schemes, see [[43], Section 4] and references therein for equations with locally Lipschitz . For finite element methods there is a number of results, including results for fast diffusions (), Stefan problems, convergence for strong and weak solutions, discontinuous Galerkin methods, see e.g. [[57], [40], [41], [39], [65], [55], [52]]. Note that the latter paper considers the general form of (1.1) with and provides a convergence analysis in using nonlinear semi-group theory. A number of results on finite difference methods for degenerate convection-diffusion equations also yield results for (1.1) in special cases, see e.g. [[42], [8], [50], [49]]. In particular the results of [[42], [50]] imply our convergence results for a particular scheme when is locally Lipschitz, , and solutions have a certain additional BV regularity. Finally, we mention very general results on so-called gradient schemes [[35], [36]] for doubly or triply degenerate parabolic equations. This class of equations include local porous medium type equations as a special case.

In the nonlocal case, the literature is more recent and not so extensive. For the linear case we refer somewhat arbitrarily to [[19], [45], [46], [54]] and references therein. Here we also mention [[21]] and its novel finite element plus semigroup subordination approach to discretizing . Some early results for nonlocal problems came for finite difference quadrature schemes for Bellman equations and fractional conservation laws, see [[48], [12], [6]] and [[34]]. For the latter case discontinuous Galerkin and spectral methods were later studied in [[18], [16], [64]]. The first results that include nonlinear nonlocal versions of (1.1) was probably given in [[15]]. Here convergence of finite difference quadrature schemes was proven for a convection-diffusion equation. This result is extended to more general equations and error estimates in [[17]] and a higher order discretization in [[38]]. In some cases our convergence results follow from these results (for two particular schemes, , and locally Lipschitz). However, the analysis there is different and more complicated since it involves entropy solutions and Kružkov doubling of variables arguments.

In the purely parabolic case (1.1), the behaviour of the solutions and the underlying theory is different from the convection-diffusion case (especially so in the nonlocal case, see e.g. [[23], [24], [61], [22], [62]] and [[37], [13], [1], [15], [2], [47]]). It is therefore important to develop numerical methods and analysis that are specific for this setting. The first (nonlocal) results in this direction seems to be [[27], [32]]. These papers are based on the extension method [[10]], and introduce and analyze finite difference and finite elements methods for the Fractional Porous Medium Equation. The present paper is another step in this direction, possibly the first not to use the extension method.

Outline. The assumptions, numerical schemes, and main results are given in Section 2. In Section 3 we provide many concrete examples of schemes that satisfy the assumptions of Section 2. We also show some numerical results for a nonlocal Stefan problem with non-smooth solutions. The proofs of the main results are given in Section 4, while some auxiliary results are proven in our final section, Section 5.

In the companion paper [[31]] there is a more complete discussion of the family of numerical methods. It includes more discretizations of the operator , more schemes, and many numerical examples. There we also provide proofs and explanations for why the different schemes satisfy the (technical) assumptions of this paper.

2. Main results

The main results of this paper are presented in this section. They include the definition of the numerical schemes, their consistency, monotonicity, stability, and convergence of numerical solutions towards distributional solutions of the porous medium type equation (1.1).

2.1. Assumptions and preliminaries

The assumptions on (1.1) are


Sometimes we will need stronger assumptions than (($\textup{A}_{\varphi}$)) and (2.1):

Remark 2.1.
  1. Without loss of generality, we can assume (replace by ), and when (($\textup{Lip}_{\varphi}$)) holds, that is globally Lipschitz (since is bounded). In the latter case we let denote the Lipschitz constant.

  2. Under assumption (2.1), for any and any ,

  3. Assumption (($\textup{A}_{f}$)) is equivalent to requiring , an iterated -space as in e.g. [[4]]. Note that .

Definition 2.1 (Distributional solution).

Let and . Then is a distributional (or very weak) solution of (1.1) if for all , and


Note that if e.g.  and continuous. Distributional solutions are unique in .

Theorem 2.2 (Theorem 3.1 [[28]]).

Assume (($\textup{A}_{\varphi}$)), (($\textup{A}_{f}$)), (($\textup{A}_{u_{0}}$)), and (2.1). Then there is at most one distributional solution of (1.1) such that .

2.2. Numerical schemes without spatial grids

Let be a nonuniform grid in time such that . Let , and denote time steps by


For , , and , we define


and we define our time discretized scheme as


where, formally, , , and

Typically , but when the is not Lipschitz, we have to approximate it by a Lipschitz to get a monotone explicit method [[31]]. Let . Depending on the choice of and , we can then get many different schemes:

  1. Discretizing separately the different parts of the operator

    e.g. the local, singular nonlocal, and bounded nonlocal parts, corresponds to different choices for and . Typically choices here are finite difference and numerical quadrature methods, see Section 3 for several examples.

  2. Explicit schemes (), implicit schemes (), or combinations like Crank-Nicholson (), follow by the choices

  3. Combinations of type (1) and (2) schemes, e.g. implicit discretization of the unbounded part of and explicit discretization of the bounded part.

Finally, we mention that our schemes and results may easily be extended to handle any finite number of and .

Definition 2.2 (Consistency).

We say that the scheme (2.5) is consistent if

  1. in as for all ,

  2. locally uniformly as ,

  3. in for ,

for satisfying (($\textup{A}_{\varphi}$)), (2.1), and , , of the form (1.2)–(1.4).

We will focus on discrete operators , in the following class:

Definition 2.3.

An operator is said to be

  1. in the class (2.1) if for a measure satisfying (2.1); and

  2. discrete if

    for and such that .

  3. is called the stencil and the weights of the discretization.

All operators in the class (2.1) are nonpositive operators, in particular they are integral or quadrature operators with positive weights. The results presented in this section hold for any operator in the class (2.1). However, in practice, when dealing with numerical schemes, the operators will additionally be discrete. Moreover, when the scheme (2.5) has an explicit part, that is, , we need to assume that satisfies (($\textup{Lip}_{\varphi}$)) and impose the following CFL-type condition to have a monotone scheme:


where we recall that is the Lipschitz constant of (see Remark 2.1). Note that this condition is always satisfied for an implicit method where . The typical assumptions on the scheme (2.5) are then:

Theorem 2.3 (Existence and uniqueness).

Assume (($\textup{A}_{\textup{NS}}$)), (($\textup{A}_{f}$)), and (($\textup{A}_{u_{0}}$)). Then there exists a unique a.e.-solution of the scheme (2.5).

Theorem 2.4 (A priori estimates).

Assume (($\textup{A}_{\textup{NS}}$)), (($\textup{A}_{f}$)), and (($\textup{A}_{u_{0}}$)). Let be solutions of the scheme (2.5) with data and . Then:

  1. (Monotone) If and , then .

  2. (-stable) .

  3. (-stable) .

  4. (Conservative) If additionally satisfies (($\textup{Lip}_{\varphi}$)),

Remark 2.5.

By (b), (c), and interpolation, the scheme is -stable for .

The scheme is also -contractive and equicontinuous in time. Combined, these two results imply time-space equicontinuity and compactness of the scheme, a key step in our proof of convergence.

Theorem 2.6 (-contractive).

Under the assumptions of Theorem 2.4,

For the equicontinuity in time we need a modulus of continuity:



for some constant , , compact with Lebesgue measure , and . In view of (2.6), we also need to assume a uniform Lévy condition on the approximations,

Remark 2.7.

Condition (($\textup{A}_{\nu^{h}}$)) is in general very easy to check. For example it follows from pointwise consistency of as we will see in [[31]].

Theorem 2.8 (Equicontinuity in time).

Assume (($\textup{A}_{f}$)) and (($\textup{A}_{u_{0}}$)), and let (2.5) be a consistent scheme satisfying (($\textup{A}_{\textup{NS}}$)) and (($\textup{A}_{\nu^{h}}$)). Then, for all such that and all compact sets ,

The main result regarding convergence of numerical schemes without spatial grids will be presented in a continuous in time and space framework. For that reason, let us define the time interpolant as

Theorem 2.9 (Convergence).

Assume (($\textup{A}_{f}$)), (($\textup{A}_{u_{0}}$)), , and for all , let be the solution of a consistent scheme (2.5) satisfying (($\textup{A}_{\textup{NS}}$)) and (($\textup{A}_{\nu^{h}}$)). Then there exists a unique distributional solution of (1.1) and

Convergence of subsequences follows from compactness and full convergence follows from stability and uniqueness of the limit problem (1.1). The detailed proofs of Theorems 2.3, 2.4, and 2.62.9 can be found in Sections 4.14.3.

2.3. Numerical schemes on uniform spatial grids

To get computable schemes, we need to introduce spatial grids. For simplicity we restrict to uniform grids. Since our discrete operators have weights and stencils not depending on the position , all results then become direct consequences of the results in Section 2.2.

Let , , and be the uniform spatial grid

Note that any discrete (2.1) class operator with stencil , has a well-defined restriction defined by

and all . Using such restricted discrete operators, we get the following well-defined numerical discretization of (1.1) on the space-time grid ,


where and are the cell averages of the - functions and :


The function and the solution are functions on , and we define their piecewise constant interpolations in space as


The next proposition shows that solutions of the scheme (2.5) with piecewise constant initial data are solutions of the fully discrete scheme (2.8) and vice versa.

Proposition 2.10.

Assume (($\textup{A}_{f}$)), (($\textup{A}_{u_{0}}$)), let , be defined by (2.9) and , by (2.10), and let , be class (2.1) discrete operators with stencils .

  1. If is an a.e. solution of (2.5) with data and , then (a version of) is constant on the cells for all , and is a solution of (2.8) with data and .

  2. If is a solution of (2.8) with data and , then defined in (2.10) is a piecewise constant solution of (2.5) with data and .

In view of this result, the scheme on the spatial grid (2.8) will inherit the results for the scheme (2.5) given in Theorems 2.3, 2.4, 2.62.9.

Theorem 2.11.

Assume (($\textup{A}_{\textup{NS}}$)), (($\textup{A}_{f}$)), (($\textup{A}_{u_{0}}$)), and the stencils .

  1. (Exists/unique) There exists a unique solution of (2.8) such that

Let be solutions of the scheme (2.8) with data and respectively.

  1. (Monotone) If and , then .

  2. (-stable) .

  3. (-stable) .

  4. (Conservative) If satisfy (($\textup{Lip}_{\varphi}$)),

  5. (-contractive)

  6. (Equicontinuity in time) If (($\textup{A}_{\nu^{h}}$)) hold, then for all compact sets ,

Assume in addition that , and for all , let be the solution of a consistent scheme (2.8) satisfying (($\textup{A}_{\textup{NS}}$)) and (($\textup{A}_{\nu^{h}}$)).

  1. (Convergence) There exists a unique distributional solution of (1.1) such that for all compact sets ,

Remark 2.12.

All these results can be formulated in terms of the space interpolant , e.g. the -contraction in part (f) then becomes

Moreover, convergence in (h) can be stated in terms of space-time interpolants as

since for any function on .

The proofs of the above results can be found in Section 4.4.

2.4. Well-posedness for bounded distributional solutions

Theorem 2.9 implies the existence of bounded distributional solutions solutions of (1.1), and uniqueness has been proved in [[28]]:

Theorem 2.13 (Existence and uniqueness).

Assume (($\textup{A}_{\varphi}$)), (($\textup{A}_{f}$)), (($\textup{A}_{u_{0}}$)), and (2.1). Then there exists a unique distributional solution of (1.1) such that

Another consequence of Theorem 2.9 is that most of the a priori results in Theorems 2.4, 2.6, 2.8 will be inherited by the solution of (1.1).

Proposition 2.14 (A priori estimates).

Assume (($\textup{A}_{\varphi}$)) and (2.1). Let be the distributional solutions of (1.1) corresponding to and satisfying (($\textup{A}_{u_{0}}$)) and (($\textup{A}_{f}$)) respectively. Then, for every :

  1. (Comparison) If and , then .

  2. (-bound) .

  3. (-bound) .

  4. (-contraction)

  5. (Time regularity) For every and every compact set ,