A Discrete Analogue of Energy Integral for a Difference Scheme for Quasilinear Hyperbolic Systems

R. D. Aloev^1*, Z. K. Eshkuvatov^2,3*, M. U. Khudayberganov¹, N. M. A. Nik Long^3,4

¹Faculty of Mathematics and Mechanics, National University of Uzbekistan (NUUz), Tashkent, Uzbekistan.
²Faculty of Science and Technology, Universiti Sains Islam Malaysia (USIM), Nilai, Malaysia.
³Institute for Mathematical Research (INSPEM), UPM, Kuala Lumpur, Malaysia.
⁴Department of Mathematics, Faculty of Science, Universiti Putra Malaysia (UPM), Kuala Lumpur, Malaysia.

DOI: 10.4236/am.2018.97055                

Abstract

The class of three-dimensional quasilinear hyperbolic systems is studied. The initial boundary value problem for this class of quasilinear hyperbolic systems is given. By constructing the energy’s integral, a priori estimate for the solution of the initial boundary value problem is obtained. Difference scheme is constructed and an a priori estimate for its solution is obtained. Numerical example exhibits the efficiency and accuracy of the method.

Keywords

Hyperbolic Systems, Stability, Quasilinear Equations

1. Introduction

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 [1] . 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 [2] - [8] .

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 [6] [9] .

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 [2] [10] [11] [12] , conservative difference schemes [2] [4] [11] , total variation diminishing (TVD) schemes [4] [6] [7] , finite volume schemes [13] and compact difference schemes [3] [5] [14] . 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 [15] and McCormack method [16] , as well as the third-order Rusanov scheme [17] , Burshtein-Mirin scheme [18] , Balakin scheme [19] , Warming-Kutler-Lomax schemes [20] 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 [6] . 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 [21] 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. [22] 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. [22] . In [23] , 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 [24] [25] . However, it is known [26] [27] 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 [26] [27] [28] [29] [30] , 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 [26] .

Definition 1. (Godunov [26] ) Let $A,B,C,D$ be symmetric matrices and matrix A be positive defined. Assume that all elements $A,B,C,D$ and f are sufficiently smooth functions of $t,x,y,z,u$ , then system of equation

$\begin{array}{l}A\left(t,x,y,z,u\right)\frac{\partial u}{\partial t}+B\left(t,x,y,z,u\right)\frac{\partial u}{\partial x}+C\left(t,x,y,z,u\right)\frac{\partial u}{\partial y}+D\left(t,x,y,z,u\right)\frac{\partial u}{\partial z}\\ =f\left(t,x,y,z,u\right)\end{array}$

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

$A\frac{\partial u}{\partial t}+B\frac{\partial u}{\partial x}+C\frac{\partial u}{\partial y}+D\frac{\partial u}{\partial z}=0,$ (1)

where

$\begin{array}{l}A=A\left(t,x,y,z,u\right),B=B\left(t,x,y,z,u\right),\\ C=C\left(t,x,y,z,u\right),D=D\left(t,x,y,z,u\right)\end{array}$

and

$A=A^{*},B=B^{*},C=C^{*},D=D^{*}.$

Suppose that the system allows the following form

$\frac{\partial }{\partial t}\left(A^{*}u\right)+\frac{\partial }{\partial x}\left(B^{*}u\right)+\frac{\partial }{\partial y}\left(C^{*}u\right)+\frac{\partial }{\partial z}\left(D^{*}u\right)=0.$ (2)

Since the matrices $A,B,C,D$ are symmetry, Equation (2) can be written as

$\frac{\partial }{\partial t}\left(Au\right)+\frac{\partial }{\partial x}\left(Bu\right)+\frac{\partial }{\partial y}\left(Cu\right)+\frac{\partial }{\partial z}\left(Du\right)=0.$ (2*)

Applying the dot product to the system (1)-(2) with vector u yields

$\begin{array}{l}\left(A\frac{\partial u}{\partial t},u\right)+\left(\frac{\partial }{\partial t}\left[Au\right],u\right)+\left(B\frac{\partial u}{\partial x},u\right)+\left(\frac{\partial }{\partial x}\left[Bu\right],u\right)\\ +\left(C\frac{\partial u}{\partial y},u\right)+\left(\frac{\partial }{\partial y}\left[Cu\right],u\right)+\left(D\frac{\partial u}{\partial z},u\right)+\left(\frac{\partial }{\partial z}\left[Du\right],u\right)=0.\end{array}$

Using the fact that $A=A^{*},B=B^{*},C=C^{*},D=D^{*}$ , we transform each terms as follows

$\begin{array}{l}\left(A\frac{\partial u}{\partial t},u\right)+\left(\frac{\partial }{\partial t}\left[Au\right],u\right)=\left(\frac{\partial u}{\partial t},Au\right)+\left(\frac{\partial }{\partial t}\left[Au\right],u\right)\\ =\left(Au,\frac{\partial u}{\partial t}\right)+\left(\frac{\partial }{\partial t}\left[Au\right],u\right)=\frac{\partial }{\partial t}\left(Au,u\right).\end{array}$ (2**)

Similarly

$\begin{array}{l}\left(B\frac{\partial u}{\partial x},u\right)+\left(\frac{\partial }{\partial x}\left[Bu\right],u\right)=\frac{\partial }{\partial x}\left(Bu,u\right),\\ \left(C\frac{\partial u}{\partial y},u\right)+\left(\frac{\partial }{\partial y}\left[Cu\right],u\right)=\frac{\partial }{\partial y}\left(Cu,u\right),\\ \left(D\frac{\partial u}{\partial z},u\right)+\left(\frac{\partial }{\partial z}\left[Du\right],u\right)=\frac{\partial }{\partial z}\left(Du,u\right).\end{array}$

Then (2**) becomes

$\frac{\partial }{\partial t}\left(Au,u\right)+\frac{\partial }{\partial x}\left(Bu,u\right)+\frac{\partial }{\partial z}\left(Du,u\right)+\frac{\partial }{\partial y}\left(Cu,u\right)=0.$ (3)

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

