Spectral Stability Conditions for an Explicit Three-Level Finite-Difference Scheme for a Multidimensional Transport Equation with Perturbations

We study difference schemes associated with a simplified linearized multidimensional hyperbolic quasi-gasdynamic system of differential equations. It is shown that an explicit two-level vector difference scheme with flux relaxation for a second-order hyperbolic equation with variable coefficients that is a perturbation of the transport equation with a parameter multiplying the highest derivatives can be reduced to an explicit three-level difference scheme. In the case of constant coefficients, the spectral condition for the time-uniform stability of this explicit three-level difference scheme is analyzed, and both sufficient and necessary conditions for this condition to hold are derived, in particular, in the form of Courant type conditions on the ratio of temporal and spatial steps.


INTRODUCTION
The hyperbolic quasi-gasdynamic system is a specially perturbed system of gasdynamic equations in which the terms containing second derivatives with respect to the space and time variables are multiplied by a small parameter τ > 0. This system is used for constructing a family of three-level and two-level vector difference schemes for the numerical solution of various problems in gas dynamics [1]. Note that parabolic quasi-gasdynamic systems are presented in [2,3]. The papers [4,5] establish a number of mathematical properties of the quasi-gasdynamic system, including timeuniform estimates for the corresponding system linearized on a constant solution. A second-order hyperbolic perturbation with the parameter τ of the parabolic initial-boundary value problem without convective terms [6,7] and with them [8] was also studied, and the stability of the corresponding three-level weighted and two-level vector numerical methods was analyzed [9]. The time-uniform stability of implicit three-level weighted and two-level vector difference schemes for the linearized quasi-gasdynamic system was proved in [10].
In practice, explicit difference schemes for the quasi-gasdynamic system are of most interest. However, the linearizations of these schemes have not yet been proved to be stable uniformly in τ . So far, the energy method has not been successfully applied to prove this by analogy with [4,10], because the linearized schemes have a rather complex structure; in particular, they involve approximations to terms containing first and second derivatives with respect to the spatial and temporal variables, the second derivatives being multiplied by the parameter τ .
The present paper is intended to make some progress in this direction. To this end, instead of the quasi-gasdynamic system, we consider the simplified case of one second-order hyperbolic equation with variable coefficients that is a perturbation of the transport equation with the parameter τ multiplying the second derivatives. First, we show that an explicit two-level vector difference fluxrelaxation scheme of the type in [11] for this equation can be reduced to an explicit three-level difference scheme. In this three-level scheme, the parameter τ is multiplied by some factor greater than unity in front of the approximation for the second t-derivative. Second, and most importantly, we analyze the spectral condition for the time-uniform stability of this explicit three-level scheme in the case of constant coefficients. Both sufficient and close-to-them necessary conditions for the validity of the spectral condition are obtained, in particular, in the form of Courant type conditions on the ratio of temporal and spatial steps. These conditions are independent of τ .
The necessary conditions include a condition for the prevalence of the coefficients of viscous terms over the transport coefficients. The derivation is based on analyzing the location of the roots of the corresponding characteristic polynomial with complex coefficients on the complex plane using the generalized Routh-Hurwitz criterion [12,Ch. XV,Sec. 18]. Criteria for a similar spectral stability condition for the differential equation itself and a difference scheme with zero parameter τ are preliminarily analyzed as well. Such a criterion for this equation is a condition for a kind of prevalence of the matrix of viscous terms over the vector of transport coefficients; conditions of this kind played an essential role in [4,10]. For the scheme indicated, the criterion is a Courant type condition. The results are rather encouraging for the subsequent analysis of the corresponding properties of difference schemes for the quasi-gasdynamic system.

