## Abstract

In an earlier study (Zhang & Shu 2010*b* *J. Comput. Phys.* **229**, 3091–3120 (doi:10.1016/j.jcp.2009.12.030)), genuinely high-order accurate finite volume and discontinuous Galerkin schemes satisfying a strict maximum principle for scalar conservation laws were developed. The main advantages of such schemes are their provable high-order accuracy and their easiness for generalization to multi-dimensions for arbitrarily high-order schemes on structured and unstructured meshes. The same idea can be used to construct high-order schemes preserving the positivity of certain physical quantities, such as density and pressure for compressible Euler equations, water height for shallow water equations and density for Vlasov–Boltzmann transport equations. These schemes have been applied in computational fluid dynamics, computational astronomy and astrophysics, plasma simulation, population models and traffic flow models. In this paper, we first review the main ideas of these maximum-principle-satisfying and positivity-preserving high-order schemes, then present a simpler implementation which will result in a significant reduction of computational cost especially for weighted essentially non-oscillatory finite-volume schemes.

## 1. Introduction

An important property of the unique entropy solution to the scalar conservation law
1.1is that it satisfies a strict maximum principle, namely, if , , then *u*(**x**,*t*)∈[*m*,*M*] for any **x** and *t*. This property is also naturally desired for numerical schemes solving equation (1.1), since numerical solutions outside of [*m*,*M*] often are meaningless physically, such as negative density, or negative percentage or percentage larger than one for a component in a multi-component mixture.

One of the main difficulties in solving equation (1.1) is that the solution may contain discontinuities even if the initial condition is smooth. Moreover, the weak solutions of equation (1.1) may not be unique. Therefore, the nonlinear stability and convergence to the unique entropy solution must be considered for the numerical schemes. Total-variation (TV) stable functions form a compact space, so a conservative TV-stable scheme will produce a subsequence converging to a weak solution by the Lax–Wendroff Theorem. The E-schemes including Godunov, Lax–Friedrichs and Engquist–Osher methods satisfy an entropy inequality and are total-variation-diminishing (TVD) thus satisfying maximum principle. However, E-schemes are at most first-order accurate. In fact, any TVD scheme in the sense of measuring the variation of grid point values or cell averages will be at most first-order accurate around smooth extrema, see Osher & Chakravarthy (1984), although TVD schemes can be designed for any formal order of accuracy for smooth monotone solutions, e.g. the high-resolution schemes.

For conventional maximum-principle-satisfying finite-difference schemes, the solution is at most second-order accurate; for instance, only the second-order central scheme was proved to satisfy the maximum principle in Jiang & Tadmor (1998). This fact has a simple proof owing to Ami Harten. For simplicity, we consider a finite difference scheme, namely is the numerical solution approximating the point values *u*(*x*_{j},*t*^{n}) of the exact solution, where *n* is the time step and *j* denotes the spatial grid index. Assume the scheme satisfies the maximum principle
1.2Consider the linear convection equation with periodic boundary conditions. Set the grid as *x*_{j}=(*j*−1/2)Δ*x* where Δ*x*=1/*N* and *N* is a multiple of 4. The numerical initial value is . Without loss of generality, assume Δ*t*=1/2Δ*x*. At the grid point *j*=*N*/4+1 and *t*=Δ*t*, the exact solution is and the numerical solution is by inequality (1.2). The error of the scheme at the grid point *j*=*N*/4+1 after one time step is equal to . That is, even after one time step, the scheme is already at most second-order accurate. A similar proof also works for finite volume schemes where the numerical solution approximates cell averages of the exact solution.

The simple derivation above implies that equation (1.2) is too restrictive for the scheme to be higher than second-order accurate. A heuristic point of view to understand the restriction is, some high-order information of the exact solution is lost as we only measure the total variation or the maximum at the grid points or in cell averages. To overcome this difficulty, Sanders proposed to measure the total variation of the reconstructed polynomials and he succeeded in designing a third-order TVD scheme for one-dimensional scalar conservation laws in Sanders (1988), which has been extended to higher order in Zhang & Shu (2010*a*). But it is very difficult to generalize Sanders’ scheme to higher space dimension. By measuring the maximum of the reconstructed polynomial, Liu & Osher (1996) constructed a third-order non-oscillatory scheme, which could be generalized to two space dimensions. However, it could be proved maximum-principle-satisfying only for the linear equation. The key step of maximum-principle-satisfying high-order schemes above is a high-order accurate time evolution, which preserves the maximum principle. The exact time evolution satisfies this property and it was used in Sanders (1988), Zhang & Shu (2010*a*) and Liu & Osher (1996). Unfortunately, it is very difficult, if not impossible, to implement such exact time evolution for multi-dimensional nonlinear scalar equations or systems of conservation laws.

Successful high-order numerical schemes for hyperbolic conservation laws include, among others, the Runge–Kutta discontinuous Galerkin (RKDG) method with a total-variation bounded (TVB) limiter, e.g. Cockburn & Shu (1989), the essentially non-oscillatory (ENO) finite volume and finite difference schemes, e.g. Harten *et al.* (1987) and Shu & Osher (1988), and the weighted ENO (WENO) finite volume and finite difference schemes, e.g. Liu *et al.* (1994) and Jiang & Shu (1996). Although these schemes are nonlinearly stable in numerical experiments and some of them can be proved to be TV stable, they do not, in general, satisfy a strict maximum principle. In Zhang & Shu (2010*b*), we proved a sufficient condition for the cell averages of the numerical solutions in a high-order finite volume or a discontinuous Galerkin (DG) scheme with the strong stability preserving (SSP) time discretization, e.g. Shu & Osher (1988) and Shu (1988), to be bounded in [*m*,*M*] for equation (1.1). We have also proved that, with a simple scaling limiter introduced in Liu & Osher (1996), this sufficient condition can be enforced and not only the cell averages but also the numerical solution itself can be guaranteed to stay in [*m*,*M*] without destroying accuracy for smooth solutions. In other words, we have constructed a high-order scheme by adding a simple limiter to a finite volume WENO/ENO or RKDG scheme and it can be proved to be high-order accurate and maximum-principle-satisfying. This was the first time that genuinely high-order schemes are obtained satisfying a strict maximum principle especially for multi-dimensional nonlinear problems.