${\displaystyle \int {\displaystyle \int {\displaystyle \int {\displaystyle {\int }_G\left[\frac{\partial }{\partial t}\left(Au,u\right)+\frac{\partial }{\partial x}\left(Bu,u\right)+\frac{\partial }{\partial y}\left(Cu,u\right)+\frac{\partial }{\partial z}\left(Du,u\right)\right]\text{d}t\text{d}x\text{d}y\text{d}z}}}}=0$ .

The integral on the left can be transformed into a surface integral by the multidimensional divergence theorem

${\displaystyle {∰}_S\left[\tau \left(Au,u\right)+\xi \left(Bu,u\right)+\eta \left(Cu,u\right)+\zeta \left(Du,u\right)\right]\text{d}s}=0,$

where $\left(\tau ,\xi ,\eta ,\zeta \right)$ the unit vector of the outer normal to the surface S. The integral identity (Vorozhtsov [1] )

${\displaystyle {∰}_S\left(\left[\tau A+\xi B+\eta C+\zeta D\right]u,u\right)\text{d}s}=0$ ,

is called the energy integral for a symmetric system.

Problem 1. Consider the initial-boundary value problem for system (1)-(2) in the region $G=\left\{\left(t,x,y,z\right):0<t\le T;0<x<l,\left|y\right|<\infty ,\left|z\right|<\infty \right\}$ with periodic boundary conditions

$u\left(t,0,y,z\right)=u\left(t,l,y,z\right)$ , ${B|}_{x=0}={B|}_{x=l}$ (4)

where $T,l$ are positive real numbers, and for any fixed $t,x$ ,

$\left(Сu,u\right)\underset{\left|y\right|,\left|z\right|\to +\infty }{\to }0,\left(Du,u\right)\underset{\left|y\right|,\left|z\right|\to +\infty }{\to }0$ (5)

with the initial data t = 0,

$u\left(0,x,y,z\right)=u_0\left(x,y,z\right)\text{,0}\le x\le l\text{,}\left|y\right|<\infty \text{,}\left|z\right|<\infty .$ (6)

and $u_0\left(x,y,z\right)$ is a given function such that

${\displaystyle {\int }_0^l{\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }\left(A\left(0,x,y,z,u_0\left(x,y,z\right)\right)u_0\left(x,y,z\right),u_0\left(x,y,z\right)\right)\text{d}z\text{d}y\text{d}x}}}<+\infty$ .

Let us consider the solution of the symmetric t-hyperbolic system (1)-(2) in the domain G, and obtain the energy integrals

${\displaystyle {∰}_S\left(\left[\tau A+\xi B+\eta C+\zeta D\right]u,u\right)\text{d}s}=0.$

It is convenient to apply this identity not to the whole region G, but only $t_1<t<t_2$ , bounded by the planes $t=t_1,t=t_2$ then we obtain

$\begin{array}{l}I\left(t_2\right)-I\left(t_1\right)+{\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\left[-{\left(Bu,u\right)|}_{x=0}+{\left(Bu,u\right)|}_{x=l}\right]\text{d}t\text{d}y\text{d}z}}}\\ +{\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{0}{\overset{l}{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left(Сu,u\right)|}_{-\infty }^{+\infty }\text{d}t\text{d}x\text{d}z}}}+{\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{0}{\overset{l}{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left(Du,u\right)|}_{-\infty }^{+\infty }\text{d}t\text{d}x\text{d}y=0}}},\end{array}$

where $I\left(t\right)={\displaystyle \underset{t=const}{\iiint }\left(Au,u\right)\text{d}x\text{d}y\text{d}z}$ .

From this equality, using the periodicity of the boundary conditions (4) and (5), we deduce the following equality

$I\left(t_2\right)=I\left(t_1\right),\forall t_1,t_2,0\le t_1<t<t_2\le T.$ (7)

Suppose that the matrix B has a special canonical form:

$B=diag\left(k_1,k_2,\cdots ,k_{n^{+}},-k_{n^{+}+1},\cdots ,-k_{n^{+}+n^{-}}\right),n^{+}+n^{-}=n;k_i>0,i=1,\cdots ,n.$

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 $x=0$

$u^I=S{u}^{II},0<t\le T,\left|y\right|<\infty ,\left|z\right|<\infty$ (8)

and for $x=l$

$u^{II}=R{u}^I,0<t\le T,\left|y\right|<\infty ,\left|z\right|<\infty$ (9)

and with the initial data (6). It is assumed that the condition (5) is satisfied. Note that

$\begin{array}{l}u={\left(u^I,u^{II}\right)}^{\text{T}}\text{,}{u}^I={\left(u_1\text{,}{u}_2\text{,}\cdots \text{,}{u}_{n^{+}}\right)}^{\text{T}}\text{,}{u}^{II}={\left(u_{n^{+}+1}\text{,}{u}_2\text{,}\cdots \text{,}{u}_n\right)}^{\text{T}},\\ S=S\left(u,t,y,z\right),R=R\left(u,t,y,z\right),\end{array}$

are the rectangular matrices of dimension $n^{+}\times n^{-}$ and $n^{-}\times n^{+}$ , respectively.

The boundary conditions (8) given on the boundary $x=0$ are said to be dissipative at the points of this boundary if vector function $\left(u_1\text{,}{u}_2\text{,}\cdots \text{,}{u}_n\right)$ satisfies the following inequality [26]

$-{\left(Bu,u\right)|}_{x=0}\ge 0.$

The boundary conditions (9) given on the boundary $x=l$ are said to be dissipative if at the points of this boundary the vector function $\left(u_1\text{,}{u}_2\text{,}\cdots \text{,}{u}_n\right)$ satisfies the following inequality [26]

${\left(Bu,u\right)|}_{x=l}\ge 0.$

For the solution of the initial boundary value problems (1), (2), (6), (8), (9), we obtain the following equality

$\begin{array}{l}I\left(t_2\right)-I\left(t_1\right)+{\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\left[-{\left(Bu,u\right)|}_{x=0}+{\left(Bu,u\right)|}_{x=l}\right]\text{d}t\text{d}y\text{d}z}}}\\ +{\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{0}{\overset{l}{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left(Сu,u\right)|}_{-\infty }^{+\infty }\text{d}t\text{d}x\text{d}z}}}+{\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{0}{\overset{l}{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left(Du,u\right)|}_{-\infty }^{+\infty }\text{d}t\text{d}x\text{d}y=0}}}.\end{array}$

