Discrete Kinetic Schemes for Multidimensional Systems of Conservation Laws.
We present here some numerical schemes for general multidimensional systems of conservation laws based on a class of  discrete kinetic approximations, which includes the relaxation schemes by S. Jin and Z. Xin. These schemes have a simple formulation even in the multidimensional case and do not need the solution of the local Riemann problems. For these approximations we give a suitable multidimensional generalization of the Whitham's stability subcharacteristic condition. In the scalar multidimensional  case we establish the rigorous  convergence of the approximated solutions to the unique entropy solution of the equilibrium Cauchy problem.
Introduction
. In this paper we present a new class of numerical schemes
based on a discrete kinetic approximation for multidimensional hyperbolic systems of
conservation laws. Consider a weak solution u
K to the Cauchy
problem
where the system is hyperbolic (symmetrizable) and the flux functions A d are locally
Lipschitz continuous on R K with values in R K . We approximate problem (1.1), (1.2)
by a sequence of semilinear systems
with Cauchy data
Here # is a positive number, # d are real diagonal L - L matrices, P is a real constant
coe#cients K - L matrix, and M is a Lipschitz continuous function defined on R K
# Received by the editors August 5, 1998; accepted for publication (in revised form) November
15, 1999; published electronically May 23, 2000. This work was partially supported by TMR project
http://www.siam.org/journals/sinum/37-6/34307.html
Math-ematiques Appliqu-ees de Bordeaux, Universit-e de Bordeaux 1, 351 cours de la Lib-eration,
Talence, France (aregba@math.u-bordeaux.fr).
# Istituto per le Applicazioni del Calcolo "M. Picone," Consiglio Nazionale delle Ricerche, Viale
del Policlinico 137, I-00161 Rome, Italy (natalini@iac.rm.cnr.it).
1974 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
with values in R L . Moreover we suppose the following relations are satisfied for all
fixed
rectangle# in R K :
It is easy to see that, if f # converges in some strong topology to a limit f and if
0 converges to u 0 , then Pf is a solution of problem (1.1), (1.2). In fact system
(1.3) is just a BGK approximation for (1.1); see [5, 12] and references therein. The
interaction term on the right-hand side is given by the di#erence between a nonlinear
function, which describes the equilibria of the system, in our case M(Pf ), and the
unknown f . Our purpose here is to construct numerical schemes for system (1.3) in
order to obtain a numerical approximation of (1.1) in the relaxed limit
As well known for general relaxation problems, see [71, 44, 55], approximation
needs a suitable stability condition to produce the correct limits. In the frame-work
of general 2 - 2 quasi-linear hyperbolic relaxation problems, this condition is
known as the subcharacteristic condition. In section 2, we shall argue in the spirit
of the Chapman-Enskog analysis [71, 44, 14], to find the following stability condition
for
for all # 1 , . , # D ) # (R K ) D and every u belonging to some fixed
rectangle# R K .
Actually, in [54] and for the scalar case convergence of Pf # towards the
Kruzkov entropy solution of (1.1), (1.2) has been obtained under a slightly stronger
version of condition (1.6): every component of the Maxwellian function M is monotone
nondecreasing on the interval I. The main tool in that case is the fact that under
this condition the right-hand side in system (1.3) is quasi-monotone in the sense of
[25] and this implies special comparison and stability properties on the corresponding
system. Unfortunately, similar properties are not verified for nontrivial examples in
the general case K > 1. Therefore, for the systems, we shall use just condition (1.6).
Note that for certain families of kinetic approximations, the results of [6] show that
(1.6) is also a necessary condition for (1.3) to be compatible with the entropies of
(see section 3).
The continuous kinetic approximation of systems of conservation laws in gas dynamics
is classical. In particular, Euler equations can be formally obtained as the
fluid dynamical limit of the Boltzmann equation; see [12, 13]. The rigorous theory of
kinetic approximations for solutions with shocks is recent and the main results were
obtained only in the scalar case. The first result of convergence of a fractional step
BGK approximation with continuous velocities, with an entropy condition for the
limit (weak) solution, was proven in [7] (see also [22]). Another convergence result
was given later in [60], using a continuous velocities BGK model. An important related
kinetic formulation can be found in [41]. Other results have been established for
special systems or partially kinetic approximations [42, 29, 9, 40]. Related numerical
schemes can be found in [19, 59]; for a general overview and many other references
see [24]. Discrete velocities models and their fluid dynamical limits have also been
considered by many people; see the review paper [61]. In particular we mention the
studies on the Broadwell model [11, 72]. Convergence for various relaxation models
was also investigated in [15, 14, 17, 47, 68, 74, 32, 33, 67]. The analysis of the stability
of various nonlinear waves for relaxation models, and in particular for the Jin-Xin
relaxation approximation, can be found in [44, 16, 50, 43, 46, 49]. A general survey of
recent results on relaxation hyperbolic problems is given in [55]. Let us also point out
some numerical references related to our approach. A lot of computational work has
been done in the last ten years in the very closed framework of lattice Boltzmann and
BGK models; see [21, 62] and references therein. Let us also mention the monotone
schemes of [8], which are an example of numerical (relaxed, i.e., discretization
of our construction in the scalar case. Other numerical investigations for
hyperbolic problems with relaxation can be found in [37, 57, 58, 73, 4, 3, 18, 70, 31].
Our numerical schemes are constructed by splitting (1.3) into a homogeneous linear
part and an ordinary di#erential system, which is exactly solved thanks to the
particular structure of the source term. In the scalar case this construction allows us
to preserve the monotonicity properties of (1.3) and to prove convergence results. Our
approximation framework generalizes to systems the construction presented in [54] for
the scalar case, and shares most of the advantages of the relaxation approximation as
proposed in [30] (see also [53, 2, 70]): simple formulation even for general multidimensional
systems of conservation laws and easy numerical implementation, hyperbolicity,
regular approximating solutions. Actually the main advantage, especially in the multidimensional
case, of both the approximations, seems to be the possibility of avoiding
the resolution of local Riemann problems in the design of numerical schemes. Moreover
our framework presents some special properties:
- the scalar and the system cases are treated in the same way at the numerical
- all the approximating problems are in diagonal form, which is very likely for
numerical and theoretical purposes;
- we can easily change the number and the geometry of the velocities involved
in our construction to improve the accuracy of the method.
In this sense our work shares most of the spirit of [56, 34, 45], where very flexible
and simple schemes, which do not need Riemann solvers in their construction, were
proposed to approximate general multidimensional systems of conservation laws. Let
us also observe that the presented algorithms are surely not optimal, but they just
illustrate how to construct an e#cient and simple approximation even for very complicated
systems. This could be useful, for example, in the numerical investigation of
large systems like those arising in the extended thermodynamics and other generalized
moment closures hierarchies for kinetic theories [52, 38, 1]. Further investigations will
be addressed to the construction of high order schemes.
The plan of the article is as follows: In section 2 we establish the stability condition
and define the monotone Maxwellian functions for the scalar case. In section
3 we propose some examples of stable approximations in the class (1.3). The issue of
entropy is discussed. In section 4 we set the numerical schemes and section 5 is devoted
to the convergence results in the scalar multidimensional case. Some numerical
experiments are given in section 6.
After the completion of this work we received a preprint from Serre [64], where he
proves, by using the methods of compensated compactness, the convergence for the
Jin-Xin relaxation approximation and some of the discrete kinetic approximations
contained in the present paper to (one-dimensional) genuinely nonlinear hyperbolic
systems of conservation laws having a positively invariant domain. The convergence
of related first-order numerical schemes has been proved by Lattanzio and Serre in
[35]. The stability conditions are in both cases strongly related to our condition (1.6).
1976 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
2. Chapman-Enskog analysis and monotone Maxwellian functions. In
this section we discuss the stability conditions for the discrete kinetic approximation
(1.3). Since the local equilibrium for that system is given by the hyperbolic system
(1.1), it is natural to seek for a dissipative first-order approximation to (1.3), which
is the analogue of the compressible Navier-Stokes equations in the classical kinetic
theory. In principle we could try to use the theory developed in a more general context
in [14]. Unfortunately it is easy to realize that their main assumption, namely the
existence of a strictly convex dissipative entropy for the relaxing system (1.3), which
verifies in particular the requirement (iii) of Definition 2.1 of [14], is not satisfied in
the present case and we need it for a di#erent construction.
Let f # be a sequence of solutions to (1.3)-(1.4) parametrized by #, for a fixed initial
data f 0 , which for simplicity we can choose as a local equilibrium, i.e., f 0
for some u 0 # L # (R D , R K ). Set
Then, from (1.3) and the compatibility assumptions (1.5), we have
Consider a formal expansion of f # in the form
Then
Reporting in (2.1) yields
Now we have
Then, up to the higher order terms in (2.4), we obtain
where
is a K-K matrix. We can state now our stability condition.
Proposition 2.1. The first-order approximation to system (1.3) takes the form
and it is dissipative provided that the following condition is verified:
for all # 1 # R K , . , # D # R K and every u belonging to some fixed
rectangle# R K .
As we shall see in the next section parabolicity of (2.6) is, at least in some cases,
necessary for the compatibility of (1.3) with the entropies of (1.1). But let us recall
that, even in the scalar case, the expansion (2.6) cannot be considered in any way
as a rigorous asymptotic description of system (1.3). Actually, to prove our rigorous
convergence results, we need a slightly stronger version of condition (2.8).
Definition 2.2. Take I # R be a fixed interval. A Lipschitz
continuous function R L is a monotone Maxwellian function (MMF)
for (1.1) and with respect to the interval I if conditions (1.5) are verified and if,
moreover,
M i is a monotone (nondecreasing) function on I, for every i # {1, . , L}.
This condition was used in [54] to show the convergence of approximation (1.3),
(1.4) at the continuous level and in the multidimensional scalar case. In the following
section we present some examples of di#erent approximations according to the choices
of the matrices of velocities # j and the local Maxwellian function M . We discuss the
issue of entropy and we investigate both stability and (only for K=1) monotonicity
conditions.
3. Examples of discrete kinetic approximations. In order to construct the
system (1.3) one must find P , M and # such that the consistency relations (1.5) are
satisfied. The first three examples presented here own a block structure. Keeping
notations of the introduction we take blocks
I K , the identity matrix in R K . Each matrix # d is constituted of N diagonal blocks of
size K- K:
I K , # (d)
With this formalism, (1.3) can be written as
and f #
Considering that A # setting
1978 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
we have the following equivalent expression for (2.8):
In the case when the M #
are symmetric we obtain
System (3.1) enters the framework proposed by Bouchut [6]. Suppose that there
exists at least one smooth strictly convex entropy for (1.1), and that the Maxwellian
functions are of form
nd A d (u),
where a n , b nd are scalar. Denote #(M #
n (u)) the set of eigenvalues of M #
under some technical assumptions, it is shown in [6] that the M #
n (u) are diagonalizable.
Moreover if
then (3.1) is compatible with any convex entropy # of (1.1): there exists a kinetic
entropy for (3.1) associated with # and Lax entropy inequalities are satisfied in the
hydrodynamic limit # 0. As well known [36, 10] these inequalities characterize the
admissible weak solutions of (1.1). Moreover, in this case, (2.6) is parabolic. More
general results for Maxwellian functions not in the form (3.4) can be found in [51].
Remark that for a general hyperbolic system of conservation laws neither (3.5) nor
imply convergence of the kinetic approximation. In the scalar case it is always
possible to write (1.3) under the form (3.1), and (3.5) is the monotonicity condition
(2.9). Under this condition convergence holds [54]. Moreover one can use the fact
that # N
and the discrete Jensen inequality to prove directly the following
proposition [54].
Proposition 3.1. Let suppose that M is an MMF. Then (2.8) is
In the general case, denote B the K - D-K - D matrix defined by the blocks B dj ,
D. The stability condition means
positively defined
for all u #
where# is some fixed rectangle of R K .
Example 1. The diagonal relaxation method (DRM). We first consider the minimal
case 1. The system (1.5) is then a squared linear one. We take # > 0
and
Let us denote by e j the canonical unit vector of R D . The characteristic speeds of
system (1.3) are -#e 1 , . , -#eD , # D
. The Maxwellian function is given by
means of the following
A d (u) # /(D
Therefore the result of [6] applies to this model as described above.
For a one-dimensional system of conservation laws this formulation coincides with
the relaxation approximation of [30]. In fact we can set
Then we have
We recall that in this case and for convergence to the unique entropy solution
was proved in [53], and convergence of the associated numerical relaxed schemes was
done in [2]. For other convergence and error estimate results on this model, see also
[33, 67]. However, in several space dimensions, there is no diagonal form for the
formulation of [30], as already pointed out in [54].
For a one-dimensional system the stability condition (2.8) is here:
for all # R K ,
which, if A # (u) is symmetric with spectral ray #(u), gives
#(u).
This coincides with Bouchut's condition (3.5) so that here both (3.5) and (2.8) are
equivalent.
In two space dimensions, (2.8) becomes
# .
In fact in the scalar case we are able to prove that (2.8) and (2.9) coincide.
Proposition 3.2. Let us suppose that and that the choice of P ,
and M is as above. Then the stability condition (2.8) and the monotonicity condition
coincide and can be written as
1980 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
Proof. Here 3B is symmetric and its characteristic polynomial is of form
where A and C are defined by the following relations:
Both eigenvalues of B are positive if and only if A and C are positive.
A is positive if and only if #
8 .
Now we remark that
Six cases have to be studied for the values of A # 1 and A # 2 . As A # 1 and A # 2 play symmetric
roles in the formulas, we can suppose that A # 1 # A # 2 and study the three following
cases:
1.
2.
3.
Let us study the first case. C is positive if and only if
It remains to determine the position of 0: two cases are under consideration.
and one can see that
Consequently the stability condition is satisfied if and only if
which is easily seen to be the monotonicity condition (2.9).
so that we again recover the
same condition.
Cases 2 and 3 are similar and we omit the proof.
Example 2. Flux decomposition method (FDM). In this example, in view of a
more accurate approximation, we take a greater number of equations,
and, following an idea due to Brenier [8], we decompose the Jacobians of the fluxes
in positive and negative part. Denoting by B d the diagonal matrix of eigenvalues of
A # d and by Q d the associated matrix of the right eigenvectors we set
- .
Then we can define
with A d
for some # d > 0 and
In the scalar convex case and for an appropriate (first-order) discretization, this choice
corresponds, in the relaxation limit # 0, to the Engquist-Osher numerical scheme
[8]. For a one-dimensional system the stability condition (2.8) is here:
for all # R K , (|A #I K - |A # 0
so that in the one-dimensional scalar case it coincides with the monotonicity condition.
The two-dimensional case condition (2.8) reads
(R K
and
# .
For K=1 this matrix is positive if and only if
and this is exactly the monotonicity condition.
Proposition 3.3. Let us suppose that K=1 and that the choice of P ,
and M is as above. Then the stability condition (2.8) and the monotonicity condition
coincide and can be written as (3.14).
Concerning the entropy properties of this model we refer to [51].
Example 3. Orthogonal velocities method (OVM). This example works with any
number of blocks N # D+ 1. We take the velocities such that
and
1982 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
This means that we choose an orthogonal family of D vectors (# 1d , . , #Nd
1, . , D) in the orthogonal space of the vector . The corresponding
Maxwellian function is now given by
A d (u)
a 2
d
where a 2
id . Here again the Maxwellian functions are of type (3.4) and the
result of [6] applies as described in the preceding section.
The one-dimensional case. Let us examine the stability condition in one space
dimension for this approximation. We take # > 0 and
# I K # 1#n#N
diag
Then
We obtain the following expression for B:
I K -
For instance if is symmetric with spectral ray # 2 (u), the stability condition
reads
However Bouchut's condition (3.5) now reads
This approximation therefore gives an example where, even in the scalar case, the
monotonicity condition is strictly stronger than the stability one.
The two-dimensional case. In two space dimensions we take length and direction
varying velocities. Fix J, N # 1, # > 0. Set
diag # n diag # cos( i#
diag # n diag # sin( i#
Here we have J . It is easy to see that (3.15) and (3.16) are satisfied and that
a 2
For the Maxwellian functions are
Proposition 3.4. Let us suppose that and that the choice of P , # 1 , # 2 and
M is as above. We denote # the argument of
sin #. If
then the monotonicity condition (2.9) is satisfied. The stability condition (2.8) is
satisfied as soon as
Proof. The first part of the proposition is immediate. Let us examine the stability
condition (2.8). We remark first that in addition to (3.15), (3.16) we also have
sin 2 i#
cos 3 i#
sin 3 i#
Consequently
for all j, d # {1, 2}, P# d
jd
a 2
d
N .
Denoting
-A
-A #2 2
# .
This matrix is positive if and only if
which ends the proof.
Example 4. Suliciu's method. Let us conclude this presentation by giving an
example which is a bit di#erent from those presented above. For K=2, D=1, let us
consider a system of the form
Take
and
# .
1984 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
For A 2 we have a p-system, and in this case the present approximation was
first introduced by Suliciu [65, 66] to study instability problems in phase transitions
described by elastic or viscoelastic constitutive equations. In this case the Chapman-Enskog
analysis gives the stability condition
Some preliminary investigations on this model can be found in [26, 48]. More recent
results have been found in [28, 27, 39, 69]. In particular see [69] for an entropy
condition.
In fact, this model enters the preceding block structure (3.1) by equivalence with
the following one:
# .
The Maxwellian functions M are of type (3.4) so that the result of [6] applies. For
an easy calculation shows that both conditions (2.8) and (3.5) coincide
with (3.26).
4. Discrete kinetic schemes. In this section we construct numerical schemes
for the relaxing semilinear problem (1.3), (1.4) associated with (1.1), (1.2). Of course
a lot of numerical schemes are available for this problem including those presented in
[30]. We present here a finite volume scheme on structured mesh based on a splitting
method.
The space time domain R D
discretized by a rectangular grid:
I # , [0, T
d be the canonical d th vector in R D .
As usual we denote by x# the center of I #x #,d the length of I # in the direction
d, #t
#x #,d ,
#x #,d
#t n .
Finally we set
f #,n
f #,n
f #,n
If
then f #
0 is approximated by
f #,0
We use for u the same notation as for f here above. Let us recall that
and
System (1.3) is split into a linear diagonal hyperbolic part and an ordinary di#erential
system. For a given f #,n
# , the function f #,n+1/2
# is an exact or approximate solution at
time t n+1 of the problem
# .
As the system is diagonal, we may consider each equation separately. We suppose
that the scheme can be put in conservation form:
f #,n+1/2
#x #,d
#,n
ed - #,n
ed # ,
where
#,n
ed
#-k1+ed ,1 , . , f #,n
#-kL+ed ,L , . , f #,n
#,n
ed ,l # 1#l#L
Here k l # Z D and
# l,d (g, . ,
In the following, the scheme on the linear part will be referred to as homogeneous
scheme (HS) and the associated evolution operator will be denoted by H# :
f #,n+1/2
# .
To take into account the contribution of the singular perturbation term on the
right-hand side, we solve on [t n , t n+1 ] the ordinary di#erential system
with initial data
1986 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
for all # Z D . Using (1.5) we obtain
so that the solution of (4.9) with data F (t n at t n can be explicitly obtained as
Hence
f #,n+1
where u is defined by
Note that
Therefore we have constructed a wide family of numerical schemes for the semilinear
system (1.3), which di#er by the choice of the HS. In the scalar case (K=1), thanks to
the monotonicity properties of the interaction, we show in the following sections that
in fact the properties of each scheme are roughly speaking the same as those of HS and
the estimates are uniform in #. We shall often refer to these properties, which are now
classical; see [23]: preservation of extrema, monotonicity, total variation diminishing
(TVD), and L 1 -contraction. In particular recall that monotonicity implies all the
other properties and the TVD property implies the preservation of extrema for initial
data in L 1 (R D , R) # L # (R D , R) #BV(R D , R). In the following, the numerical scheme
given by (4.1), (4.2), (4.7), (4.13) will be referred to as discrete kinetic scheme (DKS).
When # 0 in DKS, we obtain the relaxed limit of the scheme:
# .
This last scheme can be written in conservation form:
#x #,d
ed -A n
ed # .
The numerical fluxes are defined by
A n
ed
ed
with
ed
#-k l +ed ), . , M l (u n
#+k l
In the following section, for the scalar case K=1, we specify in what sense this relaxed
scheme is the limit of DKS. Moreover all the estimates for DKS are uniform and pass
to the limit, so that we obtain strong convergence for the limit scheme.
5. Convergence of discrete kinetic schemes for multidimensional scalar
conservation laws. In this section we use monotonicity to prove rigorous convergence
results for DKS and the associated relaxed scheme. We consider a scalar conservation
law (K=1) and {u #
a family of initial data for (1.1) such that
(R D , R) # L # (R D , R) # BV(R D , R) and
ld be such that (1.5) is satisfied.
Throughout this section we suppose without loss of generality that
and for sake of simplicity we consider a uniform mesh
for all # Z D #x #,d = #x d .
Moreover, as shown in section 3, M and the # d can be found satisfying
is an MMF on [-# ].
The two following relations will be useful in that which follows:
for all
and
for all u,
Therefore we have for all # > 0
#f #,0
#f #,0
#f #,0
#,l # ,
TV(f #,0
TV(f #,0
#.
5.1. The supremum norm bound.
Proposition 5.1. Suppose (H1), (H2) are satisfied and that HS preserves the
extrema. Then the scheme DKS is L # stable and there holds
for all # > 0 for all n # 0, #f #,n
# .
Moreover, if HS is monotone, then the same is true for DKS.
Proof.
Using (1.5) and (4.12), we have
S(t,
1988 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
By (H2) we have also
for all # Z D , for all l # {1, . , L}, M l (-# f #,0
Suppose that
for all # Z D , for all l # {1, . , L}, M l (-# f #,n
Since HS preserves the extrema, we have
l (-# f #,n+1/2
for every # Z D and l # {1, . , L}. Then
for all # Z D , -# u #,n+1
# .
Using expression (4.12) we obtain
l (-# S l # t, t n , f #,n+1/2
for every t # [t n , t n+1 ], # Z D and l # {1, . , L}. Recalling (5.2) we obtain
(5.4). To conclude we remark that if M is an MMF on the interval [-# ],
then the RHS of (4.9) is a quasi-monotone application on
for all F # [M(-#
We can then apply the results of [25]: the flow given by (4.9) preserves the order.
5.2. BV and L 1 estimates. The following lemma follows easily from the monotonicity
properties of the interaction.
Lemma 5.2. Let F 0 and G 0 be two initial conditions for system (4.9) with corresponding
solutions F (t), G(t). If F 0 , G
for all t # 0
|F l (t) -G l (t)| #
|F 0,l -G 0,l | .
As a consequence we obtain another important result.
Proposition 5.3. Suppose that (H1), (H2) are satisfied for {u #
(1) If HS is TVD then DKS is TVD.
(2) If HS is L 1 contracting then DKS is L 1 contracting and there holds
#f #,n+1
In order to prove the equicontinuity property and the boundedness of the inter-action
we need to estimate the RHS of (1.3).
Lemma 5.4. Suppose that (H1), (H2) are satisfied. Suppose HS is TVD and
can be put in conservation form (4.7) with a Lipschitz continuous flux #. Then there
exists a constant C not depending on # such that
exp #t n
Proof. From (4.13) and the fact that u #,n+
# we have
exp #t n
Moreover, by (4.7),
ld
#x d
#,n
ed ,l - #,n
ed ,l # .
Hence for 1 # l # L
|M l (u #,n+ 1# ) -f #,n+ 1#,l | # |M l (u #,n
#,l | +C
#x d
sup l |# ld ||#,n
ed ,l -#,n
ed ,l | .
Using the fact that the flux is Lipschitz continuous and that for
we have
# |w #+ed - w# | ,
we obtain the desired inequality.
Proposition 5.5. Suppose that (H1), (H2) are satisfied. Suppose HS is TVD.
Then there exists a constant C not depending on # such that
for all t # [0, T
for all t, t # [0, T
Proof. Inequality (5.8) is immediate by (5.7), thanks to a recursive argument.
Let us now prove (5.9). As the numerical flux of HS is Lipschitz continuous, as in the
proof of Lemma 5.4, we obtain
#f #,n+ 1# - f #,n
and
#f #,n+1
As M(u #,0
# we have, by a recursive argument,
#f #,n+1
The above estimates allow us to prove the convergence of our numerical schemes.
First we have convergence towards the unique solution of (1.3), (1.4) when # is fixed.
Theorem 5.6. Let T > 0, # > 0, suppose that (H1), (H2) are satisfied for
Suppose that HS
1990 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
is TVD. For any T > 0, let #t
#x be constant. As #t # 0 the sequence f #
# converges in
to the unique solution f # of (1.3), (1.4), f #
sup
sup
As a corollary we obtain global existence of the solution of (1.3), (1.4) (see [54]
for a direct proof).
Proof of Theorem 5.6. We follow exactly the method of [20] and just give a sketch
of the proof: for all t # 0, {f #
is bounded in L 1
# L # . Moreover, by
Proposition 5.3, we can apply Fr-echet-Kolmogorov theorem and obtain a relatively
compact set of L 1
loc . Using the equicontinuity property (5.9) as in the proof of Ascoli-
Arzela's theorem we obtain convergence in L # (0, almost everywhere,
and the limit is a solution of the problem by Lax-Wendro#'s theorem. Estimates
(5.10), (5.11), (5.12) follow from (5.6), (5.8), (5.9).
5.3. The relaxed limit of the scheme. In this section we are interested in
the behavior of our numerical schemes as the parameter # tends to zero. The above
estimates imply boundedness and TVD properties for the relaxed discrete kinetic
schemes. We need one further assumption on the convergence of the initial data:
(H3) the sequence {u #
converges towards a function u 0 in L 1 .
As a consequence
and f #,0 converges towards f As above we set and by (1.5)
this notation is compatible with what happens at
Theorem 5.7. Let T > 0 and suppose that (H1), (H2), and (H3) are satisfied
and HS is TVD. As # 0, f #
# converges in L # ((0, T ); L 1
loc (R) L ) to a limit f # which
satisfies
#x #,d
ed -A n
ed # ,
where the numerical fluxes are defined by
A n
ed
ld # n
ed ,l
(5.
and
ed
#-k l +ed ), . , M l (u n
#+k l
The resulting numerical scheme is TVD and converges to a weak solution of (1.1),
(1.2). Moreover, if HS is monotone, the limit scheme is also monotone and converging
to the unique entropy solution of (1.1), (1.2).
Remark 5.8. (5.17)-(5.19) are the formulas (4.16)-(4.18) obtained formally at
the end of section 4.
Proof. L 1 stability and the TVD property imply that for n # {0, . , N}, {f #,n
0} is bounded in L 1
# BV, so there exists a sequence # k # 0 such that
in L 1
loc and f # k ,n
# for every #. Consequently
in L # (0,
loc ). Estimates (5.14), (5.15), (5.16), monotonicity, and TVD properties
are immediate consequences of Propositions 5.1-5.5. Moreover the resulting scheme
can be written as (5.16)-(5.19). Thus f # is unique and the whole sequence converges.
The consistency with the conservation law (1.1) is the consequence of the consistency
of HS with the linear part of (1.3).
Therefore we can use the same arguments as in the proof of Theorem 5.6 to obtain
convergence and, when HS is monotone, the monotonicity property ensures that the
limit is the unique entropy solution of (1.1), (1.2).
6. Numerical experiments.
6.1. The numerical schemes. System (1.3) is split into a linear diagonal hyperbolic
part and an ordinary di#erential system. For given f #,n
# is an
approximate solution at time t n+1 of the problem
ld # xd f
# .
As the system is diagonal, we may consider each equation separately so that we have
to approximate the scalar problem
where the # d are real and w n
# is a piecewise constant function given in L 1 (R D
(R D ) #BV (R D ). We present here two methods. Both are constructed on a cartesian
grid; see notations in section 4. The first is a straightforward generalization of the
upwind scheme: one solves exactly (6.3) on [t n , t n+1 ], obtaining a piecewise constant
1992 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
function -
w# . w n+1
# is then calculated by taking the average of -
w# on each cell. We
obtain the following explicit formulations:
For D=1,
where
For D=2, we set
.
For D=3, we set
i-u,j-v,k-w
i-u,j-v,k
i-u,j,k .
A second order MUSCL type method can also be applied, generalizing the one-dimensional
schemes, by the following steps:
(1) given the piecewise constant function w n
# we construct a piecewise linear
function -
#,d for x #]x #- 1, x #+ 1[ ;
(2) we solve on [t n , t n+1 ] the linear system (6.2) exactly with initial data:
(3) we compute cell average of the resulting solution to obtain w n+1 .
The method depends on the choice of # n
# . For example we can choose
#x#
#+ed - w n
#-ed
d (w n
#-ed , w n
#+ed
where the minmod function is defined by
and e
This choice corresponds to a linear interpolation of the piecewise constant function
# on the neighboring cells of I # with slope limiter. The resulting formulas for
are, respectively, (6.4), (6.5) to which one adds the following correction
terms:
For D=1,
For D=2 and
A similar formula holds for D=3.
Let us write down the one-dimensional numerical flux for the relaxed scheme
A n
l,# l <0
l,# l >0
Note that each Maxwellian function appears individually in the computation of the
slopes #.
6.2. Application to some models. Let us apply the above formulas to some
of the approximations given in section 3. One has to write down (6.4), (6.8), (6.5),
for each of the N blocks of K equations
and to apply (5.17)-(5.19).
For the approximation DRM in one space dimension, actually the Jin-Xin ap-
proximation, (6.4) gives the following scheme:
#x
#t
#x
As already pointed out in [8] and [2], in the scalar case, if #t
Lax-Friedrichs scheme.
Actually, at first order the approximation OVM of Example 3 in section 3 may
be thought of as a generalization of the relaxation approximation. A straightforward
calculation gives the following scheme associated with (6.4):
#x
#t
#x
i-1,
1994 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI
where
if N is even, and
if N is odd. The viscosity coe#cient is #N #t
#x and it is easy to verify that for a scalar
linear conservation law the monotonicity condition, as well as the stability condition,
ensure L 2 stability.
The model OVM is an interesting case where the monotonicity condition (3.19) is
strictly stronger than the stability condition (3.18) issued from the Chapman-Enskog
analysis. Actually if we just impose (3.18) the TVD property is lost. Nevertheless
the numerical experiments are satisfying (see also Remark 6.1 below). Another point
is that (3.19) is not the minimal condition for (6.12) to be monotone. It is su#cient
to impose the intermediate condition
if N # 3 is odd. With this condition (6.12) can be also viewed as scheme (6.11) where
# is replaced by #N # |A # (u)|.
Refer now to the second example of section 3, the FDM. The numerical scheme
corresponding to the discretization (6.4) is given by
Remark that # does not appear explicitly but is needed to ensure that the CFL condition
is satisfied. As already pointed out in [8], in the convex scalar (one-dimensional)
case we recover the Engquist-Osher's scheme. Our construction improves this scheme
(in the general case) by the correction terms (6.8) of the MUSCL discretization.
As far as one is concerned with first-order approximations our method retrieves
some known schemes. This is not true for the second-order MUSCL discretizations:
formula (6.10) shows that the computation of the slopes # involves the Maxwellian
functions. Therefore we do not recover in a direct way any known scheme. However
our work shares the spirit of the central, Riemann solver free schemes of [56, 34] and
a closer comparison should be useful.
We write the two-dimensional schemes only for the FDM, which has been found
to be the most e#cient. Set We have
#t
second order correction terms.
Here again # does not appear explicitly.
6.3. The numerical tests. In this section we perform some numerical tests
with our schemes. The systems considered are a simple p-system, the one-dimensional
Euler equations, a two-dimensional scalar conservation law, and the two-dimensional
Euler equations.
As we have observed that the relaxed numerical solution is in all cases
better than the numerical solution of (1.3), we present here only computations performed
with The HS is always chosen to be the MUSCL type one exposed
above.
(1) The p-system. We consider the system
In all our computations we take Exact solutions are known
for this system; see, for example, [63]. We compute the 1-shock, 2-rarefaction solution
of the Riemann problem with data:
0.4 for x > 0,
Here we discretize the di#erent approximations: DRM, 4 velocities, and FDM, 6
velocities, and we compute the L 1 error between exact and calculated solutions. We
have also computed the solution for Example 4 in section 3, Suliciu's method, but
there is almost no di#erence with that of DRM. For the same system we also test the
approximation OVM, with 16 and 26 velocities (N=8 and N=13, respectively).
Let us recall that this system has the eigenvalues - # (u). At each time step
we take the minimal value provided by the stability condition
(2.8). Here we have
so that the stability condition (3.18) for the OVM reads
The computation has been performed on a space interval [-2, 2] and with a maximal
1. The space step has been kept constant and equal to 0.01 and the
ratio #t/#x varies between 1 and 0.1. The results are given in Table 6.1.
Remark 6.1. In view of (6.12) one may think that for #t small it is possible
to take # smaller than # (u) # 3(N-1)
N+1 . We tried to perform the computation with
# (u). For ratios # > 0.5 we observed oscillations around discontinuities and
rarefactions, but then when # 0.5 the results are correct and even better than those
obtained with the right condition.
1996 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI

Table
exact and numerical solution,
1.0 8.35E-03 2.984E-02 6.302E-03 6.122E-03
0.7 7.34E-03 7.262E-03 6.045E-03 6.035E-03
0.4 6.85E-03 5.281E-03 6.019E-03 5.925E-03
Then we observe the evolution of the L 1 error when the space length step varies
and the ratio #t/#x is kept constant and equal to 0.5. Here the convergence
speed is defined by
#k
where e k is the L 1 error. The results are the following.

Table
exact and numerical solution.
2.00E-02 1.226E-02 .87 9.854E-03 .94 1.058E-02 1.031E-02
5.00E-03 4.035E-03 .84 3.177E-03 .84 3.314E-03 3.287E-03
2.50E-03 2.237E-03 .85 1.769E-03 .84 1.777E-03 1.760E-03

Tables

6.1 and 6.2 show that the relaxation model is the worst. The other three
seem to be comparable but the computation is faster for the FDM, since there are
fewer equations in the model and also because the value of # is smaller, allowing a
greater time step to reach the same value of t. Actually for the same value of #t the
results are better for the FDM. We complete this results by the graphic representation
of the computed and exact solutions (Figures 6.1-6.2). we represent the v
component of the solution, the u component being similar. It appears that the OVM
with 26 velocities is nearly the same as the relaxation approximation. The scheme
FDM improves the approximation of the one-shock but not that of the rarefaction,
while the OVM improves both slightly.
(2) One-dimensional Euler system. We now consider the one-dimensional Euler
system
where #, v, p and E are, respectively, the density, velocity, pressure, and total energy
of a perfect gas:
(DRM)
exact solution
Fig. 6.1. p-system: diagonal relaxation method (DRM) and orthogonal velocities method (OVM).0.450.550.65-1 -0.5 0 0.5 1
Fig. 6.2. p-system: flux decomposition method (FDM) and orthogonal velocities method (OVM).
We have tested our schemes with a Sod shock tube. Here the solution is constituted
by three constant states connected by a one-rarefaction, a two-contact discontinuity,
and a three-shock. We have analyzed the L 1 error for the same three models as for
the p-system. In Table 6.3 #t
is defined as above.

Table
exact and numerical solution.
2.50E-03
exact solution0.20.61-0.4 -0.3 -0.2 -0.1 0
exact solution
Fig. 6.3. One-dimensional Euler system: density and velocity with flux decomposition method
exact solution
Fig. 6.4. One-dimensional Euler system: pressure with flux decomposition method (FDM).
In view of this table, it appears that the approximation FDM and the 39 velocities
model of the OVM give similar results. Actually the graphic representation shows an
oscillation around the three-shock for the 39 velocities model while the first is stable
(see

Figures

6.3, 6.4, 6.5(a)). On Figure 6.5(b) we have suppressed these oscillations
by taking #t
All the computations show that the contact discontinuity is
not well approximated. This is a general feature of kinetic schemes (see, for example,
[24]). A smaller time step does not improve the contact discontinuity very much but
a smaller space step does, as shown in Figure 6.6.
(3) A two-dimensional scalar conservation law. In order to compare exact solutions
with numerical computations we also consider a one-dimensional conservation
law with one-dimensional initial data v 0 . This problem is solved on a two-dimensional
cartesian grid not parallel to the axis. We solve
u(x, y,
where # [0, #/2[. Putting u(t, x, R is the rotation of angle
# we have
v(X, Y,
exact solution0.20.61-0.4 -0.3 -0.2 -0.1 0
(a) (b)
Fig. 6.5. One-dimensional Euler system: density with orthogonal velocities method (OVM),
-0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4
exact solution
Fig. 6.6. One-dimensional Euler system: density with flux decomposition method (FDM), two
di#erent space steps,
Here we take #/5 and a rectangular mesh such that #y = 0.8#x. The CFL
number is 0.8 for x and y. We consider Burgers' equation:
a shock wave: v
As the relaxation model has been proved to be worse, we concentrate our attention
on the second-order model FDM, here with five velocities, and on the model
OVM. Recall that for these last models one-dimensional tests show that it is not very
e#cient to take a large number J of modulus for the velocities. The two-dimensional
computations confirm this analysis. Figure 6.7 represents the solution for at
In the first we compare the
cases 1. In the second we compare the
cases It appears that the best choice
for this model is to take the minimal values N

Figure

6.8 represents the same solution with, respectively,
but here we compare the model FDM to the one with
Here again the five velocities model is the best one. We want to point out that the
computation durations are nearly the same in both cases.
2000 DENISE AREGBA-DRIOLLET AND ROBERTO NATALINI0.20.61
exact solution0.20.61
Fig. 6.7. 2D scalar conservation law at
Fig. 6.8. Two-dimensional scalar conservation law at
Two-dimensional Euler system. We end this section with the two-dimensional
Euler system
where #, v, p and E are, respectively, the density, velocity, pressure, and total energy
of a perfect gas:
The initial data are chosen in order to represent a "double Sod tube":
We have chosen here the approximation FDM, here with 20 velocities. The calculation
has been performed with a rectangular mesh where and the CFL
Fig. 6.9. Two-dimensional Euler system: density and velocity x-component.
Fig. 6.10. Two-dimensional Euler system: velocity y-component and energy.
condition has been fixed to 0.4 in x and y. Figures 6.9 and 6.10 represent the isolines
of, respectively, the density, the first and second component of the velocity, and the
energy at time 0.16.

Acknowledgment

. The authors would like to thank Vuk Milisic for performing
the numerical tests for the two-dimensional Euler system and for many useful
discussions.



--R

An extended thermodynamic framework for the hydrodynamic modeling of semicon- ductors
Convergence of relaxation schemes for conservation laws


A model for collision processes in gases.
Construction of BGK models with a family of kinetic entropies for a given system of conservation laws


A kinetic formulation for multi-branch entropy solutions of scalar conservation laws
The unique limit of the Glimm scheme
The fluid-dynamical limit of nonlinear model of Boltzmann equations
The Boltzmann Equation and Its Applications
The Mathematical Theory of Dilute Gases
Hyperbolic conservation laws with sti
Zero relaxation and dissipation limits for hyperbolic conservation laws

Convergence of the relaxation approximation to a scalar nonlinear hyperbolic equation arising in chromatography.
Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics
Numerical passage from kinetic to fluid equations
Monotone di
Special issue on lattice gas methods for PDEs: Theory
A kinetic construction of global solutions of first order quasilinear equations
Hyperbolic Systems of Conservation Laws
Numerical Approximation of Hyperbolic Systems of Conservation Laws
Weakly coupled systems of quasilinear hyperbolic equations
Stability of traveling wave solutions for a rate-type viscoelastic system
Zero relaxation limit to centered rarefaction waves for a rate-type viscoelastic system
Nonlinear stability of rarefaction waves for a rate-type viscoelastic system
Kinetic formulation of the chromatography and some other hyperbolic systems
The relaxation schemes for systems of conservation laws in arbitrary space dimensions
Convergence and error estimates of relaxation schemes for multidimensional conservation laws
Contractive relaxation systems and the scalar multidimensional conservation law

New high-resolution central schemes for nonlinear conservation laws and convection-di#usion equations
Convergence of a relaxation scheme for n-n hyperbolic systems of conservation laws
Shock waves and entropy
A linear hyperbolic system with a sti
Moment closure hierarchies for kinetic theories
Zero Relaxation Limit to Piecewise Smooth Solutions for a Rate-Type Viscoelastic System in the Presence of Shocks
Existence and stability of entropy solutions for the hyperbolic systems of isentropic gas dynamics in Eulerian and Lagrangian coordinates
A kinetic formulation of multidimensional scalar conservation laws and related equations
Kinetic formulation of the isentropic gas dynamics and p-systems
Stability of a relaxation model with a nonconvex flux
Hyperbolic conservation laws with relaxation
Positive schemes for solving multi-dimensional hyperbolic systems of conservation laws
Asymptotic stability of planar rarefaction waves for the relaxation approximation of conservation laws in several dimensions
BV solutions and relaxation limit for a model in viscoelasticity
Linear stability of shock profiles for a rate-type viscoelastic system with relaxation
Nonlinear stability of shock fronts for a relaxation system in several space dimensions
stability of travelling waves for a hyperbolic system with relaxation

Springer Tracts Nat.
Convergence to equilibrium for the relaxation approximations of conservation laws
A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws
Recent results on hyperbolic relaxation problems
Nonoscillatory central di
Numerical methods for hyperbolic conservation laws with sti
Numerical methods for hyperbolic conservation laws with sti

A kinetic equation with kinetic entropy functions for scalar conservation laws
Discrete velocity models of the Boltzmann equation: A survey on the mathematical aspects of the theory
Recent advances in lattice Boltzmann computing

Relaxation semi-lin-eaire et cin-etique des syst-emes de lois de conservation
On modelling phase transitions by means of rate-type constitutive equations
Some stability-instability problems in phase transitions modelled by piecewise linear elastic or viscoelastic constitutive equations
Pointwise error estimates for relaxation approximations to conservation laws
On the rate of convergence to equilibrium for a system of conservation laws with a relaxation term
Viscosity and relaxation approximations for hyperbolic systems of conservation laws
Convergence of relaxing schemes for conservation laws
Linear and Nonlinear Waves
The fluid dynamical limit of the Broadwell model of the nonlinear Boltzmann equation in the presence of shocks
Numerical Analysis of Relaxation Schemes for Scalar Conservation Laws

--TR

--CTR
Mapundi K. Banda, Variants of relaxed schemes and two-dimensional gas dynamics, Proceedings of the international conference on Computational methods in sciences and engineering, p.56-59, September 12-16, 2003, Kastoria, Greece
Mapundi Kondwani Banda, Variants of relaxed schemes and two-dimensional gas dynamics, Journal of Computational and Applied Mathematics, v.175 n.1, p.41-62, 1 March 2005
Ansgar Jngel , Shaoqiang Tang, Numerical approximation of the viscous quantum hydrodynamic model for semiconductors, Applied Numerical Mathematics, v.56 n.7, p.899-915, July 2006
D. Aregba-Driollet , R. Natalini , S. Tang, Explicit diffusive kinetic schemes for nonlinear degenerate parabolic systems, Mathematics of Computation, v.73 n.245, p.63-94, January 2004
Ansgar Jngel , Shaoqiang Tang, A relaxation scheme for the hydrodynamic equations for semiconductors, Applied Numerical Mathematics, v.43 n.3, p.229-252, November 2002
conservation laws;hyperbolic systems;kinetic schemes;BGK models;numerical convergence