On conditions for weak conservativeness of regularized explicit finite-difference schemes for 1D barotropic gas dynamics equations
A. Zlotnik T. Lomonosov


On conditions for weak conservativeness of regularized explicit finite-difference schemes for 1D barotropic gas dynamics equations
On conditions for weak conservativeness of explicit finite-difference schemes

A. Zlotnik National Research University Higher School of Economics, Myasnitskaya 20, 101000 Moscow, Russia azlotnik@hse.ru
T. Lomonosov National Research University Higher School of Economics, Myasnitskaya 20, 101000 Moscow, Russia tlomonosov@hse.ru

*
    A. Zlotnik and T. Lomonosov
    December 30, 2023
===============================



We consider explicit two-level three-point in space finite-difference sche­mes for solving 1D barotropic gas dynamics equations.
The schemes are based on special quasi-gasdynamic and quasi-hydrodynamic regularizations of the system.
We linearize the schemes on a constant solution and
derive the von Neumann type necessary condition and a CFL type criterion (necessary and sufficient condition) for weak conservativeness in L^2 for the corresponding initial-value problem on the whole line.
The criterion is essentially narrower than the necessary condition and wider than a sufficient one obtained recently in a particular case; moreover, it corresponds most well to numerical results for the original gas dynamics system.

Keywords: gas dynamics, barotropic quasi-gas dynamics system of equations, explicit finite-difference schemes, stability criterion, weak conservativeness



§ INTRODUCTION


The stability theory for finite-difference schemes for model problems in gas dynamics is well presented in the literature <cit.>.
In this paper we consider some finite-difference schemes for solving 1D barotropic gas dynamics equations.
The schemes are explicit, two-level in time and use a symmetric three-point stencil in space.
Their construction is based on special quasi-gasdynamic and quasi-hydrodynamic <cit.> regularizations of the original equations (without a regularization, the schemes are unstable).
The schemes of this kind were successfully applied in numerous and various practical applications, in particular, see
<cit.>, but their theory is not developed so well.
We linearize the schemes on a constant solution and derive both the von Neumann type necessary condition and a CFL type criterion for weak conservativeness in L^2 for the corresponding initial-value problem on the whole line.
The weak conservativeness in L^2 means the uniform in time bound for the norm of scaled solution by the norm of initial data instead of the energy conservation law for the linearized original system, i.e., the acoustics system of equations.
Our numerical experience show that validity of the weak conservativeness property is important since it prevents numerical solutions from the well-known possible spurious oscillations.
The property guarantees the uniform in time L^2 stability with respect to initial data.
Since in practice necessary conditions are often in use (a derivation of sufficient conditions is much more complicated in general), it is important to know to what extent this is lawful to do.
The criterion turns out to be essentially narrower than the necessary condition and at the same time wider than a sufficient condition obtained recently in a particular case in <cit.>.
Moreover, namely the criterion corresponds most well to results of numerical experiments for the original gas dynamics system.
Therefore the criterion (but not the necessary condition or sufficient one) is most adequate and useful for practical purposes.



§ SYSTEMS OF EQUATIONS, FINITE-DIFFERENCE SCHEMES AND THEIR LINEARIZATION


The 1D barotropic gas dynamics (Euler) system of equations consists in the mass and momentum balance equations

    ∂_tρ+∂_x(ρ u) = 0, ∂_t(ρ u)+∂_xp(ρ)=0,