HYPERBOLIC SECOND-ORDER DIFFERENTIAL EQUATION AND REDUCTION OF AN EXPLICIT TWO-LEVEL VECTOR DIFFERENCE FLUX RELAXATION SCHEME FOR IT TO A THREE-LEVEL SCHEME
Consider the hyperbolic second-order differential equation with a parameter τ > 0 multiplying the highest derivatives. The sought function u = u(x, t) is defined and the equation is considered for x = (x 1 , . . . , x n ) ∈ R n , t ≥ 0, and n ≥ 1. Hereafter, . . , a n (x)} is a diagonal matrix, the operators div and ∇ are taken with respect to x, and summation from 1 to n over repeated indices i and j is assumed. Let us introduce the uniform mesh ω kh with the nodes x kl = lh k , l ∈ Z, and step h k > 0; the mesh shifted by h k /2 with the nodes x k(l−1/2) = (l − (1/2))h k , l ∈ Z; and the mesh operators where v l = v(x kl ) and w l−1/2 = w(x k(l−1/2) ), 1 ≤ k ≤ n. Let us also set h = (h 1 , . . . , h n ) and introduce the mesh ω h := ω 1h × . . . × ω nh and the operators s = (s 1 , . . . , s n ) and ∇ h = (δ 1 , . . . , δ n ). Let us introduce the uniform mesh ω ht with the nodes t m = mh t , m ≥ 0, and step h t > 0 and the difference operators where y m = y(t m ),y m = y m−1 , and y m = y m+1 . Let ω ht = ω ht \{0}. For Eq. (1), consider the following explicit vector difference scheme with flux relaxation of the type in [11] that is symmetric and three-point with respect to x 1 , . . . , x n and two-level with respect to t: We seek a function v ≈ u defined on the mesh ω h × ω ht and an auxiliary vector function ϕ = (ϕ 1 , . . . , ϕ n ) T with components ϕ k defined on the mesh shifted by h k /2 along x k with respect to ω h × ω ht . Here B(x) = diag {b 1 (x), . . . , b n (x)} is a diagonal matrix, div h ϕ = δ * 1 ϕ 1 + . . . + δ * n ϕ n , and the equations of the scheme are written on the same meshes as the sought functions themselves.
Let us reduce the two-level vector scheme to a three-level scheme of the standard type for v.
Theorem 1. The function v satisfies the time-explicit symmetric difference scheme that is three-point in all variables x 1 , . . . , x n and t (three-level in t), where α 0 (ζ) := ζ coth ζ and the free term has the formf Proof. We apply the operator div h to Eq. (2) to obtain the equation Let us eliminate the function ϕ from this equation by expressing div h ϕ and div h ϕ =δ t v −f by virtue of Eq. (3) as We symmetrize the notation in the last equation by multiplying it by e ht/(2τ ) : Further, we use the simple formulas and proceed to the equation We divide the last equation by 2 sinh(h t /(2τ )) and finally obtain the explicit difference scheme (4), (5) three-level in time for v. The proof of the theorem is complete.
) as ζ → +∞; further, |α 0 (ζ) − ζ| < 0.075, 0.015, 0.0027 already for ζ ≥ 2, 3, 4, respectively. Note that, in the limit as h t /(2τ ) → +∞, the expression τ α 0 (h t /(2τ )) is replaced by the h t /2, and the three-level difference scheme (4) passes into the explicit two-level scheme Note that, in the case of constant coefficients, the multiplication of the difference scheme (4) by τ and the replacement of the steps (h, h t ) → (h/τ, h t /τ ) lead to the special case of this scheme with τ = 1 (including the argument of the coefficient α 0 (h t /(2τ ))). Therefore, some properties of this scheme are independent of τ . This is also true for the spectral condition for time-uniform stability studied below.
It is important that the above method of reducing a two-level scheme to a three-level one can be applied not only to the scalar case in question but also to the quasi-gasdynamic system itself. The reduction of other two-level vector methods for second-order hyperbolic equations to threelevel methods was carried out, in particular, in [10,13,14]. We also note that numerous two-level difference schemes for the transport equation can be found, e.g., in [15].