From this and using dissipative boundary conditions together with

$-{\left(Bu,u\right)|}_{x=0}\ge 0,{\left(Bu,u\right)|}_{x=l}\ge 0,$

we obtain the inequality

$I\left(t_2\right)\le I\left(t_1\right),\forall t_1,t_2:0\le t_1<t<t_2\le T.$

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 [27] , 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

$\begin{array}{l}{t}_m=m\cdot \Delta t,x_i=i\cdot \Delta x,y_j=j\cdot \Delta y,z_k=k\cdot \Delta z,\\ m=0,\cdots ,M,i=0,1,\cdots ,N,\left|j\right|,\left|k\right|=0,1,\cdots ,+\infty ,\\ M\cdot \Delta t=T,N\cdot \Delta x=L,\end{array}$

where $\Delta t,\Delta x,\Delta y,\Delta z$ are step size of the difference grid. For the value of the solution at the difference grid points, we introduce the following notation

$u\left(m\cdot \Delta t,i\cdot \Delta x,j\cdot \Delta y,k\cdot \Delta z\right)=u_{ijk}^m=u.$

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

$\begin{array}{l}\left[U^{m+1}+U\right]\left(u^{m+1}-u\right)+r_{x}U\left(u_{i-1/2}^L\right)B^{+}\left(u_{i-1/2}^L\right)\left(u-u_{i-1}\right)\\ +r_{x}U\left\{\left[B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L\right]-\left[B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L\right]\right\}\\ +r_{x}U\left(u_{i+1/2}^R\right)B^{-}\left(u_{i+1/2}^R\right)\left(u_{i+1}-u\right)\\ +r_{x}U\left\{\left[B^T{}^{-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R\right]-\left[B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R\right]\right\}\\ +r_{y}U\left(u_{j-1/2}^L\right)C^{+}\left(u_{j-1/2}^L\right)\left(u-u_{j-1}\right)\end{array}$

$\begin{array}{l}+r_{y}U\left\{\left[C^T{}^{+}\left(u_{j+1/2}^L\right)u_{j+1/2}^L\right]-\left[C^T{}^{+}\left(u_{j-1/2}^L\right)u_{j-1/2}^L\right]\right\}\\ +r_{y}U\left(u_{j+1/2}^R\right)C^{-}\left(u_{j+1/2}^R\right)\left(u_{j+1}-u\right)\\ +r_{y}U\left\{\left[C^T{}^{-}\left(u_{j+1/2}^R\right)u_{j+1/2}^R\right]-\left[C^T{}^{-}\left(u_{j-1/2}^R\right)u_{j-1/2}^R\right]\right\}\\ +r_{z}U\left(u_{k-1/2}^L\right)D^{+}\left(u_{k-1/2}^L\right)\left(u-u_{k-1}\right)\\ +r_{z}U\left\{\left[D^T{}^{+}\left(u_{k+1/2}^L\right)u_{k+1/2}^L\right]-\left[D^T{}^{+}\left(u_{k-1/2}^L\right)u_{k-1/2}^L\right]\right\}\end{array}$

$\begin{array}{l}+r_{z}U\left(u_{k+1/2}^R\right)D^{-}\left(u_{k+1/2}^R\right)\left(u_{k+1}-u\right)\\ +r_{z}U\left\{\left[D^T{}^{-}\left(u_{k+1/2}^R\right)u_{k+1/2}^R\right]-\left[D^T{}^{-}\left(u_{k-1/2}^R\right)u_{k-1/2}^R\right]\right\}=0;\\ m=0,1,\cdots ,M-1;i=0,1,\cdots ,N-1;\left|j\right|,\left|k\right|\text{=}\left\{0,1,\cdots +\infty \right\}\text{,}\end{array}$ (10)

and

$\begin{array}{l}{u}_{0jk}^m=u_{Njk}^m\text{,}{u}_{1jk}^m=u_{N-1jk}^m\text{, }m=0,1,\cdots ,M,\left|j\right|,\left|k\right|=\left\{0,1,\cdots ,+\infty \right\}\text{, }\\ \left(B^T{}^{+}\left(u_{N-1/2}^L\right)u_{N-1/2}^L,u_{N-1}\right)+\left(B^{T-}\left(u_{N-1/2}^R\right)u_{N-1/2}^R,u_N\right)\\ =\left(B^T{}^{+}\left(u_{1/2}^L\right)u_{1/2}^L,u_0\right)+\left(B^T{}^{-}\left(u_{1/2}^R\right)u_{1/2}^R,u_1\right),\\ \left(C^T{}^{+}\left(u_{j+1/2}^L\right)u_{j+1/2}^L,u\right)+\left(C^T{}^{-}\left(u_{j+1/2}^R\right)u_{j+1/2}^R,u_{j+1}\right)\underset{\left|j\right|\to +\infty }{\to }0,\forall i,k,\\ \left(D^T{}^{+}\left(u_{k+1/2}^L\right)u_{k+1/2}^L,u\right)+\left(D^T{}^{-}\left(u_{k+1/2}^R\right)u_{k+1/2}^R,u_{k+1}\right)\underset{\left|j\right|\to +\infty }{\to }0,\forall i,k,\end{array}$ (11)

$\begin{array}{l}{u}_{ijk}^0=u_0\left(i\cdot \Delta x,j\cdot \Delta y,k\cdot \Delta z\right),i=0,1,\cdots ,N,\\ \left|j\right|,\left|k\right|=\text{0,1,}\cdots ,\text{+}\infty ,U=\text{diag}\left(u_1,u_2,\cdots ,u_n\right)\text{,}\end{array}$ (12)

where

$\begin{array}{l}B\left(u\right)=B^{+}\left(u\right)+B^{-}(u),B^{+}\left(u\right)\ge 0,B^{-}\left(u\right)\le 0,\\ C\left(u\right)=C^{+}\left(u\right)+C^{-}\left(u\right),C^{+}\left(u\right)\ge 0,C^{-}\left(u\right)\le 0,\\ D\left(u\right)=D^{+}\left(u\right)+D^{-}\left(u\right),D^{+}\left(u\right)\ge 0,D^{-}\left(u\right)\le 0,\\ \forall u\in {ℝ}^n,r_x=\Delta t/\Delta x,r_y=\Delta t/\Delta y,r_z=\Delta t/\Delta z,\\ B=B\left(u_{ijk}^m,m\cdot \Delta t,i\cdot \Delta x,j\cdot \Delta y,k\cdot \Delta z\right),\end{array}$

with $B^{T+},B^{T-},C^{T+},C^{T-},D^{T+},D^{T-}$ are the corresponding transpose matrices. Here we have omitted the indices which is not having shifts relative to the points $\left(m\cdot \Delta t,i\cdot \Delta x,j\cdot \Delta y,k\cdot \Delta z\right)$

$u=u_{ijk}^m,u^{m+1}=u_{ijk}^{m+1},u_{i\pm 1}=u_{i\pm 1jk}^m,u_{j\pm 1}=u_{ij\pm 1k}^m,u_{k\pm 1}=u_{ijk\pm 1}^m.$

We consider the following reconstruction

$\begin{array}{l}{u}_{i+1/2}^L=u+\frac{1}{2}\Psi \left(R_i\right)\left(u-u_{i-1}\right),u_{i-1/2}^R=u-\frac{1}{2}\Psi \left(\frac{1}{R_i}\right)\left(u_{i+1}-u\right),\\ u_{j+1/2}^L=u+\frac{1}{2}\Psi \left(R_j\right)\left(u-u_{j-1}\right),u_{j-1/2}^R=u-\frac{1}{2}\Psi \left(\frac{1}{R_j}\right)\left(u_{j+1}-u\right),\\ u_{k+1/2}^L=u+\frac{1}{2}\Psi \left(R_k\right)\left(u-u_{k-1}\right),u_{k-1/2}^R=u-\frac{1}{2}\Psi \left(\frac{1}{R_k}\right)\left(u_{k+1}-u\right),\end{array}$

$\begin{array}{l}\Psi \left(R\right)=\text{diag}\left(\psi \left(R_1\right),\psi \left(R_2\right),\cdots ,\psi \left(R_n\right)\right),\\ \Psi \left(\frac{1}{R}\right)=\text{diag}\left(\psi \left(\frac{1}{R_1}\right),\psi \left(\frac{1}{R_2}\right),\cdots ,\psi \left(\frac{1}{R_n}\right)\right),\\ {\left(R_q\right)}_i=\frac{{\left(u_q\right)}_{i+1}-{\left(u_q\right)}_i}{{\left(u_q\right)}_i-{\left(u_q\right)}_{i-1}},{\left(R_q\right)}_j=\frac{{\left(u_q\right)}_{j+1}-{\left(u_q\right)}_j}{{\left(u_q\right)}_j-{\left(u_q\right)}_{j-1}},{\left(R_q\right)}_k=\frac{{\left(u_q\right)}_{k+1}-{\left(u_q\right)}_k}{{\left(u_q\right)}_k-{\left(u_q\right)}_{k-1}},\end{array}$

where $\psi :R\to R$ is a continuous function which is called limiter.

Theorem: The solutions of the finite-difference scheme (10)-(12) is stable in energetic norm $\sqrt{I_m}$ , where $I_m=\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}\left(u_{ijk}^m,u_{ijk}^m\right)}}}$ .

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 $E={\left(1,1,\cdots ,1\right)}^{\text{T}}$