For hyperbolic conservation law systems, the entropy solutions in general do not satisfy the maximum principle. We consider the positivity of some important quantities instead. For instance, density and pressure in compressible Euler equations, and water height in shallow water equations should be non-negative physically. In practice, failure of preserving positivity of such quantities may cause blow-ups of the computation because the linearized system may become ill-posed. From the point of view of stability, it is highly desired to design schemes, which can be proved to be positivity-preserving. Most commonly used high-order numerical schemes for solving hyperbolic conservation law systems do not in general satisfy such properties automatically. It is very difficult to design a conservative high-order accurate scheme preserving the positivity. In Zhang & Shu (2010*c*, 2011) and Zhang *et al.* (in press), we have generalized the maximum-principle-satisfying techniques to construct conservative positivity-preserving high-order finite volume and DG schemes for compressible Euler equations, which could be regarded as an extension of the positivity-preserving schemes in Perthame & Shu (1996).

In this paper, we first review the general framework to construct maximum-principle-satisfying and positivity-preserving schemes of arbitrarily high-order accuracy. In §2, we illustrate the main ideas in the context of scalar conservation laws. We then discuss generalizations of this idea to other equations and systems in §§3 and 4. In §5, we propose a more efficient implementation of the framework for WENO finite volume schemes, and provide numerical examples to demonstrate their performance. Concluding remarks are given in §6.

## 2. Maximum-principle-satisfying high-order schemes for scalar conservation laws

### (a) One-dimensional scalar conservation laws

We consider the one-dimensional version of equation (1.1) in this section: 2.1

#### (i) The first-order schemes

It is well known that a first-order monotone scheme solving equation (2.1) satisfies the strict maximum principle. A first-order monotone scheme has the form:
2.2where *λ*=Δ*t*/Δ*x* with Δ*t* and Δ*x* being the temporal and spatial mesh sizes (we assume uniform mesh size for the structured mesh cases in this paper for simplicity in presentation, however, the methodology does not have a uniform or smooth mesh restriction), and is a monotone flux, namely it is Lipschitz continuous in both arguments, non-decreasing (henceforth referred to as increasing with a slight abuse of the terminology) in the first argument and non-increasing (henceforth referred to as decreasing) in the second argument, and consistent . Under suitable Courant–Friedrichs–Lewy (CFL) conditions, typically of the form:
2.3for example Lax–Friedrichs scheme and Godunov scheme, one can prove that the function *H*_{λ}(*a*,*b*,*c*) is increasing in all three arguments, and consistency implies *H*_{λ}(*a*,*a*,*a*)=*a*. We, therefore, immediately have the strict maximum principle
provided .

#### (ii) High-order spatial discretization

Now consider high-order finite volume or DG methods, for example, the WENO finite volume method in Liu *et al.* (1994) and the DG method in Cockburn & Shu (1989) solving equation (2.1). We only discuss the Euler forward temporal discretization in this subsection and leave higher order temporal discretization to section §2*a*(iv). The finite volume method or the scheme satisfied by the cell averages in the DG method discretization can be written as
2.4where is the approximation to the cell averages of *u*(*x*,*t*) in the cell *I*_{j}=[*x*_{j−1/2},*x*_{j+1/2}] at time level *n*, is again a monotone flux, and *u*^{−}_{j+1/2}, *u*^{+}_{j+1/2} are the high-order approximations of the nodal values *u*(*x*_{j+1/2},*t*^{n}) within the cells *I*_{j} and *I*_{j+1}, respectively. These values are either reconstructed from the cell averages in a finite volume method or read directly from the evolved polynomials in a DG method. We assume that there is a polynomial *p*_{j}(*x*) (either reconstructed in a finite volume method or evolved in a DG method) with degree *k*, where *k*≥1, defined on *I*_{j} such that is the cell average of *p*_{j}(*x*) on *I*_{j}, *u*^{+}_{j−1/2}=*p*_{j}(*x*_{j−1/2}) and *u*^{−}_{j+1/2}=*p*_{j}(*x*_{j+1/2}).

Given a scheme in the form of equation (2.4), assuming for all *j*, we would like to derive some sufficient conditions to ensure . A very natural first attempt is to see if there is a restriction on *λ* such that, if all five arguments of *G* are in [*m*,*M*]
then we could prove . Unfortunately, one can easily build counter examples to show that this cannot be always true. The problem is that the function *G*_{λ}(*a*,*b*,*c*,*d*,*e*) in equation (2.4) is only monotonically increasing in the first, third and fourth arguments and is monotonically decreasing in the other two arguments. Hence the strategy to prove maximum principle for first-order monotone schemes cannot be repeated here. In the literature, many attempts have been made to further limit the four arguments *u*^{−}_{j+1/2},*u*^{+}_{j+1/2},*u*^{−}_{j−1/2},*u*^{+}_{j−1/2} (remember the cell average cannot be changed owing to conservation) in the arguments of *G*_{λ} in equation (2.4) to guarantee that . However, these limiters always kill accuracy near smooth extrema.

Our approach follows a different strategy. We consider an *N*-point Legendre Gauss–Lobatto quadrature rule on the interval *I*_{j}=[*x*_{j−1/2},*x*_{j+1/2}], which is exact for the integral of polynomials of degree up to 2*N*−3. We denote these quadrature points on *I*_{j} as
2.5Let be the quadrature weights for the interval such that . Choose *N* to be the smallest integer satisfying 2*N*−3≥*k*, then
2.6We then have the following theorem. We assume that the monotone flux corresponds to a monotone scheme (2.2) under the CFL condition (2.3).