MULTIDIMENSIONAL TRANSPORT EQUATION WITH PERTURBATIONS AND AN EXPLICIT THREE-LEVEL DIFFERENCE SCHEME FOR IT
Consider the multidimensional homogeneous transport equation perturbed by terms containing the second derivatives of the sought function with respect to t and x i and multiplied by parameters: where b i and a ij are constant coefficients and α 0 > 0 and τ > 0 are parameters. Recall that repeated indices i and j assume summation from 1 to n. The introduction of the parameter α 0 is due to the form of the scheme (4). (We can readily eliminate this parameter by the changes of variables α 0 τ → τ and a ij → a ij /α 0 ; however, we preserve the parameter for clarity.) The sought function u(x, t) is defined and Eq. (6) is considered on R n ×[0, +∞). The functions u| t=0 and ∂ t u| t=0 are assumed to be given. Without loss of generality, we assume that the matrix A = (a ij ) n i,j=1 is symmetric.
First, we draw attention to the fact that the differentiation of the transport equation Thus, the solution of the transport equation also satisfies the second-order equation of the type (6) where i is the imaginary unit and the symbol "·" stands for the inner product of vectors. Substituting these solutions into Eq. (6), we obtain the ordinary differential equation for w. Let us study the uniform boundedness of the solution to this equation with respect to t. Let Ker A and Im A be the kernel and image of the matrix A (i.e., of the corresponding linear operator acting in R n ). Recall that since A = A T , it follows that R n can be represented in the form of the orthogonal direct sum of the subspaces Ker A and Im A.
be satisfied for any w(ξ, 0) and ∂ t w(ξ, 0) is that A ≥ 0. Along with the condition A ≥ 0, the following conditions serve as criteria (i.e., necessary and sufficient conditions) for the estimate (8) to hold: is the restriction of A to Im A; otherwise, b = 0 for Im A = {0} (i.e., A = 0). In particular, if A = diag {a 1 , . . . , a n }, then the criterion is given by the set of conditions a k ≥ 0, and if a k = 0, then b k = 0 for 1 ≤ k ≤ n; α 0 Proof. The characteristic equation for the ordinary differential equation (7) has the form Property (8) is equivalent to the property that both roots of the polynomial r 2 (z) lie in the closed left half-planeP l := {z ∈ C: z R ≤ 0} and there exist no multiple roots on the imaginary axis (i.e., z = 0 is not a multiple root, which is exactly the case). Hereafter, we use the standard notation of complex numbers in the form z = z R + iz I . To analyze the last property, we use the generalized Routh-Hurwitz criterion for polynomials with complex coefficients [12, Ch. XV, Sec. 18]. To this end, we write the auxiliary polynomial −ir 2 (iz) = z + θ I + i(α 0 τ z 2 − θ R ) and the fourth-order determinant that is the resultant (up to the sign) of the arising two polynomials with real coefficients, Its leading principal minor of the second order V 2 = α 0 τ is positive, and the determinant itself is equal to Then the criterion for the roots r 2 (z) to lie in the open left half-plane P l := {z ∈ C: z R < 0} is the condition V 4 > 0. For V 4 = 0, as is easy to check, there exists a pure imaginary root −iθ I , and therefore, the second root is iθ I − 1/(α 0 τ ) ∈ P l . Thus, the criterion for both roots to lie inP l is the condition V 4 ≥ 0; i.e., It can be written in the form α 0 b ⊗ b ≤ A, and this is, in a sense, a condition of prevalence of A over b. This condition does not contain the parameter τ ; this is natural in accordance with the preceding. Conditions of this kind play an essential role in the papers [4,10]. It follows from the criterion (10) that A ≥ 0. If A > 0, then it can be written in the form Introducing the vector η := A 1/2 ξ, we obtain the equivalent explicit criterion in Item 1 of the theorem. If, however Ker A = {0}, then b ⊥ Ker A by virtue of the criterion (10), hence b ∈ Im A, and the criterion itself acquires the form α 0 (b · ξ) 2 ≤ Aξ · ξ = Aξ · ξ for all ξ ∈ Im A; this reduces Item 2 of the theorem to part 1 in which Im A serves as R n for Im A = {0}. The case of Im A = {0} is obvious. The proof of the theorem is complete. Now for Eq. (6) we consider an explicit symmetric difference scheme that is three-point in all variables x 1 , . . . , x n and t (three-level with respect to t): where the sought function v is defined on ω h ×ω ht . The values v(x k , 0) and v(x k , h t ) at the initial and first levels in time, where x k = (k 1 h 1 , . . . , k n h n ), k = (k 1 , . . . , k n ) ∈ Z n , are assumed to be given. Here, for clarity, we take the case of A = diag {a 1 , . . . , a n } with a i > 0 for 1 ≤ i ≤ n.
Since div h (Bs) = b iδi and div h (A∇ h ) = a i Λ i for constant coefficients b i and a i , we see that this scheme is obtained from the scheme (4). Consider complex particular solutions of the difference scheme (11) of the form v(x k , t m ) = e ik·ξ w(ξ, t m ), k ∈ Z n , m ≥ 0, ξ = (ξ 1 , . . . , ξ n ) T ∈ [−π, π] n .
The function e ik·ξ is an eigenfunction of the operatorsδ l and Λ l with the respective eigenvalues ih −1 l sin ξ l and −4h −2 l sin 2 (ξ l /2). Therefore, w satisfies the three-level difference equation (actually, a difference scheme for the ordinary differential equation (7)) By the definitions of the operators Λ t andδ t , after the multiplication by h 2 t /(α 0 τ ), the last equation can be written in the recurrent form where The characteristic equation for the difference equation (14) has the form The validity of the spectral stability condition (12) is equivalent to the property that both roots of this equation lie in the closed unit disk |q| ≤ 1 in C and there exist no multiple roots on the boundary of the disk. For comparison, let us first study the three-level symmetric scheme which is the special case of the scheme (11) for τ = 0.
Theorem 3. For the difference scheme (16), a criterion for the spectral condition of time-uniform stability to be satisfied has the form of the Courant type condition Proof. For τ = 0, the scheme (13) implies the recurrent equation which is different from (14). The corresponding characteristic equation q 2 + 2ih t γq − 1 = 0 has the roots In the case of h t |γ| > 1, the maximum of |q 1,2 | is equal to h t |γ| + h 2 t γ 2 − 1 > 1. In the case of h t |γ| ≤ 1, we have |q 1,2 | = 1, but for h t |γ| = 1 the root becomes multiple. This leads to the condition h t |γ| < 1 for all ξ ∈ [−π, π] n and finally (for ξ = ((π/2) sgn b 1 , . . . , (π/2) sgn b n )) to condition (17). The proof of the theorem is complete. Now let us proceed to analyzing the difference scheme (11).
The validity of the spectral condition for the time-uniform stability of the difference scheme (11) is independent of the parameter τ > 0.
For the validity of this condition, it is necessary that the following two Courant type conditions on h t hold: the second of these conditions almost coincides with condition (17). The following condition is also necessary: (9)).
The following condition is sufficient for the indicated spectral condition to hold: For n = 1, it acquires the form and serves as a criterion.
Remark 1. Let us give some comments on the conditions in this theorem. The second condition in (18) follows from the first condition and condition (19) by virtue of the Cauchy-Schwarz inequality, As will be proved, the sharper condition whose right-hand side, by virtue of the first condition in (18), is at most 1, is necessary as well.
Obviously, the inequality serves as a sufficient condition coarser than (20) (cf. the first condition in (18) and condition (19)).
Remark 2. One can use the 1-norms of the coefficient vectors a := (a 1 , . . . , a n ) and b and, for a = 0 and b = 0, introduce special mean space steps h a > 0 and h b > 0 by the formulas Proof. If the roots of Eq. (15) are multiple, then their modulus is equal to |1 − d|/(1 + d) and is less than 1. If θ = 0, then the polynomial p 2 (q) has 1 and (1 − d)/(1 + d) as its roots (the latter is obviously less than 1 in modulus). Consequently, in what follows we can assume that θ = 0, and hence p 2 (1) = 2θ = 0, and that these roots are simple.
Here r 2 (1) = 2(1 + d) > 0. By z 1 and z 2 we denote the roots of the polynomial r 2 (z) and obtain conditions for them to belong to the half-planeP l .
First, let θ I = 0. As is well known, for θ R = 0 one has z 1 , z 2 ∈P l under the following conditions on the coefficients of the polynomial: d/θ R ≥ 0 and (2 − θ R )/θ R ≥ 0. Since d > 0, this leads to the condition 0 < θ R ≤ 2 (with z 1 , z 2 ∈ P l for θ R < 2 or z 1 = −d < 0 and z 2 = 0 for θ R = 2). Now let θ I = 0 and hence ξ = 0. We again write the auxiliary polynomial and the fourth-order determinant that is the resultant (up to the sign) of the two polynomials with real coefficients on the right-hand side, Its leading principal second-order minor V 2 = 2dθ R is positive for ξ = 0. Let us multiply the second and fourth rows in V 4 by θ I /θ R , subtract them from the first and third rows, respectively, and expand the resulting determinant by the first column: (One can readily see that the final formula also holds for θ R = 0.) Recall that for V 4 = 0 the criterion for the roots r 2 (z) to lie in P l is given by the condition V 4 > 0. For V 4 = 0, there exists a pure imaginary root z 1 = iy 1 . It is a common root of the two real polynomials on the right-hand side in (23); it readily follows that y 1 = θ I /(dθ R ). Then z 2 = −iy −1 1 ((2/θ) − 1), and hence z 2,R = −d/(θ R |θ| 2 ) < 0; i.e., z 2 ∈ P l . Consequently, the criterion for both roots to lie inP l and for the imaginary axis not to contain multiple roots is given by the condition V 4 ≥ 0; i.e., In accordance with the preceding analysis, this condition also holds for θ I = 0. Geometrically, the resulting criterion means that θ belongs to an ellipse on C, Since θ I /d = b i (h t /h i ) sin ξ i , we conclude that the parameter τ does not occur in this criterion. Obviously, the inequalities 0 ≤ θ R ≤ 2, |θ I |/d ≤ 1 for all ξ ∈ [−π, π] n (the condition for θ to belong to a rectangle enclosing the ellipse) are a necessary condition for the criterion (25) to hold. The maximum of θ R ≥ 0 is obviously attained at ξ = (π, . . . , π), and the maximum of |θ I |/d is attained at ξ = ((π/2) sgn b 1 , . . . , (π/2) sgn b n ); this leads to conditions (18). Setting ξ = ((π/2) sgn b 1 , . . . , (π/2) sgn b n ) in (24), we arrive at the inequality i.e., at the necessary condition (22). Let us also take ξ k = εb k h k /(a k h t ), k = 1, . . . , n, in the criterion (24) and expand both sides in powers of ε → 0; we obtain Dividing both sides by ε 2 and passing to the limit as ε → 0, we derive the necessary condition (19).
On the other hand, we have | sin ξ k | = 2 σ k (1 − σ k ), where σ k := sin 2 (ξ k /2) runs through the segment [0, 1], and by the Cauchy-Schwarz inequality one has the estimate Therefore, canceling by the factor 2θ R = 0 (which is equivalent to σ := (σ 1 , . . . , σ n ) = 0), we conclude that the criterion (24) follows from the inequality Its left-hand side is an affine function σ of simple form, and therefore, as is easily seen, this inequality is equivalent to inequality (20), the latter thereby being a sufficient condition for the validity of the criterion. For n = 1, after dividing by 2θ R , the criterion (24) acquires the form (26), and therefore, conditions (21) become a criterion.
Remark 3. In addition, note that in the case of coefficients a i of arbitrary signs, the situation in Theorem 4 becomes similar to the one indicated in the criterion (9). First, it was shown in the proof of Theorem 4 that the condition θ R ≥ 0 is necessary in the case θ I = 0. For θ I = 0 and θ R < 0, we have V 2 < 0 and V 4 < 0; hence one of the roots r 2 (z) lies in the right half-plane z R > 0. Therefore, in the case θ I = 0 it is also necessary that θ R ≥ 0 for all ξ ∈ [−π, π] n . The latter implies the necessity of the condition a i ≥ 0, 1 ≤ i ≤ n. Second, if a k = 0 for some k, then we take ξ i = 0 for i = k and sin ξ k = 0. Then θ R = 0 and θ I /d = b k (h t /h k ) sin ξ k . For b k = 0, the roots of r 2 (z) = i(θ I z 2 − 2idz − θ I − 2i) are such that z 1 + z 2 = 2id/θ I ; it follows that z 1R + z 2R = 0, and for z 1 , z 2 ∈P l we have z 1R = z 2R = 0. However, the product z 1 z 2 = −1 − (2i/θ I ) should then be real, which is not the case. Hence the condition b k = 0 is necessary in the case a k = 0.
Therefore, Theorem 4 also remains valid for a i ≥ 0, 1 ≤ i ≤ n, with the only distinction that the sums in conditions (19) and (20) should be taken over all 1 ≤ k ≤ n such that a k > 0 and one should bear in mind the property that if a k = 0, then b k = 0.
In the particular case of b = 0, we have θ I = 0, and the first condition in (18) in Theorem 4 becomes a criterion. It is similar to the stability condition for the three-level difference scheme with parameter τ in [9] in the case of the weight σ = 0. In this case, one should take B h = α 0 I, B 1h = I (here I is the identity operator), and A h = −τ a i Λ i in this scheme, with the stability condition having the form h t (α h /2) ≤ ε 1 with an arbitrary 0 < ε 1 < 1, where N k being the number of partitioning segments along x k (unlike the present paper, the paper [9] considers an initial-boundary value problem), with sin 2 (ζ k /2) → 1 and the last inequality becoming the equality as N k → +∞, 1 ≤ k ≤ n. The case of ε 1 = 1 is also possible (with some relaxation of the norms used); in this connection, see [17,Remark 1].
FUNDING This work was supported by the Russian Science Foundation, project no. 19-11-00104.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.