$\begin{array}{l}\left(\left[U^{m+1}+U\right]\left(u^{m+1}-u\right),E\right)\\ +r_x\left(U\left(u_{i-1/2}^L\right)B^{+}\left(u_{i-1/2}^L\right)\left(u-u_{i-1}\right),E\right)\\ +r_x\left(U\left\{\left[B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L\right]-\left[B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L\right]\right\},E\right)\\ +r_x\left(U\left(u_{i+1/2}^R\right)B^{-}\left(u_{i+1/2}^R\right)\left(u_{i+1}-u\right),E\right)\\ +r_x\left(U\left\{\left[B^T{}^{-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R\right]-\left[B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R\right]\right\},E\right)\end{array}$

$\begin{array}{l}+r_y\left(U\left(u_{j-1/2}^L\right)C^{+}\left(u_{j-1/2}^L\right)\left(u-u_{j-1}\right),E\right)\\ +r_y\left(U\left\{\left[C^T{}^{+}\left(u_{j+1/2}^L\right)u_{j+1/2}^L\right]-\left[C^T{}^{+}\left(u_{j-1/2}^L\right)u_{j-1/2}^L\right]\right\},E\right)\\ +r_y\left(U\left(u_{j+1/2}^R\right)C^{-}\left(u_{j+1/2}^R\right)\left(u_{j+1}-u\right),E\right)\\ +r_y\left(U\left\{\left[C^T{}^{-}\left(u_{j+1/2}^R\right)u_{j+1/2}^R\right]-\left[C^T{}^{-}\left(u_{j-1/2}^R\right)u_{j-1/2}^R\right]\right\},E\right)\\ +r_z\left(U\left(u_{k-1/2}^L\right)D^{+}\left(u_{k-1/2}^L\right)\left(u-u_{k-1}\right),E\right)\end{array}$

$\begin{array}{l}+r_z\left(U\left\{\left[D^T{}^{+}\left(u_{k+1/2}^L\right)u_{k+1/2}^L\right]-\left[D^T{}^{+}\left(u_{k-1/2}^L\right)u_{k-1/2}^L\right]\right\},E\right)\\ +r_z\left(U\left(u_{k+1/2}^R\right)D^{-}\left(u_{k+1/2}^R\right)\left(u_{k+1}-u\right),E\right)\\ +r_z\left(U\left\{\left[D^T{}^{-}\left(u_{k+1/2}^R\right)u_{k+1/2}^R\right]-\left[D^T{}^{-}\left(u_{k-1/2}^R\right)u_{k-1/2}^R\right]\right\},E\right)=0;\\ m=0,1,\cdots ,M-1;i=0,1,\cdots ,N-1;\left|j\right|,\left|k\right|=\left\{0,1,\cdots ,+\infty \right\}.\end{array}$

We transform each summand to give

$1)\left(\left[U^{m+1}+U\right]\left(u^{m+1}-u\right),E\right)=\left(\left[u^{m+1}-u\right],\left[u^{m+1}+u\right]\right)=\left(u^{m+1},u^{m+1}\right)-\left(u,u\right);$