### Theorem 2.1

*Consider a finite volume scheme or the scheme satisfied by the cell averages of the DG method (2.4), associated with the approximation polynomials p*_{j}*(x) of degree k (either reconstruction or DG polynomials) in the sense that* *u*^{+}_{j−1/2}*=p*_{j}*(x*_{j−1/2}*) and u*^{−}_{j+1/2}*=p*_{j}*(x*_{j+1/2}*). If u*^{−}_{j−1/2}*, u*^{+}_{j+1/2} *and* *(α=1,…,N) are all in the range [m,M], then* *under the CFL condition
*2.7

### Proof.

With equation (2.6), by adding and subtracting , the scheme (2.4) can be rewritten as
2.8Noticing that and are monotone under the CFL condition (2.7), we can see from equation (2.8) that is a monotonically increasing function of all the arguments involved, namely *u*^{−}_{j−1/2}, *u*^{+}_{j+1/2} and for 1≤*j*≤*N*. The same proof for the first-order monotone scheme now applies to imply . ■

### Remark 2.2

We recall that the CFL condition for linear stability for the DG scheme using polynomial of degree *k* is *λa*≤1/(2*k*+1) in Cockburn & Shu (1989), which is close to the CFL condition (2.7).

#### (iii) The linear scaling limiter

Theorem 2.1 tells us that for the scheme (2.4), we need to modify *p*_{j}(*x*) such that *p*_{j}(*x*)∈[*m*,*M*] for all *x*∈*S*_{j}, where *S*_{j} is defined in equation (2.5). For all *j*, assume , we use the modified polynomial by the limiter introduced in Liu & Osher (1996), i.e.
2.9with
2.10Let and . We get the revised scheme of equation (2.4):
2.11The scheme (2.11) satisfies the sufficient condition in theorem 2.1. We will show in the next lemma that this limiter does not destroy the uniform high order of accuracy.

### Lemma 2.3

*Assume* *then equations* (2.9) *and* (2.10) *give a* (*k*+1)-*th order accurate limiter*.

### Proof.

We need to show for any *x*∈*I*_{j}. We only prove the case that *p*_{j}(*x*) is not a constant and the other cases being similar. Since and we have . Therefore,
By the definition of *θ* in equation (2.9), implies that , i.e. there is an overshoot *M*_{j}>*M*, and the overshoot *M*_{j}−*M*=*O*(Δ*x*^{k+1}) since *p*_{j}(*x*) is an approximation with error *O*(Δ*x*^{k+1}). Thus, we only need to prove that , where *C*_{k} is a constant depending only on the polynomial degree *k*. In Liu & Osher (1996), *C*_{2}=3 is proved. We now prove the existence of *C*_{k} for any *k*. Assume *p*_{j}(*x*)=*a*_{0}+*a*_{1}[(*x*−*x*_{j})/Δ*x*]+⋯+*a*_{k}[(*x*−*x*_{j})/Δ*x*]^{k} and *p*(*x*)=*a*_{0}+*a*_{1}*x*+⋯+*a*_{k}*x*^{k}, then the cell average of *p*(*x*) on is and So we have
Let , then it suffices to prove the existence of *C*_{k} such that
It is easy to check that and are both norms on the finite dimensional linear space consisting of all polynomials of degree *k* whose averages on the interval *I* are zero. Any two norms on this finite dimensional space are equivalent, hence their ratio is bounded by a constant *C*_{k}. ■

Notice that in equation (2.10), we need to evaluate the maximum/minimum of a polynomial. We prefer to avoid evaluating the extrema of a polynomial, especially since we will extend the method to two dimensions. Since we only need to control the values at several points, we could replace equation (2.10) by 2.12and the limiters (2.9) and (2.12) are sufficient to enforce . As to the accuracy, limiter (2.12) is a less-restrictive limiter than equation (2.10), so the accuracy will not be destroyed. Also, it is a conservative limiter because it does not change the cell average of the polynomial.

For the conservative maximum-principle-satisfying scheme (2.11), it is straightforward to prove the following stability result.

### Theorem 2.4

*Assuming periodic or zero boundary conditions, then the numerical solution of equation (2.11) satisfies
*

### Proof.

Taking the sum of equation (2.11) over *j*, we obtain . Since the numerical solutions are maximum-principle-satisfying, namely, , we have
The other equality follows similarly. ■

### Remark 2.5

As an easy corollary, if the solution is non-negative, namely if *m*≥0, then we have the *L*^{1} stability .

#### (iv) High-order temporal discretization

We use SSP high-order time discretizations. For more details, see Shu & Osher (1988) and Shu (1988). For example, the third-order SSP Runge–Kutta method in Shu & Osher (1988) (with the CFL coefficient *c*=1) is
where *F*(*u*) is the spatial operator, and the third-order SSP multi-step method in Shu (1988) (with the CFL coefficient *c*=1/3) is
Here, the CFL coefficient *c* for an SSP time discretization refers to the fact that, if we assume the Euler forward time discretization for solving the equation *u*_{t}=*F*(*u*) is stable in a norm or a semi-norm under a time-step restriction Δ*t*≤Δ*t*_{0}, then the high-order SSP time discretization is also stable in the same norm or semi-norm under the time-step restriction Δ*t*≤*cΔt*_{0}.

Since an SSP high-order time discretization is a convex combination of Euler forward, the full scheme with a high-order SSP time discretization will still satisfy the maximum principle. The limiters (2.9) and (2.12) should be used for each stage in a Runge–Kutta method or each step in a multi-step method. For details of the implementation, see Zhang & Shu (2010*b*).

### (b) Two-dimensional extensions

Consider the two-dimensional scalar conservation laws *u*_{t} + *f*(*u*)_{x} + *g*(*u*)_{y} = 0, *u*(*x*,*y*,0)=*u*_{0}(*x*,*y*) with . We only discuss the DG method with the Euler forward time discretization in this section, but all the results also hold for the finite volume scheme (e.g. ENO and WENO).

