Numer. Math. (2004) 97: 441–471 Digital Object Identifier (DOI) 10.1007/s00211-003-0502-9

Numerische Mathematik

Numerical approximation of entropy solutions for hyperbolic integro-differential equations Andreas Dedner, Christian Rohde Albert–Ludwigs–Universit¨at Freiburg, Institut f¨ur Angewandte Mathematik, Hermann–Herder Str. 10, 79104 Freiburg i. Brsg, Germany; e-mail: {dedner, chris}@mathematik.uni-freiburg.de Received May 2, 2002 / Revised version received July 12, 2003 / c Springer-Verlag 2004 Published online February 25, 2004 –


Summary. Scalar hyperbolic integro-differential equations arise as models for e.g. radiating or self-gravitating fluid flow. We present finite volume schemes on unstructured grids applied to the Cauchy problem for such equations. For a rather general class of integral operators we show convergence of the approximate solutions to a possibly discontinuous entropy solution of the problem. For a specific model problem in radiative hydrodynamics we introduce a convergent fully discrete finite volume scheme. Under the assumption of sufficiently fast spatial decay of the entropy solution we can even establish the convergence rate h1/4 | ln(h)| where h denotes the grid parameter. The convergence proofs rely on appropriate variants of the classical Kruzhkov method for local balance laws together with a truncation technique to cope with the nonlocal character of the integral operator. Mathematics Subject Classification (2000): 35L65, 35Q35, 65M15

1 Introduction Let T > 0 be given. We consider the following inhomogeneous Cauchy problem for the unknown function u = u(x, t) : R2 × [0, T ) → R. ˆ (1) ∂t u(x, t) + divf (u(x, t)) = T[u(., t)](x), (2) u(x, 0) = u0 (x), Correspondence to: C. Rhode

x ∈ R2 , t ∈ (0, T ), x ∈ R2 .

442

A. Dedner, C. Rohde

Here f : R → R2 denotes the flux function and u0 is a scalar function on R2 . For functions B : R → R, q : R → R, and a kernel function ˆ takes the form k : R2 × R2 → R, the balance term T (3)

ˆ T[w] =
T[w] + q(w), T[w](x) =

w : R2 → R,

k(x, y)B(w(y)) dy, R2

x ∈ R2 .

Precise assumptions on the data will be presented later on. Here we proceed with saying that the problem (1),(2) is general enough to represent a scalar model problem for the dynamics of compressible fluids described by the Euler equations of gas dynamics undergoing solution-dependent external forces. For B = 0 and k = 0, the term T reflects the global influence of heat sources or gravitation fields. In particular, integral operators arise in radiative hydrodynamics ([15, 20]) and in the context of self-gravitating fluids modeled by the Euler-Poisson system ([21]). Under certain physical assumptions the balance law (1) can be derived as a simple model problem in these cases. For instance, it plays the same rˆole for the equations of radiative hydrodynamics or the Euler-Poisson system as the scalar hyperbolic conservation law for pure hydrodynamics. In passing we note that equations of type (1), respectively a slight variant of it, have been derived also in traffic simulation ([18]). Let us return to the scalar problem (1),(2). In Section 2.1 we will specify the assumptions on the data of the Cauchy problem and present a general ˆ under consideration called admissible operators. Due to class of operators T the global character of the integral operator solutions of (1),(2) – when they exist– may not be compactly supported for positive time even if the initial datum is compactly supported. With other words: we loose the finite speed of propagation, a basic property in the theory of hyperbolic conservation laws. Infinitely fast spreading of information is a property of parabolically regularized versions of (1),(2), i.e., ∂t u(x, t) + divf (u(x, t)) = u(x, t). These allow for globally smooth solutions and one might then conjecture that (1),(2) could also have smooth solutions. However, this is not true in general. In the case of special but physically relevant examples it has been shown recently that classical solutions do not exist globally in time for smooth initial datum u0 with big amplitude (or big amplitude derivative). For initial datum below some critical threshold the solutions stay classical for all times. We refer to [3,6,11]. Partial regularization is typical for spatial integral operators as considered here. As a consequence we cannot use a notion of classical solutions. In the sequel we shall search for weak solutions, namely for entropy solutions in the sense of Kruzhkov which turn out to provide the correct notion

Numerical approximation of entropy solutions

443

of solution ([3,14]). Section 2.1 concludes with the definition of entropy solutions and a well-posedness assumption on the existence, uniqueness, and regularity of entropy solutions. The main part of this paper is devoted to the numerical analysis of a formally first order finite volume scheme for (1),(2). Finite volume type schemes on unstructured meshes are used successfully in many examples of systems of hyperbolic balance laws (cf. [8, 12], for instance), in particular also in radiative hydrodynamics. In Section 2.2 we introduce the unstructured grid, the fully discrete monotone finite volume scheme, and the approximate soluˆ ,h of tion defined by the scheme. In particular, possible discretizations T ˆT are considered. Here h denotes the grid parameter and  the computational domain. We present discrete admissible operators as counterparts of the admissible operators. In Section 2.3 we derive a general a-priori error estimate (Theorem 1) for the L1 -error between the approximate and the exact entropy solution on R2 in terms of h and the size of . For a regular family of grids the estimate will hold provided ˆ ˆ ,h is an admissible discrete operator for T, (a) the operator T ∞ (b) the approximate solution is bounded uniformly in L () with respect to h. Coupling the size of  with h in an appropriate way leads to a convergent algorithm. We obtain a positive rate of convergence for h → 0 if (among other more technical conditions) (c) the exact entropy solution decays sufficiently fast at ±∞. The proof relies on the seminal work of Kruzhkov and Kuznetsov which has been used to prove convergence(rates) in the homogeneous case ([2,7,13, 16,19]) and in the inhomogeneous case with local source terms ([1,10,17]). The first main additional difficulty and novelty in our case is the global character of the integral operator T. As noted above we loose the basic property of the exact solutions of the equation (1) in the homogeneous case: the finite speed of propagation. This problem is overcome by some cut-off technique and careful tracking of the size of the computational domain. The second main difficulty is related to the numerical approximation of T itself. To control the associated data error some compactness of T is necessary. This need fits exactly with the fact that integral operators in relevant physical situations are weakly singular in many cases. Finally in Section 3 we will apply the general convergence Theorem 1 to a model problem motivated by the system of radiative hydrodynamics. The evolution of a radiative plasma can be described by the system of gas dynamics where the energy transport through radiation is modeled by a source term Qrad in the balance law for the total energy density. The source term describes

444

A. Dedner, C. Rohde

the emission and absorption of the radiation intensity I = I (x, t, µ) at the point (x, t) ∈ R2 × [0, T ] in direction µ ∈ S , S being the unit sphere. It is given by
Qrad (x, t) = χ (x, t) (I (x, t, µ) − B(x, t)) dµ. S

The functions χ and B are given functions depending in fact on the temperature and the density of the plasma. In the case where the plasma velocities are very small compared to the speed of light one can assume an instantaneous radiation equilibrium. In this case the radiation intensity in a fixed direction µ from S satisfies µ · ∇x I (x, t, µ) + χ(x, t)I (x, t, µ) = χ (x, t)B(x, t). The assumption of instantaneous radiation equilibrium leads to the global character of the source term Qrad . In the model problem we assume that χ is constant and by formally solving the transport equation for I along characteristics we arrive at an expression for Qrad of the form (3). For the derivation of this model and analytical results we refer to [3,4]. In the more restrictive case of unidirectional radiation the model problem (1),(2) has been considered in [9,11]. We will introduce a fully discrete FV-scheme, in particular, an admissible ˆ ,h such that the conditions (a) and (b) hold true. From the discrete operator T main Theorem 1 we conclude the convergence of the scheme in Theorem 4. Exploiting decay estimates for solutions with initially compact support in the sense of (c) we obtain an error of order h1/4 | ln(h)|. This is a logarithmically small loss in order when compared to the up-to-date results for local scalar balance laws. We recall that numerical experiments give an error estimate of order h1/2 but the proof is an open problem up to the knowledge of the authors.

2 The general convergence theory In this section we present the necessary analytical background for problem (1),(2) and introduce a numerical scheme for the approximate solution of (1),(2). We will prove a general convergence theorem for this approximate solution depending on the mesh size h and the size of the computational domain . The theorem applies if a number of conditions on the solutions of problem (1),(2) and on the discretization are satisfied. In the second part of the paper –Section 3– we will show that these conditions are satisfied for the solutions of a physically relevant model problem in radiative fluid dynamics and its naturally discretized version.

Numerical approximation of entropy solutions

445