$\begin{array}{l}2)\left(U\left(u_{i-1/2}^L\right)B^{+}\left(u_{i-1/2}^L\right)\left(u-u_{i-1}\right),E\right)\\ +\left(U\left\{\left[B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L\right]-\left[B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L\right]\right\},E\right)\\ =\left(B^{+}\left(u_{i-1/2}^L\right)\left(u-u_{i-1}\right),u_{i-1/2}^L\right)\\ +\left(\left\{\left[B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L\right]-\left[B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L\right]\right\},u\right)\end{array}$

$\begin{array}{l}=\left(\left[u-u_{i-1}\right],B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L\right)\\ +\left(u,\left\{\left[B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L\right]-\left[B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L\right]\right\}\right)\\ =\left(u,B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L\right)-\left(u_{i-1},B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L\right)\\ =\left(B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L,u\right)-\left(B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L,u_{i-1}\right).\end{array}$

$\begin{array}{l}3)\left(U\left(u_{i+1/2}^R\right)B^{-}\left(u_{i+1/2}^R\right)\left(u_{i+1}-u\right),E\right)\\ +\left(U\left\{\left[B^T{}^{-}(u_{i+1/2}^R)u_{i+1/2}^R\right]-\left[B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R\right]\right\},E\right)\\ =\left(B^{-}\left(u_{i+1/2}^R\right)\left(u_{i+1}-u\right),u_{i+1/2}^R\right)\\ +\left(\left\{\left[B^T{}^{-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R\right]-\left[B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R\right]\right\},u\right)\end{array}$

$\begin{array}{l}=\left(\left[u_{i+1}-u\right],B^{T-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R\right)\\ +\left(\left\{\left[B^T{}^{-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R\right]-\left[B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R\right]\right\},u\right)\\ =\left(u_{i+1},B^{T-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R\right)-\left(B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R,u\right)\\ =\left(B^{T-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R,u_{i+1}\right)-\left(B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R,u\right).\end{array}$

$\begin{array}{l}4)\left(U\left(u_{j-1/2}^L\right)C^{+}\left(u_{j-1/2}^L\right)\left(u-u_{j-1}\right),E\right)\\ +\left(U\left\{\left[C^T{}^{+}\left(u_{j+1/2}^L\right)u_{j+1/2}^L\right]-\left[C^T{}^{+}\left(u_{j-1/2}^L\right)u_{j-1/2}^L\right]\right\},E\right)\\ =\left(C^T{}^{+}\left(u_{j+1/2}^L\right)u_{j+1/2}^L,u\right)-\left(C^T{}^{+}\left(u_{j-1/2}^L\right)u_{j-1/2}^L,u_{j-1}\right);\end{array}$

$\begin{array}{l}5)\left(U\left(u_{j+1/2}^R\right)C^{-}\left(u_{j+1/2}^R\right)\left(u_{j+1}-u\right),E\right)\\ +\left(U\left\{\left[C^T{}^{-}\left(u_{j+1/2}^R\right)u_{j+1/2}^R\right]-\left[C^T{}^{-}\left(u_{j-1/2}^R\right)u_{j-1/2}^R\right]\right\},E\right)\\ =\left(C^T{}^{-}\left(u_{j+1/2}^R\right)u_{j+1/2}^R,u_{j+1}\right)-\left(C^T{}^{-}\left(u_{j-1/2}^R\right)u_{j-1/2}^R,u\right);\end{array}$

$\begin{array}{l}6)\left(U\left(u_{k-1/2}^L\right)D^{+}\left(u_{k-1/2}^L\right)\left(u-u_{k-1}\right),E\right)\\ +\left(U\left\{\left[D^T{}^{+}\left(u_{k+1/2}^L\right)u_{k+1/2}^L\right]-\left[D^T{}^{+}\left(u_{k-1/2}^L\right)u_{k-1/2}^L\right]\right\},E\right)\\ =\left(D^T{}^{+}\left(u_{k+1/2}^L\right)u_{k+1/2}^L,u\right)-\left(D^T{}^{+}\left(u_{k-1/2}^L\right)u_{k-1/2}^L,u_{k-1}\right);\end{array}$

$\begin{array}{l}7)\left(U\left(u_{k+1/2}^R\right)D^{-}\left(u_{k+1/2}^R\right)\left(u_{k+1}-u\right),E\right)\\ +\left(U\left\{\left[D^T{}^{-}\left(u_{k+1/2}^R\right)u_{k+1/2}^R\right]-\left[D^T{}^{-}\left(u_{k-1/2}^R\right)u_{k-1/2}^R\right]\right\},E\right)\\ =\left(D^T{}^{-}\left(u_{k+1/2}^R\right)u_{k+1/2}^R,u_{k+1}\right)-\left(D^T{}^{-}\left(u_{k-1/2}^R\right)u_{k-1/2}^R,u\right);\end{array}$

Taking all these transformations into account, we obtain

$\begin{array}{l}\left(u^{m+1},u^{m+1}\right)-\left(u,u\right)+r_x[\left(B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L,u\right)-\left(B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L,u_{i-1}\right)\\ +\left(B^{T-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R,u_{i+1}\right)-\left(B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R,u\right)]+r_y[\left(C^T{}^{+}\left(u_{j+1/2}^L\right)u_{j+1/2}^L,u\right)\\ -\left(C^T{}^{+}\left(u_{j-1/2}^L\right)u_{j-1/2}^L,u_{j-1}\right)+\left(C^T{}^{-}\left(u_{j+1/2}^R\right)u_{j+1/2}^R,u_{j+1}\right)-\left(C^T{}^{-}\left(u_{j-1/2}^R\right)u_{j-1/2}^R,u\right)]\\ +r_z[\left(D^T{}^{+}\left(u_{k+1/2}^L\right)u_{k+1/2}^L,u\right)-\left(D^T{}^{+}\left(u_{k-1/2}^L\right)u_{k-1/2}^L,u_{k-1}\right)\\ +\left(D^T{}^{-}\left(u_{k+1/2}^R\right)u_{k+1/2}^R,u_{k+1}\right)-\left(D^T{}^{-}\left(u_{k-1/2}^R\right)u_{k-1/2}^R,u\right)]=0.\end{array}$ (13)

We multiply both sides of the Equation (13) by $\Delta x\Delta y\Delta z$ and sum up over $i=1,\cdots ,N-1$ , over j from $-\infty$ to $+\infty$ and over k from $-\infty$ to $+\infty$

$\begin{array}{l}\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}\left(u^{m+1},u^{m+1}\right)}}}-\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}\left(u,u\right)}}}\\ +\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}{r}_x[\left(B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L,u\right)-\left(B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L,u_{i-1}\right)}}}\\ +\left(B^{T-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R,u_{i+1}\right)-\left(B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R,u\right)]\end{array}$