#### (i) Rectangular meshes

For simplicity, we assume we have a uniform rectangular mesh. At time level *n*, we have an approximation polynomial *p*_{ij}(*x*,*y*) of degree *k* with the cell average on the (*i*,*j*) cell [*x*_{i−1/2},*x*_{i+1/2}]×[*y*_{j−1/2},*y*_{j+1/2}]. Let *u*^{+}_{i−1/2,j}(*y*),*u*^{−}_{i+1/2,j}(*y*), *u*^{+}_{i,j−1/2}(*x*), *u*^{−}_{i,j+1/2}(*x*) denote the traces of *p*_{ij}(*x*,*y*) on the four edges, respectively. A finite volume scheme or the scheme satisfied by the cell averages of a DG method on a rectangular mesh can be written as
where , are one-dimensional monotone fluxes. The integrals can be approximated by quadratures with sufficient accuracy. Let us assume that we use a Gauss quadrature with *L* points, which is exact for single variable polynomials of degree *k*. We assume denote the Gauss quadrature points on [*x*_{i−1/2},*x*_{i+1/2}], and denote the Gauss quadrature points on [*y*_{j−1/2},*y*_{j+1/2}]. For instance, (*β*=1,…,*L*) are the Gauss quadrature points on the left edge of the (*i*,*j*) cell. The subscript *β* will denote the values at the Gauss quadrature points, for instance, . Also, *w*_{β} denotes the corresponding quadrature weight on interval , so that . We will still need to use the *N*-point Gauss–Lobatto quadrature rule where *N* is the smallest integer satisfying 2*N*−3≥*k*, and we distinguish the two quadrature rules by adding hats to the Gauss–Lobatto points, i.e. will denote the Gauss–Lobatto quadrature points on [*x*_{i−1/2},*x*_{i+1/2}], and will denote the Gauss–Lobatto quadrature points on [*y*_{j−1/2},*y*_{j+1/2}]. Subscripts or superscripts *β* will be used only for Gauss quadrature points and *α* only for Gauss–Lobatto points.

Let *λ*_{1}=Δ*t*/Δ*x* and *λ*_{2}=Δ*t*/Δ*y*, then the scheme becomes
2.13

We want to find a sufficient condition for the scheme (2.13) to satisfy . We use ⊗ to denote the tensor product, for instance, . Define the set *S*_{ij} as
2.14See figure 1*a* for an illustration for *k*=2. For simplicity, let *μ*_{1}=*λ*_{1}*a*_{1}/(*λ*_{1}*a*_{1}+*λ*_{2}*a*_{2}) and *μ*_{2}=*λ*_{2}*a*_{2}/(*λ*_{1}*a*_{1}+*λ*_{2}*a*_{2}) where and . Notice that , we have
2.15

### Theorem 2.6

*Consider a two-dimensional finite volume scheme or the scheme satisfied by the cell averages of the DG method on rectangular meshes (2.13), associated with the approximation polynomials p*_{ij}*(x,y) of degree k (either reconstruction or DG polynomials). If* *and p*_{ij}*(x,y)∈[m,M] (for any (x,y)∈S*_{ij}*), then* *under the CFL condition
*2.16

### Proof.

Plugging equation (2.15) in, equation (2.13) can be written as Following the same arguments as in theorem 2.1, it is easy to check that the formulation above for is a monotonically increasing function with respect to all the arguments , , and . ■

To enforce the condition in theorem 2.6, we can use the following scaling limiter similar to the one-dimensional case. For all *i* and *j*, assuming the cell averages , we use the modified polynomial instead of *p*_{ij}(*x*,*y*), i.e.
2.17with
2.18It is also straightforward to prove the high-order accuracy of this limiter following the proof of lemma 2.3.

#### (ii) Triangular meshes

For each triangle *K*, we denote by (*i*=1,2,3), the length of its three edges (*i*=1,2,3), with outward unit normal vector *ν*^{i}(*i*=1,2,3). *K*(*i*) denotes the neighbouring triangle along and |*K*| is the area of the triangle *K*. Let be a one-dimensional monotone flux in the *ν* direction (e.g. Lax–Friedrichs flux), namely is an increasing function of the first argument and a decreasing function of the second argument. It satisfies (conservativity), and (consistency), with **F**(*u*)=〈*f*(*u*),*g*(*u*)〉. The first-order monotone scheme can be written as
Then *H*(⋅,⋅,⋅,⋅) is a monotonically increasing function with respect to each argument under the CFL condition where .

A high-order finite volume scheme or a scheme satisfied by the cell averages of a DG method, with first-order Euler forward time discretization, can be written as
where is the cell average over *K* of the numerical solution, and *u*^{int(K)}_{i} and *u*^{ext(K)}_{i} are the approximations to the values on the edge obtained from the interior and the exterior of *K*. Assume the DG polynomial on the triangle *K* is *p*_{K}(*x*,*y*) of degree *k*, then in the DG method, the edge integral should be approximated by the (*k*+1)-point Gauss quadrature. The scheme becomes
2.19where *w*_{β} denote the (*k*+1)-point Gauss quadrature weights on the interval , so that , and *u*^{int(K)}_{i,β} and *u*^{ext(K)}_{i,β} denote the values of *u* evaluated at the *β*-th Gauss quadrature point on the *i*-th edge from the interior and exterior of the element *K*, respectively.

Motivated by the derivation in the previous subsection, to find a sufficient condition for the scheme (2.19) to satisfy , we need to decompose the cell average by a quadrature rule, which includes all the Gauss quadrature points for each edge with all the quadrature weights being positive. Such a quadrature can be constructed by mapping the Gauss tensor Gauss–Lobatto points on a rectangle to a triangle. Details of the mapping can be found in Zhang *et al.* (in press). In the barycentric coordinates, the set of quadrature points for polynomials of degree *k* on a triangle *K* can be written as
2.20where *u*^{α} (*α*=1,…,*N*) and *v*^{β} (*β*=1,…,*k*+1) are the Gauss–Lobatto and Gauss quadrature points on the interval , respectively. See figure 1*b* for an illustration of .