Concerning notations Lp (S), p ∈ [1, ∞], denotes the usual Lebesgue space on S ⊂ Rn and C k (S), k ∈ N0 , the space of k-times continuously differentiable functions on S. We need the space BV of functions of bounded variation:    BV (S) = w ∈ L1 (S)  |w|BV (S)   
:= sup w(x)div(φ(x)) dx < ∞ . φ∈C0∞ (S), φ∞ ≤1

S

Furthermore we make the subsequent convention for the corresponding norms of functions u ∈ Lp (R2 × [0, T ]) and w ∈ Lp (R2 ), p ∈ [1, ∞) ∪ {∞}: up = uLp (R2 ×[0,T ]) ,

wp = wLp (R2 ) .

The expressions .C 1 and |.|BV have to be understood in the same manner. Norms of all spaces of functions acting on other subsets of R2 × [0, T ] or R2 will be given explicitly. The symbol C denotes a generic positive constant and the continuous function CTˆ : R≥0 → R≥0 is used to give estimates on ˆ and its discretization. We define the space W (0, T ) by the operator T W (0, T ) = L∞ (0, T ; L1 (R2 ) ∩ L∞ (R2 )). 2.1 Problem statement and entropy solutions We give the precise assumptions on the data of the Cauchy problem (1),(2). A truncated problem and the concept of entropy solutions for (1),(2) will be introduced. Concerning the initial datum u0 we will assume u0 ∈ L∞ (R2 ),

(4)

supp (u0 ) ⊂ BR0 (0),

