CC BY
26
Вычислительные технологии
Том 3, № 1, 1998
THE STABILITY OF NUMERICAL BOUNDARY TREATMENT FOR FINITE-DIFFERENCE SPLITTING SCHEME FOR THE ACOUSTICS EQUATIONS SYSTEM
A.M. Blokhin, D.A. Krymskikh Institute of Mathematics SB RAS, Novosibirsk, Russia e-mail: [email protected]
На примере начальной и начально-краевой задач для системы уравнений акустики изучаются вопросы влияния краевых условий на устойчивость конечно-разностной схемы расщепления. Эта схема широко используется для численной аппроксимации решений задач аэродинамики [1]. Показано, что устойчивость схемы зависит не только от типа задачи (начальная или начально-краевая), но и от ее размерности.
1. Introduction
Nowadays, researchers often use finite-difference schemes to find approximate solutions to equations of mathematical physics. A wide range of finite-difference schemes is applied to numerical calculations of external and internal flows of a gas with shock waves in problems of supersonic aerodynamics; among them there are the Godunov scheme [2], the MacCormack scheme [3] and others. The paradox of situation is that a researcher seeks numerical approximations to solutions of a mixed (i. e., initial-boundary value) problem having not a notion if the solution really exists.
Necessity of simultaneous study of the differential problem and its finite-difference analog has long been discussed in the theory of differential equations. We refer the reader to the fundamental paper [4] and monographs [5-7] where this approach was developed and applied to problems of mathematical physics.
As in [7], construction of a difference scheme to the gas dynamics equations and study of its stability were based on the requirement of adequacy between the difference model and the differential problem. The author of [7] considers the difference scheme adequate if it provides the existence theorem for the differential problem. The latter is of great importance for calculation practice since with the existence theorem in hands a researcher can be sure that approximate solutions really converge to the solution of the differential problem as grid steps tend to zero.
While considering numerical approximations to solutions of a mixed problem, a researcher naturally arrives at the question on stability of the difference scheme [8, 9].
© A.M. Blokhin, D.A. Krymskikh, 1998.
In order to answer this question in regard to the problems of gas dynamics considered in [7] the following program was realized. First, well-posedness of the linear differential mixed problem was studied with the help of dissipative energy integrals. Then the difference scheme which admits existence of a difference analog of the dissipative energy integral was constructed. Existence of such analog gives possibility to obtain an energy estimate and then to prove stability of the suggested finite-difference scheme (see also [8]). As the result, adequacy of the difference model to the linear differential problem was established.
It is seen that such program requires careful theoretical study of the differential problem. Apparently, in the framework of existing techniques it can not be done in every case even for linear problems, not to mention nonlinear mixed problems. By this reason researchers usually restrict themselves to the spectrum analysis of the difference initial value problem (see [8, 9] on spectrum analysis of stability of difference schemes) and then directly use the obtained restrictions on grid steps in the case of the difference mixed problem.
There exist examples of finite-difference scheme, both explicit and implicit, for which it is strictly proved that stability of the difference mixed problem does not follow from stability of the difference initial value problem [10]. It means absurdity of attempts to use recommendations obtained for the difference initial value problem in the case of the difference mixed problem.
In this paper, using the well-known finite-difference splitting scheme [1], we give another confirmation to this conclusion. As well as in [10], the linear mixed problem on stability of shock waves is considered.
In § 2 we describe differential problems. In § 3 we construct finite-difference analogs of these differential problems on the basis of the splitting scheme from [1]. This algorithm leads to additional boundary conditions which are used in the difference analogs of the mixed problems. This section also contains results of spectrum analysis of difference initial-value problem. The spectrum analysis of the one-dimensional difference mixed problem is carried out in § 4. The difference analog of the Hadamard-type ill-posedness example which proves instability of the two-dimensional difference mixed problem is given in § 5 (concerning the Hadamard-type ill-posedness examples see [6]). Concluding remarks form the last section.
2. Preliminaries
Following [7, 10], we consider the linear mixed problem on stability of shock wave. Nonstationary equations of gas dynamics written in the Cartesian co-ordinates x, y and linearized with respect to the main constant solution behind the oblique shock with the equation x = 0 under certain scaling form the well-known acoustics equations system. The Rankine — Hugoniot relations on the shock wave are linearized in a similar way, and then the function which describes a small perturbation of the front is excluded from the obtained relations. As the results we arrive at the following linear mixed problem on stability of the shock wave (for detail, see [7]): in the domain R++ the solution W = W(t,x,y) to the acoustics equations system
Wt + AjWx + AWy = 0 , which satisfies boundary conditions at x = 0, (t,y) £ R+:
u + dp = 0 , vt + uvy — Xpy = 0 and initial data at t = 0, x > 0, y £ R1
W(0, x, y) = $(x,y)
(2.1a)
y
(2.1b)
(2.1c)
is sought. Here
W
W1 W2 W3
P
Mu Mv
Ai
1 k 0 w 0 k
k 1 0 , A2 = 0 w 0
0 0 1 1 k 0 w
M, k, w are constants:
0 < M < 1 , k = 1/M > 1 , w > 0;
d, A are constants defined in [7]; p, u, v are small perturbations of the pressure and the components of the velocity vector;
R++ = {(t,x,y) | t,x > 0 , y e R1} , R+ = {(t,y) 1t > 0 , y e R1} .
Along with problem (2.1) we consider its one-dimensional variant: in the domain R++ the solution W = W(t,x) to the acoustics equations system
Wt + AWx = 0 : which satisfies boundary condition at t > 0, x = 0:
u + dp = 0
and initial data at t = 0, x > 0
W(0, x) = $(x)
(2.2a)
(2.2b) (2.2c)
is sought. Here
W
w1
W2
P \ A 1 1 K
Mu
k1
R++ = {(t,x) | t,x > 0} .
Finally, we consider also the initial value problem for the acoustics equations system: in the
domain R+ the solution W = W(t, x, y) to the system
Wt + AiWx + A2Wy = 0 , which satisfies initial data at t = 0, (x,y) e R2:
W(0, x, y) = $(x,y)
is sought. Here
R+ = {(t, x, y) 11> 0 , (x,y) e R2}
(2.3a)
(2.3b)
3. Difference model. Spectrum analysis of difference initial value problem
We introduce the grid
= {(nA,khx,lhy) | n,k, |/| = 0,1,...}
bx^lOy
where A, hx, hy are the grid steps, in the domain
R++ = {(t, x, y) | t, x > 0 , y e R1}.
The system (2.1a) is approximated by the splitting scheme over the physical processes and spatial variables from [1]. By its linearity, it turns into
W = — (rx/af+ ryuhv + rx(£i£ + B^) + ry+ BJij)) Wnkl, (3.1a)
(l + Af)Wki4 = W , (3.1b)
(l+ An) WH2 = Wf , (3.1c)
(/3 + &(Bi£ + B$)W^ = Wkf , (3.1d)
(/a + P4(B2V + B2^)W^ = Wkf , (3.1e)
TWki = W« , n, |1| =0,1,..., k = 1, 2,.... (3.1f)
/ 000 \ / 000" Bi = к 00 , B2 = 000 \ 0 0 0 / \ к 00
B* is the transpose to Bi; /3 is the identity matrix of order 3;
r = фп - 1 , С = фк - 1, C = 1 - ф— > П = (-Ф2 + 4ф1 - 3)/2 , n = (3 - 4Ф-1 + ф-2)/2
are the difference operators:
Фкwn = wn+i 1, Ф1 wn = wnI+i,
^ = wn+1 ./,-1wn = wn '.-1л*7"П _ ЛЛТП
Here
= wni+1, ^wn = wn_i 1
^ W = Wi-i, Wki = W(nA , khx , lhy) is the numerical approximation to W. Functions Wk are grid vector one,
ßl = fx^i , ß2 = ryW«2 , ßs = ,
ß4 = fya , fx = A/hx , fy = A/hy; a1, a2, a3, a4 > 0 are constants.
We begin with construction of the difference analog of mixed problem (2.1). Algorithm (3.1) requires some additional boundary conditions for grid vector functions Wk1. Following the algorithm from [1], we assign the following boundary conditions at points of the grid with k = 0:
W0/4 = 0, (3.2a)
(1+ ^Wjf = 0 , (3.2b)
(w2)w4 = 0 , K&fWO/4 = Mi/2, (w3)w4 = (w3)0/2 , (3.2c)
/3 + A(B2n + B2n)) W01 = Wf , (3.2d)
tW = W^, m = 0,1,.... (3.2e)
Besides, we assign the boundary conditions of "tending to zero" and "periodicity" for Wk so that the solution W^ ^ 0 at k ^ to and W^ = Wni+m where m is a positive integer number. In what follows we will seek solutions to the difference mixed problem with these properties.
Now we give a brief description how the splitting scheme (3.1) is realized (its detailed description can be found in [1]). Let the solution to this scheme on the n-th time layer be known, it tends to zero by k and is m-periodic by /. First from (3.1a) by the explicit formulae we determine the grid vector function Wk1, k =1, 2,... , |1| =0,1,.... With the help of (3.1b)
~ 1/4 _
we determine the grid vector function Wk1 at k = 1, K0 — 1, |/| = 0,1,...
W1/4 = W ki+AW k-4i i k1 1 + A '
~ 1/4
assuming for the sake of simplicity that Wk1 = 0 at k > K0, |1| = 0,1,... . Here K0 is sufficiently large positive integer number. Consequently, in a view of (3.2a), the grid vector
~ 1/4 _
function Wk1 is determined at k = 1,K0 — 1 by the running count. Then from (3.1c) for every fixed index k = 0, 1, . . . we obtain a vector difference equation and pose the periodicity boundary conditions:
1/2 - 1/2 - 1/2 - 1/4
A2Wfc;_2 — 4^2Wki_1 + (2 + 3&)Wkl = 2Wki , l = 2, m +1,
-1/2 ~ 1/2 ~ 1/2 ~ 1/2 Wk0 = Wkm , Wk1 = Wkm+1 .
The solution to this vector difference equation is found by the scalar sweep method which is well-posed for sufficiently small time step A. Now from (3.1d) it follows that (w-3)k/4 = (w-3)k/2, for
the remained components of the grid vector function Wk/4 we obtain the three-point difference equations at each fixed index I = 0, ±1,...
—K2A32(w-1)k/411 + (1 + 2K2A2)(w-1)k/4 — k2A2(*-1 )k+411 = (w-1)k/2 — «^wif,
k = l,Ki - 1,
-K2^32(wi2)k/-11 + (1 + 2K2^32)(W2)k/4 - «^(W1 = " - (¿1)^
k = 1,K2 - 1
We take first two relations from (3.2c) as the boundary conditions for these difference equations, and we assume that (w1)k/4 = 0 at k > K1 and (w2)k/4 = 0 at k > K2 where K1,
K2 are sufficiently large positive integer numbers. Solutions of these difference equations are determined by the sweep method which is well-posed at sufficiently small A. Then it follows from (3.1e) that (wi2)k1 = (w2)k/4, and it is possible to obtain five-point difference equations for the grid functions (ww1)k1, (w«3)ki, the boundary conditions are the conditions of periodicity. Solutions to these equations also exist for sufficiently small A. Finally the solution on the (n + 1)-th time layer is determined from (3.1f). The obtained solution W]n+1 tends to zero by k and is m-periodic by l. Thus, given the initial data W^ = such that — 0 at k —> to, = $k1+m, after realization of the algorithm from above we come to the solution Wn, n = 1, 2,... to (3.1) which tends to zero by k and is m-periodic by l. Boundary conditions (2.1b) are approximated in the following way
+ <i+1 = 0 , ryAnpS+1 — (t + ry^ K = 0 . (33)
Algorithm from [1] provides subjection of the solution on the (n + 1)-th time layer to the boundary conditions (33), i.e., the so-called determining component, for example (w1)0+1, is chosen and the rest of components of the vector W0+1 are redetermined with the help of boundary conditions (33).
In order to carry out theoretical investigation, in particular, the spectrum analysis, it is convenient to exclude the grid vector functions Wk1 from (3.1). As the result we arrive at the implicit scheme:
(1 + A£)(1 + An) (/3 + A (B1e + Bio) (/3 + A(B2n + Bin)) rW™ + + (rx/af + ryW/3/7 + rx(B1e + BiO + ry(B2n + B^)) W = 0 . (34)
In order to stay within the limits of the grid it is necessary to take k = 2, 3,...
in the first equation of scheme (34) and k = 1, 2, . . . in the second and third equations of scheme (34). Thus, with four boundary conditions on the left boundary of the grid and
the "periodicity" and "tending-to-zero" conditions, similar to the ones for (3.1), scheme (34) becomes solvable, and its solutions tend to zero by the index k and are m-periodic by l if the initial data $kl — 0 at k — to, = $k1+m. As we have only two boundary conditions
(33) on the left boundary, we, consequently, have to find two additional boundary conditions. We derive these conditions taking into account (3.1). By this algorithm, one component of the vector W0+1 (let it be (w1)0+1) is determined immediately after realization splitting scheme (3.1). Excluding corresponding components of grid vector functions W01 in (3.2), we obtain the relation for (w1)ni+1:
(1 + an)(1 — k2aV? )£r № = 0,
which is considered as the first additional boundary condition. Besides, the first equation in
(34) at points of the grid with k =1 should be realized in a special way. Excluding Wk1
with k =1 from (3.1) and taking into account that (w2)3/4 = 0 (see the first condition in (3.2c)), we have the first equation of (34) at points of the grid with k =1
«1(^1 + «2(^2)^ + a3n (w3)n = 0 .
Here
«1 = br + rx° + ry, «2 = ft(A&r + rx°) ,
«3 = K(Abr + ry), b = (1 + AO)(1 + An) •
This relation serves as the second additional boundary condition.
We now formulate a difference analog of mixed problem (2.1): we look for the solution Wk of the scheme
a!(w1)n + a^w)£ + a^W)£ = 0 , k = 2, 3,... ,
«2^(w1)n + a1(w2)ni + a4(w3)n = 0 ,
«3^1)^ + «1(^3)^ = 0 , k = 1, 3,...,
n, |1| = 0,1,..., (3.5a)
satisfying the boundary conditions
d(w1)ni+1 + kMS+1 = 0 ,
ryAnM0+1 - K(T + ry= 0 ,
(1 + ^2n)(1 - (W1)S = 0 ,
«1(^1)^ + «2(^2)^1 + 03/7(^3)^ = 0 ,
n, |1| = 0,1,...
Wki = , k, |/| = 0,1,.... (3.5b)
and the initial data Here
«2 = + fx) , ß4 = K2ß3ß4&£nr •
Scheme (3.5a) is the component-wise presentation of scheme (34). Boundary conditions (3.5b) consist of boundary conditions (33) and two additional boundary conditions. We now construct a difference analog of mixed problem (2.2). In the domain
—2
R++ = {(t,x) | t,x > 0}
we introduce the grid
= {(nA, kh) | n,k = 0,1... } , where A, h are grid steps. The one-dimensional variant of splitting scheme (3.1) is given below
Wk = —r((/2 + B*)£ + wn ,
1+ W1/2 = W k, (/2 + + B t)) W k = W k/2,
t wn = w k,
n = 0,1,... , k = 1, 2,... . (36)
Here
B1= (0 0
\ k 0
/2 is the identity matrix of order 2, Wn = W(nA , kh) is the numerical approximation to W, Wk are grid vector functions,
A1 = ra1, A2 = ra2, r = A/h;
ai, a2 > 0 are constants.
It is easy to show that, given
W,
1/2
0 0,
(W2)1 = o
T Wn
Wn
(37)
and "tending-to-zero" (W^" ^ 0 at k ^ to) boundary conditions, at sufficient small time step A the solution to (3.6) exists, and it tends to zero by k, if the initial data $k ^ 0 at k ^ to. Excluding Wk from (3.6), we formulate a difference analog to (2.2): we seek the solution
Wn of the scheme
a1(w1)n + «2^ (w2)n «2^(w1)n + «1(w2)n = 0 , k
which satisfies the boundary conditions
0, k
: 1, 3,.
2, 3,
n
0, 1,
(3.8a)
d(wi)n+1 + K(w2)n+1 = 0
ct (wi)n = 0, «1(w1)n + a2(w2)n = 0, n = 0,1,
and the initial data
Here
W0
, k = 0,1,
(3.8b) (3.8c)
«1 = (1 + )t + rê, «2 = k(&(1 + A£)t + r) ,
02 = «(&(!+ A£> + rê) .
It should be noticed that the first condition from (3.8b) is obtained by approximation of boundary condition (2.2b), the second and third ones are determined by scheme (3.6) in view of boundary conditions (37).
Finally we construct a difference analog of initial value problem (2.3). We consider scheme (3.1) on the grid
= |(nA,khx,/hy ) | n, |k|, |1| =0,1,...} .
It is easy to show that, given the periodicity boundary condition (Wn+TO1l = W^, Wni+TO2 = W£i, m1, m2 are positive integers), at a sufficiently small time step A, the solution to (3.1) exists, it is m1-periodic by k and m2-periodic by I if the initial data are periodic: i = ^fci,
k i+m2 ^ ki.
We now formulate a difference analog of initial value problem (2.3): at the grid we
seek the solution W^ to (34), satisfying initial data (3.5b).
Further, we carry out the spectrum analysis of the difference initial value problem, i. e. we check out the fulfilment of the necessary von Neumann stability condition (see also [8, 9]). Under assumption on periodicity of initial data (3.5b) $ki = eik^Wo, we look for the solution to (34) of the form
Wn = qneikVl^Wo , (39)
where W0 is a constant vector; ^ are arbitrary real values, i = \J -1, q = is a
complex value. Values q for the initial-value problem form the spectrum of the operator of
scheme (34) while passing from one time layer to another one (see [9]). The necessary von Neumann stability condition for (34) is in fulfilment of the inequality
|q(^)|< 1, (310)
for every ^ that means that the spectrum of the difference initial-value problem lies in the unit circle.
As inequality (310) is difficult to be cheeked out for arbitrary ^ we restrict ourselves to a particular case. Assumed ^ = we determine conditions when (310) becomes true. With this purpose in a view we substitute (39) into (34). As the result we obtain a linear algebraic system to determine components of the vector W0
rWo = 0. (311)
Here
(71 72^ 73m 72 a 71 74 73m 0 71
Y1 = a + rxa + rywm , 72 = (Aa + rx)re, 73 = (Aa + ry)k, Y4 = K2^3^4aam , a = (1 + Aa)(1 + A2m)S , S = q — 1, a = e^ — 1, a =1 — e-^, m = (—e^ + 4e^ — 3)/2 , m = (3 — 4e-i^ + e-i2^)/2 .
If (39) is the nontrivial solution to (34), then (311) must also have the nontrivial solution, consequently, the determinant of the matrix r must be equal to zero. The equality detr = 0 is treated as an equation to determine q = q(^). It is difficult to obtain general explicit formulae for q(^). By this reason, we use the technique of almost eigenfunctions from [9]. Let consider the case of "long waves", i.e., we consider partial solutions (39) with the wave length much more then grid steps hx, hy, i.e., ^ = o(1). At ^ = 0 detr = 0 if q =1. Consequently, if = e^', where e > 0 is a small parameter, is real, then q = 1 + eq' + e2q'' + O(e3). In this case condition (310) transforms into: Re q' < 0; if Re q' = 0 then 2Re q''+(Im q')2 < 0.
Put = e^' into detr = 0, then find three different roots
qj = 1 + eqj + e2qj + O(e3), j = 1, 2, 3 , q1 = —ib^',
q'i' = — ( (A + A)b1 + rx/2)(^')2, q2,3 = —i(b1 ± v^V, q2',3 = — ((A + A)(b1 ± Vb2) + rx/2 + bs(1 ± b1/(^')2,
where
b1 = rx + ryw, b2 = K2(rX + ry)
63 = K2(Arx + Ary) .
Here and in the sequel, a square root means its principal value. Consequently, condition (310) is fulfilled for any if
(61 - 2(A + &))&! - rx < 0 ,
(61 ±^/62 - 2(^1 + &))(&! ± v/62) - rx - 263(1 ± 61/^62) < 0 .
The last inequalities are fulfilled, in particular, for arbitrary rx, ry if
a1 = a2 = a3 = a4 > 1/2.
Thus, the spectrum analysis of the difference initial value problem has not revealed any restrictions to the time step A.
It is seen, however, that influence of the boundary conditions, including additional conditions, on stability of (34) are not considered by the spectrum analysis carried out above. So, it is expedient to analyze the spectrum of difference mixed problems (3.5), (3.8).
4. Spectrum analysis of mixed difference scheme in one-dimensional case
In this section we analyze the spectrum of difference mixed problem (3.8). We assume that the initial data are of the form
3
** = E xkwj), (41)
j=1
where Wj) are constant vector, Xj are different complex numbers, |Xj| < 1, j = 1, 3. The solution to (3.8) with initial data (41) under certain choice of Xj = Xj (q), j = 1, 3 have the form
3
Wn = XkWj) , (42)
fc=1
where q is a complex number. Thus, for (3.8) to be stable, the inequality
|q| < 1 (43)
must be fulfilled for all partial solutions (42) with parameters |Xj| < 1, j = 1, 3.
In the sequel, instead of testing of this inequality, we will show that, under certain circumstances, difference mixed problem (3.8) does not have solutions of the type (42) with |q| > 1. First we find Xj = Xj (q), j = 1, 3 as functions of q which should be taken in (42). To do this we look for partial solutions Wn = qnXkW0 to (3.8a) where X is a complex number. Substituting this solution into (3.8a), we obtain a linear algebraic system to determine the vector W0
rW0 = 0 . (44)
Here
r f YX ^ V w Y /
Y = (x + + ra ,
u
= (ß2(X + ß1a)S + rX)Ka S = q - 1 , a = X - 1.
As Wn = qnXkW0 is the nontrivial solution to (3.8a), the determinant of the matrix r must be equal to zero. Considering detr = 0 as the equation to determine x, we derive four roots
Xj = Xj (q), j = 1, 4. The component w2j) corresponding to the root Xj ) are connected by the relations
of the vectors W0j) (W0" are solutions to (44)
,(1) + w21) =0,
,(2) - w22) = 0,
,(3) - w23) = 0,
,(4) + w24) = 0.
Explicit formulae for the roots Xj are written down for the particular case q £ > 0 is a small parameter, q' is a complex number:
M2 + O(£),
■ 2 X2 = £t rw* i
(45)
1 + eq', where
X1
^ 2 + o(£3)
X3 = 1 - e
q'
r(K + 1)
X4 = 1 +
+ O(e2),
+ O(e2)
r(K — 1)
It is apparent that |xi|, |x2| < 1 for any q'. If Req' > 0 (in this case condition (43) is broken down), then |XsI < 1 and |x4| > 1.
Thus, Xi, X2, X3 are the roots which must be taken in (42). It remains to choose q such that (42) satisfies boundary conditions (3.8b). Substituting (42) into (3.8b) and taking into account (45), we arrive at a linear algebraic system to determine the vector W0
r W n
0,
(46)
where
/ d K d K d K
ai 0 a2 0 a3 0
Yi Ui Y2 U 2 Y3 U3
v/X1 1 0 0 0 0
0 0 ^X2 -1 0 0
00
0
4".
0 VX3 -1 /
u,
a,
Wo = (w1;, ••• , wr)*,
(Xj + A^j+ Kraj ,
'3 = Xj — 1, Yj = Y(Xj), j = 173 • System (46) must have the nontrivial solution, consequently, its determinant must be equal to zero. Expending detr into series by powers of £, we have
detr = —Kr(1 — M2)M(d + k) + O(£) •
It is apparent that if d = — k then detr = 0 (it is easy to show that for the so-called perfect gas d = —k, see [7]). Thus, difference mixed problem (3.8) does not have solutions of the form (42) with parameters |q| > 0, |xjI < 1, j = 1, 3 at least for q = 1 + £q'.
2
q
5. Instability of difference mixed problem in two-dimensional case
In this section we prove instability of difference mixed problem (3.5). Let initial data (3.5b) be of the form
4 j=i
= e^J] xkW0j) , (51)
where W0j) are constant vectors, Xj are different complex numbers, |xj| < 1, j = 1 , 4; ^ i real, i = V — 1. Solution of difference mixed problem (3.5) with initial data (51) under certain choice of Xj(q, j = 1, 4 have the form:
4
W = xkW0j) , (52)
j=i
where q = q(^) is complex. Problem (3.5) is stable if all its partial solutions (52) with parameters |Xj(q, <^)| < 1, j = 1, 4 satisfy the following condition for any
|q(¥OI< 1 • (53)
Further we show that for (3.5) there exists a difference analog of the Hadamard-type ill-posedness example (see [6, 7]), i.e., a partial solution of the form (52) with |xj(q,^)| < 1, j = 1,4 at |q(<^)| > 1. In this case, on one hand solution (52) at n = 0 tends to zero at k ^ to, on the other hand it indefinitely grows up at n ^ to, what means instability of (3.5) in any grid norm (see [6, 7]).
To prove that we look for partial solutions to (3.5a) of the form W^ = qnXke^W0, where X is complex. Substituting this solution into (3.5a), we obtain the linear algebraic system to determine the vector W0
rWo = 0 , (54)
where
/ YiX Y2 Y3m \ r= Y2 Yi Y4 , Yi = via + YX,
\ Y3m 0 Yi / Y2 = + v3X), Y3 = KX((4a + v4X),
Y4 = K2Cima(Aa + X), a = X — 1, m = (—ei2^ + 4e^ — 3)/2 , m = (3 — 4e-^ + e-i2^)/2 , Y = h8 + ry wm , 8 = q — 1,
h = 1 + Am , vi = Ah8 + rx ,
V3 = Ah8 + rx , V4 = Ah8 + ry , Zi = AAh8 , (3 = AAh8 , (4 = AAh8 •
Apparently, the determinant of the matrix r must be equal to zero. The equality detr = 0 serves as the equation to determine X = X(q, We study (3.5) on "long" waves, i.e., under
is
assumption that ^ = e^', where e > 0 is a small parameter, is real. In this case q, x are sought in the form q = 1 + eq' + O(e2), x = X0 + ex' + O(e2). With this remark in a view we derive the roots of detr = 0
Xi = M2 + O(e) 2
'3H
X2=e ) +O(e)
X3 = 1 - + O(e2)
X4 = 1 + e
X5 = 1 + e
Y' — y a
fx(K2 — 1)
+ O(e2)
Y ^ + O(e2)
fx(K2 — 1)
where
y' = 4 + iry
a = k2(y')2 + (k2 - 1)K2r^')2 .
It is obvious that |xi| < 1, |x2| < 1 for any q', If Re q' > 0 (in this case |q| > 1), then
(j) „„(j)
i
w.
2
Ix3| < 1) |X41 < 1) iXs| > 1- Thus, in (52) we should take Xi, X2, X3, X4- Components w
of the vector W0j) (it is the solution to (54) corresponding to the root Xj, j = 1, 4) must be connected by relations
(j) (j) - (j) Yij Xj w1 + Y2j w2 + Y3j mw3
0,
(j) (j) Y3j mwi + Yij w3
0
j = 1, 4 , (55)
where Yij = Yi(Xj), Y2j = Y2(Xj), Y3j = Y3(Xj)• Substituting solutions of the form (52) into boundary conditions (3.5b) and taking into account relations (55), we obtain the linear algebraic system to determine the vector Wo
r W 0 = 0 ,
(56)
where
r
d K 0 d K 0 d K 0 d K 0
K i 0 K 3 KYi 0 K3 K1 0 K3 K1 0 K3
0 0 ^2 0 0 ^3 0 0 ^4 0 0
Yii Y2i Y3i Yi2 Y22 Y32 Yi3 Y23 Y33 Yi4 Y24 Y34
YiiXi Y2i Y3im 0 0 0 0 0 0 0 0 0
Y3im 0 Yii 0 0 0 0 0 0 0 0 0
0 0 0 Yi2X2 Y22 Y32m 0 0 0 0 0 0
0 0 0 Y32m 0 Yi2 0 0 0 0 0 0
0 0 0 0 0 0 Yi3 X3 Y23 Y33m 0 0 0
0 0 0 0 0 0 Y32m 0 Yi2 0 0 0
0 0 0 0 0 0 0 0 0 Yi4X4 Y24 Y34"
0 0 0 0 0 0 0 0 0 Y34m 0 Yi4
W 0 =( wSi), . (4)\ * .. , w3 ;) .
72j = ((3 + v3) , 73j = Km((4aj + v4Xj) ,
= Xj — 1, j = M , KKi = ryAmq , KK3 = —k(8 + rywm) .
System (56) has the nontrivial solution if detr = 0. This equality is treated as the equation to determine q. Expending detr into series by powers of e, we finally come to
detr = e8 ■ iK3rXry^AWfM(1 — M2)x x ^ k2 y' — VK ) ^y' — ) x
x(Ar^')2 — d(Y')2 — yVS) + O(e9) = 0 .
The expression at e8 in this relation turns into zero if, for example, k2y' — \F& = 0. The last inequality is fulfilled if y' = ryThus, if > 0 then Re q' = ry> 0. Consequently, at q' = ry— iryw^', > 0 (52) is the solution to (3.5) with the parameters
|q| > 1, iXil < 1, |X2| < 1, iXai < 1, |X41 < 1.
That means existence of the Hadamard-type ill-posedness example and, consequently, instability of (3.5).
6. Conclusive remarks
Here we comment the obtained results. First, we showed that in the case of difference mixed problem it makes very little sense to follow directly the recommendations on choice of the grid steps, obtained for the difference initial value problem. Second, it was established that stability of a problem can depend on its dimension. In principle, these conclusions are not of great surprise for researchers since they have been already discussed in scientific literature (see [9, 12]). But, in our opinion, these facts must be necessarily reminded since researchers restrict themselves to very cursory theoretical analysis of the applied difference schemes. As the result of such nonstrict analysis, usual calculation instability is often treated as a new physical phenomenon.
Results of linear analysis of the finite-difference splitting scheme which was carried out in this paper does not mean that this scheme can not be used in calculation practice to find approximate solutions to nonlinear mixed problems. It has been only shown that the recommendation obtained for the linear case can not be directly applied to nonlinear problems. This fact also has been noted in scientific literature (see, for example, [13]).
The finite-difference splitting scheme which was discussed in this paper is often used for numerical approximation to solutions of nonlinear aerodynamics problems, in particular, the problem on supersonic regimes of flowing around blunted solids with application of steady-state calculations [1]. From our viewpoint, instability of the finite-difference splitting scheme established in this paper at the linear level appears in the way that steady-state calculations are very sensitive to the choice of initial data [1].
Calculation instability of numerical algorithm based on the considered finite-difference splitting scheme was discovered in [14] while calculating transsonic and hypersonic stationary regimes for flowing around blunted solids. Thus, instability of the finite-difference splitting scheme established at the linear level appears in one or other form in nonlinear problems too.
References
[1] KovENYA V. M., Tarnavskii G. A., Chernyj S. G. Method of splitting in problems of aerodynamics. Nauka, Novosibirsk, 1990 (in Russian).
[2] Godunov S. K., ZABRODIN A. B., IVANOV M.J. ET.AL. Numerical calculation of multidimensional problems of gas dynamics. Nauka, Moscow, 1976 (in Russian).
[3] MAC CORMACK R. W. The effect of viscosity in hypervelocity impact cratering. AIAA Paper, 69-354, 1969.
[4] Courant R., Friedrichs K. O., Levy H. Uber die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann., 100, 1928, 32.
[5] LADYZHENSKAJA O. A. Boundary problems of mathematical physics. Nauka, Moscow, 1973 (in Russian).
[6] Godunov S. K. Mathematical physics equations. Nauka, Moscow, 1979 (in Russian).
[7] Blokhin A. M. Energy integrals and their application to gas dynamics problems. Nauka, Novosibirsk, 1986 (in Russian).
[8] Richtmyer R. D., Morton K.W. Difference methods for initial-value problems. Interscience publishers, 1967.
[9] Godunov S.K., Ryaben'kij V. S. Difference schemes. Nauka, Moscow, 1973 (in Russian).
[10] Blokhin A. M., Alajev R. D. Energy integrals and their application to study of stability of difference scheme. Novosibirsk University Publishing House, 1993 (in Russian).
[11] Shokin Yu. I., Janenko N. N. Method of differential approximation. Applications to gas dynamics. Nauka, Novosibirsk, 1985 (in Russian).
[12] Krause E. Aeron. J. 764, 1974, 337-354.
[13] Samarskii A. A., Gulin A. V. Stability of difference schemes. Nauka, Moscow, 1973 (in Russian).
[14] TARNAvSKII G.A. On existence of domains of instability for flows with shock waves in perfect gas. Part two. Prepr. Inst. Theor. and Appl. Mech., USSR Acad. of Sciences, Siberian Branch, Novosibirsk, 1985 (in Russian).
Received for publication November 18, 1997