$\begin{array}{l}+\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}{r}_y[\left(C^T{}^{+}\left(u_{j+1/2}^L\right)u_{j+1/2}^L,u\right)-\left(C^T{}^{+}\left(u_{j-1/2}^L\right)u_{j-1/2}^L,u_{j-1}\right)}}}\\ +\left(C^T{}^{-}\left(u_{j+1/2}^R\right)u_{j+1/2}^R,u_{j+1}\right)-\left(C^T{}^{-}\left(u_{j-1/2}^R\right)u_{j-1/2}^R,u\right)]\\ +\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}{r}_z[\left(D^T{}^{+}\left(u_{k+1/2}^L\right)u_{k+1/2}^L,u\right)-\left(D^T{}^{+}\left(u_{k-1/2}^L\right)u_{k-1/2}^L,u_{k-1}\right)}}}\\ +\left(D^T{}^{-}\left(u_{k+1/2}^R\right)u_{k+1/2}^R,u_{k+1}\right)-\left(D^T{}^{-}\left(u_{k-1/2}^R\right)u_{k-1/2}^R,u\right)].\end{array}$

Then, using the following relations

$\begin{array}{l}1)\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}{r}_x[\left(B^T{}^{+}\left(u_{i+1/2}^L\right)u_{i+1/2}^L,u\right)-\left(B^T{}^{+}\left(u_{i-1/2}^L\right)u_{i-1/2}^L,u_{i-1}\right)}}}\\ +\left(B^{T-}\left(u_{i+1/2}^R\right)u_{i+1/2}^R,u_{i+1}\right)-\left(B^T{}^{-}\left(u_{i-1/2}^R\right)u_{i-1/2}^R,u\right)]\\ =\Delta x\Delta y\Delta z{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}{r}_x[\left(B^T{}^{+}\left(u_{N-1/2}^L\right)u_{N-1/2}^L,u_{N-1}\right)-\left(B^T{}^{+}\left(u_{1/2}^L\right)u_{1/2}^L,u_0\right)}}\\ +\left(B^{T-}\left(u_{N-1/2}^R\right)u_{N-1/2}^R,u_N\right)-\left(B^T{}^{-}\left(u_{1/2}^R\right)u_{1/2}^R,u_1\right)]=0\end{array}$

$\begin{array}{l}2)\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}{r}_y[\left(C^T{}^{+}\left(u_{j+1/2}^L\right)u_{j+1/2}^L,u\right)-\left(C^T{}^{+}\left(u_{j-1/2}^L\right)u_{j-1/2}^L,u_{j-1}\right)}}}\\ +\left(C^T{}^{-}\left(u_{j+1/2}^R\right)u_{j+1/2}^R,u_{j+1}\right)-\left(C^T{}^{-}\left(u_{j-1/2}^R\right)u_{j-1/2}^R,u\right)]=0;\end{array}$

$\begin{array}{l}3)\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}{r}_z[\left(D^T{}^{+}\left(u_{k+1/2}^L\right)u_{k+1/2}^L,u\right)-\left(D^T{}^{+}\left(u_{k-1/2}^L\right)u_{k-1/2}^L,u_{k-1}\right)}}}\\ +\left(D^T{}^{-}\left(u_{k+1/2}^R\right)u_{k+1/2}^R,u_{k+1}\right)-\left(D^T{}^{-}\left(u_{k-1/2}^R\right)u_{k-1/2}^R,u\right)]=0;\end{array}$

we have

$\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}\left(u^{m+1},u^{m+1}\right)}}}=\Delta x\Delta y\Delta z{\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\displaystyle \underset{k=-\infty }{\overset{+\infty }{\sum }}\left(u,u\right)}}},m=0,\cdots ,M-1.$

Hence

$I_{m+1}=I_m,m=0,\cdots ,M-1,$ (14)

which means the stability of the difference scheme (10)-(12) in the energetic norm $\sqrt{I_m}$ is established.

4. Numerical Example

Example 1. As an example, we consider the Burger’s hyperbolic equation given by:

$u_t+u{u}_x=0$ .

We introduce the following notations

$\begin{array}{l}f\left(u\right)=\frac{u^2}{2},f\left(u\right)=f^{+}\left(u\right)+f^{-}\left(u\right),\frac{\text{d}{f}^{+}}{\text{d}u}\ge 0,\frac{\text{d}{f}^{-}}{\text{d}u}\le 0,\left|u\right|\underset{\left|x\right|\to +\infty }{\to }0,u\in R,\\ u_i^{m+1}=u_i^m-\frac{\Delta t}{\Delta x}\left[f^{+}\left(u_{i+1/2}^L\right)-f^{+}\left(u_{i-1/2}^L\right)+f^{-}\left(u_{i+1/2}^R\right)-f^{-}\left(u_{i-1/2}^R\right)\right]\end{array}$

Consider the following reconstruction:

$\begin{array}{l}{u}_{i+1/2}^L=u_i^m+\frac{1}{2}\psi \left(R_i^m\right)\left(u_i^m-u_{i-1}^m\right),\\ u_{i-1/2}^R=u_i^m-\frac{1}{2}\psi \left(\frac{1}{R_i^m}\right)\left(u_{i+1}^m-u_i^m\right),R_i^m=\frac{u_{i+1}^m-u_i^m}{u_i^m-u_{i-1}^m},\end{array}$

where $u_i^k\approx \frac{1}{\Delta x}{\displaystyle {\int }_{x_{i-1/2}}^{x_{i+1/2}}u\left(\xi ,k\Delta t\right)\text{d}\xi }$ .

