Two a posteriori error estimates for one-dimensional scalar conservation laws1 Laurent Gosse2 and Charalambos Makridakis3 Abstract
In this paper, we propose a-posteriori local error estimates for numerical schemes in the context of one dimensional scalar conservation laws. The main tool to derive them is a synthetic version of Kruzkov's estimates recently introduced by Bouchut and Perthame, further developed by Katsoulakis, Kossioris and Makridakis. We rst consider the schemes for which a consistent in-cell entropy inequality can be derived. Then, we extend this result to second order schemes written in viscous form satisfying weak entropy inequalities. As an illustration, we show several numerical tests on the classical Burgers equation and we propose an adaptive algorithm for the selection of the mesh.
Key words: conservation law, entropy dissipation, a posteriori error estimates, adaptive nite di erence
method.
AMS subject classi cation: 65M15, 65M50. 1 Work partially supported by TMR project HCL N . ERBFMRXCT960033 2 Foundation for Research and Technology-Hellas, Institute of Applied and Computational Mathematics - P.O. Box 1527 - 71110 Heraklion - Crete, Greece:
[email protected] 3 Department of Mathematics, University of Crete - 71409 Heraklion-Crete, Greece and Foundation for Research o
and Technology Hellas, Institute of Applied and Computational Mathematics - P.O. Box 1527 - 71110 Heraklion - Crete, Greece:
[email protected]
1
1 Introduction We consider one-dimensional nite di erence schemes for the approximation of the scalar conservation law: ut + f (u)x = 0; x 2 R; t > 0 with initial data in the space of functions of bounded variation. In this paper, we propose a simple framework which is used to derive explicit a posteriori error estimates in the local L1 norm for this class of numerical schemes. It is mainly based on a recent result of Bouchut and Perthame [3], further developed by Katsoulakis, Kossioris and Makridakis [14]. In these papers, a synthetic version of Kruzkov-type error estimates is derived starting from an approximate entropy inequality the approximate solution is supposed to satisfy. From a numerical point of view, this provides a way to convert any in-cell discrete entropy inequality into an error estimate. In the literature, these estimates are used to derive convergence rates mainly for monotone schemes [8] or E-schemes [26, 28]; for the higher-order numerical schemes, there are usually no consistent discrete entropy inequalities available, with the exception of [2]. In this work we derive local a posteriori error estimates of the form Z
K
ju(x; t) ? vh(x; t)j dx E (vh; t; K )
(1.1)
where vh denotes the approximation of u obtained by an E-scheme or a (possibly) high-order scheme, K is a typical cell of the space partition, and E (vh ; t; K ); is a computable quantity that depends on the data of the problem and on values of vh in appropriate extended cones of dependence, cf., Theorems 3, 4. Error estimates of the form (1.1) are proposed for schemes satisfying two types of entropy inequalities. First we consider schemes that satisfy \strong" (consistent) entropy inequalities, Theorem 3. Then starting from general (high-order) schemes written in viscous form, we are able to derive a very weak entropy inequality which is not consistent in the sense of Lax and Wendro [19]. In the nal estimate this lack of consistency gives rise to an additional error term in the error estimate of Theorem 4. It seems though, that this term is in fact of higher order, if compared with the other terms in the estimate, cf., Remarks 1,3. Indeed it is shown that under stronger CFL conditions that this term is of higher order in the case of convex uxes, Remark 3. In the general case the numerical tests con rm that the contribution of this term is actually very small. Note that estimates of the form (1.1) are meaningful even for schemes for which we do not know if they are convergent or not. Indeed if our computed approximations are such that E ; given tolerance then the error in (1.1) is at most : In particular these estimates can be used in computations with, e.g., MUSCL schemes even when the ux is not convex. Estimates of the form (1.1) can be used to design adaptive algorithms. As an illustration, we present some numerical runs on the modi ed Lax-Friedrichs rst-order scheme introduced by Tadmor [32] and on its second-order extension obtained by means of the piecewise constant viscosity modi cation proposed by Osher and Tadmor [28]. These tests involve some Riemann problems and the computation of a stationary shock for the inviscid Burgers equation. In 2
addition, we present some preliminary computations based on an adaptive algorithm using the estimate proposed in theorem 4. The idea of obtaining a posteriori bounds based on the estimates in [16, 17] was originated in [6]. In [13] a posteriori estimates in a weighted space-time L2 norm for nite element schemes approximating the "-viscous perturbation of one-dimensional systems of conservation laws are obtained. Estimates in a dual Lip+ norm for convex uxes were proved in [23]. See also [20, 21, 1, 29] and the their references for adaptive schemes for conservation laws. The results of this paper were announced in [11]. In [11] a posteriori estimates for monotone nite volume schemes were also derived, cf. also [15] for similar results in the space-time L1 norm. This paper is organized as follows: in the sequel, we introduce the schemes and our notation. We consider the coecients of Harten [12] which are of great use in this work, and we present a variation of the error theorem from [3, 14]. In the section 2, we derive a local a posteriori error estimate for entropy consistent schemes using a strong in-cell discrete entropy inequality. In the section 3, we extend this result to more general schemes written in viscous form by means of a weaker discrete entropy inequality. To prove the main result of this section, some lemmas are needed and their proofs are forwarded in the appendix at the end of the paper. Then, in the section 4, we present some numerical computations using the modi ed Lax-Friedrichs scheme and its second order extension. This provides a way to compare the estimated and the exact error for these numerical approximations. As a nal illustration, we display some approximate solutions generated by an adaptive second order algorithm based on our second error estimate. Preliminaries. We are interested in the Cauchy problem for the following partial di erential equation: ut + f (u)x = 0 with (x; t) 2 R R+ (1.2) u(:; 0) = u0 2 BV (R) where f 2 C 2 (R) and BV (R) denotes the space of functions of bounded variation, [9]. Problem 1.2 admits a unique entropy solution, [16], i.e. for all positive 2 C01 (R R+ );
Z
RR+
[U (u) t + F (u) x ]dxdt +
Z
R
U (u0 )(x) (x; 0)dx 0;
where U 2 C 2 (R) is strictly convex and F 0 = U 0 f 0 : In the case of strictly convex uxes, f 00 a; a > 0; the entropy solution u satis es, Oleinik [25], the following one-sided Lipschitz estimate: 1 u ( x; t ) ? u ( y; t ) def + 2 = k u ( :; t ) k 8(t; s) 2 (R ) ; t s; sup + (R) Lip 1 x?y x6=y + ku(:;s)k + (R) + a(t ? s) (1.3) For each a 2 R, the notation (a)+ stands for: max(a; 0). Numerical schemes and approximate solutions. Following for instance [9], we recall here the basic properties of the most usual nite di erence schemes for (1.2). We introduce consequently a space-time discretization de ned by a uniform time-step t > 0 and a countable collection of grid points xj + 21 for j 2 Z. Let also Lip
h ; h = sup hj hj = xj+ 21 ? xj? 21 ; 0 < h = jinf 2Z j j 2Z
3
(1.4)
Then, we denote xj the middle-points of the cells Kj def = [xj + 21 ; xj ? 12 [: xj = xj + 21 ? (hj =2) = xj? 21 + (hj =2). Similarly, we de ne time intervals of the type: I n def = [tn ; tn+1 [ where tn = nt for n 2 N . In all this work, we will consider piecewise constant numerical approximations vh of u generated by conservative algorithms in the sense of Lax and Wendro [19].
v h : R R+ ! R (x; t) 7! vh (x; t) = vjn if x 2 Kj ; t 2 I n
(1.5)
R
Let us de ne the following discretization for u0 : 8j 2 Z; vj0 = h1 K u0 (x)dx: Then we consider schemes that are constructed using an consistent \essentially 3-point" numerical ux function g: In particular let g; g : R2p ! R be a Lipschitz continuous function such that: j
j
= g(vjn?p+1 ; :::; vjn+p ) (vjn?p+1; :::; vjn+p ) ! 7 gjn+ 21 def g(v; :::; v) = f (v) : The explicit schemes we will consider in this section are written as: vn+1 = vn ? t (gn ? gn ) j
j
hj
j + 12
j ? 12
(1.6)
(1.7)
Following [28], we write these (possibly high-order) numerical schemes in the equivalent form vjn+1 = vjn ? 2ht (fjn+1 ? fjn?1 ) + 2ht (Qnj+ 21 (vjn+1 ? vjn) ? Qnj? 21 (vjn ? vjn?1 )) (1.8) j j where the quantities fjn and Qnj+ 1 are respectively called modi ed ux and numerical viscosity 2 for the scheme. We assume that the modi ed ux is of the form fjn = f (vjn) + f~jn (1.9) f~jn depends on vjn?1 ; vjn ; vjn+1 In [28] was shown that these schemes are SOR (second order accurate in the smooth regions) as soon as the following requirement is ensured: 8 > < > :
~
~
+f + O(jvjn+1 ? vjn j) Qnj+ 21 = (anj+ 21 )2 + fv +1 +1 ?v
anj+ 21 = f (vv+1+1)??fv (v ) n j n j
n j
n j
n j n j
n j n j
(1.10)
More precisely, second order accuracy is given up in the neighborhood of nonsonic extremal points. We mainly refer to [28] for detailed results. In order to study the convergence properties of the schemes (1.7), (1.8), it is convenient to introduce their discrete spatial total variation:
TV (vh )(:; tn ) = 4
X
j 2Z
jvjn+1 ? vjnj
Moreover, we will use the following notation: def = t=h. Following Harten [12], we de ne now the coecients mnj+ 1 ; pnj? 1 which are going to be widely used throughout this paper: 2
2
8
x 2 K j ) R < A = R(RR+ K B = R(RR+ )2 (x; t)j x (x; t; y; s)jdxdtdyds > : C= (RR+ )2 (x; t)j t (x; t; y; s)jdxdtdyds j
j
We precise now the choice of the functions and : for any positive constants ; , we de ne: (x; t) = x(x) t (t), where x; t satisfy R
t (t)dt = 1; R x (x)dx = 1 R t (t) = 1t (t=)= and supp(1t ) ] ? 1; 0[ x(x) = x (x=)= and supp( x ) ] ? 1=4; 1=4[ R
1
1
Then j1t (t)j (1 + m), j1x (x)j 2(1 + m) with m > 0. Similarly, weR de ne a regularized Heaviside function Y such that Y (?1) = 0, Y0 (t) = Y10 (t=)= and Y10 (t)dt = 1. With another parameter " > 0, we set 0 (t) = Y"(t) ? Y"(t ? T ) 1 2 C 1 (]0; T + "[) and, nally, we get the expression of 2 C01 (R R+ ): (x; t) = (t):[1 ? Y (jxj ? R ? =2 ? M (T ? t)] def = (t) (x; t) 1 which remains smooth as long as M" R + =2. We also deduce the corresponding bounds for the derivatives of : jx (x; t)j max2R jY10 ( )j=, jt (x; t)j j 0 (t)j + M: max2R jY10 ( )j= and j 0 (t)j 2 max2R jY10 ( )="j where we have once again jY10 (t)j (1 + m) with m > 0. Because of the Lipschitz property coming from the ux function regularity and the uniform bound on the amplitude of the solutions, we get:
?
Z
(RR+ )2
ju(x; t) ? vh(y; s)j 0 (t) (x; t) (x ? y; t ? s)dxdtdyds A + B + C
Taking " ! 0 in the above relation and using the bounds R
A jxjR+M (T ?t)+ (x; t)dxdt R 0 B kY1 k 1 + 2 k1k 1 jxjR+M (T ?t)+ (x; t)dxdt R 0 C M kY1 k 1 n+ 2 k1 k 1 jxjR+M (T ?t)+ (x; t)dxdto R +2kY10 kL1 sup0tT + jxjR+M (T ?t)+ (x; t)dx x
L
L
L
t
L
we derive the desired bound following the arguments in [3, 14]. 6
2 The case of strongly entropy consistent numerical schemes In this section, we consider the particular case of numerical schemes which are endowed with a strong consistency relation with the continuous problem (1.2). More precisely, this means that for all entropy functions U , there exists a consistent and Lipschitz continuous numerical entropy
ux G such that the following inequality holds: U (vn+1 ) ? U (vn ) + t (Gn ? Gn ) 0 (2.1) j
j
hj
j + 21
j ? 21
with the same notation conventions as before. This requirement (2.1) is met for example in the case of E-schemes [26, 28] or monotone schemes [8] (which are special E-schemes). We also know that this class of algorithms is most of the times limited to (formal) rst-order accuracy, so we can restrict ourselves to numerical approximations generated by a 3-points scheme corresponding to the value p = 1 in (1.6): (2.2) vn+1 = vn ? t [g(vn ; vn ) ? g(vn ; vn )] j
j
j
hj
j +1
j ?1 j
Thanks to the regularity requirements, it is possible to rewrite them under a unique viscous form determined by the viscosity coecients denoted Qnj+ 1 : 2
vjn+1 = vjn ? 2ht (f (vjn+1 ) ? f (vjn?1)) + 2ht (Qnj+ 21 (vjn+1 ? vjn ) ? Qnj? 21 (vjn ? vjn?1 )) j
(2.3)
j
In this case, any modi ed ux used in (1.8) boils down to the continuous one. We can also simply de ne in a unique way Qnj+ 1 = [f (vjn ) + f (vjn+1 ) ? 2gjn+ 1 ]=(vjn+1 ? vjn ) as a function 2 2 of the two adjacent values (vjn ; vjn+1 ). It is possible to derive an analytical expression of these coecients at least for the formal rst-order schemes like Godunov, Lax-Friedrichs, MurmanRoe, Engquist-Osher [28]. Of course, provided we use a suciently small time-step, we are in position to use the theorem 1 to bound the resulting numerical approximation vh for which one gets a somehow simpli ed form for the CFL restriction and the incremental coecients: f (v n ) j +1 N; vjn+1
8(j; n) 2 Z
? f (vjn) n 1 and Qj + 1 n 2 2 ? vj
8 < :
mnj+ 1 = Qnj+ 1 ? f (vv+1+1)??fv (v ) 2 2 n n pj? 1 = Qj? 1 + f (vv )??fv (?v 1?1 ) 2
2
n j n j n j n j
n j
n j n j
n j
Derivation of the error estimate. The main purpose of this paragraph is to prove the
following result: Theorem 3. Assume that the restrictions of theorem 1 are met for the scheme (2.2) which is supposed to satisfy (2.1) for all Lipschitz continuous entropies U . Then we have for all T > 0: x + 1 +MT + 2 jvh(x; 0) ? u(x; 0)jdx K x ? 1 ?MT ? 2 N 1 X p p + 2C (3 2 + 1) TV (u0 ) tXVn 2 + C sup fXVn g n=0;:::;N +1 n=0
Z
jvh (x; T ) ? u(x; T )jdx
j
Z
j
j
7
(2.4)
where
XVn =
X
`12J
p
n j;
h` Qn`+ 12 jv`n+1 ? v`nj with Jj;n = f` 1 : x` 2 B (xj ; hj =2 + M (T ? t) + )g
= 3 MTh=2; C = 2 M = supfjf 0 ()j with jj ku0 kL1 (R)g; N = T=t
(2.5) Proof. We have to compute the following quantity for all k 2 R and all positive test-functions ' in C01 (R R+ ):
I=?
Z
(RR+ )
[jvh ? kj't + sgn(vh ? k)(f (vh ) ? f (k))'x ](x; t)dxdt
(2.6)
It is therefore sucient to only consider the one-parameter family of Kruzkov's entropies: U k (u) = ju ? kj, F k (u) = sgn(u ? k)(f (u) ? f (k)). Since vh is a piecewise constant function, we can rewrite I as:
I=
X
j;n
[U k (vjn+1 ) ? U k (vjn )]
Z |
K
j
'(x; tn+1 )dx +[F k (vjn+1 ) ? F k (vjn)] {z
' +1
}
Z
|I
n Kj
n
'(xj+ 21 ; t):dt {z
}
' + 1 2 In
(2.7)
j
Thanks to the strong entropy requirement (2.1), there exists a consistent numerical entropy ux Gk ; (Gk (vjn ; vjn ) = F k (vjn ) ) associated to each U k : Thus one can show that
I
n+1 n k n n k n n I k n+1 k n j;n[U (vj ) ? U (vj )]('K ? 'j =t) + [G (vj ; vj +1 ) ? G (vj ; vj )]('j + 12 +[Gk (vjn ; vjn ) ? Gk (vjn?1 ; vjn )]('Ij ? 1 ? 'nj =hj ) 2 P
n
j
n
? 'nj =hj ) (2.8)
To estimate the right hand side of (2.8) we proceed as follows jU k (vjn+1 ) ? U k (vjn )j jvjn+1 ? vjnj " " n ) ? f (vn ) # n ) ? f (vn ) # f ( v f ( v j jvn ? vn j + (=2) Qn + j j ?1 jvn ? vn j; (=2) Qnj+ 21 ? vj+1 j +1 j j j ?1 n ? vn j ? 21 vn ? vn j +1
{z
|
m + 1 0 2
j
}
n
j
|
{z
p ? 1 0 2 n
j
j ?1
}
j
(2.9)
and
jGk (vjn; vjn+1 ) ? Gk (vjn ; vjn )j + jGk (vjn; vjn ) ? Gk (vjn?1; vjn)j M [jvjn+1 ? vjnj + jvjn ? vjn?1j] R R Moreover, j'nK+1 ? 'nj =tj K I j't (x; t)jdxdt and j'Ij 1 ? 'nj =hj j K I j'x (x; t)jdxdt : 2 At this point, we can de ne the coecients ; ; used in Theorem 2 for (x; t) 2 Kj I n : (x; t) 0 (x; t) = M [jvjn+1 ? vjnj + jvjn ? vjn?1 j]
(x; t) = (=2)[mnj+ 1 jvjn+1 ? vjnj + pnj? 1 jvjn ? vjn?1 j]: 2 2 n
j
j
n
j
8
n
According to this theorem, we have the following error estimate for all T > 0: R
R
h jxjR+MT + ju(s; 0)n? vh(s; 0)j:ds jxjR ju(s; T ) ? v (s; T )j:ds o ? R R +2C ( + M)TV (u0 ) + 0T B(x ;h =2+M (T ?t)+) (s; t) + 4C (s;t) + C 2M + 1 (s; t) ds:dt R +C sup0tT + B(x ;h =2+M (T ?t)+) (s; t)ds j
j
j
j
with B (xj ; hj =2+ M (T ? t)+) = [xj ? 21 ? M (T ? t) ? ; xj + 12 + M (T ? t)+] for all t 2 [0; T ]. Since the functions ; are piecewise constant, the right-hand side can be rewritten as a discrete summation. In fact, we denote by N the value T=t and we consider the following two terms: n
P
o
P
E1 = 2C TV (u0) + M Nn=0 t `2J (4 + :Qn`+ 21 ) h` jv`n+1 ? v`nj P P E2 = 2C MTV (u0 ) + C Nn=0 t `2J Qn`+ 1 h` jv`n+1 ? v`n j 2 The best values for the ; constants are implicitly given by: n j;
n j;
q
nP
o1
P
= TVM(u0 ) Nn=0 t `2J (4 + :Qn`+ 1 ) h` jv`n+1 ? v`n j 2 nP o1 q P N 1 n = 2MTV (u0 ) n=0 t `2J :Q`+ 1 h` jv`n+1 ? v`n j 2 n j;
n j;
2
2
According to the CFL condition (1.12) which implies the TVD property, an upper bound for is given by: s
"
N X
# 21
1 jv0 ? v0 j 3 MT h 12 TVMh t 4 + (u0 ) n=0 `2Z 2 `+1 ` 2 So, for this choice of ; ; the terms E1 and E2 are bounded by X
i1 2 9 n n n n=0 t `2J 2 Q`+ 12 h` jv`+1 ? v` j ;
hP N
p
P
4C TV (u0 ) and i1 hP p P 2C TV (u0 ) Nn=0 t `2J Qn`+ 1 h` jv`n+1 ? v`n j 2 2 n j;
n j;
1
2 respectively (where in Jj;n ; is replaced by its upper bound 3 MT 2 h ) It remains only to add the rst-order time term and thus to complete the proof.
3 Error estimates for Second Order Resolution (SOR) schemes In this section, we will derive estimates for Second Order Resolution (SOR) schemes of the form (1.7). Given any \essentially 3-point" conservative discretization of the type (1.7), one is always able to de ne a viscous form (1.8) considering a modi ed ux fjn and a modi ed numerical viscosity Qnj+ 1 = mnj+ 1 + pnj+ 1 for which the following analogue of the rst-order formula holds 2
2
2
Qnj+ 21
fjn + fjn+1 ? 2gjn+ 12 = vjn+1 ? vjn 9
We will use also a discrete one-sided Lipschitz semi-norm of the numerical approximations as another indicator in the error estimate:
kvh (:; tn)klip+(R) = sup j 2Z
vjn+1 ? vjn hj +
(3.1)
A numerical scheme is said to be lip+ -stable if the quantity (3.1) remains bounded as h ! 0 for all tn > 0 (see e.g. [23, 24, 31, 33] for details). A weaker in-cell entropy inequality. In this paragraph we derive a weak in-cell entropy inequality for general numerical schemes written in viscous form (1.8). This will be useful to use once again theorem 2. Lemma 1. Assume that the scheme (1.8) satis es the stability conditions of theorem 1, then, for all entropies U k (v) = jv ? kj; k 2 R and all (j; n) 2 Z N , we have the following in-cell inequality: U k (vjn+1 ) ? U k (vjn ) ? 2ht [mnj+ 21 (U k (vjn+1 ) ? U k (vjn)) ? pnj? 21 (U k (vjn ) ? U k (vjn?1 ))] 0 (3.2) j Proof. Multiplying the viscous form (1.8) by sgn(vjn+1 ? k), one gets for every k 2 R:
jvjn+1 ? kj ? sgn(vjn+1 ? k)(vjn ? k) ? 2ht mj+ 21 sgn(vjn+1 ? k)(vjn+1 ? vjn) + 2ht pj? 21 sgn(vjn+1 ? k)(vjn ? vjn?1) = 0 j
j
Since sgn(vjn+1 ? k)(vjn ? k) jvjn ? kj, using the CFL condition of the theorem 3, we get: t n +1 n jvj ? kj ? jvj ? kj 1 ? 2h (mj+ 21 + pj? 21 ) ? 2ht mj+ 21 jvjn+1 ? kj ? 2ht pj? 12 jvjn?1 ? kj 0 j j j This yields the desired result. An approximate entropy inequality for vh: As in the previous section, we want now to bound the following quantity for all k 2 R and all positive test-functions ' in C01 (R R+ ). The main point will be to compensate the lack of consistency of (3.2). scheme.
I=?
Z
[j|vh{z? k}j 't + sgn( vh ? k)({z f (vh ) ? f (k))} 'x ](x; t)dxdt + | (RR ) U (v ) F (v ) k
h
k
h
Since the approximate solution vh is piecewise constant,
f +1 ?f n+1 k n+1 k n I R v +1 j;n 'K [U (vj ) ? Uh (vj )] + 'j + 21 vi sgn(v ? k) v +1 ?v R ?f +'Ij + 1 vv +1 sgn(v ? k) f 0 (v) ? vf +1 ?v :dv +1 2
P
I =
n
j
n
def
n j n j
n j n j
= I1 + I2 + I3
R
R
n j n j
where 'nK+1 = K '(x; tn+1 )dx, 'Ij + 1 = I '(xj + 12 ; t)dt . 2 j
n
j
n
We have now the following estimate 10
n j n j
n j n j
n j n j
:dv
(3.3)
Lemma 2. I1 + I2
X1h
i
mnj+ 21 jvjn+1 ? vjnj + pnj? 21 jvjn ? vjn?1j 2 j;n
Z
K I j
n
j't (x; t)j + j'x (x; t)jdxdt
!
(3.4)
Proof. We have:
I1 + I2 =
n+1 k n+1 k n j;n"'K [U (vj ) ? U (vj )] # Z v Z v n ? fn n ? fn +1 f f 1 j +1 j j j ? 1 + 2 'Ij + 21 sgn(v ? k) vn ? vn :dv + 'Ij ? 21 sgn(v ? k) vn ? vn :dv v v ?1 j +1 j j j ?1
P
j
n j
n
n j
n
n j
n j
R
R
Using the weak entropy inequality (3.2) multiplied by 'nj = I K '(x; t)dxdt 0 we nally obtain n+1 n k n+1 k n I1 + I2 PP j;n('K ? 'j =t)[U (vj ) ? U (vj )] + j;n 'nj =(2hj )[Qnj+ 1 (U k (vjn+1 ) ? U k (vjn )) ? Qnj? 1 (U k (vjn ) ? U k (vjn?1 ))] 2 R 2 P ? + 12 j;n('Ij + 1 ? 'nj =hj ) vv +1 sgn(v ? k) v +1 :dv +1 ?v 2 R P v + 21 j;n('Ij ? 1 ? 'nj =hj ) v ?1 sgn(v ? k) v ??v ??11 :dv n
j
j
n j n j n j n j
n
n
2
fjn
n j fjn n j
fjn
fjn
n j
n j
To estimate the space terms of the right-hand-side of this inequality, rst we observe that: +1 ? +1 ? k n k n v sgn(v ? k) v +1 ?v :dv = v +1 ?v (U (vj +1 ) ? U (vj )). So, the quantity
R vjn+1
fjn
n j
n j
fjn
fjn
n j
n j
fjn
n j
R = 12 Pj;n 'nj=hj [Qnj+ 21 (U k (vjn+1) ? U k (vjn)) ? Qnj? 21 (U k (vjn ) ? U k (vjn?1))] P ? [U k (vjn+1 ) ? U k (vjn )] + 12 j;n ('Ij + 1 ? 'nj =hj ) v +1 +1 ?v 2 P + 21 j;n ('Ij ? 1 ? 'nj =hj ) v ??v ??11 [U k (vjn ) ? U k (vjn?1 )] 2 fjn
n
n j fjn n j
n
fjn
fjn
n j
n j
becomes after a summation by parts X R = 12 Qnj? 21 [U k (vjn) ? U k (vjn?1)][('Ij? 12 ? 'nj=hj ) ? ('Ij? 12 ? 'nj?1=hj?1 )] j;n X fn ? fn 1 + 2 vjn ? vjn?1 [U k (vjn ) ? U k (vjn?1 )][('Ij ? 12 ? 'nj =hj ) + ('Ij ? 12 ? 'nj?1 =hj ?1 )] : j ?1 j;n j n
n
n
n
Therefore,
R = 12
X
j;n
h
i
[U k (vjn ) ? U k (vjn?1 )] pnj? 21 ('Ij ? 12 ? 'nj =hj ) ? mnj? 21 ('Ij ? 21 ? 'nj?1 =hj ?1 ) : n
n
And we have: I1 + I2 Pj;n('nK+1 ? 'nj =t)[U k (vjn+1h) ? U k (vjn )] i P + 21 j;n[U k (vjn ) ? U k (vjn?1 )] pnj? 1 ('Ij ? 1 ? 'nj =hj ) ? mnj? 1 ('Ij ? 1 ? 'nj?1 =hj ?1 ) : j
n
2
11
n
2
2
2
The proof of the lemma is complete by observing that X ('nK+1 j
j;n
? 'nj =t)[U k (vjn+1 ) ? U k (vjn)] iZ Xh n n n n n n 2 mj+ 21 jvj+1 ? vj j + pj? 21 jvj ? vj?1j
K I
j;n
and
1 X[U k (vn ) j 2 j;n 1 Xh
h
j
? U k (vjn?1)] pnj? 21 ('Ij? 12 ? 'nj=hj ) ? mnj? 21 ('Ij? 12 ?
2
j;n
n
n
iZ
mnj+ 21 jvjn+1 ? vjnj + pnj? 21 jvjn ? vjn?1j
K I j
n
n
j't (x; t)jdxdt
i 'nj?1 =hj?1 )
j'x (x; t)jdxdt :
The proofs of the next two lemmas are forwarded to the appendix. Lemma 3. For every (a; b) 2 R2 and 2 [0; 1], let us de ne the following quantity:
A(a; b; ) def = 2 (1 ? )f (a) + f (b) ? f (1 ? )a + b
?
Then, we have:
I3
X
Z
+ jn+ 21 ) sup {z } K I x0 2K q + 12
(!jn+ 21
j;n |
j
n
j
'(x0 ; t) dxdt hj
(3.5)
n j
where:
!jn+ 21
= ( 1v v +1 ? n j
n j
jn+ 21 = jf~jn+1 ? f~jnj
1v v +1 ) max A(vn ; vn ; ) and 2[0;1] j j +1 + n j
n j
Remark 1. Since the ux function f is smooth, we have !jn+ 1 = O(jvjn+1 ? vjn j2 ). In 2 addition, the di erence of two adjacent correction terms is known [28] to be such that jn+ 1 = 2 O(jvjn+1 ? vjnj2 ). Consequently, the contribution of I3 will be actually very small in the areas of smoothness of the entropy solution of (1.2). Moreover, one notices that Oleinik's entropy condition [25] implies that if the values vjn ; vjn+1 are connected by a single shock, then !jn+ 1 = 0. 2
Assuming that f is convex, one can re ne the expression of !jn+ 1 . 2
Lemma 4. In the case where f is convex there holds f 00 ()jvjn+1 ? vjn j:(vjn+1 ? vjn)+ 8j 2 Z; !jn+ 21 2[max v ;v ] n j
+1
n j
12
(3.6)
Summarizing, under the assumptions of Theorem 1, vh satis es the following approximate entropy inequality
?
Z
(RR+ )
[jvh ?kj't + sgn(vh ? k)(f (vh ) ? f (k))'x ](x; t)dxdt
+ +
X
j;n
X
j;n
qnj+ 21 h1
Z
'(x0 ; t)dxdt sup 0
j K I x 2K j
n
j
mn jvn ? vnj + pn jvn ? vn j Z j't (x; t)jdxdt j ?1 j ? 21 j 2 j + 21 j +1 j K I
X 1
mnj+ 21 jvjn+1 ? vjn j + pnj? 21 jvjn ? vjn?1j 2 j;n
Z
j
K I j
(3.7)
n
n
j'x (x; t)jdxdt:
Now, we can use the framework provided by theorem 2 starting from the approximate inequality (3.7). Theorem 4. Assume that the restrictions of theorem 1 are met for the scheme (1.8) which satis es (3.2) for all Lipschitz continuous entropies U k , k 2 R. Then we have for all T > 0: x + 1 +MT + 2 jvh(x; 0) ? u(x; 0)jdx K x ? 1 ?MT ? 2 N N o1=2 X n X p p + t XEn + C sup fXVn g tXVn + 2C ( 10 + 1) TV (u0 ) n=0;:::;N +1 n=0 n=0
Z
jvh (x; T ) ? u(x; T )jdx
j
where
Z
j
j
XVn = P`2J h` Qn`+ 21 jv`n+1 ? v`nj XEn = P`2J qn`+ 21 with qn`+ 21 = !`n+ 21 + `n+ 12 q n n Jj; = f` 1 : x` 2 B (xj ; hj =2 + M (T ? t ) + )g; = 54Th ; 0 M = supfjf ()j with jj ku0 kL1 (R)g; C = 2; N = T=t:
(3.8)
n j;
n j;
(3.9)
Remark 2. Note that the error bound in (3.8) is at best O(h1=2 ): This is also observed computationally, cf. Section 4. In any case there are not known asymptotic estimates of higher order for high resolution schemes. Remark 3. It seems that the new terms XEn which measure the lack of the entropy consistency of the scheme are of higher order if compared to the other terms in the estimate. Indeed, rst note that in the case where fjn = f (vjn ) and the ux is convex then this term is O(h) while the other terms in the estimate are O(h1=2 ): This is a consequence of the lip+ stability of these schemes proved by Tadmor. In the general case this is supported by heuristic arguments based on Remark 1 and on the fact that the correction term f~jn should be \close" to h(Qnj+ 1 ? (anj+ 1 )2 )=2, cf. [9]. We will see the existence of numerical evidence of the nice 2 2 behavior of this indicator in the tests displayed in the forthcoming section.
13
Proof. Starting from (3.7), we de ne the following piecewise constant functions for (x; t) 2 Kj I n:
(x; t) = qnj+ 21 =hj (x; t) = 12 mnj+ 1 jvjn+1 ? vjn j + pnj? 1 jvjn ? vjn?1 j 2 2
(x; t) = 2 mnj+ 1 jvjn+1 ? vjnj + pnj? 1 jvjn ? vjn?1 j 2
(3.10)
2
Then, Theorem 2 implies, R R h jxjR+MT + ju(ns; 0) ? vh(s; 0)j:ds jxjR ju(s; T ) ? v (s; T )j:ds o ? R R +2C ( + M)TV (u0 ) + 0T jxjR+M (T ?t)+ (s; t) + 4C (s;t) + C 2M + 1 (s; t) ds:dt R +C sup0tT + jxjR+M (T ?t)+ (s; t)ds
with B (xj ; hj =2+ M (T ? t)+) = [xj ? 21 ? M (T ? t) ? ; xj + 12 + M (T ? t)+] for all t 2 [0; T ]. By the de nition of Jj;n we have Z
B(x ;h =2+M (T ?t)+) j
(s; t)ds:dt
j
N X n=0
t
X
`2J n j;
qn`+ 21 :
(3.11)
Similar bounds hold for the other terms. These terms are bounded by E1 + E2 where h
i
E1 = C 2TV (u0 ) + 1 Nn=0 t `2J (4 + 2M)Qn`+ 1 h` jv`n+1 ? v`n j 2 i h P P E2 = C 2MTV (u0 ) + 1 Nn=0 t `2J Qn`+ 1 h` jv`n+1 ? v`nj P
P
n j;
n j;
2
Using in this expression the CFL condition 2M 1 yields the optimal value of : N hX i1=2 X t 5Qn`+ 21 h` jv`n+1 ? v`n j : = p 1 2TV (u0 ) n=0 `2J n j;
The TVD property and supj;n Qnj+ 1 12 , imply that 2 following value for
=p
N hX
1
2M TV (u0 ) n=0
t
X
`2J n j;
q
5 4 Th. We choose also the i1=2
Qn`+ 21 h` jv`n+1 ? v`nj
? Then, replacing in Jj;n by its upper bound 45 Th 1=2 we conclude that E1 and E2 are bounded respectively by p
hP
i1=2
P
2C TV (u0 ) Nn=0 t `2J 10Qn`+ 1 h` jv`n+1 ? v`n j 2 and i1=2 hP p P 2C TV (u0 ) Nn=0 t `2J Qn`+ 1 h` jv`n+1 ? v`n j 2 and thus the proof is complete. n j;
n j;
14
4 Numerical results In this section, we display some numerical experiments on the inviscid Burgers equation: (
ut + u22 x = 0 u0 2 BV (R)
(4.1)
We will restrict ourselves to suciently simple initial data in order to be able to derive the exact solution and to compare the two estimates (2.4), (3.8) with the real error of the numerical schemes.
4.1 Two Riemann problems for the modi ed Lax-Friedrichs scheme Here, we consider the modi ed Lax-Friedrichs scheme introduced by Tadmor [32]. It corresponds to the upper bound for the numerical viscosity coecient in theorem 1: Qnj+ 1 21 . We only 2 used regular grids for which hj x for all j 2 Z. That means in particular that we have simply = t=x. We performed a numerical test on two Riemann problems for (4.1):
u0 (x)
=
1 for x < 1 ; u (x) = 0 for x < 1 0 1 for x > 1 : 0 for x > 1
(4.2)
Since the numerical scheme is proved to be endowed with a strong in-cell entropy inequality, we are in position to use both of the estimates (2.4), (3.8). We compare also these values with the real global L1 error. The results for T = 0:4 are displayed in gure 1 for the shock and for the expansive wave. We used = 21 and sixteen values for x between x = 0:05 and x = 0:000276. The computational domain is the interval [0; 2]. 10
10 "errors3.num" "errors2.num" "errors1.num"
"errors3.num" "errors2.num" "errors1.num"
1 1
0.1 0.1 0.01
0.01 0.001
0.0001 0.0001
0.001
0.01
0.1
0.001 0.0001
0.001
0.01
0.1
Figure 1: Comparison between the error estimates (2.4), (3.8) (doted lines) and the real error. Left: shock wave - Right: rarefaction wave. The solid line corresponds to the real error of the numerical scheme on the whole computational domain. It shows a nearly rst-order convergence in both cases regardless to the smoothness of the exact solution. The two dotted lines display the estimates given by (2.4), (3.8) on the same interval.
15
4.2 An experiment with its second order (SOR) extension
We present now some results obtained by a second order extension of the former Lax-Friedrichs scheme. This SOR algorithm will be derived by means of the piecewise constant viscosity modi cation introduced by Osher and Tadmor. To carry this out, we have to choose convenient values for f~jn and Qnj+ 1 in order to use them in (1.8). One possible choice suggested in [28] is 2 given by: n o f~jn = sgn(v +1 ?v )+4 sgn(v ?v ?1 ) min (Qn 1 ? (an 1 )2 )jvjn+1 ? vjn j; (Qn 1 ? (an 1 )2 )jvjn ? vjn?1 j ; n j
n j
n j
n 1 = Qn 1 + f~jnn+1 ?f~jnn j + 2 vj+1 ?vj j+ 2
Q
n j
j+ 2
with Qnj+ 1 = 2
j+ 2 1 n 2 and aj + 21
=
f (v +1 )?f (v ) v +1 ?v n j n j
n j
n j
j? 2
j? 2
h These values satisfy the second order requirement (1.10). Moreover, the numerical solution v f (v +1 )?f (v ) will be endowed with the TVD property under the following CFL condition: v +1 ?v 32 Consequently, we will consider the following initial data for (4.1): u0 (x) = sin(x) for x 2 [0; 2] (4.3) n j n j
n j
n j
Its exact solution develops a standing shock at x = 1 and can be computed accurately using the method of characteristics. This provides a way to compare the real global error of the numerical scheme with the estimate (3.8). We used once again a uniform grid and also the same range of parameters as in the preceding paragraph. The results a time T = 0:4 are displayed in the gure 2 where the solid line represents the exact errors and the dotted one the estimated ones. For the sake of completeness, we present also some results before the formation of the stationary shock at time T = 0:15. The second order accuracy for this smooth solution is clearly noticeable. On the contrary, the estimated errors for the SOR modi cation are extremely close to the ones obtained for the former rst order modi ed Lax-Friedrichs scheme, cf. Section 3, Remark 2. In addition the contribution of the XEn terms in the estimate is indeed very small in this case. 100
100 "errors3.num" "errors2.num"
"errors3.num" "errors2.num" 10
10
1 1 0.1 0.1 0.01 0.01 0.001
0.001
0.0001 0.0001
0.0001
0.001
0.01
0.1
1e-05 0.0001
0.001
0.01
0.1
Figure 2: Comparison between the error estimate (3.8) and the real error. Left: standing shock at T = 0:4 - Right: smooth solution at T = 0:15.
4.3 An example of adaptive algorithm built on these estimates
As a last illustration we present preliminary computational experiments on the inviscid Burgers equation (4.1) with an SOR adaptive numerical scheme using as a building block the former 16
1 "solLXF.num" "solEX.num" 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 3: Numerical and exact solutions at time T = 0:8. 1 "grille.num"
0.1
0.01
0.001 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 4: Computational grid at time T = 0:8. extension of the Lax Friedrichs scheme. The algorithm consists in solving this P.D.E. on a variable grid whose cells are denoted (Kj )j 2Z. At each time tn , the local estimate (3.8) is used to compute an upper bound of the error eK (tn ) in each mesh cell. Given a certain tolerance (tn ) > 0, we adjust the computational grid in such a way that the following requirement is ful lled: n (4.4) 8j; (20t ) eK (tn) (tn) This is achieved by re ning and de-re ning the grid for each time step. The test problem we considered consists in solving the equation (4.1) associated with the following Lipschitz continuous initial data: 1 for x < 0 u0 (x) = (1 ? x) for 0 x 2 ?1 for x > 2 j
j
17
Its exact solution develops a stationary shock at x = 1 for t 1. Everywhere else, the solution is constant: it is therefore relevant to try to optimize the computational grid on such a problem since it is merely in the transition zone that a higher number of points are needed to describe accurately the solution. Consequentely, we display in the gures 3 and 4, the numerical results and the space discretization obtained at time T = 0:8 before the formation of the shock. The CFL number is xed 0:45: Then, the grid is adjusted to ful ll the requirement (4.4) for each time step 0 < tn T . 1 "solLXF.num" "solEX.num" 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 5: Numerical and exact solutions at time T = 1:2. 1 "grille.num"
0.1
0.01
0.001
0.0001 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 6: Computational grid at time T = 1:2. We turn now to the time T = 1:2 for which the stationary shock is clearly noticeable on the gure 5. The computational grid ts quite well with the structure of this steady-state solution, see the gure 6. This run was carried out using the same parameters as the preceding one.
18
A Appendix: Estimates on I3 Proof of lemma 3 We want to estimate for k 2 R: ! ! Z v Z v X +1 +1 f (vjn+1 ) ? f (vjn) f~jn+1 ? f~jn I 0 I3 = 'j+ 21 sgn(v ? k) f (v) ? vn ? vn :dv ? sgn(v ? k) vn ? vn :dv v v j +1 j +1 j j j;n n j
n
n j
n j
n j
R If k 62]a; b[, the rst integral is zero. So, the desired result follows by observing that vv +1 sgn(v ? ~ ?f~ ~n dv jfj+1 ? f~jn j. On the contrary, if k lies in ]a; b[, we have: k) fv +1 +1 ?v (i) b < k < a We can write k = b + (1 ? )a for 2]0; 1[: n j n j
n j n j
n j n j
f (b)?f (a) dv = 2f (k) ? f (a) ? f (b) ? (2k ? (a + b)) f (b)?f (a) 0 b?a a sgn(v ? k) f (v) ? b?a = 2f (k) ? f (a) ? f (b) + (1 ? 2)(f (b) ? f (a))
Rb
(ii) a < k < b In this case, we can write k = a + (1 ? )b for 2]0; 1[:
f (b)?f (a) dv = f (a) + f (b) ? 2f (k) ? ((a + b) ? 2k) f (b)?f (a) 0 a sgn(v ? k) f (v) ? b?a b?a = ?2f (k) + f (a) + f (b) + (1 ? 2)(f (b) ? f (a))
Rb
Proof of lemma 4 The goal is to get an estimate on the following term, provided that f is convex and k 2 R: ! Z v +1 f (vjn+1 ) ? f (vjn) 0 n sgn(v ? k) f (v) ? vn ? vn dv !j+ 21 = v j +1 j n j
n j
First of all, it is useful to nd the sign of ab sgn(v ? k) f 0 (v) ? f (bb)??fa(a) dv for all (a; b; k) 2 R3 . If k 62]a; b[, this integral is zero. So, we look at the two other cases: (i) b < k < a We can write k = b + (1 ? )a for 2]0; 1[: R
f (b)?f (a) dv = 2f (k) ? f (a) ? f (b) ? (2k ? (a + b)) f (b)?f (a) 0 b?a a sgn(v ? k) f (v) ? b?a = 2f (k) ? f (a) ? f (b) + (1 ? 2)(f (b) ? f (a)) ?(1 ? 2)f (b) + (2 ? 2)f (a) ? f (a) + (1 ? 2)(f (b) ? f (a)) 0 We used Jensen's inequality: f (b + (1 ? )a) f (b) + (1 ? )f (a). (ii) a < k < b In this case, we can write k = a + (1 ? )b for 2]0; 1[: Rb
f (b)?f (a) dv = f (a) + f (b) ? 2f (k) ? ((a + b) ? 2k) f (b)?f (a) 0 b?a a sgn(v ? k) f (v) ? b?a = ?2f (k) + f (a) + f (b) + (1 ? 2)(f (b) ? f (a)) ?(1 ? 2)[f (b) ? f (a)] + (1 ? 2)(f (b) ? f (a)) 0
Rb
19
So, we already have that: Z
b a
sgn(v?k) f 0 (v) ? f (b) ? f (a)
b?a
f (a) dv 1a
