Hyperbolic systems of conservation laws describe in a non-viscous approximation the phenomena that arise when flowing around aerodynamic forms, in rocket nozzles, gas jets, propagation of polluting gases in the atmosphere and nuclear explosions  . To date, various methods have been developed in Eulerian coordinates for the numerical solutions of these systems. Description of the most common methods can be found in  -  .
It is possible to divide the existing numerical methods for solving quasilinear partial differential equations of the hyperbolic type into two large groups:
・ First group is the methods that essentially use the solution of the Riemann problem but do not use the approximate solution.
・ Second group is called Riemann solvers (RS method). The most complete description of RS methods for solving hyperbolic problems of various dimensions is available in   .
Numerical methods for solving hyperbolic equations that do not use the solution of the Riemann problem are called Non-Riemann-Solvers (NRS). Some known NRS are homogeneous difference schemes with artificial viscosity     , conservative difference schemes    , total variation diminishing (TVD) schemes    , finite volume schemes  and compact difference schemes    . In numerical solutions obtained by both RS and NRS methods, shock waves are smeared on several intervals of the spatial computational grid, while the thickness of the transition zone remains approximately constant in time. Early methods of second-order accuracy such as the Lax-Wendrof method  and McCormack method  , as well as the third-order Rusanov scheme  , Burshtein-Mirin scheme  , Balakin scheme  , Warming-Kutler-Lomax schemes  were obtained by expanding the grid functions into Taylor series. However, on the other hand, early schemes of high accuracy are characterized by the presence of parasitic oscillations of the numerical solution in the surrounding area of strong discontinuities. During the last 30 years, some methods have been developed for reducing the amplitude of these oscillations. Descriptions of some of these methods can be found in  . Such monotonic and quasi-monotonic schemes of high accuracy are highly in accurate as compared with the first-order schemes for numerical simulation of multidimensional problems with many interacting shock waves and contact discontinuities.
TVD schemes are proposed in  and they are used as the main tool of calculators working in the field of supersonic aerodynamics. The main advantages of TVD-schemes are the absence of the unphysical oscillations on the discontinuities and the fulfillment of the condition of non-decreasing entropy. As known, in TVD-schemes the transition to first-order accuracy schemes is carried out with the aim of monotonic numerical solution, but as a result of this, intense smearing of the discontinuous occurs. Bona et al.  proposed the total variation bounded (TVB) schemes with third and fifth orders accuracy. However, the replacement of the TVD condition by the TVB condition led to the appearance of significant parasitic oscillations of the numerical solution in the surrounding area of strong discontinuities, as it was clearly shown in Bona et al.  . In  , comparisons were made between the existing RS and NRS methods for solving the Euler equations on a large number of one and two dimensional test problems. It was found that the accuracy of both RS and NRS methods were comparable.
It should be noted that all of the above-mentioned schemes are basically for solving the Cauchy problems   . However, it is known   that in spite of the stability of the difference Cauchy problem, the initial boundary value difference problem may not be stable. Therefore it is essential to take into account the boundary conditions. As a fundamental factor in the construction and investigation of difference schemes      , we shall require the adequacy of the difference model to the original differential problem. The difference model for a hyperbolic system was constructed in such a way that it eventually allowed the derivation of difference analogs of the a priori estimate of the solution of the original differential problem. The latter circumstance seems to be an extremely important fact, since in numerical calculations approximate solution tends to be the solution of the original differential problems.
In the present paper, we select a class of quasi-linear hyperbolic systems that allow the construction of the energy integrals. The aim of our work is to construct the difference schemes for which the discrete analogue of the energy integrals is valid. We obtain a global―a priori estimate of the solutions, and construct the corresponding difference schemes. It is an attempt to systematically expound the technique of constructing the difference analogue of the energy integrals and its application in the study of the stability of difference schemes in computational practice.
2. Differential Problem
In this section we give a differential statement of the problem, for which then in the next section we construct a difference scheme. In addition, we investigate some properties of the solution of the differential problem. We start with the definition of a hyperbolic system according to  .
Definition 1. (Godunov  ) Let
be symmetric matrices and matrix A be positive defined. Assume that all elements
and f are sufficiently smooth functions of
, then system of equation
is called a symmetric t-hyperbolic system.
Symmetric hyperbolic systems might have important relationship with their solutions. These relations that generalize the law of conservation of energy for the solution of acoustics or Maxwell’s equations, is called energy integrals. It plays a fundamental role in the construction of the whole theory of symmetric systems.
Consider the quasilinear systems of the form
Suppose that the system allows the following form
Since the matrices
are symmetry, Equation (2) can be written as
Applying the dot product to the system (1)-(2) with vector u yields
Using the fact that
, we transform each terms as follows
Then (2**) becomes
We consider some region G lying inside the domain of existence of the solution u bounded by a piecewise smooth surface S. Integrating Equation (3) over G, gives
The integral on the left can be transformed into a surface integral by the multidimensional divergence theorem
the unit vector of the outer normal to the surface S. The integral identity (Vorozhtsov  )
is called the energy integral for a symmetric system.
Problem 1. Consider the initial-boundary value problem for system (1)-(2) in the region
with periodic boundary conditions
are positive real numbers, and for any fixed
with the initial data t = 0,
is a given function such that
Let us consider the solution of the symmetric t-hyperbolic system (1)-(2) in the domain G, and obtain the energy integrals
It is convenient to apply this identity not to the whole region G, but only
, bounded by the planes
then we obtain
From this equality, using the periodicity of the boundary conditions (4) and (5), we deduce the following equality
Suppose that the matrix B has a special canonical form:
Such a diagonal form of the matrix B is essentially used when specifying the boundary conditions for the following problems.
Problem 2. Let us consider the initial-boundary value problem for system (1)-(2) in the region G with boundary conditions for
and with the initial data (6). It is assumed that the condition (5) is satisfied. Note that
are the rectangular matrices of dimension
The boundary conditions (8) given on the boundary
are said to be dissipative at the points of this boundary if vector function
satisfies the following inequality 
The boundary conditions (9) given on the boundary
are said to be dissipative if at the points of this boundary the vector function
satisfies the following inequality 
For the solution of the initial boundary value problems (1), (2), (6), (8), (9), we obtain the following equality
From this and using dissipative boundary conditions together with
we obtain the inequality
Remark: As an example, we consider a system of equations describing the three-dimensional motion of a gas, under the assumption that the gas is in viscid, not thermally conductive, and is in a state of local thermodynamic equilibrium. In  , it is shown that for the defined conditions the three-dimensional system of equations of gas dynamics could be represented in the form (1) (2).
3. Difference Scheme with Limiter
In this section, we describe a difference scheme, which can be used to approximately solve a dissipative boundary value problem. Then we obtain the estimates for the solutions of difference equations which is analog to estimates of energy integrals.
Let us consider the difference grid points
are step size of the difference grid. For the value of the solution at the difference grid points, we introduce the following notation
We assume that А is a unit matrix. The difference model for problem (1) (2) with initial-boundary value (6), (8), (9) is formulated as follows
are the corresponding transpose matrices. Here we have omitted the indices which is not having shifts relative to the points
We consider the following reconstruction
is a continuous function which is called limiter.
Theorem: The solutions of the finite-difference scheme (10)-(12) is stable in energetic norm
Proof: Now we prove that the difference model (10)-(12) admits the presence of difference analog of the dissipative energy integral. This gives us the possibility to get the energetic estimation (the difference analog a priori estimation (7)), from which follows stability of the difference Scheme (10)-(12). We multiply the system (10) by the vector
We transform each summand to give
Taking all these transformations into account, we obtain
We multiply both sides of the Equation (13) by
and sum up over
, over j from
and over k from
Then, using the following relations
which means the stability of the difference scheme (10)-(12) in the energetic norm
4. Numerical Example
Example 1. As an example, we consider the Burger’s hyperbolic equation given by:
We introduce the following notations
Consider the following reconstruction:
corresponds to the scheme of the first order,
is the one sided upwind scheme of second order. Consider the following modification of the obtained scheme
We now prove that the difference scheme (15) admits the availability of difference analogue of the energy integral. Multiplying the system (15) by
Using formulas of difference differentiation, we obtain the following identity
Taking into account all of these transformations, we have
Multiply both sides of the (16) by
and sum up over i from
, and noting that the function u tends to zero at infinity and denoting the quantity
we obtain the equality
From this it is easily follows the energetic estimate
which means stability of the difference scheme in the norm
Example 2. Consider the following Cauchy problem
Rewrite the difference scheme (15) in the following form
In solving the difference scheme (19) we apply the iteration method with respect to the nonlinear coefficients
We have carried out a posteriori error analysis of the proposed scheme. Table 1 presents the exact and numerical solutions at points
for the Cauchy problem (18).
Figures 1-4 exhibits the graphical results of the exact and numerical solutions. Comparisons between numerical and exact solutions are presented in Figures 1-4. It can be concluded that the scheme with limiter well modulates the jump. All results are obtained in Mathcad.
Table 1. Exact and numerical solution for problem (18).
Figure 1. Exact solution of the problem of (18).
Figure 2. Numerical solution by scheme (19).
Figure 3. Red circle and solid blue lines are the exact and numerical solutions respectively for t = 0.4.
Figure 4. Red circle and solid blue lines are the exact and numerical solutions respectively for t = 1.4.
The class of three-dimensional quasilinear hyperbolic systems is studied. The formulation of initial boundary value problem for this class of quasilinear hyperbolic systems in two variants is given. A priori estimate of the solution of initial boundary value problem is obtained by constructing an energy integral. Difference scheme with limiter is constructed and a priori estimate for its solution is obtained. Numerical results for the developed schemes show agreement with exact solution.
This work was supported by Research Management Center (RMC), Universiti Sains Islam Malaysia (USIM) under Research Grand PPP/USG-0216/FST/30/15316.
This research was supported in part under R. A. Welch Foundation Grant E-0608.