When $\psi =0$ corresponds to the scheme of the first order, $\psi =1$ is the one sided upwind scheme of second order. Consider the following modification of the obtained scheme

$\begin{array}{l}{u}^{m+1}-u+\frac{4}{3\left(u^{m+1}+u\right)}r\times [\left(f^{+}\left(u_{i-1/2}^L\right)\cdot \left(u_i-u_{i-1}\right)\right)+\left(\left[f^{+}\left(u_{i+1/2}^L\right)-f^{+}\left(u_{i-1/2}^L\right)\right]\cdot u\right)\\ +\left(f^{-}\left(u_{i+1/2}^R\right)\cdot \left(u_{i+1}-u_i\right)\right)+\left(\left[f^{-}\left(u_{i+1/2}^R\right)-f^{-}\left(u_{i-1/2}^R\right)\right]\cdot u\right)]=0,r=\Delta t/\Delta x.\end{array}$ (15)

We now prove that the difference scheme (15) admits the availability of difference analogue of the energy integral. Multiplying the system (15) by $\left[u^{m+1}+u\right]$ :

$\begin{array}{l}\left[u^{m+1}-u\right]\cdot \left[u^{m+1}+u\right]+\frac{4}{3}r\cdot [\left(f^{+}\left(u_{i-1/2}^L\right)\cdot \left(u-u_{i-1}\right)\right)+\left(\left[f^{+}\left(u_{i+1/2}^L\right)-f^{+}\left(u_{i-1/2}^L\right)\right]\cdot u\right)\\ +\left(f^{-}\left(u_{i+1/2}^R\right)\cdot \left(u_{i+1}-u\right)\right)+\left(\left[f^{-}\left(u_{i+1/2}^R\right)-f^{-}\left(u_{i-1/2}^R\right)\right]\cdot u\right)]=0.\end{array}$

Using formulas of difference differentiation, we obtain the following identity