where ρ>0, u and p are the gas density and velocity (the sought functions) and pressure.
We assume that p'(ρ)>0 and consider the equations for x∈ℝ and t>0.
The 1D barotropic quasi-gas dynamics (QGD) system of equations consists in the regularized mass and momentum balance equations

    ∂_tρ+∂_xj = 0, ∂_t(ρ u)+∂_x(ju+p(ρ)-Π)=0,
    
    
    j=ρ(u-w), w = τ/ρu∂_x(ρ u)+ŵ, ŵ = τ/ρ[ρ u∂_xu + p'(ρ)],
    
    Π = Π_NS + ρ uŵ + τ p'(ρ)∂_x(ρ u),  Π_NS=μ(ρ)∂_xu.

Here j and Π are the regularized mass flux and stress,
w and ŵ are the regularizing velocities,
τ=τ(ρ)>0 is a regularization parameter
and Π_NS is the Navier-Stokes viscous stress with μ(ρ)≥ 0 being proportional to the viscosity coefficient.
In the barotropic case, quasi-gasdynamic and quasi-hydrodynamic systems were introduced and investigated (in multidimensional case) in <cit.>.

The QGD system is simplified into the original system (<ref>) for τ=μ=0 and the Navier-Stokes system of equations for viscous compressible barotropic gas flow for τ=0 and μ>0.

System (<ref>) can be linearized on a constant solution ρ_*≡const>0 and u_*=0.
Substituting the solution in the form ρ=ρ_*+Δρ and u=u_*+Δ u in the equations and neglecting the terms having the second order of smallness with respect to Δρ and Δ u and their derivatives leads us to the following system of equations:

    ∂_tΔρ+ρ_*∂_xΔ u=0, ρ_*∂_tΔ u+p'(ρ_*)∂_xΔρ u=0.

For the dimensionless unknowns ρ̃=Δρ/ρ_* and ũ=Δ u/√(p'(ρ_*)) we gain the acoustics system of equations:

    ∂_tρ̃+c_*∂_xũ=0, ∂_tũ+c_*∂_xρ̃=0.

Hereafter c_*=√(p'(ρ_*)) is the background velocity of sound.
Given the initial data ρ̃|_t=0=ρ̃_0 and ũ|_t=0=ũ_0 (that one can consider complex-valued), for the solution to the last system the following energy conservation law holds

    ρ̃(·,t)_L^2(ℝ)^2+ũ(·,t)_L^2(ℝ)^2
     =ρ̃_0_L^2(ℝ)^2+ũ_0_L^2(ℝ)^2  for  t≥ 0.


Now we pass to discretization.
Let ω_h be a uniform mesh on ℝ with the nodes x_k=kh, k∈ℤ, and step h=X/N.
Let ω^*_h be an auxiliary mesh with the nodes x_k+1/2 = (k+0.5)h, k∈ℤ.
Define a uniform mesh in t with the nodes t_m=mΔ t, m≥ 0, and step Δ t>0.
We define the shift, averaging and difference quotient operators

    v_±,k=v_k± 1, (sv)_k-1/2=v_k+v_k+1/2,  (δ v)_k-1/2=v_k-v_k-1/h,
    (δ^∗ y)_k = y_k+1/2-y_k-1/2/h,  δ_tv=v^+-v/Δ t,  v^+,m=v^m+1.


We first consider a standard explicit two-level in time and three-point symmetric in space discretization of the QGD equations (<ref>)-(<ref>):

    δ_t ρ + δ^∗ j = 0, δ_t(ρ u) + δ^∗(jsu +p(sρ) - Π)=0  on  ω_h,
    
    
     j =(sρ)su-(sρ)w, (sρ)w =(sτ)[δ(ρ u)]su + (sρ)ŵ,
    
    
     (sρ)ŵ = (sτ)[(sρ)(su)δ u + δ p(ρ)],
    
    Π = μδ u + (su)(sρ)ŵ + (sτ)[p'(sρ)]δ(ρ u).

The main unknown functions ρ>0, u and the parameter τ are defined on ω_h whereas j,w,ŵ,Π and μ are defined on ω_h^*.

In <cit.> two non-standard spatial discretizations of the QGD equations (<ref>)-(<ref>) were constructed which are weakly conservative in energy (see their generalization to a multidimensional case in <cit.>).
One of them has the “enthalpy” form

    δ_t ρ + δ^∗ j = 0, δ_t(ρ u) + δ^∗(jsu-Π) + s^∗[(sρ)δh(ρ)]=0,
    
    
     j =sρ· su-sρ· w,  sρ· w = [(τ∂_x)_h(ρ u)]su + (sρ)ŵ, 
    ŵ = (sτ)[(su)δ u + δh(ρ)], Π =μδ u + (su)(sρ)ŵ + p'(sρ)(τ∂_x)_h(ρ u),
    
    
     (τ∂_x)_h(ρ u) = (sτ/h'(ρ)){[δh(ρ)]su + p'(sρ)δ u},

where h(ρ)=∫_r_0^ρp'(r)/r dr, with some r_0>0, is the gas enthalpy and thus h'(ρ)=p'(ρ)/ρ.
In the isentropic case p(ρ)=p_1ρ^γ with γ>1, one can take r_0=0 and then h(ρ)=γ/γ-1p(ρ)/ρ and h'(ρ)=γp(ρ)/ρ^2.
Notice the non-standard h(ρ)-dependent discretizations of ∂_xp(ρ) in (<ref>) and (<ref>) and τ∂_x(ρ u) in (<ref>)-(<ref>), see (<ref>).

We linearize scheme (<ref>)-(<ref>) on a constant solution ρ_*≡const>0 and u_*=0.
To do that, we write its solution in the form ρ=ρ_*+Δρ and u=u_*+Δ u, neglect terms having the second order of smallness with respect to Δρ and Δ u and obtain

    δ_tΔρ+ρ_*δ^*sΔ u-τ(ρ_*)p'(ρ_*)δ^*δΔρ=0,
    
    ρ_*δ_tΔ u+p'(ρ_*)δ^*sΔρ-[μ(ρ_*)+τ(ρ_*)ρ_*p'(ρ_*)]δ^*δΔ u=0.

For the dimensionless unknowns ρ̃=Δρ/ρ_* and ũ=Δ u/c_* we get equations

    δ_tρ̃+c_*δ^*sũ-τ(ρ_*)c_*^2δ^*δρ̃=0,
    
    δ_tũ+c_*δ^*sρ̃
     -[μ(ρ_*)/ρ_*+τ(ρ_*)c_*^2]δ^*δũ=0

(cf. systems (<ref>) and (<ref>)).
The linearization of scheme (<ref>)-(<ref>) is the same.
Notice that since u_*=0, the linearization result remains the same if it would be w=ŵ and the terms dependent on u were omitted in the definition of these variables, i.e., for example, w=ŵ=τ̃/sρδ p(ρ) instead of formulas in (<ref>)-(<ref>).

We assume that the regularization parameter and viscosity coefficient are given by usual QGD-formulas

    τ(ρ)= α h/√(p'(ρ)), μ(ρ) = α_sτ(ρ)ρ p'(ρ),

where α>0 and α_s≥ 0 are parameters.
Then omitting tildes above ρ and u, equations (<ref>)-(<ref>) can be rewritten in the following recurrent form

    ρ^+=ρ-β/2(u_+-u_-)+αβ(ρ_–2ρ+ρ_-),
    
    
     u^+=u-β/2(ρ_+-ρ_-)+ϰαβ(u_+-2u+ u_-)

with three parameters α, β:=c_*Δ t/h and ϰ:=α_s+1≥ 1.
The functions ρ^0 and u^0 are given, i.e., we consider the initial-value problem for the scheme.
Below it is convenient to consider ρ and u as complex-valued mesh functions.



§ WEAK CONSERVATIVENESS ANALYSIS


Let 𝐲^m=(ρ^m u^m)^T, m≥ 0, be a column-vector function on ω_h
and the linearized difference scheme (<ref>)-(<ref>) be rewritten in a matrix form

    𝐲^+
    = [  αβ β/2; β/2 ϰαβ ]𝐲_-
    +[  1-2αβ      0;      0 1-2ϰαβ ]𝐲
    +
    [   αβ -β/2; -β/2  ϰαβ ]𝐲_+.

Let H be a Hilbert space of complex valued square-summable on ω_h vector functions, i.e. having a finite norm

    𝐲_H=(h∑_k=-∞^∞|𝐲_k|^2)^1/2.

For 𝐲^0=(ρ^0 u^0)^T∈ H we have that 𝐲^m∈ H for all m≥ 1.
We define a weak conservativeness of scheme (<ref>) as validity of the bound

    sup_m≥ 0𝐲^m_H≤𝐲^0_H  ∀𝐲^0∈ H.

This definition is motivated by the energy conservation law (<ref>) for the acoustics system of equation (<ref>).
It is essential to notice that for the linearized QGD-system (<ref>)-(<ref>) namely the corresponding inequality holds in place of equality (<ref>) so that it is natural to study the bound for schemes based on such a system.
Of course, estimate (<ref>) guarantees the uniform in time stability in H with respect to initial data.
We first substitute a partial solution in the form 𝐲_k^m=e^𝐢kξ𝐯^m(ξ), k∈ℤ, m≥ 0, where 𝐢 is the imaginary unit and 0≤ξ≤ 2π is a parameter, into (<ref>) and obtain

    𝐯^+(ξ)=G(ξ)𝐯(ξ), G(ξ)
     =[  1-ω_1  -𝐢ω_2;  -𝐢ω_2 1-ϰω_1 ],

where we denote
ω_1=4αβθ, θ=sin^2ξ/2∈[0,1] and ω_2=βsinξ for brevity.
Below it is important that ω_2^2=4β^2θ(1-θ).

It is known  (see similar formulas in <cit.>) that if 𝐲^0=(ρ^0 u^0)^T∈ H, then there exists a function
𝐯^0∈ L^2(0,2π) such that

    𝐯^0(ξ)=1/√(2π)∑_k=-∞^∞𝐲_k^0 e^-𝐢kξ,

and we can write the solution to scheme (<ref>) in an integral form

    𝐲_k^m=1/√(2π)∫_0^2π𝐯^m(ξ)e^𝐢kξ dξ,  k∈ℤ,

where 𝐯^m∈ L^2(0,2π) due to (<ref>).
The following Parseval identity also holds

    𝐲^m_H=√(h) 𝐯^m_L^2(0,2π),  m≥ 0.

The von Neumann type spectral condition

    max_0≤ξ≤ 2πmax_l|λ_l(G(ξ))|≤ 1

is known to be a necessary condition for property (<ref>) to hold (see similar result in  <cit.>).
Hereafter λ_l(A) are eigenvalues of a matrix A.
Let us determine the spectral form of property (<ref>).


Validity of the spectral bound

    max_0≤ξ≤ 2πmax_lλ_l((G^*G)(ξ))≤ 1

is necessary and sufficient for the weak conservativeness property (<ref>) to hold.


Due to the Parseval identity (<ref>) and formula (<ref>) we have

    h^-1ŷ_H^2=𝐯̂_L^2(0,2π)^2=G𝐯_L^2(0,2π)^2
     =(G^*G𝐯,𝐯)_L^2(0,2π).

Since (G^*G)(ξ)≥ 0 is a Hermitian matrix, it has a spectral decomposition (G^*G)(ξ)=U^*(ξ)Λ(ξ)U(ξ), where U(ξ) is a unitary matrix, and Λ(ξ) is a diagonal matrix with numbers
λ_l((G^*G)(ξ))≥ 0 forming its diagonal.
Hence for z(ξ):=U(ξ)𝐯(ξ) we have

    (G^*G𝐯,𝐯)_L^2(0,2π)=(Λ𝐳,𝐳)_L^2(0,2π)=Λ^1/2𝐳_L^2(0,2π)^2.

Thus 𝐲^m_H^2=hΛ^m/2𝐳^0_L^2(0,2π)^2 for m≥ 0, and bound (<ref>) is equivalent to the following one

    sup_m≥ 0Λ^m/2𝐳^0_L^2(0,2π)^2≤𝐳^0_L^2(0,2π)^2  ∀𝐳^0∈ L^2(0,2π).

It holds if and only if the spectral bound (<ref>) holds.


Under validity of the spectral bound (<ref>), the norm 𝐲^m_H is actually non-increasing in m≥ 0 that serves as a stronger property than (<ref>).

In the proof of this lemma, the specific form and dimension of the matrix G are clearly inessential, and actually it holds in general case.
In our case the matrix G^*G has the form

    G^*G=
    [    (1-ω_1)^2+ω_2^2      -𝐢(1-ϰ)ω_1ω_2;       𝐢(1-ϰ)ω_1ω_2 (1-ϰω_1)^2 + ω_2^2 ].

Note that G^∗G=[(1-ω_1)^2+ω_2^2]I in the simplest case ϰ=1, where I is a unit matrix.



The necessary spectral condition (<ref>) holds if and only if

    β≤min{(ϰ+1)α,1/2ϰα}.



The characteristic polynomial for the matrix G has the following form

    q_1(λ)=λ^2-( G)λ+ G
     =λ^2+ [(ϰ+1)ω_1-2]λ+
    +[ϰω_1^2+ω_2^2+1-(ϰ + 1)ω_1].

We set a_0:=q_1(1)=1-( G)+ G=ϰω_1^2+ω_2^2≥ 0 and notice that a_0=q_1(1)=0 if and only if G=I.
We transform the unit circle {|λ|≤ 1} with a punctured point (1,0) on ℂ into the closed left half-plane  { z≤ 0} and put

    q̂_1(z):=(z-1)^2q_1(z+1/z-1)
     =a_0z^2+2a_1z+a_2,

where a_1=1- G, a_2=1+ G+ G.
It is well known that for a_0>0 the roots q̂_1(z) lie in { z≤ 0} under the conditions a_1≥ 0 and a_2≥ 0, i.e.

    ϰω_1^2+ω_2^2-(ϰ+1)ω_1≤ 0,
     ϰω_1^2+ω_2^2-2(ϰ+1)ω_1+4≥ 0.

We rewrite these conditions as

    β(4ϰα^2θ+1-θ)-(ϰ+1)α≤ 0  for  0≤θ≤ 1,
    
    
     r(θ):=β^2(4ϰα^2-1)θ^2-β(2(ϰ+1)α-β)θ+1≥ 0  for  0≤θ≤ 1.

The left-hand side of (<ref>) is linear in θ, thus it suffices to test it for θ=0,1 that leads us to the condition

    β≤min{(ϰ+1)α,ϰ+1/4ϰα}.

Next we analyze condition (<ref>).
Notice that r(0)=1 and due to (<ref>) we have 2(ϰ+1)α-β>0.
For a=4ϰα^2-1≠ 0 the vertex of the parabola r(θ) is given by

    θ_v=2(ϰ+1)α-β/2β(4ϰα^2-1).

For 4ϰα^2-1>0 the property θ_v>1 means that
β<(ϰ+1)α/4ϰα^2-0.5,
and it holds due to (<ref>).
Hence condition (<ref>) reduces to r(1)≥ 0, i.e.

    4ϰα^2β^2-2(ϰ+1)αβ+1
     =4ϰα^2(β-1/2ϰα)(β-1/2α)≥ 0.

For 4ϰα^2-1<0 we have θ_v<0, so that condition (<ref>) reduces again to r(1)≥ 0.
For  4ϰα^2-1=0 the condition also reduces to r(1)≥ 0 (since r(0)=1).
Since ϰ≥ 1, inequality (<ref>) means that either of the conditions

    β≤1/2ϰα,  β≥1/2α

holds.
Combining them with (<ref>), we obtain (<ref>).


Now we turn to the spectral criterion (<ref>).


The spectral criterion (<ref>) holds if and only if

    β≤min{2α,1/2ϰα}.



The characteristic polynomial of G^*G has the following form

    q_2(λ)=λ^2-(G^*G)λ+( G)^2.

Since λ_l(G^*G)≥ 0,
the property |λ_l(G^*G)|≤ 1 means validity of the conditions

    1/2(G^*G)≤ 1,  q_2(1)=1-(G^*G)+( G)^2≥ 0.


The first of them has the form

    ϰ^2+1/2 ω_1^2+ω_2^2-(ϰ+1)ω_1≤ 0

and can be specified as

    8α^2β^2(ϰ^2+1)θ^2 + 4β^2θ(1-θ)-4αβ(ϰ + 1)θ≤ 0 for  0≤θ≤ 1.

After dividing by 4βθ we get that it suffices to confine ourselves with the values θ=0,1 that leads to the condition

    β≤min{(ϰ + 1)α,ϰ + 1/2(ϰ^2+1)α}.


In order to transform the second condition (<ref>) we notice that

    (G^*G)=2(b+1)+(ϰ-1)^2ω_1^2,  ( G)^2=(b+1)^2,

where b:=ϰω_1^2+ω_2^2-(ϰ+1)ω_1, see (<ref>).
Hence the following factorization holds

    q_2(1)=b^2-(ϰ-1)^2ω_1^2=(b-(ϰ-1)ω_1)(b+(ϰ-1)ω_1),

that is decisive for the simplicity of our analysis.
Since ϰ≥ 1, the condition q_2(1)≥ 0  is equivalent to validity of either of the conditions

    ϰω_1^2+ω_2^2-2ω_1≤ 0,  ϰω_1^2+ω_2^2-2ϰω_1≥ 0,

i.e., more specifically, to validity of either of the conditions

    β(4ϰα^2θ+1-θ)-2α≤ 0  for  0≤θ≤ 1,
    
    β(4ϰα^2θ+1-θ)-2ϰα≥ 0  for  0≤θ≤ 1.

As above they respectively mean that the inequalities

    β≤min{2α,1/2ϰα}, β≥max{2ϰα,1/2α}

hold. Combining them with  (<ref>) we complete the proof.


It is essential that the function on the right-hand side of condition (<ref>) reaches its maximal value at
α=α_∗:=1/2√(ϰ)≤1/2
and the maximal value equals 1/√(ϰ)≤ 1.
Hence the criterion coincides with the standard CFL stability condition β≤ 1 if and only if α=α_* and ϰ=1.
The criterion gives an important information on the optimal choice of α since in practice for the original non-linear problem α is normally sought experimentally.
Note also that criterion (<ref>) and the necessary condition (<ref>) coincide only in the case α≥α_*.
We call attention to a paradoxical moment: criterion (<ref>) becomes stronger as the coefficient of  “effective viscosity” ϰ increases (it is harder to say that about the necessary condition (<ref>)).
Therefore the best choice in the present bounds is α_s=0, i.e. ϰ=1.
But this conclusion is not universal in practice and it is known that in some situations α_s>0 has to be taken (see, for example, <cit.>).
In Fig. <ref> we compare the necessary condition, the criterion of stability and the sufficient condition as well as the results of numerical experiments for the original system (<ref>) for p(ρ)=ρ^2 (the scaled case of the shallow water equations)
and ϰ=7/3.
The sufficient condition was obtained in <cit.> only for these p(ρ) and ϰ by the energy method and has the form

    β≤min{2α/1+6α+4α^2,4α/1+6α+16α^2}.

Notice that the first fraction in it is less than the second one for 0<α<3+√(17)/8≈ 0.890.
The corresponding graph is almost flat for 0.3≤α≤ 0.9 in contrast to the cases of necessary condition and criterion.
The computations are accomplished for 0≤ t≤ 0.5 for the Riemann problem with the discontinuous initial data

    ρ_0(x) =
    
     1,     x < 0
    
     0.1,     x > 0
    , u_0(x) =
    
     0.1,     x < 0,
    
     0,     x > 0

for both schemes  (<ref>)-(<ref>) and (<ref>)-(<ref>) with h=1/125.



We observe a good correspondence of the obtained criterion with the experimental results, and that the sufficient condition underestimates the criterion up to several times in the most interesting region α≈α_*.
Also the results for the “enthalpy” scheme (<ref>)-(<ref>) are clearly different from and better than for the standard one (<ref>)-(<ref>) though the above linearized analysis gives the same results for them.

We have identified non-conservative computations by noticeable well-known oscillations of the numerical solutions (some of  computations have not even been completed due to overflow).
In Fig. <ref> we give an example of conservative and non-conservative solutions ρ and u
for the “enthalpy” scheme (at time t=0.5) for α =0.4 and two neighboring values of β from Fig. <ref> (b).





§ THE CASE OF THE SCHEMES BASED ON A SIMPLIFIED REGULARIZATION


We also consider a simplified (quasi-hydrodynamic <cit.>) regularization (<ref>)-(<ref>), where the terms with ∂_x(ρ u) are omitted, in particularly, it becomes w=ŵ.
Correspondingly in schemes (<ref>)-(<ref>) and (<ref>)-(<ref>) we have to omit both terms with respectively δ(ρ u) and (τ∂_x)_h(ρ u).
In the linearized scheme the term τ(ρ_*)c_*^2 disappears from equation (<ref>), hence now ϰ=α_s.
Notice that usually 0<α_s≤ 1 though in specific cases α_s>1 can be also taken.

For the simplified scheme based on the quasi-hydrodynamic regularization the following results are valid:
(1) in the case 0≤α_s≤ 1
the necessary condition (<ref>) and criterion (<ref>) hold if and only if respectively

    β≤min{(α_s+1)α,1/2α},
    
    β≤min{2α_sα,1/2α};


(2) in the case α_s≥ 1 the results of Theorems <ref> and <ref> remain valid with ϰ=α_s.


The above given analysis holds true except for some changes in the case 0≤ϰ=α_s≤ 1.
In this case, inequalities (<ref>) are replaced by the following ones

    β≤1/2α,  β≥1/2ϰα;

they being combined with (<ref>) lead to (<ref>).
Also inequalities (<ref>) are replaced by the following ones

    β≤min{2ϰα,1/2α}, β≥max{2α,1/2ϰα},

they being combined with (<ref>) lead to (<ref>).

The maximal value of the function on the right-hand side of criterion (<ref>) is reached at
α=α_*:=1/2√(α_s)≥1/2
and equals √(α_s)≤ 1.
We notice that the necessary condition (<ref>) is especially rough compared to criterion (<ref>) for α_s≈ 0, including the case α_s=0 when actually the stability is absent at all.


The study was partially supported by the RFBR, project nos. 16-01-00048 and 18-01-00587.


99
B73
Bakhvalov, N.S., Zhidkov, N.P., Kobelkov, G.M.: Numerical methods. Binom, Moscow, (2011) [in Russian]

BZS16
Balashov, V., Zlotnik, A., Savenkov, E.: Analysis of a regularized model for the isothermal two-component mixture with the diffuse interface. Russ. J. Numer. Anal. Math. Model. 32 (6) (2017), 347–358.
See also: Keldysh Inst. Appl. Math. Preprint 89 1–26 (2016).
http://library.keldysh.ru/preprint.asp?id=2016-89 [in Russian]

Ch04
Chetverushkin, B.N.:Kinetic schemes and quasi-gas dynamic system of equations, CIMNE, Barcelona (2008)

Cou13
Coulombel, J.-F.: Stability of finite difference schemes for hyperbolic initial boundary value problems. In: HCDTE Lecture Notes. Part I. Dispersive and Transport Equations, G. Alberti, F. Ancona, S. Bianchini et al., eds,
Amer. Inst. Math. Sci., Springfield, 97–226 (2013)

GV96
Ganzha, V.G., Vorozhtsov, E.V.: Computer-aided analysis of difference schemes for partial differential equations, Wiley, New York (1996)

GR73
Godunov, S.K., Ryabenkii, V.S.: Difference schemes. Studies in Mathematics and its Applications 19, North Holland, Amsterdam (1987)

GKO95
Gustafsson, B., Kreiss, H.-O., Oliger, J.: Time dependent problems and difference methods. Wiley, New York (1995)

E07
Elizarova, T.G.: Quasi-gas dynamic equations. Springer, Dordrecht (2009)

EB11
Elizarova, T.G., Bulatov, O.V.: Regularized shallow water equations and a new method of simulation of the open channel flows.
Comput. Fluids 46 206–211 (2011)

EZI17
Elizarova, T.G., Zlotnik, A.A, Istomina, M.A.: Hydrodynamic aspects of spiral-vortex structure formation in rotating gas discs.
Astron. Reports 62 (1) (2018), 9–18.
See also: Keldysh Inst. Appl. Math. Preprint 1 1–30 (2017).
http://library.keldysh.ru/preprint.asp?id=2017-1 [in Russian]

LeV04
LeVeque, R.J.: Finite volume methods for hyperbolic problems. Cambridge University Press,  Cambridge (2004)

RM72
Richtmyer, R.D., Morton, K.W.: Difference methods for initial-value problems. 2nd ed., Wiley-Interscience (1967)

Sh09
Sheretov, Yu.V.: Continuum dynamics under spatiotemporal averaging.  RKhD, Moscow-Izhevsk (2009) [in Russian])

SSh13
Suhomozgii, A.A., Sheretov, Yu.V.: Stability analysis of a finite-difference scheme for solving the Saint-Venant equations in the shallow water theory. In: Appl. of functional anal. in approx. theory, Tver State Univ. (2013) 48–60 [in Russian]

Z08
Zlotnik, A.A.: Parabolicity of a quasihydrodynamic system of equations and the stability of its small perturbations.
Math. Notes 83 (5) 610-623 (2008)

ZCMMP10
Zlotnik, A.A:
Energy equalities and estimates for barotropic quasi-gasdynamic and quasi-hydrodynamic systems of equations.
Comput. Math. Math. Phys. 50 (2) 310-321  (2010)

ZMM12a
Zlotnik, A.A:
On construction of quasi-gasdynamic systems of equations and the barotropic system with the potential body force.
Math. Modeling 24 (4) 65-79 (2012) [in Russian]

ZMM12b
Zlotnik, A.A.: Spatial discretization of the one-dimensional barotropic quasi-gasdynamic system of equations
and the energy balance equation, Math. Modeling 24 (10) 51–64 (2012) [in Russian].

ZCMMP16
Zlotnik, A.A.: On conservative spatial discretizations of the barotropic quasi-gasdynamic system of equations with a potential body force. Comput. Math. Math. Phys. 56 (2) 303–319 (2016)

ZCh08
Zlotnik, A.A., Chetverushkin, B.N.: Parabolicity of the quasi-gasdynamic system of equations,
its hyperbolic second-order modification, and the stability of small perturbations for them.
Comput. Math. Math. Phys. 48 (3) 420–446 (3) (2008)

ZG16
Zlotnik, A., Gavrilin, V.: On a conservative finite-difference method for 1D shallow water flows based on regularized equations. In: Math. Problems in Meteorological Model. Math. in Industry, A. Bátkai, P. Csomós, A. Horányi, etc., eds. 24 16–31 (2016)