where R0 > 0. For the flux f we require f ∈ C 1 (R, R2 ). ˆ Let B, q ∈ C 1 (R, R) with q(0) = 0 Definition 1 (Admissible operator T) 2 2 ˆ from (3) is called admisand k : R × R → R be given. The operator T sible if there is a continuous function CTˆ : R≥0 → R≥0 such that for all w, w˜ ∈ L∞ (R2 ) ∩ L1 (R2 ) with w∞ , w ˜ ∞ ≤ M, M > 0 we have ˆ : L∞ (R2 ) ∩ L1 (R2 ) → L∞ (R2 ) ∩ L1 (R2 ), T ˆ T[w] ˆ (M), ∞ ≤ CT ˆ T[w]1 ≤ CTˆ (M)w1 ,

    ˆ  w(x) − w(x) ˆ w](x) (8) (iv) ˜  dx. ˜ T[w](x) − T[  dx ≤ CTˆ (M) (5) (i) (6) (ii) (7) (iii)

R2

The function CTˆ can depend on B, q, and k.

R2

446

A. Dedner, C. Rohde

Note 1 (Existence of admissible operators) Let us comment on the validity of the estimates (6), (7), and (8) as far as the integral operator T in (3) is concerned. Clearly operators with smooth bounded kernel function k with k(x, .) ∈ L1 (R2 ), x ∈ R2 , e.g. convolution-type operators, are included. More interestingly, the conditions can be verified for (at least all linear) weakly singular operators, i.e. operators T with B ≡ I d and kernels satisfying (9) 0 ≤ k(x, y) ≤ M|x − y|α−2 ,

x, y ∈ R2 ,

x = y,

α ∈ (0, 2].

In Section 3 we will consider a weakly singular operator from radiative hydrodynamics in detail. ˆ be an admissible operator. For some compact subset  of Now, let T R , consider additionally the following truncated Cauchy problem for u : R2 × [0, T ) → R. 2

(10) (11)

∂t u (x, t) ˆ  [u (., t)](x), (x, t) ∈ R2 × (0, T ), +divf (u (x, t)) = T x ∈ R2 . u(x, 0) = u0 (x),

ˆ  is defined by Here the truncated operator T (12)

ˆ  [w] = T [w] + q(w), T
T [w](x) := χ (x)

w ∈ L∞ (R2 ) ∩ L1 (R2 ),

k(x, y)B(w(y)) dy, R2

x ∈ R2 .

ˆ is admissible by assumption. ˆ  is admissible since T Note that the operator T We have CTˆ  = CTˆ . Due to the definition of T [w] and q(0) = 0 by Definition 1 any distributional solution of (10),(11) has compact support if u0 is compactly supported. Recall that a distributional solution of (1),(2) is not compactly supported in general. Let a function u ∈ W (0, T ), κ ∈ R, and φ ∈ C0∞ (R2 × [0, T )) be given. We introduce the form ETˆ = ETˆ (u, κ, φ) by
ETˆ (u, κ, φ) =

T


R2

0

|u(x, t) − κ|φt (x, t) + (f (u(x, t) κ)

−f (u(x, t)⊥κ)) · ∇φ(x, t) dxdt
T
ˆ sgn (u(x, t) − κ)T[u(., t)](x)φ(x, t) dxdt. + 0

R2

Here we put a b = max{a, b} and a⊥b = min{a, b} for a, b ∈ R. Following the discussion on regularity of solutions for (1),(2) we formulate

Numerical approximation of entropy solutions

447

Definition 2 (Entropy solution) Consider (1),(2) such that (4) holds. Let ˆ be given. an admissible operator T A function u ∈ W (0, T ) is called an entropy solution (or Kruzhkov solution) of (1),(2) if we have for all κ ∈ R and φ ∈ C0∞ (R2 × [0, T ), R≥0 ) the inequality
|u0 (x) − κ|φ(x, 0) dx. ETˆ (u, κ, φ) ≥ − R2

For the rest of the paper we make the following assumption on the existence, uniqueness and regularity of entropy solutions for (1),(2). ˆ be an Assumption 1 (Existence/uniqueness of entropy solutions) Let T admissible operator. Consider the Cauchy problems (1),(2) and (10),(11) such that (4) holds for u0 . (i) There exists a unique entropy solution u ∈ W (0, T ) of (1),(2) and a unique entropy solution u ∈ W (0, T ) of (10),(11). (ii) If u0 satisfies additionally u0 ∈ BV (R2 ) we have u, u ∈ BV (R2 × [0, T ]), and in particular for almost all t ≥ 0 u(., t), u (., t) ∈ BV (R2 ). Note 2 (Existence results) Assumption 1 holds true at least for small times T > 0. In fact the existence of an entropy solution u ∈ L∞ (R2 × [0, T ]) ∩ C(0, T ; L1 (R2 )∩L∞ (R2 )) can be shown. Existence for arbitrary times T > 0 can be proven for special choices of T and q, (for instance in the case of radiating fluids to be considered later on in Section 3). For both results we refer to [3]. Our convergence results will depend on the (spatial) modulus of contiˆ ˆ applied to the entropy nuity of T[u(., t)] for t ∈ [0, T ], i.e. the operator T solution u. For a function w ∈ L∞ (R2 ),  > 0, and a subset S of R2 , the (integrated) modulus of continuity ε(, S, w) is given by  
|w(x + x) − w(x)| dx . ε(, S, w) = sup |x|≤

S

ˆ and the entropy solution u from Assumption 1 For an admissible operator T we define the function γε : R≥0 → R≥0 by   ˆ (13) t)]) , γε () = ess sup ε(, R2 , T[u(., t∈[0,T ]

for  > 0 and γε (0) = 0. Since we have u(., t) ∈ L∞ (R2 ) Definition 1(i) ensures by classical results on L∞ -functions that γε satisfies γε () = o(1) ˆ ([14]). We note that T[u(., t)] ∈ BV (R2 ), t ∈ [0, T ], would lead to γε () = O().

448

A. Dedner, C. Rohde

2.2 The numerical scheme 2.2.1 The unstructured grid and the time step. Let h0 ∈ (0, 1). For h ∈ (0, h0 ], we consider a family of unstructured grids {Th } on R2 . Let a natural number K ≥ 3 be given. We assume that for each grid parameter h there exists an index set I h such that the grid Th satisfies h (i) T h = {Tj | j ∈ 2I , Tj closed k-polygon, k ∈ {3, . . . , K}}, (ii) j ∈I h Tj = R , (iii) the intersection of any two polygons of Th is either a joint edge, a vertex, or empty, (iv) h = supj ∈I h {diam(Tj )}.

We introduce the following notations for a grid Th . |Tj | N (j ) Ej l |Ej l | ωj Ih νj l

: : : : : : :

area of Tj , j ∈ I h , indices of elements that have a joint edge with Tj , joint edge of Tj , j ∈ I h and Tl , l ∈ N (j ), length of Ej l , center of gravity of Tj , index set {j | Tj ∩  = ∅} for some set  ⊂ R2 , normal of Ej l (pointing from Tj to Tj l ) with length |Ej l |, j ∈ I h , l ∈ N (j ).

The time step of the numerical method will be denoted by t > 0 and t n = nt for n ∈ N0 . For T > 0, the number NT ∈ N is the smallest natural number such that tNT > T . We suppose that the family of unstructured grids {Th } and t satisfy the subsequent uniformity condition. Assumption 2 (Grid and time step) Let {Th } be a family of unstructured grids on R2 and t > 0. There exist constants cG , c1 > 0 such that for all j ∈ I h and l ∈ N (j ) (14)

h ≤ c1 , t

cG h2 ≤ |Tj |,

cG |Ej l | ≤ h.

2.2.2 The Finite Volume Scheme. We start with the standard definition of monotone numerical flux functions. Definition 3 (Monotone numerical flux) Let a family {Th } of unstructured grids be given. For each h ∈ (0, h0 ], the family Fh = {Fj l ∈ C 2 (R2 ) | j ∈ I h , l ∈ N (j )} is called a family of monotone numerical fluxes for (f, Th ) if we have for all j ∈ I h , l ∈ N (j ) (i) Fj l is consistent with f , i.e. for all u ∈ R Fj l (u, u) = νj l · f (u),

Numerical approximation of entropy solutions

449

(ii) Fj l is conservative, i.e. for all u, v ∈ R Fj l (u, v) = −Flj (v, u), (iii) Fj l is Lipschitz continuous, i.e. for each M > 0 there exists a constant LF (M) > 0 such that for all u1 , u2 , v1 , v2 ∈ BM (0) we have |Fj l (u1 , v1 ) − Fj l (u2 , v2 )| ≤ |Ej l |LF (M)(|u1 − u2 | + |v1 − v2 |), (iv) Fj l is monotone non-increasing in the second variable, i.e. we have for all (u, v) ∈ R2 ∂ Fj l (u, v) ≤ 0. ∂v Examples of monotone numerical flux functions can be provided by versions of the Lax-Friedrichs or the Engquist-Osher flux [12]. For each h ∈ (0, h0 ], let Vh denote the set of functions on R2 which are piecewise constant on the cell volumes Tj , j ∈ I h . ˆ ,h ) Consider an admissible Definition 4 (Admissible discrete operator T ˆ ˆ and, for  ⊂ R2 compact, the associated truncated operator T operator T from (12). Let a family of unstructured grids {Th } be given. ˆ ,h : Vh → Vh is called an admissible discrete For h ∈ (0, h0 ], an operator T ˆ  on Th if there exist continuous functions CTˆ , γapp : R≥0 → operator for T R≥0 such that for all wh , w˜ h ∈ Vh ∩L∞ (R2 )∩L1 (R2 ) with wh ∞ , w˜ h ∞ ≤ M, M > 0, we have (15) (16) (17)

ˆ ,h [wh ] ≤ CTˆ (M), (i) T ∞ ˆ ,h [wh ] ≤ CTˆ (M)wh 1 , (ii) T 1 ˆ ,h [wh ]) ⊂ {x ∈ R2 | dist(x, ) ≤ h}, (iii) supp (T    ˆ  ˆ |Tj |T ˜ h ](ωj ) (iv) ,h [wh ](ωj ) − T,h [w j ∈I h

≤ CTˆ (M)

(18)



    |Tj |(wh − w˜ h )(ωj ),

j ∈I h




  T ˆ ,h [wh ](x) dx ˆ  [wh ](x) − T

(v) Tj

(19)

≤ |Tj |γapp (h) (j ∈ Ih ).

The function γapp can depend on M (which is omitted in the notation) but not on . It satisfies for  > 0 (20)

γapp () > 0,

γapp (0) = 0.

As the next step we define the finite volume scheme.

450

A. Dedner, C. Rohde

Definition 5 (Finite volume scheme) Consider problem (1),(2) in R2 × [0, T ) for T > 0. Let a family of unstructured grids {Th }, h ∈ (0, h0 ], and t > 0 be given. Assume that Fh is a family of numerical fluxes for (f, Th ) ˆ  on Th . ˆ ,h is an admissible discrete operator for T and that T h For n = 1, . . . , NT − 1 and j ∈ I , we define iteratively (21)

= unj − un+1 j

t  ˆ ,h [u,h (., t n )](ωj ), Fj l (unj , unl ) + t T |Tj | l∈N (j )

where u0j is given by (22)

u0j =

1 |Tj |


u0 (x) dx. Tj

The approximate solution u,h : R2 × [0, NT t] → R is defined by u,h (x, 0) = u0j for x ∈ R2 , j ∈ I h and by (23)

for (x, t) ∈ Tj × (t n , t n+1 ] u,h (x, t) = un+1 j

for n ∈ {0, . . . , NT − 1}, j ∈ I h . Note 3 (Support of uh ) For some set S ⊂ R2 and δ > 0, let B (S, δ) = {x ∈ R2 | dist(x, S) ≤ δ}. Consider an approximate solution u,h given by Definition 5. ˆ ,h [u,h (., 0)]) ⊂ B (, h) Since the scheme (21) is explicit and supp (T by (17) we get the inclusion supp u,h (., t 1 ) ⊂ BR0 +h (0) ∪ B (, 2h). This implies inductively the inclusion supp u,h ⊂ BR0 +NT h ∪ B (, (NT +1)h)× [0, T ]. From (14) we obtain supp u,h ⊂ B(R0 +(T +t)/c1 ) ∪ B (, (T + 2t)/c1 ) × [0, T ]. We define RT () > 0 to be the smallest number such that for all t ∈ [0, T ] (24)

supp (u,h (., t)) ∪ supp (u (., t)) ⊂ BRT () (0),

Note 4 (Explicit/implicit scheme) In this paper we consider only the fully explicit scheme (21). Stiffness of the (local or global) source terms may lead one to favour an implicit discretization with less restrictive time step conditions. Convergence results similar to the one presented here for the explicit case can then be derived along the same lines.

Numerical approximation of entropy solutions

451

2.3 Error analysis for the finite volume scheme Let Assumption 1 hold with entropy solution u ∈ W (0, T ) for (1),(2) and entropy solution u ∈ W (0, T ) for (10),(11). Furthermore, for h ∈ (0, h0 ], let the family of approximations {uh } from Definition 5 be given. In this subsection we present the principal steps towards the estimate of the error u − u,h 1 . The main result is the following theorem. ˆ be an admissible operator. Theorem 1 Let u0 ∈ BV (R2 ) satisfy (4) and T Consider the Cauchy problem (1),(2) with the entropy solution u ∈ W (0, T ). For h ∈ (0, h0 ] and a compact set  ⊂ R2 with || ≥ 1, let an unstructured grid Th satisfying Assumption 2, an associated family of monotone numerical ˆ ,h for T ˆ  on Th be given. fluxes Fh and an admissible discrete operator T We suppose that there is a constant MT > 0 (independent of h and ) such that we have for u and the function u,h given by Definition 5 u∞ , u,h ∞ ≤ MT .

(25)

If the time step t > 0 satisfies for some ξ ∈ (0, 1) (26)

2 (1 − ξ )cG t ≤ h 2LF (MT )

(j ∈ I h ),

then there exists a constant C > 0 such that the estimate
T
|u(x, t) − u,h (x, t)| dxdt 2 R 0 ≤ C (1 + RT ()2 )h1/4 + γε (h1/4 ) + RT ()2 γapp (h)


(27) |T[u(., t)](x)| dxdt + ess sup t∈[0,T ]

R2 \

holds. The constant C > 0 depends on the quantities CTˆ (MT ), cG , ξ, T , u0 but not on h or . For the definition of RT we refer to (24). Proof. The proof follows from Theorems 2 and 3 below which treat the errors u − u 1 and u − u,h 1 according to (28)

u − u,h 1 ≤ u − u 1 + u − u,h 1 .

 

Discussion of Theorem 1: Theorem 1 gives an a-priori error estimate in terms of h and the size of . To obtain convergence for a numerical scheme in terms of the discretization parameter h alone we couple the up to now arbitrary set  with h.

452

A. Dedner, C. Rohde

Corollary 1 Let the assumptions of Theorem 1 be satisfied for  = h , h ∈ (0, h0 ] where the set h is given by a ball of radius min{h−1/16 , (γapp (h))−1/4 } around 0. Then we have lim u − uh ,h 1 = 0.

h→0

Proof. From the definition of h and RT in (24) we conclude that there is a constant C > 0 that does not depend on h such that (29) RT (h )2 (h1/4 + γapp (h)) ≤ C(h1/8 + γapp (h)1/2 ) → 0

for h → 0.

Furthermore we have limh→0 |h | → ∞ and therefore
lim ess sup (30) |T[u(., t)](x)| = 0 h→0 t∈[0,T ]

R2 \h

for the L1 -function T[u]. Combining (27),(29) and (30) gives the convergence statement since γε = o(1).   Of course the choice for h in the corollary might be not optimal for the convergence speed. To obtain even an algebraic rate of convergence for the error ˆ we need sufficiently fast decay of T[u](x) for |x| → ∞ and, in particular, the existence of a constant α > 0 such that (31)

γε (h1/4 ) + γapp (h) = O(hα ).

ˆ Let us consider the most basic choice to discretize T
ˆ ,h [wh ](x) := 1 ˆ  [wh ](y) dy (32) (x ∈ Tj , wh ∈ Vh ). T T |Tj | Tj ˆ  : L1 (R2 ) ∩ L∞ (R2 ) → BV (R2 ) standard theory for Provided we have T BV -functions shows α = 1/4 ([14]). Integral operators with smooth kernel functions that decay sufficiently fast provide examples of such compact operators (e.g. convolution operators with smooth kernels). The question becomes more delicate when we consider weakly singular integral operators (cf. (9)) which in general map into a compact subset of L1 (R2 ) but not necessarily into BV (R2 ). Here one has to consider each case separately. In Section 3 we present an admissible discrete operator for an important physically relevant example that satisfies (31). The problem involves a weakly singular operator. In fact the admissible discrete operator avoids exact evaluation of integrals as in (32) but relies on (implementable) point evaluations.

Numerical approximation of entropy solutions

453

2.3.1 The error u − u,h 1 . We start to prove estimates for the error according to the splitting in (28). It is supposed that the assumptions of Theorem 1 hold. The proofs follow the guidelines as they have been developped in particular in the work of Eymard et al. ([7]) and in [17]. Discrete entropy inequalities and weak BV-estimates: Without loss of generality we assume that the bounded entropy solution u of (10),(11) satisfies the same L∞ -bound as u does: u ∞ ≤ MT .

(33)

Lemma 1 (Kruzhkov cell entropy inequality) We have for all κ ∈ R, j ∈ I h and n ∈ {0, . . . , NT − 1} the inequality |un+1 − κ| − |unj − κ| j t 1   Fj l (unj κ, unl κ) − Fj l (unj ⊥κ, unl ⊥κ) ≤− |Tj | l∈N (j )

n ˆ +sgn (un+1 j , κ)T,h [u,h (., t )](ωj ).

(34)

Proof. The proof uses only local properties of the scheme and can therefore be achieved as in the local inhomogeneous case using Definition 3 and (26) (cf. [17], for instance).   For S ⊂ R2 and n = 0, . . . , NT − 1, define ESh (n) to be set of all pairs (j, l) such that j ∈ ISh , l ∈ N (j ), and unj > unl . Lemma 2 (Weak BV-estimate) There is a constant C > 0 such that for any S ⊂ R2 , S compact and |S| ≥ 1 N T −1



t

n=0 (j,l)∈E h (n) S

(35)

−Fj l (d, d)) +

max

unl ≤c≤d≤unj

max

unl ≤c≤d≤unj

(Fj l (d, c)  1

(Fj l (d, c) − Fj l (c, c)) ≤ C|S|h− 2

and (36)

N T −1



1

|Tj ||un+1 − unj | ≤ C|S|h− 2 . j

n=0 j ∈IS

The constant C depends on MT , CTˆ (MT ), cG , u0 , ξ, T but not on h, , or S.

454

A. Dedner, C. Rohde

 NT −1 Proof. We use for the sake of brevity the notation n=0 n,(j,l) ≡  (j,l)∈ESh (n) . Weak BV-estimates of type (35),(36) for the local inhomogeneous case have been derived in [1], Proposition 3. If we proceed as in the proof of this proposition using our Assumption 2, the L∞ -estimate (25), and (26) we arrive at

 t |Ej l |

max

unl ≤c≤d≤unj

n,(j,l)

+

(Fj l (d, c) − Fj l (d, d))2 

(Fj l (d, c) − Fj l (c, c))

max

unl ≤c≤d≤unj

2

≤ C(MT , CTˆ (MT ), cG , u0 , ξ, T ) N T −1  

ˆ ,h [u,h (., t n )](ωj )| × |S| + t|Tj ||T n=0 j ∈I h S

 ≤ C(MT , CTˆ (MT ), cG , u0 , ξ, T ) |S|. For the last inequality we used the L1 (R2 )-bound on u,h which can be derived from Lemma 1, (16) and |S| ≥ 1. Note that this bound is independent of  and h. By H¨older inequality and Assumption 2 we deduce 

t

n,(j,l)

+ ≤

max

unl ≤c≤d≤unj

max

unl ≤c≤d≤unj

(Fj l (d, c) − Fj l (d, d)) 

(Fj l (d, c) − Fj l (c, c))

 2   t max (F (d, c) − F (d, d)) jl jl unl ≤c≤d≤unj |Ej l | n,(j,l)  2  21 + n max n (Fj l (d, c) − Fj l (c, c)) ul ≤c≤d≤uj

   21 × t|Ej l | n,(j,l) 1

≤ C(MT , CTˆ (MT ), cG , u0 , ξ, T )|S|h− 2 . That is statement (35). To prove (36) we use the formula (21) and obtain with (35) and once more (16)

Numerical approximation of entropy solutions N T −1



455

|Tj ||un+1 − unj | j

n=0 j ∈I h S

≤

N T −1

 

t|Fj l (unj , unl ) − Fj l (unj , unj )|

n=0 j ∈I h l∈N (j ) S

+

N T −1



ˆ ,h [u,h (., t n )](ωj )| t|Tj ||T

n=0 j ∈I h S 1

≤ C(MT , CTˆ (MT ), cG , u0 , ξ, T )|S|(h− 2 + 1).

 

Entropy residual and error estimate: Let M(R2 × [0, T ]), M(R2 ) denote the space of Radon measures on R2 × [0, T ], R2 respectively. Lemma 3 (An entropy residual for uh ) There exist measures µhxt , µhTˆ ∈ M(R2 × [0, T ]), and µh0 ∈ M(R2 ), such that we have for all φ ∈ C0∞ (R2 × [0, T )), φ ≥ 0 and κ ∈ R (37)




ETˆ ,h (uh , κ, φ) ≥ − −

T


0 T 0

For the measures µhxt 2




|φt (x, t)| + |∇φ(x, t)|) dµhxt (x, t)

R2 h φ(x, t) dµTˆ (x, t) − φ(x, 0) dµh0 (x).

R2

R2

and µhTˆ

there exists a constant C > 0 such that for each S ⊂ R , S compact, |S| ≥ 1 we have 1

(38)

µhTˆ (S × [0, T ]) ≤ CRT ()2 h 2 ,

(39)

µhxt (S × [0, T ]) ≤ C min{|S|, ||}h 2 .

1

If u0 ∈ BV (R2 ) we have (40)

µh0 (S) ≤ Ch.

The constant C depends on MT , CTˆ (MT ), cG , u0 , ξ, T but not on h, , or S. Note 5 The estimate (38) expresses the non-local influence of the integral ˆ ,h . In contrast to the local estimates (39),(40) this estimate is not operator T true when replacing RT ()2 with the eventually smaller volume of S. Proof of Lemma 3. For the smooth test function φ ≥ 0, n ∈ {0, . . . , NT −1}, and j ∈ ISh we multiply the discrete entropy inequality (34) from Lemma 1 with [t n ,t n+1 ]×Tj φ(x.t) dxdt. Adding up with respect to n and j results in (41)

Txt + TTˆ ≤ 0,

456

A. Dedner, C. Rohde

where we have Txt =

N T −1


t n+1
 1 n+1 n (|u − κ| − |uj − κ|) φ(x, t) dxdt t j tn Tj h

n=0 j ∈I S N T −1 



1  Fj l (unj κ, unl κ) |T | j n=0 j ∈I h l∈N (j ) S
n+1
 t n n φ(x, t) dxdt −Fj l (uj ⊥κ, ul ⊥κ) +

TTˆ = −

N T −1

tn



Tj

n ˆ sgn (un+1 j , κ)T,h [u,h (., t )](ωj )

n=0 j ∈I h S





t n+1 φ(x, t) dxdt. tn

Tj

Now, as in [7, Theorem 4.1], we conclude from Lemma 2 that there are measures µh0 ∈ M(R2 ) and µhxt ∈ M(R2 × [0, T ]) satisfying (39),(40) and 
T
 Txt + |uh (x, t) − κ|φt (x, t) dxdt  2 
0 T
R

 + sgn (uh (x, t), κ) (f (uh (x, t) − f (κ)) · ∇φ(x, t) dxdt  2
T 0
R


h (42) ≤ φ(x, 0) dµh0 (x). |φt (x, t)| + |∇φ(x, t)| dµxt (x, t) + 0

R2

R2

Therefore it remains to consider the global term.       T ˆ − 2 T sgn (uh (x, t), κ) T  ˆ ,h [uh (., t)](x)φ(x, t) dxdt  R 0  T N T −1     ˆ n ≤ [u (., t )](ω ) T  ,h h j  n=0 j ∈I h   t n+1 S   , κ) − sgn (u (x, t), κ) × t n Tj sgn (un+1 φ(x, t) dxdt h j  t n+1 NT −1  + n=0 j ∈ISh t n Tj |sgn (uh (x, t), κ)|  ˆ ˆ ,h [uh (., t)](x)φ(x, t) dxdt. ×T,h [uh (., t n )](ωj ) − T

The first term on the right hand side of the last inequality vanishes due to the ˆ ,h [uh (., t)](x) definition of u,h (cf. (23)). From Definition 4 we conclude T ˆ ,h [uh (., t n+1 )](ωj ) for (x, t) in (t n , t n+1 ] × Tj . We define the Radon =T measure µhTˆ for g ∈ C00 (R2 × [0, T ]) by

Numerical approximation of entropy solutions

µhTˆ , g

:=

N T −1

457

  ˆ ,h [uh (., t n )](ωj ) T

n=0 j ∈I h



t n+1  n+1 ˆ g(x, t) dxdt. −T,h [uh (., t )](ωj ) tn

Tj

ˆ ,h [uh (., t n )]− Note that the Radon measure µhTˆ is well defined since supp (T ˆ ,h [uh (., t n+1 )]) ⊂ B (, h) by (17) for n ∈ {0, . . . , NT − 1}. Altogether T we get  

T   T ˆ − ˆ ,h [uh (., t)](x)φ(x, t) dxdt  ≤ µhˆ , φ, sgn (uh (x, t), κ) T   T T R2

0

and obtain with (42) the inequality (37). To prove the estimate (38) take a family of smooth functions {gε }ε>0 such that gε ≡ 1 in S × [0, T ] for S ⊂⊂ R2 and gε L∞ ≤ 1 for ε > 0, gε → χS a.e. for ε → 0. By (18) for t ∈ [0, T ] we get µhTˆ , gε  =

N T −1

  ˆ ,h [uh (., t n )](ωj ) T

n=0 j ∈I h

  n+1 ˆ ,h [uh (., t n+1 )](ωj ) tn −T Tj gε (x, t) dxdt t N T −1    ˆ  n n+1 ˆ ≤ t|Tj |T [u (., t )](ω ) − T [u (., t )](ω ) ,h h j ,h h j  n=0 j ∈I h S N T −1

≤ CTˆ (MT )



   n t|Tj |un+1 − u j . j

n=0 j ∈I h

Now, Lemma 2 implies (38) since µhTˆ , gε  → µhTˆ (S × [0, T ]) for ε → 0.   Before we can proceed to prove an error estimate between the approximate solution u,h and the solution u of the truncated model problem (10),(2) we need some notation and an auxiliary lemma. For functions v ∈ L∞ ([0, ∞); L1 (R2 )), and S ⊂ R2 ,  > 0 we define the space-time modulus of continuity by (43) εxt (, S, u) =




T

|u(x + x, t + t) − u(x, t)| dxdt .

sup

|x|,t≤

0




S

We introduce the functions ρ ∈ C0∞ (R2 , R≥0 ), ρ¯ ∈ C0∞ (R, R≥0 ) with supp (ρ) ⊂ B1 (0), supp (ρ) ¯ ⊂ [−1, 0] and satisfying

ρ(x) dx = ρ(s) ¯ ds = 1. R2

R

458

A. Dedner, C. Rohde

Then we define the functions ρr , ρ¯r , for r ≥ 1, by (44)

ρr (x) = r 2 ρ(rx),

ρ¯r (s) = r ρ¯r (rs)

(x ∈ R2 , s ∈ R).

The proofs of the subsequent theorems rely on the following lemma which enables us to estimate the difference between two entropy solutions in an appropriate norm. ˆ 2 , the ˆ 1, T Lemma 4 (Averaged contraction) Let two admissible operators T ∞ ∞ 2 2 function u0 ∈ L (R ), and µ : C0 (R × [0, T )) → R≥0 be given. Assume that the functions u1 , u2 ∈ W (0, T ) satisfy for all  ∈ C0∞ (R2 × [0, T )),  ≥ 0 and all κ ∈ R
|u0 (x) − κ|(x, 0) dx, ETˆ 1 (u1 , κ, ) ≥ − 2 R
ETˆ 2 (u2 , κ, ) ≥ − |u0 (x) − κ|(x, 0) dx − µ[]. R2

Then we have for each ψ ∈ C0∞ (R2 × [0, T )) and r ≥ 1 the contraction estimate
T
ˆ ˆ ¯ |u1 (x, t) − u2 (x, t)|ψt (x, t) dxdt E(u1 , u2 , T1 , T2 , ψ) := 2 0
R
T + (f (u1 (x, t) u2 (x, t)) 0

R2

−f
(u1
(x, t)⊥u2 (x, t))) · ∇ψ(x, t) dxdt  T ˆ + T1 [u1 (., t)](x) R2 0  ˆ 2 [u2 (., t)](x)ψ(x, t) dxdt ≥ −R(r, ψ, u1 ). −T

(45)

For ρr , ρ¯r defined in (44) the term R is given by
T
T

ψ(x, t)ρr (x − y)ρ¯r (t − s) dµ(x, t)dyds R(r, ψ, u1 ) = 0

0

R2

R2

+ψC 1 E (r, ψ, u1 ) and

ˆ  [u1 ]) E (r, ψ, u1 ) ≤ C εxt (r, supp (ψ), u1 ) + εxt (r, supp (ψ1 ), T

1 +ε(r, supp (ψ(., 0)), u0 ) + . r

The constant C depends on u1 ∞ , u0 ∞ , but not on r and ψ. For wellˆ  [u1 ]) we posedness of the terms εxt (r, supp (ψ), u1 ) and εxt (r, supp (ψ1 ), T 2 2 extended u1 to R × R≥0 by u1 ≡ 0 in R × (T , ∞).

Numerical approximation of entropy solutions

459

Proof. The proof (of a more general form of) Lemma 4 can be found in [3]. It is an adaption of results in [7, 14, 17] to our case.   Now we can formulate and prove the main result of this section. Theorem 2 Let u0 ∈ BV (R2 ) satisfy (4). For a compact set  ⊂ R2 , || ≥ 1, consider the Cauchy problem (10),(2) with unique entropy solution u and the approximate function u,h obtained from Definition 5. Then there exists a constant C > 0 such that we have
|u (x, t) − u,h (x, t)| dxdt D

(46) ≤ C (1 + RT ()2 h1/4 + γε (h1/4 ) + RT ()2 γapp (h) . where D := supp (u,h ) ∪ supp (u ). The constant C > 0 depends on MT , CTˆ (MT ), cG , ξ , T , u0 but not on h or . Proof (of Theorem 2). We apply Lemma 4 with the choices u1 = u ,

ˆ1 = T ˆ , T

u2 = u,h ,

ˆ2 = T ˆ ,h T

ˆ , T ˆ ,h , ψ ) (cf. 45) for a partic¯  , u,h , T and consider the expression E(u ular choice of the smooth function ψ , namely ψ (x, t) = ψ1 (t)δ (|x| + Lf (t − T )),

(47)

ψ1 (t) =

(48)

e−CTˆ t − e−CTˆ T . (CTˆ + 1)

Here we put CTˆ = CTˆ (MT ) for short. The function δ ∈ C ∞ (R≥0 , R≥0 ) is supposed to satisfy (49)

δ (s) = 1 δ (s) = 0

for s ∈ [0, RT () + 1], for s ∈ [RT () + 2, ∞),

 −1 ≤ δ ≤ 0,

and Lf = max−MT ≤u≤MT {|f  (u)|}. Note that D = supp (u,h )∪supp (u ) ⊂ supp (ψ ).  ≤ 0 we deduce From the definition of ψ and δ ˆ , T ˆ ,h , ψ ) ¯  , uh , T E(u
T
≤ |u (x, t) − u,h (x, t)|ψ1 (x, t)δ (|x| + L(t − T )) dxdt R2 0
T
   ˆ ˆ + [u (., t)](x) − T [u (., t)](x) T ψ (x, t) dxdt.    ,h ,h 0

R2

460

A. Dedner, C. Rohde

We proceed by exploiting the definition of ψ1 and δ ≤ 1 ˆ , T ˆ ,h , ψ ) ¯  , uh , T (CTˆ + 1)E(u
T
≤− CTˆ |u (x, t) − u,h (x, t)|e−CTˆ t dxdt R2 0
T
  ˆ ˆ  [u,h (., t)](x) + T [u (., t)](x) − T 0

R2

0

R2

× (e−CTˆ t
T
+

− e−CTˆ T ) dxdt   ˆ ˆ ,h [u,h (., t)](x) T [u,h (., t)](x) − T

× (e−CTˆ t − e−CTˆ T ) dxdt
(50) ≤ − CTˆ |u (x, t) − u,h (x, t)|e−CTˆ T dxdt + T |RT ()2 γapp (h). D

For the last estimate we used (8) and (19). We return to the application of Lemma 4. The latter and Lemma 3 imply for  = (x, t, y, s) = ψ (x, t)ρr (x − y)ρ¯r (t − s) and r ≥ 1 the estimate ˆ T ˆ ,h , ψ ) ¯  , u,h , T E(u 
T
 ,
T
|t (x, t, y, s)| ≥− 0

R2

0

R2

+|∇(x, t, y, s)| dµhxt (x, t)dyds
T

T
+ (x, t, y, s) dµhTˆ (x, t)dyds 2 0 2 R R 0 
T

h + (x, 0, y, s) dµ0 (x)dyds − ψ C 1 E (r, ψ , u ) 0

R2

(51) =: A1 + A2 .

R2

As in [7], Lemma 5.1, we conclude from the definition of , ψ , ρr , and ρ¯r for all r ≥ 1 |A1 | ≤ ψ C 1 (r + 1)µhxt (supp (ψ ))

(52) +µhTˆ (supp (ψ )) + µh0 (supp (ψ (., 0))) . Since the initial datum u0 is in BV (R2 ) we have ε(r, supp (ψ (., 0)), u0 ) ≤ ε(r, R2 , u0 ) = O(1/r). Assumption 1 implies εxt (r, supp (ψ ), u ) ≤ εxt (r, R2 , u ) = O(1/r). We observe by (43) and (8)

ˆ  [u ]) ≤ γε 1 + CTˆ εxt (r, R2 , u ). εxt (r, supp (ψ ), T r Combining these results with ψ C 1 ≤ 1 and (13) we obtain a positive constant CBV = C(|u0 |BV , |u |BV ) such that 1

CBV (53) + γε . |A2 | ≤ r r

Numerical approximation of entropy solutions

461

If we combine (50),(51),(52),(53) we obtain  CTˆ e−CTˆ T D |u  (x, t) − u,h (x, t)| dxdt ≤ (CTˆ + 1) (r + 1)µhxt (supp (ψ )) + µhTˆ (supp (ψ ))

+µh0 (supp (ψ (., 0))) C 1

BV + γε + T RT ()2 γapp (h). + (CTˆ + 1) r r To compute µhxt (supp (ψ )), µhTˆ (supp (ψ )), and µh0 (supp (ψ (., 0))) we use Lemma 3. Then for some constant C = C(MT , CTˆ (MT ), cG , ξ, T , u0 ) >0  D |u (x, t) − u,h (x, t)| dxdt ≤ C (r + 1)RT ()2 h1/2 + RT ()2 h1/2



+h + 1r + γε 1r + RT ()2 γapp (h) . The statement of the theorem follows now by setting the free constant r equal to h−1/4 .   2.3.2 The error u − u 1 . Theorem 3 Let u0 ∈ BV (R2 ) satisfy (4). Consider the Cauchy problem (1),(2) with the entropy solution u and, for a compact set  ⊂ R2 , the Cauchy problem (10),(2) with the entropy solution u . Then there exists a constant C > 0 such that we have (54)
T
0

R2


|u(x, t) − u (x, t)| dxdt ≤ C ess sup t∈[0,T ]

R2 \

ˆ |T[u(., t)](x)| dxdt.

The constant C depends on CTˆ (MT ), MT , T but not on . Proof. We apply Lemma 4 with u1 = u,

ˆ 1 = T, ˆ T

u2 = u ,

ˆ2 = T ˆ . T

For R > 0, let ψR (x, t) = ψ1 (t)δR (|x| + L(t − T )). The function ψ1 is defined in (47) and δR ∈ C ∞ (R≥0 , R≥0 ) is a function supposed to satisfy (55) δR (s) = 1

for s ∈ [0, R],

δR (s) = 0

for s ∈ [R + 1, ∞),

δR ≤ 0.

Since u and u are entropy solutions we get in particular for r ≥ 1 by Lemma 4 (56)

ˆ T ˆ  , ψR ) ≥ ψR C 1 E (r, ψR , u). ¯ E(u, u , T,

462

A. Dedner, C. Rohde

ˆ Since u, u0 (by definition) and T[u] (by (6) and (7)) are functions in L∞ (R2 × 1 2 [0, T ]) ∩ L (R × [0, T ]) we have limr→∞ E (r, ψR , u) = 0. Then we obtain from (56)
T
|u(x, t) − u (x, t)|δR (|x| + L(t − T ))ψ1 (t) dxdt − 2 0
TR
  ˆ ˆ  [u (., t)](x)ψR (x, t) dxdt. ≤ T[u(., t)](x) − T R2

0

ˆ ˆ  [u (., t)] ∈ L1 (R2 ) t)], T Now we let R → ∞. Since u(., t), u (., t), T[u(., ∞ 2 ∩ L (R ) for almost all t ∈ [0, T ] and δR L∞ ≤ 1 we obtain from Lebesgue’s theorem
T
  u(x, t) − u (x, t)ψ  (t) dxdt − 1 0
R2
  T  ˆ ˆ ≤ T[u(., t)](x) − T [u (., t)](x)ψ1 (t) dxdt, 0

R2

and therefore with CT from (8) and the definition of the truncated operator T
T
− |u(x, t) − u (x, t)|ψ1 (t) dxdt R
2 0
T   u(x, t) − u (x, t) dxdt ≤ CT ψ1 (t) 2   R   T0  ˆ t)](x) dxdt. + 0 ψ1 (t) R2 \ T[u(., A calculation as in Theorem 2 using the definition (47) for ψ1 implies (54).   3 Application to a model problem in radiation hydrodynamics We now study a problem from fluid dynamics which models the non–local energy transport through radiation and takes the form (1), (3). To be specific we search for the nonnegative entropy solution u of the initial value problem for ˆ t)](x) ut (x, t) + divf (u(x, t)) = T[u(.,
 (57) B(u(y, t)) − B(u(x, t)) k(x, y)dy, := R2

(58)

u(., 0) = u0 ≥ 0.

Here B : R → R≥0 is given for σ > 0 by B(u) = σ u4 (u ≥ 0)

B(u) = 0 (u ≤ 0).

Numerical approximation of entropy solutions

463

It models Planck’s law of black body radiation. The kernel function k : R2 × R2 → R is given by  ˜ y)|x − y|−1 : x = y k(x, k(x, y) = (59) , 0 : x=y
π

˜ y) = exp − |x − y| dϑ (x, y ∈ R2 ). k(x, (60) sin(ϑ) 0

To see that (57) fits in the define q(u) = −4πB(u) and  framework of (1) −1 ˜ observe q(u) = −B(u) R2 k(x, y)|x − y| dy for all x ∈ R2 . The equation (57) describes the nonlinear convection of a compressible radiating fluid. The unknown u can be seen as a lumped quantity of density, momentum, and temperature. The full mathematical model for radiating plasma flow is given by the compressible Euler equations for the latter quantities, coupled with a linear Boltzmann equation for radiation intensity [20]. To derive (57) from this system under the assumption of instantaneous radiation equilibrium we refer to [3, 4]. In the sequel we will apply Theorem 1 to finite volume schemes for (57). We will do this in a more general setting that includes (57). Our analysis does not rely on the special form of B and k˜ given in (60), in fact we only suppose that k˜ : R2 × R2 → R is a smooth function and that there exists a monotone decreasing function β : R≥0 → R≥0 and a constant Ck > 0 such that
∞ β(s)ds ≤ Ck ,

(61) 0

(62) (63)

1 , |x − y| + 1 1 ˜ y)|, |∇y k(x, ˜ y)| ≤ Ck . |∇x k(x, (|x − y| + 1)2

˜ y) ≤ β(|x − y|) ≤ Ck 0 ≤ k(x,

We note that we merely assume that the integral operator is weakly singular. For the function B it suffices to suppose B ∈ C 1 (R) and B  ≥ 0. For M > 0 we define CB (M) := maxw∈[−M,M] |B  (w)|. As the first step towards the error estimate we will verify that the operˆ belongs to the class of admissible operators defined in Section 2.1, ator T Definition 1. ˆ from (57), (59) such that (61), (62), (63) Lemma 5 Consider the operator T ˜ are satisfied for the function k. ˆ is an admissible operator, i.e. there exists in particular a function Then T CTˆ : R≥0 → R≥0 such that (6), (7), and (8) hold.

464

A. Dedner, C. Rohde

Furthermore there exists a constant C > 0 such that we have for w ∈ L∞ (R2 ) ∩ BV (R2 ) and  > 0 ε(, R2 , T[w]) ≤ C.

(64)

The constant C depends only on B, Ck , w∞ , w1 , |w|BV . Proof. Straightforward computations show that (6), (7), and (8) are satisfied for each w ∈ L1 (R2 ) ∩ L∞ (R2 ) with w∞ ≤ M, M > 0 by the choice

CTˆ (M) = 4π CB (M)Ck . This implies also condition (i) of Definition 1. It remains to show the estimate on γε . Let w ∈ L∞ (R2 ) ∩ BV (R2 ). For x ∈ R2 fixed with |x| < ,  > 0 we consider
ˆ ˆ |T[w](x + x) − T[w](x)|dx R2



  k(x,   ˜ y) ≤ dydx B(w(y)) − B(w(y + x)) |x − y| R2 R2

  k(x ˜ y)    ˜ + x, y + x) − k(x, + B(w(y + x)) dydx   |x − y| R2 R2
 
 k(x ˜ y)     ˜ + x, y + x) − k(x, + B(w(x + x))  dydx |x − y| 2 2 R R
 
k(x, ˜ y)   + B(w(x + x)) − B(w(x)) dydx |x − y| R2 R2

  β(|x − y|)   ≤ 2CB (w∞ ) w(y + x) − w(y) dxdy |x − y| 2 2 R
R     +2CB (w∞ ) w(y + x)
× R2

(65)

R2

˜ + x, y + x) − k(x, ˜ y)| |k(x dxdy |x − y|

=: I1 + I2 .

By the assumption w ∈ BV (R2 ) and (61), (62) we get I1 ≤ C1 |w|BV , where C1 depends only on B, w∞ , and Ck . Due to the smoothness of k˜ we can expand k˜ around (x, y). Exploiting (63) leads to I2 ≤ C2 w1 with C2 depending solely on B, w∞ , and Ck . Both results imply ε(, R2 , w) ≤ C for some constant C > 0.  

Numerical approximation of entropy solutions

465

Now we turn to the discrete setting. Consider a time step t > 0 and, for h ∈ (0, h0 ], grids Th as defined in Section 2.2.1 such that Assumption 2 holds. Furthermore, for j ∈ I h , let ρj = ρj (h) be the radius of the biggest circle with midpoint ωj that lies in the closure of Tj . We define ρ = inf{ρj |j ∈ I h }.

(66)

Note that ρ depends on h. Additionally to Assumption 2 assume that there is a constant cˆ > 0 independent of h such that h ˆ 2 ≤ cˆ and |Tj | ≤ ch ρ

(67)

(j ∈ I h ).

For later use consider Ilh (x) := {j ∈ I h | lρ ≤ |x − ωj | < (l + 1)ρ}, l ∈ N and x ∈ R2 . Note that there exists a constant c˜ independent of h, ρ, and x ˜ such that card(Ilh (x)) ≤ cl. ˆ ,h : Vh → Vh for grid functions wh ∈ Vh Finally, consider the operator T given by   ˆ ,h [wh ](x) = χ (x) T B(wh (ωi )) i∈I h \{j }

(68) (69)

−B(wh (ωj )) |Ti |kh (ωj , ωi ) (x ∈ Tj , j ∈ I h ), ˜ y) k(x, (x, y ∈ R2 ). kh = kh (x, y) = |x − y| + h

Here  is an arbitrary compact subset of R2 . Obviously, in (68) it would have been not necessary to exclude j from the summation since the summand is 0 for j = i. However, this notation is more convenient in the sequel. From the definition of ρ we deduce for j ∈ I h     β(|ωj − ωi |) ≤ C(c, ˆ Ck ) |Ti |kh (ωj , ωi ) ≤ h2 hβ(lρ). |ωj − ωi | + h h h l∈N l∈N i∈I \{j }

i∈Il (ωj )

The that  last sum is bounded due to (61). In the same wayh one shows |Ti |kh (ωi , ωj ) is bounded. We introduce the constant Ck = Ckh (c, ˜ c, ˆ Ck ) such that   (70) |Ti |kh (ωj , ωi ), |Ti |kh (ωi , ωj ) ≤ Ckh . i∈I h \{j }

i∈I h \{j }

ˆ ,h from (68) is an admissible discrete Lemma 6 The discrete operator T operator in the sense of Definition 4. A function γapp such that (19) is satisfied for each wh ∈ Vh ∩ L∞ (R2 ) ∩ 1 L (R2 ) with wh ∞ ≤ M, M > 0, is given by (71)

γapp (h) = Ch| ln(h)| (h ∈ (0, h0 ]),

ˆ c) ˜ does not depend on h. where C = C(B, M, Ck , c,

466

A. Dedner, C. Rohde

Proof. For wh given as in the statement we put wj = wh (ωj ), j ∈ I h . To verify (15) we observe by (70) for all j ∈ I h and x ∈ Tj  ˆ ,h [wh ](x)| ≤ CB (M)M |T |Ti |kh (ωj , ωi ) i∈I h \{j }

|Ti |kh (ωi , ωj ) ≤ C(B, M, Ckh ).



+

i∈I h \{j }

To verify (16) we observe again by (70)   ˆ ,h [wh ] ≤ CB (M) |Tj | (|wi | + |wj |)|Ti |kh (ωj , ωi ) T 1 j ∈I h

i∈I h \{j }

≤ C(B, M, Ckh )wh 1 . Property (17) is clear. The estimate (18) follows similar to (16). Let also w˜ h ∈ Vh ∩ L∞ (R2 ) ∩ L1 (R2 ), w˜ h ∞ ≤ M. We get    ˆ  ˆ |Tj |T ˜ h ](ωj ) ,h [wh ](ωj ) − T,h [w j ∈I h

≤ 2CB (M)



 |Tj | 

j ∈I h

+





|wi − w˜ i ||Ti |kh (ωj , ωi )

i∈I h \{j }



|wj − w˜ j ||Ti |kh (ωj , ωi )

i∈I h \{j }

≤ C(B, M, Ckh )



|Tj ||wj − w˜ j |.

j ∈I h

It remains to prove (19). For some x ∈ Tj , j ∈ I h , consider   T ˆ ,h [wh ](x) ˆ  [wh ](x) − T   
   (B(wi ) − B(wj )) k(x, y) − kh (ωj , ωi )dy  ≤ i∈I h \{j }

Ti

 
≤ C(B, M) |k(x, y) − kh (x, y)|dy +


i∈I h \{j }

+

|kh (x, y) − kh (x, ωi )|dy

Ti




i∈I h \{j }

i∈I h \{j } Ti

 |kh (x, ωi ) − kh (ωj , ωi )|dy

Ti

=: C(B, M) (A1 (x) + A2 (x) + A3 (x)) .

Numerical approximation of entropy solutions

467

Using the definition of ρ in (66) and (62) it follows
˜ y) hk(x, dy A1 (x) ≤ |x − y|(|x − y| + h) R2 \Bρ (ωj )




≤ hCk R2 \Bρ (ωj )

1 dy |x − y|(|x − y| + h)(|x − y| + 1)

≤ C(Ck )h| ln(h)|. To estimate A2 we use Taylor expansion for the smooth function kh :
 A2 (x) ≤ ∇y kh (x, .)L∞ (Ti ) |y − ωi |dy i∈I h \{j }

≤h

(72)

3

 

Ti

∇y kh (x, .)L∞ (Ti ) .

l∈N i∈I h (x) l

Since h < 1 and due to (61), (62) we can bound |∇y kh (x, y)| for y ∈ Ti by Ck which implies (|x−y|+1)(|x−y|+h)2 ∇y kh (x, .)Ti ≤

Ck . (lρ−h+1)(lρ)2

Note that we have |x − y| ≥ |x − ωi | − |ωi − y| ≥ lρ − h for all y ∈ Ti . With (67) we deduce by elementary calculus ˆ c)h| ˜ ln(h)|. A2 (x) ≤ C(Ck , c, Similar calculations as for A2 show that we have for all i ∈ Ilh (ωi ) ∇x kh (., ωi )Tj ≤

Ck . (lρ − h + 1)(lρ)2

Therefore we arrive at the same bound for A3 as for A2 . Adding the three estimates for A1 , A2 , A3 leads us to
  T[w ˆ h [wh ](x) ≤ C(B, M, Ck , c, ˆ h ](x) − T ˆ c)|T ˜ j |h| ln(h)|.   Tj

With Lemma 6 at hand and Definition 5 we obtain a fully explicit finite volume scheme to approximate entropy solutions of (57), (58). Note that although the scheme can be used to approximate (57), (58) the computation cost for evaluating the discrete operator grows quadratically with the number of elements. Approximations which only grow linearly are available (e.g. [5] and the references therein) but any approximation results in the sense of Definition 4 known to the authors require some additional smoothness assumptions which are not fulfilled in the setting studied here.

468

A. Dedner, C. Rohde

Lemma 7 Let u0 ∈ L∞ (R2 ) be nonnegative and Ckh given by (70). If the CFL-like condition t ≤

(73)

2 h cG KLF (u0 ∞ ) + CB (u0 ∞ )Ckh

holds we have for all h ∈ (0, h0 ] and compact sets  ⊂ R2 (74) ess inf {u0 (x)} ≤ u,h (x, t) ≤ ess sup{u0 (x)} (x ∈ R2 , t ∈ [0, T ]). x∈R2

x∈R2

This implies in particular that the approximation u,h takes only values in the physically relevant set R≥0 . Proof. For all n ∈ {0, . . . , NT }, i, j ∈ I h , and l ∈ N (j ), we define Fjnl

:=

Fj l (unj , unl ) − Fj l (unj , unj ) unj − unl

,

Bjni

:=

B(uni ) − B(unj ) uni − unj

provided unj = unl , unj = uni . Otherwise set Fjnl = 0 and Bjni = 0. Rewriting formula (21) for n = 1 and the discretization (68) we obtain  t  Fj0l (u0j − u0j l ) + t Bj0i (u0i − u0j )|Ti |kh (ωj , ωi ) |Tj | l∈N (j ) i∈I h     t 0 0 0 Fj l − t Bj i |Ti |kh (ωj , ωi ) = uj 1 − |Tj | h

u1j = u0j −

+

t |Tj |



l∈N (j )

Fj0l u0j l + t

l∈N (j )

i∈I



Bj0i u0i |Ti |kh (ωj , ωi ).

i∈I h

Due to Definition 3(iv) we have Fj0l ≥ 0. Due to B  ≥ 0 for u ≥ 0, and (62) we have also Bj0i kh (ωj , ωi ) ≥ 0. This implies the lower bound in 0≤

 1  Fj0l + Bj0i |Ti |kh (ωj , ωi ) |Tj | h l∈N (j )

i∈I

KLF (u0 ∞ ) ≤ + CB (u0 ∞ )Ckh . 2 cG h The upper bound follows from Definition 3, (14), and (70). Recall that K denotes the maximal number of edges of a cell volume.

Numerical approximation of entropy solutions

469

Due to the CFL condition (73) and the compact support of u,h we have    t  1 0 0 0 uj ≤ max{uj } 1 − Fj l − t Bj i |Ti |kh (ωj , ωi ) |Tj | j ∈I h h l∈N (j )

+ max{u0j } j ∈I h

t |Tj |



j ∈I

Fj0l + max{u0j }t j ∈I h

l∈N (j )



Bj0i |Ti |kh (ωj , ωi )

j ∈I h

≤ max{u0j } ≤ ess sup{u0 (x)}. j ∈I h

x∈R2

In the same way we can prove that u1j ≥ ess inf {u0 (x)}. x∈R2

With a standard induction argument we get the estimate (74).

 

In the concluding theorem we summarize all lemmata above to obtain a convergent finite volume algorithm for the approximation of entropy solutions of (57), (58). As in Corollary 1 we let the domain  depend on h, i.e.  = h . Since in the case of this special model problem much stronger estimates than in the general case of Corollary 1 are available we obtain convergence with only mild (in fact logarithmical) growth of h for h → 0. If we assume additionally the natural decay rate for the exact entropy solution we recover (up to a logarithmically small factor) the classical convergence rate 1/4 as obtained in the scalar homogeneous case ([2,7,13,16, 19]). Theorem 4 Let u0 ∈ L∞ (R2 ) ∩ BV (R2 ) have compact support. We assume that there is a unique entropy solution u ∈ W (0, T ) of the problem (57),(58) such that u(., t) ∈ BV (R2 ) for t ∈ [0, T ]. The same is supposed for the associated truncated problem. For h ∈ (0, h0 ], let an unstructured grid Th satisfying Assumption 2 and an the associated family of monotone numerical fluxes Fh be given. Consider √ ˆ ,h from (68) with  = h = {x ∈ R2 | |x| ≤ | ln(h)|}. discrete operator T Suppose furthermore that the CFL-like conditions (26) and (73) hold. Then the sequence of approximate solutions {uh ,h } from Definition 5 satisfies (75)

lim u − uh ,h 1 = 0.

h→0

If additionally the decay estimate (76)

ˆ T[u(., t)](x) = O(exp(−|x|)) (|x| → ∞)

holds for all t ∈ [0, T ] we obtain for some positive constant C independent of h the error estimate (77)

1

u − uh ,h 1 ≤ Ch 4 | ln(h)|.

470

A. Dedner, C. Rohde

Proof. The Lemmata 5, 6, 7 show that the conditions of Theorem 1 are satisfied. Note that we can choose γε (h1/4 ) = Ch1/4 from (64). With (71) and ˆ T[u(., t)] ∈ L1 (R2 ) we get the same convergence statement (75) as in Corollary 1. Straightforward calculus using (76) establishes the estimate (77) from (27).   References 1. Chainais-Hillairet, C., Champier, S.: Finite volume schemes for nonhomogeneous scalar conservation laws. Error estimate. Numer. Math. 88(4), 607–639 (2001) 2. Cockburn, B., Coquel, F., Lefloch, P.: An error estimate for finite volume methods for multidimensional conservation laws. Math. Comput. 63(207), 77–103 (1994) 3. Dedner, A., Rohde, C.: Existence uniqueness and regularity of weak solutions for a model problem in radiative gas dynamics. In preparation 4. Dedner, A.: Solving the system of radiattion magnetohydrodynamics. Albert– Ludwigs–Universit¨at, Fakult¨at f¨ur Mathematik und Physik. Freiburg, http://www. freidok.uni-freiburg.de/volltexte/1098/ (September 2003) 5. Dedner, A., Vollm¨oller, P.: An adaptive higher order method for solving the radiation transport equation on unstructured grids. J. Comput. Phys. 178(2), 263–289 (2002) 6. Engelberg, S., Liu, H., Tadmor E.: Critical thresholds in Euler-Poisson equations. Indiana Univ. Math. J. 50, 109–157 (2001) 7. Eymard, R., Gallou¨et, T., Ghilani, M., Herbin R.: Error Estimates for the Approximate Solutions of a Nonlinear Hyperbolic Equation Given by some Finite Volume Schemes. IMA J. Numer. Anal. 18, 563–594 (1998) 8. Eymard, R., Gallou¨et, T., Herbin, R.: Finite volume methods. In: Ciarlet P.G. et al., (eds.), Handbook of numerical analysis. Vol. 7, Solution of equations in Rn (Part 3). Techniques of scientific computing (Part 3). North-Holland/Elsevier, 2000 9. Hamer, K.: Nonlinear effects on the propagation of sound waves in a radiating gas. Quart. J. Mech. Appl. Math. 24, 155–168 (1971) 10. Katsaounis, T., Makridakis, C.: Finite volume relaxation schemes for multidimensional conservation laws. Math. Comput. 70(234), 533–553 (2001) 11. Kawashima, S., Nikkuni, Y., Nishibata, S.: The initial value problem for hyperbolicelliptic coupled systems and applications to radiation hydrodynamics. Freist¨uhler H. Analysis of systems of conservation laws. Chapman & Hall/CRC, 1998 12. Kr¨oner, D.: Numerical schemes for conservation laws. Wiley, 1997 13. Kr¨oner, D., Rokyta, M.: Convergence of upwind finite volume schemes for scalar conservation laws in two dimensions. SIAM J. Numer. Anal. 31(2), 324–343 (1994) 14. Kruzhkov, S.N.: First order quasilinear equations in several independent variables. Math. USSR Sb. 10, 217–243 (1970) 15. Mihalas, D.: Radiation hydrodynamics. Steiner O. et al., Computational methods for astrophysical fluid flow. Springer, 1998 16. Noelle, S.: A note on entropy inequalities and error estimates for higher-order accurate finite volume schemes on irregular families of grids. Math. Comput. 65(215), 1155–1163 (1996) 17. Rohde, C.: Upwind finite volume schemes for weakly coupled hyperbolic systems of conservation laws in 2D. Numer. Math. 81(1), 85–123 (1998) 18. Treiber, M., Hennecke, A., Helbing, D.: Derivation Properties and Simulation of a Gas-Kinetic-Based Non-Local Traffic Model. Phys. Rev. E. 59, 239–253 (1999)

Numerical approximation of entropy solutions

471

19. Vila, J.-P.: Convergence and error estimate in finite volume schemes for general multidimensional conservation laws. Math. Model. Numer. Anal. 28, 267–285 (1994) 20. Vincenti, W.G., Krueger, C.H.: Introduction to physical gas dynamics. Wiley, 1965 21. Zel’dovich, Y.B., Rajzer, Y.P.: Physics of shock waves and high temperature hydrodynamic phenomena. Academic Press, 1966