### Theorem 2.7

*For the scheme (*2.19*) with the polynomial p*_{K}*(x,y) (either reconstruction or DG polynomial) of degree k to satisfy the maximum principle* *a sufficient condition is that each p*_{K}*(x,y) satisfies p*_{K}*(x,y)∈[m,M],* *, where* *is defined in equation (2.20), under the CFL condition* *Here,* *is still the quadrature weight of the N-point Gauss–Lobatto rule on* *for the first quadrature point.*

The proof is similar to that for the structured mesh cases, see Zhang *et al.* (in press) for the details. We can still use the same scaling limiter to enforce this sufficient condition.

## 3. Positivity-preserving high-order schemes for compressible Euler equations in gas dynamics

### (a) Ideal gas

The one-dimensional Euler system for the perfect gas is given by
3.1where **w**=(*ρ*,*m*,*E*)^{T}, **f**(**w**)=(*m*,*ρu*^{2}+*p*,(*E*+*p*)*u*)^{T}, *m*=*ρu*, *E*=(1/2)*ρu*^{2}+*ρe*, *p*=(*γ*−1)*ρe*, *ρ* is the density, *u* is the velocity, *m* is the momentum, *E* is the total energy, *p* is the pressure, *e* is the internal energy and *γ*>1 is a constant (*γ*=1.4 for the air). The speed of sound is given by and the three eigenvalues of the Jacobian **f**^{′}(**w**) are *u*−*c*, *u* and *u*+*c*.