$\begin{array}{l}[\left(f^{+}\left(u_{i-1/2}^L\right)\cdot \left(u-u_{i-1}\right)\right)+\left(\left[f^{+}\left(u_{i+1/2}^L\right)-f^{+}\left(u_{i-1/2}^L\right)\right]\cdot u\right)\\ =f^{+}\left(u_{i+1/2}^L\right)u-f^{+}\left(u_{i-1/2}^L\right)u_{i-1},\\ \left(f^{-}\left(u_{i+1/2}^R\right)\cdot \left(u_{i+1}-u\right)\right)+\left(\left[f^{-}\left(u_{i+1/2}^R\right)-f^{-}\left(u_{i-1/2}^R\right)\right]\cdot u\right)\\ =f^{-}\left(u_{i+1/2}^R\right)u_{i+1}-f^{-}\left(u_{i-1/2}^R\right)u.\end{array}$

Taking into account all of these transformations, we have

$\begin{array}{l}{\left(u^{m+1}\right)}^2-{\left(u\right)}^2+\frac{4}{3}r\cdot \{f^{+}\left(u_{i+1/2}^L\right)u-f^{+}\left(u_{i-1/2}^L\right)u_{i-1}\\ +f^{-}\left(u_{i+1/2}^R\right)u_{i+1}-f^{-}\left(u_{i-1/2}^R\right)u\}=0.\end{array}$ (16)

Multiply both sides of the (16) by $\Delta x$ and sum up over i from $-\infty$ to $+\infty$ , and noting that the function u tends to zero at infinity and denoting the quantity $\Delta x{\displaystyle \underset{i=-\infty }{\overset{\infty }{\sum }}{\left(u\right)}^2}$ by $I_m$ we obtain the equality

$I_{m+1}-I_m=0$ .

From this it is easily follows the energetic estimate

$I_m=I_0\text{,}m=1,\cdots ,M$ (17)

which means stability of the difference scheme in the norm $\sqrt{I_m}$ .

Example 2. Consider the following Cauchy problem

$\begin{array}{l}{u}_t+{\left(u^2/2\right)}_x=0,-\infty <x<\infty ,t>0,\\ \left|u\right|\underset{\left|x\right|\to +\infty }{\to }0,\\ u\left(x,0\right)=\{\begin{array}{ll}1,\hfill & x<0,\hfill \\ 0,\hfill & x>0.\hfill \end{array}\end{array}$ (18)

Rewrite the difference scheme (15) in the following form

$\begin{array}{c}{u}^{m+1}=u+\frac{2}{3\left(u^{m+1}+u\right)}r\times [\left(f^{+}\left(u_{i-1/2}^L\right)\left(u-u_{i-1}\right)\right)+\left(\left[f^{+}\left(u_{i+1/2}^L\right)-f^{+}\left(u_{i-1/2}^L\right)\right]u\right)\\ +\left(f^{-}\left(u_{i+1/2}^R\right)\left(u_{i+1}-u\right)\right)+\left(\left[f^{-}\left(u_{i+1/2}^R\right)-f^{-}\left(u_{i-1/2}^R\right)\right]u\right)].\end{array}$ (19)

In solving the difference scheme (19) we apply the iteration method with respect to the nonlinear coefficients

$\begin{array}{c}{u}^{\left(s+1\right)}=u+\frac{2}{3\left(u^{(s)}+u\right)}r\times [\left(f^{+}\left(u_{i-1/2}^L\right)\left(u_i-u_{i-1}\right)\right)+\left(\left[f^{+}\left(u_{i+1/2}^L\right)-f^{+}\left(u_{i-1/2}^L\right)\right]u\right)\\ +\left(f^{-}\left(u_{i+1/2}^R\right)\left(u_{i+1}-u_i\right)\right)+\left(\left[f^{-}\left(u_{i+1/2}^R\right)-f^{-}\left(u_{i-1/2}^R\right)\right]u\right)].\end{array}$

We have carried out a posteriori error analysis of the proposed scheme. Table 1 presents the exact and numerical solutions at points $\left(t_m,x_i\right)$ 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.

5. Conclusion

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.

Acknowledgements

This work was supported by Research Management Center (RMC), Universiti Sains Islam Malaysia (USIM) under Research Grand PPP/USG-0216/FST/30/15316.

Fund

This research was supported in part under R. A. Welch Foundation Grant E-0608.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Vorozhtsov, E.V. (2016) Construction of Third-Order Schemes Using Lagrange-Burmann Expansions for the Numerical Integration of Inviscid Gas Equations. Computational Methods and Programming, 17, 21-43. (In Russian)
[2]	Samarskii, A.A. and Popov, Y.P. (1980) Difference Schemes of Gas Dynamics Nauka, Moscow. (In Russian)
[3]	Tolstykh, A.I. (1990) Compact Difference Schemes and Their Application to Aero-Hydrodynamic Problems. Nauka, Moscow. (In Russian)
[4]	LeVeque, R.J. (1994) Numerical Methods for Conservation Laws. Birkh?user, Basel.
[5]	Tolstykh, A.I. (1994) High Accuracy Non-Centered Compact Difference Schemes for Fluid Dynamics Applications. World Scientific Publishing, Singapore. https://doi.org/10.1142/2269
[6]	Pinchukov, V.I. and Shu, C.W. (2000) High Order Numerical Methods for the Problems of Aerodynamics. Doklady Akademii Nauk, Novosibirsk. (In Russian)
[7]	Toro, E.F. (2000) Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Springer, Berlin.
[8]	Volkov, K.N., Deryugin, Y.N., Emel’yanov, V.N., Kozelkov, A.S. and Teterina, I.V. (2014) Difference Schemes in Gas Dynamics on Unstructured Grids. Fizmatlit, Moscow. (In Russian)
[9]	Boscheri, W., Balsara, D.S. and Dumbser, M. (2014) Lagrangian ADER-WENO Finite Volume Schemes on Unstructured Triangular Meshes Based on Genuinely Multidimensional HLL Riemann Solvers. Journal of Computational Physics, 267, 112-138. https://doi.org/10.1016/j.jcp.2014.02.023
[10]	Richtmyer, R. and Morton, K. (1972) Difference Methods for Initial Value Problems. Wiley, New York.
[11]	Rozhdestvenskii, B.L. and Yanenko, N.N. (1983) Systems of Quasilinear Equations and Their Application to Gas Dynamics. Nauka, Moscow.
[12]	Popov, I.V. and Fryazinov, I.V. (2014) Method of Adaptive Artificial Viscosity for Solving the Gas Dynamics Equations. Moscow. (In Russian).
[13]	LeVeque, R.J. (2004) Finite-Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge.
[14]	Vorozhtsov, E.V. (2000) Exercises for the Theory of Difference Schemes. Novosibirsk State University Press, Novosibirsk. (In Russian)
[15]	Lax, P.D. and Wendroff, B. (1960) Systems of Conservation Laws III. Communications on Pure and Applied Mathematics, 13, 217-237. https://doi.org/10.1002/cpa.3160130205
[16]	MacCormack, R. (2003) The Effect of Viscosity in Hypervelocity Impact Cratering. Journal of Spacecraft and Rockets, 40, 757-763.
[17]	Rusanov, V.V. (1968) Difference Schemes of the Third-Order Accuracy for Continuous Computation of Discontinuous Solutions. Doklady Akademii Nauk SSS, 180, 1303-1305.
[18]	Burstein, S.Z. and Mirin, A.A. (1970) Third Order Difference Methods for Hyperbolic Equations. Journal of Computational Physics, 5, 547-571. https://doi.org/10.1016/0021-9991(70)90080-X
[19]	Balakin, V.B. (1970) Methods of the Runge-Kutta Type for Gas Dynamics. USSR Computational Mathematics and Mathematical Physics, 10, 208-216. https://doi.org/10.1016/0041-5553(70)90192-8
[20]	Warming, R.F., Kutler, P. and Lomax, H. (1973) Second-and Third-Order Noncentered Difference Schemes for Nonlinear Hyperbolic Equations. AIAA Journal, 11, 189-196. https://doi.org/10.2514/3.50449
[21]	Harten, A. (1983) High Resolution Schemes for Hyperbolic Conservation Laws. Journal of Computational Physics, 49, 357-393. https://doi.org/10.1016/0021-9991(83)90136-5
[22]	Bona, C., Bona-Casas, C., and Terradas, J. (2009) Linear High-Resolution Schemes for Hyperbolic Conservation Laws: TVB Numerical Evidence. Journal of Computational Physics, 228, 2266-2281. https://doi.org/10.1016/j.jcp.2008.12.010
[23]	Liska, A. and Wendroff, B. (2003) Comparison of Several Difference Schemes on 1D and 2D Test Problems for the Euler Equations. SIAM Journal on Scientific Computing, 25, 995-1017. https://doi.org/10.1137/S1064827502402120
[24]	Roe, P.L. (1981) Approximate Riemann Solvers, Parameter Vectors and Difference Schemes. Journal of Computational Physics, 43, 357-372.
[25]	Chakravarthy, S. and Osher, S. (1985) A New Class of High Resolution TVD Schemes for Hyperbolic Conservation Laws. 23rd Aerospace Sciences Meeting, Reno, 14-17 January 1985. https://doi.org/10.2514/6.1985-363
[26]	Godunov, S.K. (1979) Equations of Mathematical Physics. Nauka, Moscow, 372.
[27]	Blokhin, A.M. and Aloev, R.D. (1993) Energy Integrals and Their Applications to the Study of the Stability of the Difference Schemes. Novosibirsk State University Press, Novosibirsk, 224.
[28]	Aloev, R.D., Eshkuvatov, Z.K., Davlatov, S.O. and Nik Long, N.M.A. (2014) Sufficient Condition of Stability of Finite Element Method for Symmetric T-Hyperbolic Systems with Constant Coefficients. Computers and Mathematics with Applications, 68, 1194-1204. https://doi.org/10.1016/j.camwa.2014.08.019
[29]	Aloev, R.D., Blokhin, A.M. and Hudayberganov, M.U. (2014) One Class of Stable Difference Schemes for Hyperbolic System. American Journal of Numerical Analysis, 2, 85-89.
[30]	Aloev, R.D., Davlatov, S.O., Eshkuvatov, Z.K. and Nik Long, N.M.A. (2016) Uniqueness Solution of the Finite Elements Scheme for Symmetric Hyperbolic Systems with Variable Coefficients. Malaysian Journal of Mathematical Sciences, 10, 49-60.