NUMERICAL ANALYSIS OF GRID-CLUSTERING RULES FOR PROBLEMS WITH POWER OF THE FIRST TYPE BOUNDARY LAYERS Текст научной статьи по специальности «Математика»

Вычислительные технологии, 2020, том 25, № 1, с. 49-65. © ИВТ СО РАН, 2020 Computational Technologies, 2020, vol. 25, no. 1, pp. 49-65. © ICT SB RAS, 2020

ISSN 1560-7534 elSSN 2313-691X

COMPUTATIONAL TECHNOLOGIES

DOI:10.25743/ICT.2020.25.1.004

Numerical analysis of grid-clustering rules for problems with power of the first type boundary layers

V. D. LlSEIKIN1'2, S. KARASULJIC3

1 Institute of Computational Technologies SB RAS, 630090, Novosibirsk, Russia 2Novosibirsk State University, 630090, Novosibirsk, Russia 3University of Tuzla, 75 000 Tuzla, Bosnia and Herzegovina Corresponding author: Liseikin Vladimir D., e-mail: [email protected]

Received June 10, 2019, revised December 2, 2019, accepted December 18, 2019

This paper demonstrates results of numerical experiments on some popular and new layer-resolving grids applied for solving one-dimensional singularly-perturbed problems having power of the first type boundary layers.

Keywords: singularly perturbed equations, small parameter, boundary and interior layers, grid generation.

Citation: Liseikin V.D., KarasuljiC S. Numerical analysis of grid-clustering rules for problems with power of the first type boundary layers. Computational Technologies. 2020; 25(1):49-65.

Introduction

The present paper describes experiments on some popular and other forms of layer-resolving grids — above and beyond those already well known and having broad acceptance, namely, those developed by Bakhvalov [1], Vulanovic [2], and Shishkin [3]. Their grids have been applied to diverse problems, but only to problems with exponential-type layers [3-5], typically represented by functions exp(-bx/ek) occurring in problems for which the solutions of reduced (e = 0) problems do not have singularities. Hereinafter k is the scale of a layer. The grids of Bakhvalov and Shishkin require knowledge of the constant b affecting the width of the exponential layer — when such knowledge is not always available, for example, for boundary layers in fluid-dynamics problems modelled by Navier — Stokes equations, or for interior layers in solutions to quasilinear nonautonomous problems. One spectacular example of the new layer-resolving grids being presented in the current paper, engendered by a function £rk/(ek + x)r, r > 0, is suitable for dealing not only with exponential layers having arbitrary widths, but with power of the first type layers occurring in problems for which the solutions of reduced problems have singularities as well. Another example of a new layer-resolving grid is aimed at dealing with logarithmic layers represented by a function ln(efc + x)/ln £k. It seems that the new layer-resolving grids described in this paper should empower and spark researchers to solve broader and more important classes of problems having not only exponential-, but power-, logarithmic-, and mixed-type boundary and interior layers.

By the application of algebraic methods or inverted Beltrami and diffusion equations in control metrics, the layer-resolving grids can be used for solving multidimensional problems [6].

1. Explicit generation of layer-damping transformations

This section gives a detailed description of basic layer-damping functions near the boundary point x0 = 0 which are applied to specify global layer-damping transformations and corresponding global layer-resolving grids on the entire interval of calculations with arbitrarily allocated layers, by the procedures of shifting, blending, scaling, inverting, composing, and matching them with themselves and polynomial mappings.

1.1. Basic layer-damping transformations

Local contraction transformations xb(£, £, b, k), XLi(C, £,b,k), i = 2, 3, 4, have the following form:

x„((,e,b,k) = -j ln(l - k> 0, b> 0,

xL2(i,£,b,k) = £k((1 -d0-1/b - 1), k> 0, b> 0, xL3(C ,£,b, k) = (ekb + d01/b - £k, k> 0, 1 >b> 0,

(1) (2) (3)

xLA(i,£,b, k) = £k((1 + £-k- 1), k> 0, b> 0,

(4)

where e E (0,1] is a small parameter. Differential equations with the small parameter e multiplying the highest-order derivative terms model viscous flows, where e is typically the reciprocal of the nondimensional Reynolds number Re; these equations describe problems of elasticity, where the parameter represents the shell thickness, or simulate flows of liquid in regions having orifices with a small diameter. As a rule, the solutions of these problems have highly localized regions (boundary and interior layers) of rapid variation.

The transformation xb(£,£,b,k), for k = 1, was introduced by Bakhvalov [1], while the transformations xLi(C,£,b,k),i = 2, 3, 4, were introduced by Liseikin [7-9]. A particular shape of the contraction mapping x^2(C,£,a, k) for a = 1, k = 1/2, having the form

xL2(C ,£, 1,1/2) = e

1/2

1 -dC

was proposed by Vulanovic [2] to generate grid nodes within some exponential layers of scale k =1/2.

The points i = 1, 2, 3, 4, such that the pth derivative of the mapping xb(£,£,b,k), xLi(C,£,b, k) on the corresponding interval [0, ] is e-uniformly bounded, and the points xb (k), xLi (£%,e,b, k), i = 2, 3, 4, which are the widths of the corresponding boundary layers, are described by the following equations:

8

1 - £k/p d ; 1 - £kfi

P

1+

m,

F kp

xB (ei,s,b,k) = — ln £-k,

xL2($,e,b,k)= £k(1-f/b) - £k, xL3(£, b, k) = (£kb + dm)1/b - £k

ln£-k - pln[ln(1 + £-k)] b ln(1 + £-k)

, xLA(Cl,£,b, k) =

ln1/p(1 + £-k)

k

1

Hence, for sufficiently small e, the widths of these boundary layers are connected by the following inequalities:

xB(iPi, £, b, k) < xL2(il,£, b, k) < xL4(t%,e, b, k) < xL3(t3,e, b, k).

In order to define a boundary-layer damping transformation x(£,£,b, k) for the target interval [0,m] through the use of the local univariate mappings xb(£,£,b, k), ,£,b, k), i = 2, 3, 4, from (1)-(4), specified on the corresponding intervals [0, which will provide adequate clustering of grid nodes near the boundary point x0 = 0, these mappings need to be extended continuously or smoothly over the interval [0,m1] to map it monotonically onto the interval [0,m]. This can be done by "gluing" these local nonuniform transformations xb(£,£,b, k), XLi(i,£,b, k) to other mappings which are more uniform, for example, polynomial functions. The glued transformation extending xb(£, s, b, k), , £, b, k) should be smooth, or at least continuous.

1.2. Local transformations eliminating singularities of high order

This section describes local coordinate transformations x(£, e) which eliminate singularities of arbitrary order in the boundary layer near the point x0 = 0 by specifying coefficients in the local functions (1)-(4). With the help of high-order approximations, such transformations are suitable for generating layer-resolving grids Xi = x(i/N, £), i = 0,1, N, providing high-order -uniform convergence and interpolations for numerical solutions of singularly-perturbed equations.

1.2.1. Transformations for exponential singularities

Power transformations. For a function u(x, e) whose derivatives up to n in the vicinity of the boundary point x0 = 0 (0 < x < m ) are estimated by an exponential function and M,

i. e.,

| u(p)(x, e) |<M[e-kp exp(-bx/ek) + 1], b> 0, 1 < p < n, 0 <x <m, (6) we have that

| u[p)(x, e) |<M, 1 < p <n, m >x >xï

k nek ~b~

ln(e-i)

(7)

while inside the interval [0,x™] the derivatives are not e-uniformly bounded.

In order to eliminate the exponential singularity (6) using a coordinate transformation, we rely on the basic contraction function (2) for the construction of nonuniformly clustering grids within both exponential and power (of the first type) boundary layers. This coordinate transformation, designated also as x^2(C,£,&, k) E Cl[0,m1], n >l> 0, has the following form:

f c£k((1 - d£)-1/a -

xl2(c ,£,a,k)

1), 0

£k(1-H/a) _ £k +

((1 - d01/a)

($)(£ - %)+

1

+ 77

H(1 -dî )1/a

( &)(£ - &)2 + ... +

(8)

1

+ 77

Z! V(1 -d£)1/a

(0

(e ?)(£ - £)' + co(£ -

& <C<m1,

k

k

where d = (1 — £k/3)/m0 > ££ > 0 (for example = mi/2); 0 < m2 < ft = a/(1 + na); a is a positive constant; n = t + 1, where t is the order of the numerical scheme under consideration; c> 0 is such as satisfies a necessary boundary condition x2(m\,e,a,k) = m; and

((r—d&r-)"'(S) = d a (1 +1)... (a+s—1) ^(n-'V(1+n'',• n > < > L

In the present paper we use the value 1 = 2, for the numerical solving boundary-values problems having the boundary layer near x = 0 we use mi = m = 1, while in the case of two boundary layers i.e. x = 0, x = 1 we use mi = m =1/2, to determinate the constant c. It was proved in [10] that the transformation (8) eliminates singularity (6) up to order n, i.e.,

dn

—u[xL2(C,£,a,k), £]

< M, 0 <C<mi. (9)

Logarithmic transformations. In the same manner it was proved in [10] that, in order to eliminate locally (in the vicinity of the boundary layer near xo = 0) the exponential singularity (6) of the function u(x, e) up to order n in a new coordinate we can use the basic logarithmic contraction function xb(£,£,a,k) in the form (1) on the corresponding interval [0, (see (5)):

k

xb(£, £,a,k) =--ln(1 — d£), 0 < £ < (10)

where d = (1 — £k/n)/^i, but with the restriction b/n2 > a > 0, and then prolongate it smoothly on the interval [0,mi]. The local transformation of this kind with k = 1 was introduced by Bakhvalov [1].

The transformation (8) is more convenient for eliminating exponential singularities than the transformation (10), since the constant a in (8) is not dependent on b in (6), so that, with an arbitrary fixed constant a > 0, this transformation alone is valid for all constants b G (0, to) in (6) for eliminating singularities of u(x, e) up to order n. Another common piecewise uniform transformation

f 2a£, 0 <£< 1/2,

xSh(Z,£, &)=< (11)

[ a + 2(1 — a)e, 1/2 <£< 1,

where a = min{0.5, (n/b)eln N}, proposed by Shishkin [3] for generating grids in exponential layers, is also dependent on constant b in (6), so that such a grid with a fixed constant will not be suitable for all b G (0, to) in (6). Compared with the grid of Bakhvalov, the grid of Shishkin provides less uniform accuracy.

1.2.2. Transformations for power singularities

Transformations for power singularities of the first type. The local power transformation (8) with a proper choice of constant a > 0 is also suitable for eliminating power singularities of the first type near x0 = 0, i.e., when solution derivatives are estimated by the following formula:

|u(p,(x, e)| < M [e kb/(£k + x)b+p + 1], 1 < p < n, 0 <x <m. (12)

Here, the boundary-layer interval, where all the derivatives up to n of u(x, e) are not uniformly bounded over e, is [0, xn], x'n = m2£kb/(b+n, ^ x'n = (kn/b)£k ln(e-i) for sufficiently small , so that the transformations (10) and (11) may not be suitable for generating layer-resolving grids for such singularities having incomparably wider layers than any exponential layer.

It can be proved as in [10] that the transformation (8), but with the following restrictions 3 = a/(1 + na) and 0 < a < b/n2, eliminates singularity (12) up to order n.

Transformations for logarithmic singularities. Solution derivatives near x0 = 0 can

also be estimated by

\u(p,(x, e)l < M [1 + 1/((ek + xf\ ln e\)], 1 < p < n, 0 <x < m.

(13)

Unfortunately, the transformation which would eliminate this singularity up to order n > 1 has not yet been found. The following transformation, based on (4), eliminates this singularity up to order 1 only:

x( , k) =

(

1+

1

0

t/to

1

£k ln(e-k) ln-i(e-k) + 2(ek + ln-i(e-k ))x 1

0 < < 0

(14)

^ x l^1 + ik^Fk))

( — 0) + 0( — 0)2

Co < C < mi.

k

2. Semilinear boundary-value problem

In this section we consider a semilinear boundary-value problem

— (e + rx)u" + a(x)u' + f(x,u) = 0, 0 <x< 1, u(0) = u0, u(1) = ui, (15)

with the following conditions:

0 < e < 1, r = 0, or r = 1, a(x) G Cn[0,1],

f(x, u) GCn'n+i([0,1] xR), fu(x,u) > c > 0 (16)

for (x,u) G [0,1] x R. This problem, for r = 0, modells qualitative behavior of viscous flows. For r = 1 it was considered in [11]. Solutions to this problem with small e may have boundary and interior layers of exponential and power types: they have power boundary layers of the first type near x = 0 when r = 0, a(0) = 0, a'(0) > 0; or when r = 1, a(0) < —1 [12]. Various rules for grid clustering are analyzed in this section for these very cases of powertype layers. The case of power boundary layers of the second type near interior point x = x0 was considered in [13] and [10]. Some numerical experiments with schemes of high order for solving (15) for r = 0 with various types of layers on the grids defined through (8) were carried out in [14].

2.1. Numerical algorithm

We use as an approximation of the singularly-perturbed problem (15) the standard upwind scheme on a nonuniform grid xi, i = 0,1,..., N, x0 = 0 < x\ < ... < xn = 1:

2(e + r Xi) hi + hi-i

Ui+1 Ui hi

Ui Ui-1

hi-i

+ a-(xi)

Ui+1 Ui hi

+ a+(xi)Ui , Ui 1 + f(xi, n't) = 0,

hi 1

i= 1, 2,

N 1,

u' = Uo,

U

n

U1 ,

where hi = xi+1 — xi, and a± = (a± |a|)/2. The nodes xi, i = 0,..., N, of the layer-resolving grid are obtained either explicitly by means of a transformation based on the layer-damping mappings xb(£,£,a,k), xlj(£,e,a,k), j = 2, 3, 4, described in Sect.2, namely,

xi = xb (i h,£,a, k), xi = x^j (i h,e,a,k), i = 0,1,...,N, h =1/N.

Calculations of problem (15) are conducted for various values of e: the results in the 1st — 4nd examples were carried out for the values 10- 6, 10 -14, 10 -20; in the 5th example we used the values 10-2, 10-8, 10-10; while for plotting the graphics we used the values 10-2, 10-4, 10-6, 10-8. For each of these values there are used sequences of grids with doubled numbers of grid steps: Nt = 2tNi, t = 0,1,..., where Nt is the number for the rough grid. Usually N' = 50, tmax = 5, i.e. the calculations are carried out on sequences of five grids with No = 50, N1 = 100, N2 = 200, N3 = 400, N5 = 800. The numerical solution at the ith node of the grid related to Nt, is designated by uNt, % = 0,1,..., Nt.

For estimating the accuracy of the numerical algorithm, the following characteristics are introduced:

rt£ = max |un — u^1 |, t = 0,1,... ' 0<i<Nt i 2i

and, in the case when the accurate solution u(x, e) is known,

Aut £ = max |u(xi, e) — unt I, t = 0,1, 0<i<Nt

Besides this, one more characteristic is introduced

u

t ,£

max IU.- + — u.11, 0<i<Nt-1 i+1

(17)

(18)

(19)

which is related to the jump of the numerical solution in the neighboring nodes. The characteristics rt, e, Aut,e are applied to estimate the order of the accuracy of the numerical solution:

¡31 = log2( rt,£/rt+1,£), P2 = log2(Aui;£/Aui+1;£), i = 0,1,..., (20)

and, consequently, dut,£ to estimate the order of the numerical solution jump in the neighboring nodes

¡3 = log2(dui;£/dui+M), i = 0,1,... (21)

Notice, if a solution to (15) has neither boundary nor interior layers, then for the numerical solution of this problem through the use of a stable scheme of order p on the uniform grid xi = ih the values [3\ and ¡2 are close to p, while ¡3 is close to 1. Recall, that we use the standard upwind finite difference scheme of order 1, and it is better for all three characteristics ([\, [2, [3) to have a value closer to 1. The aim of the present paper is to find out wether this property is valid for solving problems with power boundary layers of the first type by using popular grids and the grids defined through transformations (8) and (14).

3. Numerical experiments

In this section we present results obtained by applying the standard upwind finite difference scheme (17) on nonuniform grids. For Example 3.1-3.3 and 3.4 we assume r = 0, while for Example 3.5 we assume r = 1.

The analytical solutions of the first and second examples have a single power boundary layer of the first type and of scale k = 1/2 near the point x0 = 0, while the solutions of the third and fourth examples have two power boundary layers of the first type and scale 1/2, near the points x0 = 0 and x0 = 1, finally the solution of the last i.e. the fifth example has a single power boundary layer of the first type and scale 1 near the point x0 = 0.

The corresponding transformations for the first three grids according (10), (11), and [2], which will be used in Example 3.1, 3.2 and 3.5, have the forms given below. Modified Shishkin grid 1 is given by the transformation

{

x if^M^ ^ 0 ^ ^ 1/2, (22)

xsh1^,£ , 0) = 1 a + 2a(e - 1/2)+u(£ - 1/2)3, 1/2 ^^ 1, ( )

where a = min {0.5, (n/b)£k ln^} , and u is chosen so that hold xshi (1, £k, b) = 1. Bakhvalov grid is given by

xbfc, ^, a)=< ^):= -7lM1 - ^ 0 ^ ^ (23)

t(r) + 0 (r)(£ - t), T<£ < 1,

<f>(r) - 1

2 > 0, and the point r satisfies y (r) =

Vulanovic grid is given by

where q G (0, 0.5), b/n2 ^ a > 0, and the point r satisfies $(r)

1

it k \ i d(0 := aek--, 0 ^ £ ^ r, /0

xVui(i,£ ,a)= <j ^ q-g ^^ ' (24)

(t)(£ - T), T^^ 1,

where q G (0, 0.5), and the point r is calculated from condition xVui(1, £k,a) = 1.

Since the solutions of Example 3.3 and 3.4 have two boundary layers near the points x = 0 and x = 1 , we will use the grids given by the following formulas. Modified Shishkin grid 2 is given by

4a£, 0 ^ £ ^ 1/4,

xsh2(C, ek, b)={ a + 4a(^ - 1/4) + - 1/4)3, 1/4 < ^ 1/2, (25)

1 -xsh2(1 - £k, b), 1/2 < 1,

and now the parameter a is defined by a = min{1/4, (n/b))£k ln^}, and u is chosen from the condition xSh2(1/2, £k, b) = 1/2. Modified Bakhvalov grid is given by

xB (£, sk ,a) = <

m :=--ln(1 ), 0 r,

a V V (26)

+ # (r)(C - r), T<C ^ 1/2,

1 -xb (1 - £, £k ,a), 1/2 <£ ^ 1,

„ , <f>(r) - 1/2

where the parameter r is calculated from the condition ft (r) =- , and q = 1/4.

Modified VulanoviC grid is given by

r- 1/2

xvui(î, £k, a)

0(O := ae

k_^

o ^ e ^ r,

-

<f>(r) + 0(r)(Z - r), t < ^ 1/2,

1 -xvui(1 - £k, a), 1/2 <e ^ 1,

(27)

„ , 0(r) - 1/2

again, the parameter r is calculated from the condition ft(r) =--—, and q = 1/4.

- 1/2

For the plotting some graphics in Example 3.3, we will use modified Shihskin grid given by

o ^ £ ^ 1/4,

xSh(Ç, £ , b)={ a + 4(1/2 -a)(£ - 1/4), 1/4 <£ ^ 1/2, 1 -xSh(1 - £, £k, b), 1/2 <£ ^ 1,

here the parameter a is defined by a = min{1/4, (n/b))£k ln^}.

(28)

Remark 3.1. In the sequel,the grid given by (8) we designate by 1st Liseikin grid, and the other grid given by (14) we designate by 2nd Liseikin grid.

Modified 1st Liseikin grid is given by

f C1^((1 — d0-1/a — 1), 0 £o,

Cl

1

Xl2(Î,£,a, k) = <

gkan/(1+na) __ ^k + ^^£ka(n- 1)/(1+na) (£ __ £Q) +

if1 ( 1 '

2 a \ a

(H

+ -cP-[ -+ 1} £ka(n-2)/(1+aa)(^ - +

(29)

+c0(£-e0)3j, eo^e^ 1/2,

I 1 -xL2(1 - Î,£,a,k), 1/2 ^^ 1,

where d = (1 - £ka/(1+na))/= 1/4, a is a positive constant subject to a > m1 > 0, c0 > 0, and

1

— = 2 1

£kan/(1+na) - £k + ^£ka(n-1)/(1+na) + V\ + ^ £ka(n-2)/(1+na) (1/4)2 + Cq(1/4)

4 a 2 a\a J

Modified 2nd Liseikin grid is given by

x(£, efc )

1

(1 + ek ln(e-k))

Ç/Çq

1

1

ek ln(e-k) ln-1(e-k ) + 2(ek + ln-1(e-k ))x

0 ^ i ^ Co

(30)

x ln 1 +

1

ekln(e-k)

( 1 -x(1 -1, £k ), 1/2 1

(£-£o) + c0(£-Co)2 , 1/2

3

k

where = 1/4, c0 = 1, and

1

— = 2 Cl

-l(r-k\ I ^k , lI - ■ 1 I ■ °0

ln-1(£-k) + ^k + ln-1(£-k))+ ¡M^) +

The analytical solutions of 2nd, 4th and 5th examples are unknown, so we are going to use the characteristics [ and [3 given by the formulas (20) and (21), in order to investigate a behavior of the numerical solutions calculated on the different grids, but the analytical solutions of the 1st and 3rd examples we know and we are going to use the characteristics [2 and in our analysis.

Example 3.1. Let us consider semilinear boundary-value problem

-eu" - 8x(x - 1/2)V + /(x,u) = 0, u(0) = 0, u(1) = 2, (31)

where

_ 1 + £1/2 / 4e3/2 _ 16e 1/2x(x - 1/2)3

J(x,u) u 1 /2 I ( t 1 /2 I 2 1 /2 i + 2x

£l/2 + x\ (e1/2 + x)2 £1/2

x

.1/2

// £1/2 \

The analytical solution is u(x, e) = 2 ( 1--1/2 — 1 I 1--1/2 ^ I . Here we have that

fu(x, u) = 1 > 0, a(0) = 0, a'(0) = 1 > 0, a(1) = -1 < 0, a'(1/2) = 0, the scale of this problem is k = 1/2, and the solution of this problem has a single power boundary layer of the first type.

This problem we tested on modified Shihskin 1, Bakhvalov, Vulanovic, 1st and 2nd Li-seikin, and Shishkin grids. The first five grids were used to calculate the values of [2 and [3, beside previously mentioned grids also Shishkin grid was used in making some graphics. Other parameters we use: n = 2, a =1/8, c0 = 1, c1 is defined from the condition x2(1,£,a, 1/2) = 1 for 1st Liseikin grid; c0 = 0 and c is calculated from the condition x(C,£) = 1 for 2nd Liseikin grid; q = 0.45 for Bakhvalov and Vulanovic grids. The data corresponding this example are in Table 1, and the appropriate graphics are given in Fig. 1,2.

The graphics of two numerical solutions and two analytical solutions are presented in Fig. 1 (left). The first numerical solution is calculated by using 1st Liseikin grid with N = 100 and e = 10-2, while the second one is calculated also by using 1st Liseikin grid but N = 100 and e = 10-6. Every 2nd points of the graphics of the numerical solutions are shown. It's a well-known fact that analytical solutions of singularly-perturbed boundary problems exhibit rapid changes in the boundary and/or interior layer/layers. Usually those changes are getting faster when the perturbation parameter is smaller. From the graphics (Fig. 1 — left) we can see the creation of such a layer near the end point x = 0. The graphic (the red one), corresponding to the perturbation parameter e = 10-6 is narrower and we see a faster change of both solutions in the layer than in case of the graphic-blue one, corresponding to the perturbation parameter e = 10-2. Also on the same figure the valid points of 1st and 2nd Liseikin grids are presented, calculated by using two values of the perturbation parameter (10-2 and 10-6). It is easy to notice that the grid points are condensed in the portion corresponding to the layer, and that this condensation is higher in the case of grids with a smaller perturbation parameter value.

The graphics of the six numerical solutions obtained by using 1st, 2nd Liseikin, Vulanovic, Bakhvalov, Shishkin grids and modified Shishkin grid 1 are presented in Fig. 1 (right). The

Table 1. The results of Example 3.1

Modified

N Shishkin 1 Bakhvalov Vulanovic 1st Liseikin 2nd Liseikin

Ä Â Ä Â Ä Ä Â Ä

£ = 10- -6

50 1.1361 0.6637 0.7468 0.6459 0.8746 0.8701 0.9709 0.9680 1.0193 0.9169

100 1.1454 0.7361 0.5383 0.4209 0.9309 0.9387 0.9849 0.9840 1.0090 0.9602

200 1.0374 0.7863 0.7334 0.4599 0.9642 0.9418 0.9926 0.9917 1.0043 0.9805

400 0.8617 0.8211 0.8317 0.5465 0.9815 0.9697 0.9963 0.9959 1.0022 0.9903

800 0.8688 0.8458 0.9184 0.7495 0.9871 0.9893 0.9981 0.9979 1.0011 0.9952

£ = 10" 14

50 0.3452 0.5710 0.6457 0.3448 0.0952 0.4975 0.9529 0.9413 1.0593 0.7425

100 0.9459 0.5078 0.3271 0.2533 0.6880 0.6875 0.9729 0.9711 1.0249 0.8884

200 0.6648 0.5392 0.1485 0.2328 0.7873 0.7866 0.9873 0.9847 1.0101 0.9476

400 0.9594 0.9042 0.1158 0.2170 0.8614 0.8607 0.9938 0.9927 1.0059 0.9745

800 1.0561 0.9718 0.1043 0.2054 0.9188 0.8984 0.9968 0.9963 1.0028 0.9874

£ = 10" 20

50 0.2122 0.5295 0.4402 0.2741 0.0905 0.3967 0.9390 0.9300 1.0975 0.5790

100 0.1837 0.1829 0.2836 0.2345 0.3857 0.5857 0.9711 0.9626 1.0378 0.8264

200 0.1706 0.1625 0.2641 0.2121 0.4831 0.6831 0.9855 0.9826 1.0135 0.9217

400 0.2370 0.1691 0.2483 0.1938 0.5536 0.7535 0.9928 0.9915 1.0040 0.9626

800 0.8333 0.3607 0.2082 0.1791 0.6040 0.8040 0.9964 0.9958 1.0040 0.9816

Analytical solution for e = 10~2 Analytical solution for e = 10"6 Numerical solution for £ = 10~2, Numerical solution for e = 10~6 Points of 1st Liseikin grid for £ = 10~2 Points of 1st Liseikin grid for £ = 10~6

0

1»«»» «-« -» K- K- <t- 4»

Î.4 0.6

x - axis

v?

J

•V

*'V

iV ».* <

.V f/ if

•V J j j ■ j «

- >

/>*•;' / v ( *

r / /

r / *f *

<

•«

« Num solution on Shishkin grid

-*■ ■ Num solution on mod. Shishkin grid 1

■ Num solution on Bakhvalov grid ■f ■ Num solution on Vulanovic grid

■ Num solution on 2nd Liseikin grid * ■ Num solution on 1st Liseikin grid

■ Points of uniform grid

0 0 0 2 0 4 0 6 0 8 1 0

Ç - axis

Fig. 1. The graphics of the analytical and numerical solutions on 1st Liseikin grid N = 100, e = 10-2,10-6 (left), the graphics of the numerical solutions on all considered grids but plotted by using the uniform grid N = 100, e = 10-4 (right)

last four grids are constructed to resolve exponential layers, 1st Liseikin grid is constructed to resolve a power layer of the first type, and 2nd Liseikin grid has a goal to resolve a logarithm layer. In short the purpose of using layer grids to resolve layers is to reduce the distance between the grid points in parts of the domain where the analytical solution changes rapidly. Reasons for reducing the distances have been analyzed many times in the literature by Bakhvalov, Liseikin, Gartland, Shishkin, and others. In constructing layer-

resolving grids usually we use simple-model problems exhibiting an appropriate layer with known analytical solutions, and inverse functions corresponding these analytical solutions. Using these inverse functions we generate the grid points in the layers. Taking into the account previously written, in the case when the numerical solution was calculated by using a layer-resolving grid but plotted by using the uniform grid, we expect that the graphic of this numerical solution corresponding to the layer, be close to the graphic of the linear function.

The graphics plotted on a described way, N = 100, e = 10-4, and every 2nd points are presented in Fig. 1 (right). The grid better constructed for a specific problem, gives the graph that is closer to the graphic of linear function in the part of grid corresponding the layer. We used in this example 0.5 N of the grid points to resolve the layers for 1st, 2nd Liseikin, Shishkin grids and modified Shishkin grid 1 and slightly less than 0.5 N of the grid points in the same purpose for Bakhvalov and Vulanovic grids. From Fig. 1 (right), we see that graphic closest to the graphic of linear function, is the graphic obtained by using 1st Liseikin grid. The graphics obtained by using Shiskin, modified Shishkin 1, Bakhvalov and Vulanovic grids suddenly change in the part corresponding the layer, this kind of behavior we can explain by the fact that these four grids were constructed to resolve a "sharper" exponential layer. The graphic obtained by using 2nd Liseikin grid also deviates from the graphic of linear function but this graphic doesn't have abrupt changes.

Figure 2 (left) shows the graphics of the parts of the numerical solutions and the graphics of the part of the analytical solutions near the end point x = 0, of the problem given in Example 3.1. To calculate the numerical solutions we used N = 100, e = 10-8, and every 2nd point are plotted. Each graph of the numerical solutions (except the first one) is moved up 1 cm from the bottom one in order to better transparency. In the same way the graphic of the analytical solution is shifted and plotted, i. e. the graphics of the numerical and analytical solutions are plotted together in order to compare. As mentioned a few time before, Bakhvalov, Shishkin, Vulanovic grids and their variations were constructed to resolve exponential layers, a well-known fact is that a exponential layer comparing to a power layer of the first type is "sharper". The grid points intended for a layer are too condensed near the end point, and they resolve just a part of the layer. The rest of the layer is not divided by a sufficient number of grid points, this results in higher error values. Grouping of the grid points near the end point x = 0 is especially evident in Bakhvalov and Shishkin grids, the situation with Vulanovic and modified Shishkin 1 grids is less worse, while the points of 1st and 2nd Liseikin grids are much better distributed (Fig. 2 — left). These are the expected results because Liseikin grids were constructed for boundary-value problems having wider layers, and don't have gaps in the distribution of the graphics points, unlike the four grids mentioned previously.

Figure 2 (right) shows the graphics of the part of the numerical solutions obtained by using the uniform and 1st Liseikin grids, and the graphic of the analytical solutions. We may see that the error value is larger in case of using the uniform grid for this semilinear boundary-value problem. This example demonstrates the justification for using layer-resolving grids over uniform grids for the numerical solving of problems of these type.

In Table 1 the results of Example 3.1 are presented, i.e. the values of [2 and [3. In our experiments we use the standard upwind finite difference scheme. Let us recall that the order of accuracy of this scheme on a uniform grid for the numerical solving of boundary-value problems without boundary or interior layers is 1. We expect that ¡2 and ¡3 have a value of approximately 1. We performed calculations of the numerical solutions for

Anal t ical solution for e = 10~8 solution on 1st Liseikin grid

Num.

■ ■» ■ Num. solution on 2nd Liseikin grid ■ * • Num. solution on Vulanovic grid ■ * ■ Num. solution on Bakhvalov grid ■ * ■ Num. solution on mod. Shishkin grid 1 ■ * ■ Num. solution on Shishkin grid

0.0015 ( - axis

r

f ............

t*

/ ■ » ■ Num. solution on uniform grid ■ ■» ■ Num. solution on 1st Liseikin grid -1-I-

Fig. 2. The graphics of the parts of the numerical solutions on all considered grids N = 500, e = 10-8 (left), the graphic of the parts of the numerical solutions on 1st Liseikin and the uniform grids N = 200, e = 10-8 and the graphic of the part of the analytical solution e = 10-8 (right)

N = 50, 100, 200, 400, 800, and e = 10-6, 10-14, 10-20. From Table 1 we see that only the 1st and the 2nd Liseikin grids give satisfactory values for and 03 for all used values of N and e. VulanoviC grid gives also the satisfactory values of and for e = 10-6, but decreasing the value of the perturbation parameter e, VulanoviC grid becomes useless for the numerical solving of this boundary-value problem. Also the modified Shishkin 1 and Bakhvalov grids give unsatisfactory values of and ft3.

Example 3.2. Let us consider the following problem

-eu" - 8x(x - 1/2)3 u' + 1/2« = exp(x), 0 <x< 1, u(0) = 0, u(1) = 1. (32)

We don't know the analytical solution of this problem, and here we calculated (31 instead fi2. The results are very similar to the previous one, and here we gave only the results for the smallest value of e in Table 2.

Example 3.3. Let us consider the following problem

-eu" + 8x(x - 1/2)3(x - 1) + / (x,u) = 0, 0 <x< 1, u(0) = 2, u(1) = 0,

Table 2. The results of Example 3.2

Modified

N Shishkin 1 Bakhvalov Vulanovic 1st Liseikin 2nd Liseikin

A Â A Ä A ßi Ä ßi A

£ = 10- -20

50 0.3211 0.8908 0.1861 0.9266 0.6048 0.9266 0.8332 0.8508 0.6776 0.8444

100 0.2257 0.9085 0.1733 0.5925 0.5732 0.9526 0.9153 0.9490 0.7684 0.9601

200 0.2807 0.9554 0.1469 0.1577 0.5468 0.9802 0.9573 0.9725 0.9100 0.9534

400 0.5612 0.6085 0.1222 0.1411 0.5288 0.9898 0.9789 0.9856 0.9643 0.9587

800 0.2936 0.4968 0.1020 0.1297 0.5181 0.8107 0.9899 0.9926 0.9845 0.9775

where

f {x,u) = u -<| 1 - 2e3/2(1+ £1/2) 1

1

1

x(x - 1/2)3

(£1/2 + X)3 (£1/2 + 1 - X)3 1

+

(£1/2 + X)2 (£l/2 + 1 - X)2

+ £l/2 (1+ £i/2)

1/2\

(

- 8ei/2(1 + £i/2)x(x - 1)x 11

£i/2 + X £i/2 + 1- X

The analytical solution of this problem is

i/2

i/2

u(x, e) = 1 +

£i/2 + x £i/2 + 1 - x

i/2

1

£i/2 + 1

Here we have fu(x,u) = 1 > 0, a(0) = 0, a'(0) = 1 > 0, a(1) = 0, a'(1) = 1 > 0 and a(1/2) = 0, a'(1/2) = 0 ^ 0, the scale of this problem is k = 1/2, and the analytical solution of this problem has two power boundary layers of the first type near the points x = 0, x = 1.

In our calculation we use: n = 2, a = 1/8 and c0 = 1, the parameter c1 is calculated from the condition XL2(1/2,e,a,k) = 1/2 for 1st Liseikin grid; the parameter c is calculated from the condition x(1/2,e) = 1/2 for 2nd Liseikin grid; and q =1/4 for Bakhvalov and Vulanovic grids. The data corresponding this example is presented in Table 3.

In Figure 3 the parts of the numerical solutions as well as the parts of the analytical solution are presented. The five graphics of the parts of the numerical and analytical solutions have been moved up for better transparency. The points of the numerical solutions obtained by using mod. Shishkin, mod. Shishkin 2, mod. Bakhvalov and mod. Vulanovic grids are highly concentrated near the ends points x = 0, and x =1. We can explain this phenomenon very easy, these grids are constructed to resolving an exponential layer, not a power layer of the first type.

Table 3 shows the values of and . Because the analytical solution is known, we calculated the value of rather than ¡31. Also, here the results are better whenever values

- Analytical solution for £ = 10~8 ■ ■» ■ Num. solution on mod. 1st Liseikin grid ■ ■ Num. solution on mod. 2nd Liseikin grid ■ -f ■ Num. solution on mod. Vulanovic grid ■ ■ Num. solution on mod. Bakhvalov grid ■ ■ Num. solution on mod. Shishkin grid 2 ■ ■ Num. solution on mod. Shishkin grid

<-

>

- Analytical solution for £ = 10"8 ■ * ■ Num. solution on mod. 1st Liseikin grid ■ ■» ■ Num. solution on mod. 2nd Liseikin grid -^

■ ■» ■ Num. solution on mod. Vulanovic grid ■ ■ Num. solution on mod. Bakhvalov grid ■ f ■ Num. solution on mod. Shishkin grid 2 ■ * ■ Num. solution on mod. Shishkin grid

0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0.0040 0.9960 x - axis

0.9965 0.9970 0.9975 0.9980 0.9985 x - axis

0.9990 0.9995 1.0000

Fig. 3. The graphics of the parts of the numerical and analytical solutions N = 500, e = 10 8

Table 3. The results of Example 3.3

Modified

N Shishkin 1 Bakhvalov VulanoviC 1st Liseikin 2nd Liseikin

Ä Ä Ä Ä Ä Ä Ä Ä

£ = 10- -6

50 1.1284 0.5678 0.7682 0.6107 0.5977 0.6766 0.9358 0.9214 1.0575 0.7839

100 0.9543 0.6754 0.4918 0.4375 0.7515 0.7220 0.9642 0.9632 1.0247 0.9013

200 0.8885 0.7501 0.6869 0.4697 0.8940 0.8813 0.9830 0.9818 1.0129 0.9527

400 0.8808 0.8003 0.8288 0.5432 0.9440 0.9178 0.9917 0.9909 1.0060 0.9768

800 0.8798 0.8340 0.8991 0.7198 0.9574 0.9245 0.9959 0.9955 1.0031 0.9885

£ = 10" 14

50 0.2940 0.5677 0.5529 0.3423 0.4567 0.4111 0.8578 0.8463 1.0591 0.6651

100 0.4823 0.6553 0.3460 0.2761 0.5712 0.5712 0.9352 0.9336 1.0909 0.7098

200 0.8335 0.7501 0.2778 0.2521 0.7027 0.7026 0.9710 0.9660 1.0362 0.8485

400 0.8815 0.8003 0.1290 0.2327 0.7954 0.7953 0.9858 0.9831 1.0211 0.9284

800 0.9722 0.8340 0.1158 0.2179 0.8000 0.7942 0.9923 0.9878 1.0103 0.9651

£ = 10" 20

50 0.2872 0.5677 0.3711 0.2827 0.4421 0.4001 0.8408 0.8470 1.0128 0.6518

100 0.2013 0.6753 0.3042 0.2575 0.5671 0.5701 0.9418 0.9249 1.1255 0.7082

200 0.1729 0.3751 0.2825 0.2317 0.7002 0.6904 0.9643 0.9588 1.0639 0.7609

400 0.1719 0.1546 0.2638 0.2100 0.7532 0.7894 0.9835 0.9588 1.0289 0.8901

800 0.3322 0.1914 0.2487 0.1924 0.7986 0.7900 0.9908 0.9808 1.0156 0.9472

of [2 and [3 are closer to 1. It's also desirable that values of [2 and [3 are getting closer and closer to 1, when the number of points are increasing. Based on the presented results in Table 3 and just mentioned above, we obtained the best results of [2 and [3 by using 1st and 2nd Liseikin grids.

Example 3.4. Let us consider the following problem

-eu" + 8x(x - 1/2)3(x - 1)u' + 0.5« = sin x, 0 <x< 1, u(0) = 1, u(1) = -1. (33)

We don't know the analytical solution of this problem, and here we calculated the values of [\ instead [2. The results are very similar to the previous one, and we here gave only the results for the smallest value of £ in Table 4.

Table 4. The results of Example 3.4

Modified

N Shishkin 1 Bakhvalov VulanoviC 1st Liseikin 2nd Liseikin

ßi Ä ßi Ä ßi Ä ßi Ä ßi Ä

£ = 10" 20

50 0.3567 0.3241 0.8155 0.3031 0.7846 0.5462 0.8292 0.7621 0.8269 0.6019

100 0.3032 0.5049 0.1689 0.1848 0.5267 0.5188 0.9088 0.8666 1.1309 0.6882

200 0.2477 0.6481 0.1494 0.1637 0.5253 0.5156 0.9496 0.9350 1.2372 0.8643

400 0.3745 0.7426 0.1284 0.1479 0.5101 0.5118 0.9758 0.9652 1.0861 0.9722

800 0.5150 0.4358 0.1097 0.1361 0.5149 0.5092 0.9876 0.9829 1.0253 1.0012

Example 3.5. Let us consider the following problem

-(e + x)au" - (1.2)« + u = - sin(10x), 0 < x < 1, u(0) = 0, u(1) = 1,

where a =1. Here is fu(x,u) = 1 > 0, -a(0) = 1.2 > 1, the scale of this problem is k = 1, and the analytical solution of this problem has a single power boundary layer power of the first type near x = 0. The parameters we used have the values: n = 2, a = 1/20, c0 = 1, and c\

Table 5. The results of Example 3.5

Modified

N Shishkin 1 Bakhvalov Vulanovic 1st Liseikin 2nd Liseikin

A Â A Ä A A A A A A

£ = 10- -2

50 0.2134 0.8906 0.8762 1.0182 0.8186 1.0214 0.8190 1.0086 0.8648 1.0077

100 0.1665 0.6049 0.9383 1.0064 0.9094 1.0178 0.9069 1.0049 0.9323 1.0043

200 0.1416 0.6974 0.9696 1.0071 0.9552 1.0079 0.9534 1.0026 0.9664 1.0024

400 0.1303 0.7654 0.9850 1.0022 0.9778 1.0064 0.9767 1.0013 0.9833 1.0011

800 0.1258 0.8128 0.9926 1.0019 0.9889 1.0019 0.9884 1.0006 0.9917 1.0006

£ = 10- -8

50 -0.035 1.0127 0.0281 0.3130 0.0221 0.3239 0.6543 1.0201 0.6333 1.0228

100 -0.049 1.0079 0.0532 0.2275 0.0347 0.2492 0.6936 1.0167 0.8321 1.0157

200 -0.113 0.9575 0.0577 0.1725 -0.003 0.2166 0.7543 1.0112 0.8898 1.0067

400 -0.053 0.2885 0.0543 0.1352 -0.135 0.2274 1.0465 1.0043 1.0237 1.0019

800 0.3328 0.5200 0.0330 0.1086 -0.356 0.3082 1.0902 1.0014 1.0311 1.0006

£ = 10" 10

50 -0.028 1.0127 0.0282 0.3130 0.0277 0.3141 0.6148 1.0187 0.6145 1.0218

100 -0.010 1.0079 0.0534 0.2274 0.0523 0.2296 0.6368 0.8017 0.7957 1.0182

200 -0.008 0.9264 0.0581 0.1723 0.0555 0.1766 0.7267 1.2263 0.8268 1.0096

400 -0.030 0.1492 0.0551 0.1348 0.0481 0.1434 1.0414 1.0071 1.0742 1.0025

800 -0.105 0.1141 0.0341 0.1079 0.0253 0.1252 1.2276 1.0020 1.0791 1.0006

1.0 ■■»■ num. solution on 1st Liseikin grid

■ + ■ num. solution on 2nd Liseikin grid

■ ■» ■ points of 1st Liseikin grid

0.8 points of 2nd Liseikin grid „

«««•» N- :

»«■It«-»

■ «■»■«■«■•■■»■■«■■«■■»■■l—»■■»■•»— It—

0.4 0.6

x - axis

* num. solution on Shishkin grid

■* ■ num. solution on mod. Shishkin grid 1

■ num. solution on Bakhvalov grid num. solution on Vulanovic grid

* ■ num. solution on 2nd Liseikin grid

* ■ num. solution on 1st Liseikin grid

* ■ points of uniform grid

Â

¿7/ * *

* .* S t

f I» .

rasan^-'

tm#iMt###ii * (t#####t m * it###iiitii M ü###n##nit###!tiMMMHi##tittiM»

0.4 0.6

Ç - axis

Fig. 4. The graphics of numerical solutions N = 100, e = 10-4 (left), the graphics of numerical solutions obtained by using all considered grids, but plotted using the uniform grid (right), for both figures the parameters have the values N = 100, e = 10-4

is calculated from the condition XL2(1/2,e,a,k) = 1/2 for 1st Liseikin grid; c is calculated from the condition x(1/2,e) = 1/2 for 2nd Liseikin grid, and q = 0.45 for Bakhvalov and Vulanovic grids. The data of this example are in Table 5, and the appropriate graphics are given in Fig. 4.

Based on the data given in Table 5, we can conclude that the numerical solutions calculated on 1st and 2nd Liseikin grids are better than other numerical solutions we calculated. The scale of the layer in this example is 1, and from Fig. 4 (left) we can notice a faster change of the numerical solutions in the layer, comparing by the previous examples for the same value of the parameter e. Bearing in mind our earlier remarks, now from Fig. 4 (right) we can see that the numerical solutions calculated on 1st and 2nd Liseikin grids are better than the rest.

Conclusion

In this paper we have calculated the numerical solutions for one-dimensional singularly-perturbed problems having power boundary layers of the first type by using different layer-resolving grids. Our goal is to compare the new layer-resolving grids with well-known grids and show a benefit of using the new layer-resolving grids over to other mentioned grids. In order to do this, we tested five different examples on six layer-resolving grids. All the results are presented in the tables. The expected value for all of the examined characteristics (i. e. , @2, ) is 1. From the tables we can see that the weakest results were demonstrated by using Shishkin's grid (1 or 2), and that the best results were obtained by using the new layers-resolving grids proposed by Liseikin. These properties of grids are mainly manifested for the cases of very small values of the parameter e.

Acknowledgements. The reported study was funded by RFBR, project No. 20-01-00231. References

[1] Bakhvalov N.S. On the optimization of the methods for solving boundary value problems in the presence of a boundary layer. USSR Computational Mathematics and Mathematical Physics. 1969; 9(4):139-166.

[2] VulanoviC R. Mesh construction for numerical solution of a type of singular perturbation problems. Numer. Meth. Approx. Theory; 1984:137-142.

[3] Miller J.J.K., O'Riordan E., Shishkin G.I. Fitted numerical methods for singular perturbation problems. Singapore, New Jersey, London, Hong Kong: World Scientific; 2012: 191.

[4] Roos H.-G., Stynes M., Tobiska L. Numerical methods for singularly perturbed differential equations. Convection-Diffusion and Flow Problems. New York: Springer; 2010: 604.

[5] Linss T. Layer-adapted meshes for reaction-convection-diffusion problems. Berlin: SpringerVerlag; 2010: 320.

[6] Liseikin V.D. Grid generation methods. Third ed. Berlin: Springer; 2017: 530.

[7] Liseikin V.D., Yanenko N.N. On the numerical solution of equations with interior and exterior layers on a nonuniform mesh. BAIL III. Proc. 3th Internarional Conference on Boundary and Interior Layers-Computational and Asymptotic Methods. 1984, Dublin, Ireland: Trinity College; 1984: 68-80.

[8] Liseikin V.D. Numerical solution of equations with power boundary layer. USSR Computational Mathematics and Mathematical Physics. 1986; 26(6):133-139.

[9] Liseikin V.D. On the numerical solution of singularly perturbed equations with a turning point. J. Comput. Maths. Math. Phys. 1984; 24(6):135-139.

[10] Liseikin V.D. Grid generation for problems with boundary and interior layers. Novosibirsk: NGU; 2018: 296. (In Russ.)

[11] Polubarinova-Kochina, P.Ya. Theory of motion of phreatic water. Moscow: Nauka; 1977:665. (In Russ.)

[12] Liseikin V.D. Layer resolving grids and transformations for singular perturbation problems. Utrecht: VSP; 2001: 284.

[13] Becher S. FEM-analysis on graded meshes for turning point problems exhibiting an interior layer. 2016:1-17. Available at: https://arxiv.org/abs/1603.04653

[14] Liseikin V.D., Paasonen V.I. Compact difference schemes and layer-resolving grids for numerical modeling of problems with boundary and interior layers. Numerical Analysis and Applications. 2019; 12(1):1-17.

Вычислительные технологии, 2020, том 25, № 1, с. 49-65. © ИВТ СО РАН, 2020 ISSN 1560-7534

Computational Technologies, 2020, vol. 25, no. 1, pp. 49-65. © ICT SB RAS, 2020 eISSN 2313-691X

ВЫЧИСЛИТЕЛЬНЫЕ ТЕХНОЛОГИИ

DOI:10.25743/ICT.2020.25.1.004

Численный анализ законов сгущения сеток для задач со степенными погранслоями первого типа

В. Д. ЛИСЕЙКИН1'2, С. КАРАСУЛИЧ3

Институт вычислительных технологий СО РАН, Новосибирск, Россия

2

Новосибирский государственный университет, Новосибирск, Россия

3

Университет Тузлы, 75000 Тузла, Босния и Герцеговина Контактный автор: Лисейкин Владимир Д., e-mail: [email protected]

Поступила 10 июня 2019 г., доработана 2 декабря 2019 г., принята в печать 18 декабря 2019 г.

Аннотация

В статье приведены результаты численных расчетов обыкновенных сингулярно-возмущенных задач, решения которых имеют степенные, первого типа, пограничные слои. Расчеты проведены с использованием как известных адаптивных сеток, сгущающихся в слоях, так и новых. Численные эксперименты демонстрируют преимущество новых сеток.

Ключевые слова: уравнение с малым параметром, погранслой, адаптивная сетка. Цитирование: Лисейкин В.Д., Карасулич С. Численный анализ законов сгущения сеток для задач со степенными погранслоями первого типа. Вычислительные технологии. 2020; 25(1):49-65. (На англ. яз.)

Благодарности. Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 20-01-00231.