Let *p*(**w**)=(*γ*−1)(*E*−(1/2)(*m*^{2}/*ρ*)) be the pressure function. It can be easily verified that *p* is a concave function of **w**=(*ρ*,*m*,*E*)^{T} if *ρ*≥0. Define the set of admissible states by *G*={**w**|*ρ*>0 and *p*=(*γ*−1)(*E*−(1/2)(*m*^{2}/*ρ*))>0}, then G is a convex set. If the density or pressure becomes negative, the system (3.1) will be non-hyperbolic and thus the initial value problem will be ill-posed. In this section, we discuss only the perfect gas case, leaving the discussion for general gases to §3*b*.

We are interested in schemes for equation (3.1) producing the numerical solutions in the admissible set *G*. We start with a first-order scheme
3.2where is a numerical flux. The scheme (3.2) and its numerical flux are called positivity preserving, if the numerical solution being in the set *G* for all *j* implies the solution being also in the set *G*. This is usually achieved under a standard CFL condition
3.3where *α*_{0} is a constant related to the specific scheme. Examples of positivity-preserving fluxes include the Godunov flux, the Lax–Friedrichs flux, the Boltzmann type flux, and the Harten-Lax-van Leer flux, see Perthame & Shu (1996). In Zhang & Shu (2010*c*), we proved that the Lax–Friedrichs flux is positivity preserving with *α*_{0}=1.

In Perthame & Shu (1996), a high-order scheme preserving the positivity was proposed, but it is quite difficult to implement the method, especially in multi-dimensions. In Zhang & Shu (2010*c*, 2011) and Zhang *et al.* (in press), we generalized the ideas in the previous section to construct high-order schemes preserving the positivity of density and pressure for the Euler system.

#### (i) One-dimensional compressible Euler equations

First, we consider the first-order Euler forward time discretization. A general high-order finite volume scheme, or the scheme satisfied by the cell averages of a DG method solving equation (3.1), has the following form:
3.4where is a positivity-preserving flux under the CFL condition (3.3), is the approximation to the cell average of the exact solution **v**(*x*,*t*) in the cell *I*_{j}=[*x*_{j−1/2},*x*_{j+1/2}] at time level *n*, and **w**^{−}_{j+1/2}, **w**^{+}_{j+1/2} are the high-order approximations of the point values **v**(*x*_{j+1/2},*t*^{n}) within the cells *I*_{j} and *I*_{j+1}, respectively. These values are either reconstructed from the cell averages in a finite volume method or read directly from the evolved polynomials in a DG method. We assume that there is a polynomial vector **q**_{j}(*x*)=(*ρ*_{j}(*x*),*m*_{j}(*x*),*E*_{j}(*x*))^{T} (either reconstructed in a finite volume method or evolved in a DG method) with degree *k*, where *k*≥1, defined on *I*_{j} such that is the cell average of **q**_{j}(*x*) on *I*_{j}, **w**^{+}_{j−1/2}=**q**_{j}(*x*_{j−1/2}) and **w**^{−}_{j+1/2}=**q**_{j}(*x*_{j+1/2}). Next, we state a similar result as in the previous section.

### Theorem 3.1

*For a finite volume scheme or the scheme satisfied by the cell averages of a DG method (3.4), if* *for all j and α, then* *under the CFL condition
*

The proof is similar to that for theorem 2.1 and can be found in Zhang & Shu (2010*c*).

SSP high-order Runge–Kutta in Shu & Osher (1988) and multi-step in Shu (1988) time discretization will keep the validity of theorem 3.1 since *G* is convex. If the numerical solutions have positive density and pressure, it follows that the scheme is *L*^{1} stable for the density *ρ* and the total energy *E* owing to theorem 2.4.

#### (ii) A limiter to enforce the sufficient condition

Given the vector of approximation polynomials **q**_{j}(*x*)=(*ρ*_{j}(*x*),*m*_{j}(*x*),*E*_{j}(*x*))^{T}, with its cell average , we would like to modify **q**_{j}(*x*) into , such that it satisfies the sufficient condition in theorem 3.1 without destroying the cell averages and high-order accuracy.

Define . Then and for all *j*. Assume there exists a small number *ε*>0, such that and for all *j*. For example, we can take *ε*=10^{−13} in the computation.

The first step is to limit the density. Replace *ρ*_{j}(*x*) by
3.5Then the cell average of over *I*_{j} is still and for all *α*.

The second step is to enforce the positivity of the pressure. We need to introduce some notations. Let and denote . Define *G*^{ε}={**w**:*ρ*≥*ε*,*p*=(*γ*−1)(*E*−(1/2)(*m*^{2}/*ρ*))≥*ε*}, ∂*G*^{ε}={**w**:*ρ*≥*ε*,*p*=*ε*}, and
3.6∂*G*^{ε} is a surface and **s**^{α}(*t*) is the straight line passing through the two points and . If , then the straight line **s**^{α}(*t*) intersects with the surface ∂*G*^{ε} at one and only one point since *G*^{ε} is a convex set. If , let denote the parameter in equation (3.6) corresponding to the intersection point; otherwise let . We only need to solve a quadratic equation to find , see Zhang & Shu (2010*c*) for details. Now we define
3.7

It is easy to check that the cell average of over *I*_{j} is and for all *α*. See Zhang & Shu (2010*c*) for the proof of the accuracy.

#### (iii) Two-dimensional cases

In this section, we extend our result to finite volume or DG schemes of (*k*+1)-th order accuracy solving two-dimensional Euler equations
3.8where *m*=*ρu*,*n*=*ρv*,*E*=(1/2)*ρu*^{2}+(1/2)*ρv*^{2}+*ρe*,*p*=(*γ*−1)*ρe*, and 〈*u*,*v*〉 is the velocity. The eigenvalues of the Jacobian **f**^{′}(**w**) are *u*−*c*, *u*, *u* and *u*+*c* and the eigenvalues of the Jacobian **g**^{′}(**w**) are *v*−*c*, *v*, *v* and *v*+*c*. The pressure function *p* is still concave with respect to **w** if *ρ*≥0 and the set of admissible states *G*={**w**|*ρ*>0,*p*>0} is still convex.

With the same notions as in §2, a finite volume scheme or the scheme satisfied by the cell averages of a DG method for equation (3.8) can be written as, on a rectangular mesh, 3.9or on a triangular mesh, 3.10

Assume at time level *n* there are approximation polynomials of degree *k*, **q**_{ij}(*x*,*y*) with the cell average on the (*i*,*j*) rectangular cell, or **q**_{K}(*x*,*y*) with the cell average on the triangle *K*, let , and , then we have the following.

### Theorem 3.2

*For a finite volume scheme or the scheme satisfied by the cell averages of a DG method (3.9) on a rectangle, if* *q*_{ij}*(x,y)∈G for all i,j and (x,y)∈S*_{ij} *defined in equation (2.14), then* *under the CFL condition*

### Theorem 3.3

*For a finite volume scheme or the scheme satisfied by the cell averages of a DG method (3.10) on a triangle, if* *q*_{K}*(x,y)∈G for all K and* *defined in equation (2.20), then* *under the CFL condition*

We can construct the same type of limiters as in the previous subsection to enforce the sufficient conditions in these two theorems. See Zhang & Shu (2010*c*) and Zhang *et al.* (in press) for the proof of the theorems and implementation of limiters.

### (b) General equations of state and source terms

Now we consider the one-dimensional Euler system (3.1) with a general equation of state *E*=*ρe*(*ρ*,*p*)+(1/2)*ρu*^{2} where *e*(*ρ*,*p*) is the internal energy. As we have seen in the previous subsection, to construct high-order schemes preserving the positivity of density and pressure, there are four important steps:

Prove

*G*={**w**:*ρ*>0 and*p*>0} is a convex set.Prove the first-order scheme (3.2) preserves the positivity.

Find a sufficient condition for the Euler forward time discretization as in theorem 3.1. Then high-order SSP Runge–Kutta or multi-step will keep the positivity owing to the convexity of

*G*.Construct a limiter to enforce the sufficient condition as in equations (3.5) and (3.7).

Notice that steps 3 and 4 above both heavily depend on the convexity of *G*. Therefore, to easily generalize the previous results to general equations of state, we should not give up the convexity. In Zhang & Shu (2011), we proved steps 1 and 2 will hold for any equation of state satisfying *e*≥0⇔*p*≥0 if *ρ*≥0. Once steps 1 and 2 are valid, it is very straightforward to complete steps 3 and 4 by following the ideas in theorem 3.1, equations (3.5) and (3.7). Two-dimensional extensions are also trivial by following theorems 3.2 and 3.3.

For Euler equations with source terms, for instance, the axial symmetry, gravity, chemical reaction or cooling effect, it is still possible to construct positivity-preserving high-order schemes. It is straightforward to extend all the previous results to Euler systems with various source terms, see Zhang & Shu (2011).

## 4. Applications

### (a) Maximum-principle-satisfying high-order schemes for passive convection equations with a divergence free velocity field

We will discuss how to take advantage of maximum-principle-satisfying high-order schemes for scalar conservation laws to construct such schemes for passive convection equations with a divergence free velocity field. We will explain the main idea for the two-dimensional incompressible Euler equation.

#### (i) Two-dimensional incompressible Euler equation

The two-dimensional incompressible Euler equations in the vorticity stream-function formulation are given by
4.1and
4.2with suitable initial and boundary conditions. The definition of 〈*u*,*v*〉 in equation (4.2) gives us the divergence-free condition *u*_{x}+*v*_{y}=0, which implies equation (4.1) is equivalent to the non-conservative form:
4.3The exact solution of equation (4.3) satisfies the maximum principle *ω*(*x*,*y*,*t*)∈[*m*,*M*], for all (*x*,*y*,*t*), where and . For discontinuous solutions or solutions containing sharp gradient regions, it is preferable to solve the conservative form (4.1) rather than the non-conservative form (4.3). However, without the incompressibility condition *u*_{x}+*v*_{y}=0, the conservative form (4.1) itself does not imply the maximum principle *ω*(*x*,*y*,*t*)∈[*m*,*M*] for all (*x*,*y*,*t*). This is the main difficulty to get a maximum-principle-satisfying scheme solving the conservative form (4.1) directly.

We recall the high-order DG method solving (4.1) in Liu & Shu (2000) briefly. For simplicity, we only discuss triangular meshes here. First, solve equation (4.2) by a standard Poisson solver for the stream-function *ψ* using continuous finite elements, then take *u*=−*ψ*_{y},*v*=*ψ*_{x}. Notice that on the boundary of each cell, 〈*u*,*v*〉⋅*ν*=〈−*ψ*_{y},*ψ*_{x}〉⋅*ν*=∂*ψ*/∂*τ*, which is the tangential derivative. Thus 〈*u*,*v*〉⋅*ν* is continuous across the cell boundary since the tangential derivative of *ψ* along each edge is continuous. The cell average scheme with Euler forward in time of the DG method in Liu & Shu (2000) is equivalent to
4.4Suppose is the DG polynomial on the triangle *K*. Then we can show that the right-hand side of equation (4.4) is a monotonically increasing function of the values of evaluated at in equation (2.20). See Zhang *et al.* (in press) for the proof. Therefore, to have , we only need to show the right-hand side of equation (4.4) is consistent. Namely, it is equal to *M* if , . This fact was proved in Zhang *et al.* (in press). We, therefore, have the following theorem.

### Theorem 4.1

*For a finite volume scheme or the scheme satisfied by the cell averages of a DG method (4.4) solving equation (4.1) on a triangle, if* *defined in equation (2.20), then* *under the CFL condition*

### Remark 4.2

The same result on rectangular meshes as in theorem 2.6 also holds, see Zhang & Shu (2010*b*).

### Remark 4.3

If one chooses another method to solve the velocity field, then the result still holds as long as the quadrature rules are exact for the velocity field in the scheme. This can be easily achieved if we pre-process the divergence-free velocity field so that it is piecewise polynomial of the right degree for accuracy, continuous in the normal component across cell boundaries, and pointwise divergence-free.

#### (ii) The level-set equation with a divergence-free velocity field

Let *ϕ*(*t*,*x*,*y*,*z*)=0 define the implicit interface, then the Eulerian formulation of the interface evolution can be written as
4.5If the velocity field satisfies *u*_{x}+*v*_{y}+*w*_{z}=0, then the solution of equation (4.5) satisfies the maximum principle, i.e. *ϕ*∈[*m*,*M*], where *m* and *M* are the minimum and maximum of *ϕ*_{0}. With the same idea, it is straightforward to construct maximum-principle-satisfying high-order finite volume or DG schemes solving equation (4.5).

#### (iii) Vlasov–Poisson equations

To describe the evolution of the electron distribution function *f*(*x*,*v*,*t*) of a collisionless quasi-neutral plasma in one space and one velocity dimension, where the ions have been assumed to be stationary, the Vlasov–Poisson system is given by
4.6

The exact solution of equation (4.6) also satisfies the maximum principle, which implies that the exact solution should always be non-negative. The positivity of the numerical solution for solving equation (4.6) is very difficult to achieve without destroying the conservation and high-order accuracy, as indicated in Banks & Hittinger (2010). Since *v*_{x}=*E*_{v}=0, the equation (4.6) is the same type as equation (4.1). Thus, theorem 4.1 also applies to equation (4.6). See figure 2 for the result of the positivity-preserving fifth-order finite volume WENO schemes for the two stream instability problem. The implementation detail of positivity-preserving limiter can be found in §5. As we can see, the traditional WENO schemes will produce negative values, which was also reported in Banks & Hittinger (2010). The positivity-preserving high-order scheme guarantees non-negativity and the result is comparable to those in Banks & Hittinger (2010) and Rossmanith & Seal (submitted).

Even though theorem 4.1 is only for the Eulerian schemes solving equation (4.6), the positivity-preserving techniques can also be extended to semi-Lagrangian schemes, see Rossmanith & Seal (submitted) and Qiu & Shu (submitted).

### (b) Shallow water equations

The shallow water equation with a non-flat bottom topography has been widely used to model flows in rivers and coastal areas. The water height is supposed to be non-negative during the time evolution. If it ever becomes negative, the computation will break down quite often since the initial value problem for the linearized system will be ill-posed. The positivity-preserving techniques can be also applied to one- or two-dimensional shallow water equations. In Xing *et al.* (2010), we constructed high-order DG schemes that preserve the well-balanced property and the non-negativity of the water height.

### (c) Vlasov–Boltzmann transport equations

The Vlasov–Boltzmann transport equations describe the evolution of a probability distribution function *f*(*x*,*v*,*t*) representing the probability of finding a particle at time *t* with position at *x* and phase velocity *v*. It models a dilute or rarefied gaseous state corresponding to a probabilistic description when the transport is given by a classical Hamiltonian with accelerations component given by the action of a Lorentzian force and particle interactions taken into account as a collision operator. Following the ideas described in previous sections, a high-order positivity-preserving DG method was proposed in Cheng *et al.* (in press).

### (d) Positivity-preserving schemes for a population model

When the numerical solutions denote the density or numbers, it is desired to have non-negative solutions. In Zhang *et al.* (submitted), a positivity-preserving high-order WENO scheme was constructed for a hierarchical size-structured population model, which involves global terms through integrals in the equation and boundary conditions.

## 5. A simplified implementation of the maximum-principle-satisfying and positivity-preserving limiter for weighted essentially non-oscillatory finite volume schemes

### (a) Motivation

As described in previous sections, the maximum-principle-satisfying and positivity-preserving high-order finite volume or DG schemes are easy to implement if the approximation polynomials are available. In the DG method, these are simply the DG polynomials. In the finite volume ENO schemes, the polynomials are constructed during the reconstruction procedure. However, the WENO reconstruction returns only some point values rather than approximation polynomials. Therefore, to implement the maximum-principle-satisfying and positivity-preserving high-order WENO schemes according to the procedure described in the previous sections, one must first obtain the approximation polynomials beyond the WENO reconstructed point values, for example, by constructing interpolation polynomials as we did in Zhang & Shu (2010*b*). Thus implementation of the limiter for WENO schemes is more expensive and cumbersome, especially for multi-dimensional problems. In this section, we will propose an alternative and simpler implementation to achieve the same maximum principle or positivity without using the approximation polynomials explicitly, which results in a reduction of computational cost and complexity of the procedure for WENO schemes and even for the DG method.

Let us revisit maximum-principle-satisfying schemes for the one-dimensional scalar conservation laws in §2. To have , for all *α* is sufficient but not necessary. By the mean value theorem, there exists some *x*^{*}_{j}∈*I*_{j} such that . Then equation (2.8) can be rewritten as
5.1Therefore, we can have a much weaker sufficient condition.

### Theorem 5.1

*For the scheme (2.4), if* *then* *under the CFL condition*

To enforce this new sufficient condition, we can use the same limiter (2.9) with *M*_{j} and *m*_{j} redefined as
5.2

Equations (2.9) and (5.2) will not destroy the high-order accuracy since it is a less-restrictive limiter than limiters (2.9) and (2.10).

Equation (2.6) implies that Therefore, *θ* defined in equations (2.9) and (5.2) can be calculated without the explicit expression of the approximation polynomial *p*_{j}(*x*) or the location *x*^{*}_{j}. We only need to know the existence of such polynomials to prove the accuracy of the limiter. For WENO schemes, the existence of such approximation polynomials can be established by the interpolation, for example, Hermite interpolation for the one-dimensional case as in Zhang & Shu (2010*b*).

Extensions to two-dimensional cases are straightforward.

### Theorem 5.2

*Consider the scheme (2.13). There exists some point (x*^{*}_{i}*,y*^{*}_{j}*) in the (i,j) cell such that
*5.3*If* *, then* *under the CFL condition*

### Theorem 5.3

*Consider the scheme (2.19). There exists some point (x*^{*}_{K}*,y*^{*}_{K}*) in the triangle K such that
**If* *, then* *under the CFL condition*

### Remark 5.4

All the results above can also be easily extended to positivity-preserving schemes for compressible Euler equations.

### (b) Easy implementation for weighted essentially non-oscillatory finite volume schemes

We only state the algorithm for two-dimensional scalar conservation laws on rectangular meshes, the counterparts for the triangular meshes and compressible Euler equations are similar. For each stage in the SSP Runge–Kutta or each step in the SSP multi-step methods of the finite volume WENO schemes (2.13), the algorithm flowchart of the new limiter is

— For each rectangle, given and constructed by the WENO reconstruction, compute with equation (5.3), where

*M*_{ij}and*m*_{ij}are the max and min of .— Set and .

— Replace by the revised nodal values , , , in the scheme (2.13).

### Remark 5.5

The new algorithm is simpler and less expensive than the implementation in Zhang & Shu (2010*b*), since no extra reconstructions need to be performed for the limiter.

### Remark 5.6

The new algorithm is also cheaper for the DG method because it avoids the evaluation of the point values in *S*_{j}, *S*_{ij} and . The algorithm flowchart for the DG method is almost identical to the one described above for the finite volume method, and is therefore omitted to save space.

### (c) Numerical tests for the fifth-order weighted essentially non-oscillatory schemes

We show some numerical tests for the fifth-order finite volume WENO schemes with the simplified implementation of the limiter on rectangular meshes described above. The time discretization is the third-order SSP Runge–Kutta and the CFL is taken as equation (2.16). The algorithm for finite volume WENO schemes on rectangular meshes was described in Shu (2009) and the linear weights can be found in the appendix of Zhang & Shu (2010*b*), where the negative linear weights should be dealt with by the method in Shi *et al.* (2002). Extensive tests for scalar conservation laws were done to test the accuracy for the new limiter mentioned above. The results are similar to those in Zhang & Shu (2010*b*). We will not show the accuracy tests here to save space.

### Example 5.7

(Two stream instability for Vlasov–Poisson equations). The initial and boundary conditions are the same as in Banks & Hittinger (2010). See figure 2 for the results. The numerical solution on the top-right in figure 2 is non-negative everywhere. This can be clearly seen in the cuts along *x*=0, especially in the zoomed cuts on the bottom-right in figure 2.

### Example 5.8

(Low-density or low-pressure problems for compressible Euler equations). We consider the two-dimensional Sedov blast wave and 90^{°} shock diffraction problem in Zhang & Shu (2010*c*), where the results of the positivity-preserving third-order DG method were reported. Traditional finite volume and finite difference WENO schemes will blow up for such problems. Here, we show the results of the fifth-order finite volume WENO scheme with the new positivity-preserving limiter. See figures 3 and 4. The results are comparable to those of the DG method.

## 6. Concluding remarks

We have given a review of the recently developed maximum-principle-satisfying high-order finite volume or DG schemes for scalar conservation laws, including generalizations and applications to two-dimensional incompressible Euler equations and passively convection equations with a divergence free velocity field, and positivity-preserving schemes for compressible Euler equations, shallow water equations, Vlasov–Boltzmann transport equations and a population model. We also propose a simpler and less-expensive implementation especially for the finite volume WENO schemes, and provide several numerical examples to demonstrate their performance.

## Acknowledgements

Support by AFOSR grant FA9550-09-1-0126 and NSF grant DMS-0809086 is acknowledged.

- Received March 8, 2011.
- Accepted April 5, 2011.

- This journal is © 2011 The Royal Society