NUMERICAL METHOD FOR SYSTEM OF SPACE-FRACTIONAL EQUATIONS OF SUPERDIFFUSION TYPE WITH DELAY AND NEUMANN BOUNDARY CONDITIONS Текст научной статьи по специальности «Математика»

Izvestiya Instituta Matematiki i Informatiki Udmurtskogo Gosudarstvennogo Universiteta

2022. Volume 59. Pp. 41-54

MSC2020: 65M06, 65M12, 65M15 © M. Ibrahim, V. G. Pimenov

NUMERICAL METHOD FOR SYSTEM OF SPACE-FRACTIONAL EQUATIONS OF SUPERDIFFUSION TYPE WITH DELAY AND NEUMANN BOUNDARY CONDITIONS

We consider a system of two space-fractional superdiffUsion equations with functional general delay and Neumann boundary conditions. For this problem, an analogue of the Crank-Nicolson method is constructed, based on the shifted Griinwald-Letnikov formulas for approximating fractional Riesz derivatives with respect to a spatial variable and using piecewise linear interpolation of discrete prehistory with extrapolation by continuation to take into account the delay effect. With the help of the Gershgorin theorem, the solvability of the difference scheme and its stability are proved. The order of convergence of the method is obtained. The results of numerical experiments are presented.

Keywords: superdiffusion equations, Neumann conditions, functional delay, Riesz derivatives, Grunwald-Letnikov approximation, Crank-Nicholson method, order of convergence.

DOI: 10.35634/2226-3594-2022-59-04

Introduction

Systems of ordinary differential equations with delay effect, including distributed ones, occur in many mathematical models, for example, predator-prey [1]. Modern mathematical models use, among other things, systems of partial differential equations with delay, including those with Neumann boundary conditions [2]. In recent years systems of fractional differential equations [3, 4] have become more and more popular. These papers also use the Neumann boundary conditions.

Due to the complexity of analytical studies of fractional equations, the development of adequate numerical methods becomes relevant. Numerical algorithms for solving equations, both time-fractional and space-fractional, are the subject of a huge number of works, we note some [5-14]. But they consider the Dirichlet boundary conditions.

In this work, we use numerical algorithms developed and studied earlier for the Dirichlet problem for partial differential equations with delay [15] and for diffusion equations with a fractional derivative with respect to space without delay [17].

§ 1. Problem definition

Let us consider a system of equations of the superdiffusion type with fractional Riesz derivatives with respect to space and with a functional delay of the form

du(x, t) dau ~

o,' = + f(x,t,u(x,t),ut(x, ■),v(x,t),vt(x, •)),

dv(x,t) ^ dav - '

— = + j{x,t,u(x,t),ut{x, •)),

at d |x|a

where 0 ^ t ^ T, 0 ^ x ^ X are independent variables, u(x,t), v(x,t) are the required functions, ut(x, ■) = {u(x,t + s), t ^ s < 0} and vt(x, ■) = {v(x,t + s), t ^ s < 0} are prehistories of desired functions by the time t,T > 0 is the value of delay, 1 < a < 2, Du > 0, Dv > 0.

Riesz derivatives of order a are defined by the relations [16, p. 3]

^rfda u(x,t) da u(x,t)\ 1

- K \ n H—0 , K = —

(1.2)

d|x|a V xa d- xa J' 2 cos(an/2)'

where the left and right Riemann-Liouville partial derivatives of order a are defined respectively as

dau(x,t) _ 1 d2 fx u(£,t)

~ r(2 - a) dx2 X (x-e)""1

dau(x,t) _ 1 d2 rx u(£,t)

d_xa ~ r(2 - a) dx2 Jx (x-O"'1

^ . . dav(x,t) , dav(x,t) , „ Derivatives ———, —--and —--are defined similarly.

d|x|a d+xa 5-xa

Initial conditions are given

u(x, t) = p(x, t), 0 ^ x ^ X, -T ^ t ^ 0, v(x, t) = ^(x, t), 0 ^ x ^ X, -T ^ t ^ 0.

Homogeneous boundary conditions of the second type (Neumann conditions) are also set

= ^Uv = 0, OSiST, dx dx (1 3)

dv(x.t), dv(x,t). ^ '

dx dx

We assume that the solution to the problem (1.1)-(1.3) exists and is unique. Moreover, when proving the convergence of numerical algorithms, we will assume the necessary smoothness of the solution u(x, t), v(x, t).

Denote by Q = Q[—t, 0) the set of functions w(s), piecewise continuous to [—t, 0) with a finite number of discontinuity points of the first kind, at discontinuity points continuous on the right. We define the norm of functions on Q ||w(-)||q = supse[-r0) |w(s)|. Additionally, we will assume that the functionals /(x, t, w, v(-)) and /(x, t, w, v(-)) are defined on Q x [0, T] x R x Q and are Lipschitz in the last two arguments, i.e., there exists a constant Lf, such that for all x G Q, t G [0, T], w1 G R, w2 G R, w1^) G Q, w2(-) G Q the following inequalities are satisfied:

|/~(x,t,w1,w1(-)) — /(x,t,w2,w2(-))| ^ Lf(Iw1 — w2| + ||w1(-) — w2(-)||q), |/(x,t,w1,w1(-)) — /(x,t,w2,w2(-))| ^ Lf (Iw1 — w2| + ||w1(-) — w2(-)||q).

We introduce the vector notation

U = [UY F = [/Y D =

/' VD

then the system (1.1) can be written as

where A o U denotes the vector with coordinates ( a1U ), if A =

a2v / \ a2

§ 2. Difference scheme

1 - Mo ' ÎO a "11U 1V± - LA

Introduce the time step A = where M0 is a natural number and let M = [J]. Introduce

points tj = jA, j = -M0,...,M. Let us divide the segment [0,X] into parts with a step h = X/N, N is an integer, N > 2, by introducing the points x = ih, i = 0,... ,N.

The approximation of the vector function U(x,tj) at the grid nodes will be denoted by the

fuj;

vector Vj with coordinates i

1 w

For every fixed i = 0,..., N we introduce a discrete prehistory to the moment tm, m = 0,...,M: {Vj }m = {Vj ,m - M0 ^ j ^ m}.

The operator interpolation (with extrapolation by half a step) of a discrete prehistory is a mapping I: that associates a discrete prehistory {Vj}m with a vector function Vm(t)i defined on

[tm tm

We will say that the interpolation operator has the order of error p on the exact solution U(x,t), if there are constants C\ and C2 such that for all i, m and / G [tm — r,tm + -f-] the following inequality holds:

IIVm(t)i - U(Xi,t)\\ ^ C1 max \\Vj - U(x%,tj)|| + C2Ap.

m-M o

Here and in what follows, the norm in two-dimensional space is determined by the relation

11U || = max{|u|, |v|}.

In what follows, for the methods under consideration, we will use piecewise linear interpolation

Vm{t)i = U{tj-t)V?-1 + {t-tj-1)V?), tj-^t^tj, j^m, (2.1)

with extrapolation by continuation

Vm(t)t = i((im - t)Vr~l + (t- tm_i)VT), + (2.2)

This interpolation operator is of second order if the exact solution U(x, t) is twice continuously differentiable with respect to t [18, p. 98, 102].

To approximate the left and right fractional Riemann-Liouville derivatives, we will use the shifted Grunwald formulas

+ k=0

+ Ql, (2.4)

k=0

where the normalized weights are defined as

i / ^k(.a)(.a~ l)---(Q!-fc+l)

9a, 0 = 1, 9a,k = {-1) --, k= 1,2,3,....

Note some properties of normalized weights [17]

ro m

9a,i = -a; 1 > ga,2 > 9a, 3 > ...> 0; "^2,ga,k = 0; 0 <^ga,k < 1, m > 1.

k=0 k=0

If the exact solution U(x, t) is four times continuously differentiable with respect to x, then [16, p. 56]

||Rj|| ^ Cha, IIQj|| ^ Ch.

From (2.3)-(2.4) for the Riesz derivative we get the representation

d«u(x t ) K i+1 N-i+1

Q Z 3 = -^(Y,9«,kU(xi-k+i,tj)+ Y, 9a,kU(xi+k-i,tj)) + P-. (2.5) 1 1 fc=0 fc=0

If the exact solution U(x, t) is four times continuously differentiable with respect to x, then

l|Pj || ^ Ch. (2.6)

Let us discretize (1.4) at the nodes (x*, tm+1/2), applying a two-site approximation to the middle for the time derivative, using the shifted formulas (2.5) for the Riesz derivative with respect to space on the m-th and m + 1-th time layers and using piecewise linear interpolation (with extrapolation by half a step) of the prehistory of the discrete model, we obtain an analogue of the Crank-Nicolson scheme

Vm+1 - Vm K

rm

9a,s^ + ga,sv^ +

A 2ha .

s=0 s=0

N-i+1 N-i+1

+ £ + £ + F"

2 m + ^

(2.7)

s=0 s=0

^ = + ¿= 1, m = 0,..., M — 1

z 2

where Vm(im+i)i is the result of piecewise linear interpolation (2.1) with extrapolation by continuation (2.2) at the point tm + V"' (■)+ is the history of interpolation with extrapolation at this point.

Using the formulas for numerical differentiation by the boundary and the Neumann boundary conditions (1.3), we supplement the scheme (2.7) with the equalities

VT+1 = vr+1, V™+1 = Vmi, V0m = vr, VNm = VN-1, (2.8)

then the scheme (2.7) will take the form

vm+1 — vm K / % %

1 A * 9*,sV™s+l + ga,+1vr + £ ga,sV™tl! +

s=0 s=0

N-i N-i

+ ga,i+1 Vim+1 + £ ga,sVims-1 + ga,N-i+1 V^-i + £ ga^TS-1! +

s=0 s=0

"2

+ ^.Jv-i+iV^!1) + , i = 1,..., iV - 1, m = 0,..., M — 1. Scheme (2.9) is completed with initial conditions from (1.2)

= = ^ , i = 0,..., N, j = -M00. (2.10)

The scheme (2.9-2.10) represents, for each fixed m, two systems of linear algebraic equations of order N — 1. Let us consider the question of their solvability and stability of the method.

§ 3. Solvability and stability of the method

We introduce vectors of size 2(N - 1)

( vr \

Vn-iJ

F m

Vn -1 /

and matrices A, A, A = A + A of size (N — 1) x (N — 1), whose elements satisfy the relations

Ali

-2ga,1, S = i,

-(9a,2 + 9a,o), S = i - 1,

(ga,0 + 9a,2) , S = i +1,

9a,i-s+1, S<i - 1,

9a,s-i+1, S>i +1,

~ga,i+1,

S = 1,

Al,,

-9a,N-i+1, S = N - 1, 0, 1 < S < N - 1,

then (2.9) can be rewritten as

vm+1 - Vm K ,

Do — {- AVm -

A

2ha

+ Fm

or

K A K A

ym+l +D<> -= yra _ D -Aym + _

2 ha 2 ha

L e m m a 3.1. The eigenvalues of the matrix A have positive real parts.

(3.1)

Proof. Let us show the diagonal dominance of the matrix A = Ai + A. Its elements A, satisfy the relations

N-1

s=2

N-2

lA1,1| - > ! lA1,s| = I - 29a,1 - 9a,21 - I - (9a,0 + 9a,2)| - | - 9a,s | - | - 9a,N-1 - 9a,N | =

s=3

N-2 N

-29a,1 - 9a,2 - 9a,0 - 9a,2 - ^^ 9a,s - 9a,N-1 - 9a,N = -9a,1 - 9a,2 - ^^ 9a,s < 0;

/a,1 - 9a,2 - 9a,0 - 9a,2 - / J 9a,s - 9a,N-1 - 9a,N = -9a,1 - 9a,'2 / J ya,s

s=3 s=0

N-1 N -1

|^2,2| - | A2,11 - |A2,s | = -29a,:1 - 29a,0 - 29a,2 - 9a,3 - 9a,s-1

s=3

3

9a,N -1 =

N-1

ta,s ^ ^ 9a,s < 0; s=0 s=0

9a

N-3

s=4

N-4

^N-2,N-2 | - ^^ ^N-2,s | - |AN-2,N-1| = -29a,1 - 9a,N-1 - ^^ 9a,N-1-s

s=1

s=1

3 N-1

- 29a,0 - 29a,2 - 9a,3 = - / Qa,s - / y 9a,s < 0;

a,s

s=0 s=0

N-2

N-3

|AN-1,N-l| - ^^ |AN-1,^

s=1

s=1

— 2ga,1 — 50,2 — (ga,0 — SO^ — ^^ 1 — S^N-^ — ^a,N1 N

= — ga,1 — ga,2 — ^^ ga,s < 0;

s=0

If i _ 1,i _ 2,i = N - 2, i = N - 1, then

N—1 i-2

|Ai,1| — ^^ |A1,s| = —2ga,1 — 2(ga,0 — 2ga,2) — ^^ 50,1-s+1 — 50,1+1

s=2 N-1

a,2) i+1

s=1 N—i+1

S0,s-i+1 — ga,N—i+1 = — ^^ 50,s — ^^ 50,s < 0

s=i+2

s=0

s=0

Since there is diagonal dominance and all elements of the principal matrix A are positive, then by the Gershgorin theorem, see for example [19, p. 78], the real parts of the matrix eigenvalues are positive. □

For further study of the solvability and stability of the method (3.1), we rewrite it in the form of two systems of dimensions N — 1, typing vector

um =

/ \

\uJV-1/

vm =

( vm \

VVv-1/

f

fm+ 2 \J N-1 '/

f

i™+2 V N-1 '/

where

jm+ 2

A

Z 2

jm+ 2

A

N - 1, m = 0,

M - 1.

Then the system (3.1) can be written as

K A

K A K A

+ Du—Aum+l = um- Du——Aum + Afm,

u oZ,a u 2ha

^.....+ = ^ " D^AVm +

m+1

(3.2)

(3.3)

Theorem3.1. System (3.2)-(3.3) has a unique solution.

Proof. By Lemma 3.1, the real parts of the eigenvalues of the matrix A are positive, then the real parts of the matrix DU^A are also positive, hence the real parts of the eigenvalues of the matrix E + DU^A is greater than one, where E is the identity matrix.

But then the absolute values of all eigenvalues of the matrix E + DU^A are greater than one, which implies that the matrix E + D,,^A is nondegenerate, i.e., the system (3.2) is uniquely resolvable.

The unique solvability of the system (3.3) is shown similarly. □

By virtue of the theorem (3.1), the system (3.2)-(3.3) can be rewritten as

v

m+1 _

m+1_

u"" ^ _ S1Um + ASf 1m,

_ S2Vm + AS2 f 2m,

K A K A

S2 = (E + Dv—A)-\E-Dv—A)i

2h

2 h0

K A

51 = (E + Du—A)-\

KA

52 = (E + DV—A)-\

Theorem3.2 (stability of the method). The eigenvalues of the matrices Si and S2 are less than 1 in absolute value.

Proof. Let A = A(A) = Re A + iIm A be an eigenvalue of the matrix A and k > 0. Then

fE- kA\ _ 1 - k\(A)

\E + kA) ~ l + k\(AY

11 - kA|2 = (Re (1 - kA))2 + (Im (1 - kA))2 = 1 - 2kRe A + (kRe A)2 + (kIm A)2,

|1 + kA|2 = (Re (1 + kA))2 + (Im (1 + kA))2 = 1 + 2kRe A + (kRe A)2 + (kIm A)2,

whence, by virtue of the lemma 3.1, it follows |1 - kA|2 < |1 + kA|2, or |1 - kA| < |1 + kA|. From here, by the definition of the matrices Si and S2,, the conclusion of the theorem follows.□

§ 4. Error analysis

Denote the vector value of the error at the nodes Ej = U (x ,tj) - Vj = (ej ,ej), i = 0,... ,N, j = -M0, ...,M .By definition ej = ej = O for all i = 0,...,N,j = -M0,..., 0.

We say that the method converges with order hp + Aq, if there exists a constant C such that ||Ej || ^ C (hp + Aq) for all i = 0,...,N,j = 1,...,M.

For any m = 0,1,..., M - 1 and i = 1,..., N - 1, the residual of the method (2.9) is the value

^m _ U(Xi,tm+1) — U(Xi,tm)

A

K

K /

Do —-\^2gatSU(xi-s+1,tm) + ga,i+iU(xi,tm) +

s=0

i N-i

+ ga,sV (Xi-s+i,tm+l) + ga,i+iU (xi,tm+l) + ga,sU (xi+s-i,tm) +

s=0 s=0

N-i

+ ga,N-i+iU (xN

-i,tm) + / ga,sU (xi+s-i , tm+i ) + ga,N—i+iU (xN-i)tm+iH Fi ,

s=0

-m+i A

Fi 2 = F(xi,tm + -,U(xi,tm+i),Ut (xi,-)).

Z 2 m+!J

Lemma 4.1. If the exact solution u(x,t), v(x,t) of the problem (1.1)—(1.3) is four times continuously differentiable with respect to x and twice continuously differentiable with respect to t, then for the residual of the method the following is true:

||^m|| ^ Ca(h + A2), i = 0,..., N, m = 0,...,M - 1.

The proof is carried out using the Taylor decomposition of the residual at the nodes (xi, tm+i/2) using (2.6) and the fact that the piecewise linear interpolation is of the second order. Approximation (2.8) of homogeneous Neumann boundary conditions is also used.

Theorem 4.1. Let the conditions for the smoothness of the solution formulated in the Lemma 4.1 be satisfied, then the error of the method (2.9) has order h + A2.

The assertion of the theorem is verified by embedding method (2.9) into the general difference scheme with aftereffect [15] using the Lemma 4.1 and the Theorem 3.2.

§ 5. Numerical experiments

Example 5.1. Consider the system

' du dau

— v(x, t) + u(x, i)(| cos(x)| — \Jv?{x} t) + v2(x, t)) + u(x, t — t) cos(t) — cos(x) cos(t — T),

sin(x — n/2 + an/2)

2 cos(an/2)

— = + w(x, i) + w(x, i)(| cos(x)| — y^U2(x, t) + w2(x, ¿)) + v(x, t — t)

dt d | x |

sin(x — n/2 + an/2)

sin(t) — cos(x) sin(t — T),

2 cos(an/2)

where x G [0, n], t G [0, 4n], a = 1.8, t = 1. The Neumann boundary conditions are set

du(x, t)

du(x,t)

x=0 dx

= 0, 0 < t < T,

x=X

ôx

and initial conditions

u(x, t) = cos(x) cos(t), 0 < x < n, —T < t < 0

v(x, t) = cos(x) sin(t), 0 < x < n, —T < t < 0. The exact solution is u = cos(x) cos(t), v = cos(x) sin(t).

M = N e £

23 1.615026 1.413223

24 1.302408 7.381870 x 10"1

25 1.109274 4.852900 x 10"1

26 1.035347 3.171708 x 10"1

27 1.010894 1.935689 x 10"1

2s 1.003349 1.116905 x 10"1

29 1.001074 6.225966 x 10"2

2io 1.000374 3.405036 x 10"2

Table 1. Error of the numerical solution of system (5.1) based on the proposed scheme when

M = N.

M N e £

24 25 1.012499 4.972836 x 10"1

25 26 1.004089 3.191082 x 10"1

26 27 1.001382 1.938703 x 10"1

27 28 1.000607 1.117390 x 10"1

28 29 1.000310 6.226816 x 10"2

29 2io 1.000166 3.405202 x 10"2

Table 2. Error of the numerical solution of system (5.1) based on the proposed scheme for different values of M and N.

M N e £

25 24 1.440611 6.650061 x 10"1

26 25 1.156315 4.722424 x 10"1

27 26 1.050867 3.150544 x 10"1

28 27 1.015641 1.932433 x 10"1

29 28 1.004720 1.116410 x 10"1

2 io 29 1.001456 6.225190 x 10"2

Table 3. Error of the numerical solution of system (5.1) based on the proposed scheme for different values of M and N.

E x a m p l e 5.2. Consider the system

du dau

dt B\xY ( )

dv dav

dt d\x\a

where x G [0,1], t G [0,1], a = 1.2, t = 0.5. The Neumann boundary conditions are set

du(x, t)

dx

and initial conditions

q du(x,t)

x=0 dx

= 0, 0 < t < T,

x=X

u(x, t) = e-t(x3)(1 - x)3, 0 < x < n, -t < t < 0, v(x,t) = —e-t(x3)(1 - x)3, 0 < x < n, -t < t < 0,

The source term is given by:

-t

( r(4K ,s,„ _u.sn 3r(5)

/ = —e-i(.r3)(l - ,)■ + + (1 - xn - + (1 - xn +

+ + » " *>"> " S«»" + <> " ^e^jflV + > = - »)' + o - *>") - S<*" + o - *>"> +

+il^8+i1 - - ifl<*18+(1 - *)■++ The exact solution is: u = e-t(x3)(1 - x)3, v = -e-t(x3)(1 - x)3.

M = N e £

23 4.788462 x 10"04 4.806325 x 10"04

24 2.665764 x 10"04 2.672903 x 10"04

25 1.390116 x 10"04 1.393004 x 10"04

26 7.288451 x 10"05 7.301299 x 10"05

27 3.856764 x 10"05 3.862954 x 10"05

2s 1.099096 x 10"05 1.100758 x 10"05

29 6.395759 x 10"06 6.409708 x 10"06

2io 2.374209 x 10"06 2.381952 x 10"06

Table 4. Error of the numerical solution of system (5.2) based on the proposed scheme when

M = N.

M N e £

23 24 3.185572 x 10"04 3.194584 x 10"04

24 25 1.457558 x 10"04 1.461059 x 10"04

25 26 7.436044 x 10"05 7.450573 x 10"05

26 27 3.894708 x 10"05 3.901345 x 10"05

27 28 2.063072 x 10"05 2.066324 x 10"05

28 29 1.099096 x 10"05 1.100758 x 10"05

29 2io 6.100022 x 10"06 6.111265 x 10"06

Table 5. Error of the numerical solution of system (5.2) based on the proposed scheme for different values of M and N.

M N e £

24 23 4.773200 x 10"04 4.788415 x 10"04

25 24 2.594554 x 10"04 2.600584 x 10"04

26 25 1.366533 x 10"04 1.369104 x 10"04

27 26 7.221225 x 10"05 7.233228 x 10"05

28 27 3.838352 x 10"05 3.844318 x 10"05

29 28 1.587867 x 10"05 1.605439 x 10"05

2io 29 8.174771 x 10"06 8.194037 x 10"06

Table 6. Error of the numerical solution of system (5.2) based on the proposed scheme for different values of M and N.

Funding. The study of the second author was funded by the Russian Science Foundation, project No. 22-21-00075.

REFERENCES

1. Volterra V. Leçons sur la théorie mathématique de la lutte pour la vie, Paris: Gauthier-Villars, 1931.

2. Kayan S-, Merdan H., Yafia R., Goktepe S. Bifurcation analysis of a modified tumor-immune system interaction model involving time delay, Mathematical Modelling of Natural Phenomena, 2017, vol. 12, no. 5, pp. 120-145. https://doi.org/10.1051/mmnp/201712508

3. Pindza E., Owolabi K. M. Fourier spectral method for higher order space fractional reaction-diffusion equations, Communications in Nonlinear Science and Numerical Simulation, 2016, vol. 40, pp. 112-128. https://doi.org/10.1016/j.cnsns.2016.04.020

4. Owolabi K. M. High-dimensional spatial patterns in fractional reaction-diffusion systems arising in biology, Chaos, Solitons and Fractals, 2020, vol. 134, 109723. https://doi.org/10.1016/j.chaos.2020.109723

5. Tian W., Zhou H., Deng W. A class of second order difference approximations for solving space fractional diffusion equations, Mathematics of Computation, 2015, vol. 84, no. 294, pp. 1703-1727. https://doi.org/10.1090/S0025-5718-2015-02917-2

6. Jin X.-Q., Lin F.-R., Zhao Z. Preconditioned iterative methods for two-dimensional space-fractional diffusion equations, Communications in Computational Physics, 2015, vol. 18, no. 2, pp. 469-488. https://doi.org/10.4208/cicp.120314.230115a

7. Lin X.-L., Ng M. K., Sun H.-W. Stability and convergence analysis of finite difference schemes for time-dependent space-fractional diffusion equation with variable diffusion coefficients, Journal of Scientific Computing, 2018, vol. 75, no. 2, pp. 1102-1127. https://doi.org/10.1007/s10915-017-0581-x

8. Lin X.-L., Ng M. K. A fast solver multidimensional time-space fractional diffusion equation with variable coefficients, Computers and Mathematics with Applications, 2019, vol. 78, no. 5, pp. 1477-1489. https://doi.org/10.1016/j.camwa.2019.04.012

9. Hendy A. S., Macias-Diaz J. E. A conservative scheme with optimal error estimates for a multidimensional space-fractional Gross-Pitaevskii equation, International Journal of Applied Mathematics and Computer Science, 2019, vol. 29, no. 4, pp. 713-723. https://doi.org/10.2478/amcs-2019-0053

10. Yue X., Shu S., Xu X., Bu W., Pan K. Parallel-in-time multigrid for space-time finite element approximations of two-dimensional space-fractional diffusion equations, Computers and Mathematics with Applications, 2019, vol. 78, no. 11, pp. 3471-3484. https://doi.org/10.1016/j.camwa.2019.05.017

11. Du R., Alikhanov A. A., Sun Z.-Z. Temporal second order difference schemes for the multidimensional variable-order time fractional sub-diffusion equations, Computers and Mathematics with Applications, 2020, vol. 79, no. 10, pp. 2952-2972. https://doi.org/10.1016/j.camwa.2020.01.003

12. Zaky M. A., Machado J. T. Multi-dimensional spectral tau methods for distributed-order fractional diffusion equations, Computers and Mathematics with Applications, 2020, vol. 79, no. 2, pp. 476-488. https://doi.org/10.1016/j.camwa.2019.07.008

13. Zaky M. A., Hendy A. S., Alikhanov A. A., Pimenov V. G. Numerical analysis of multi-term time-fractional nonlinear subdiffusion equations with time delay: What could possibly go wrong?, Communication in Nonlinear Science and Numerical Simulation, 2021, vol. 96, 105672. https://doi.org/10.1016/j.cnsns.2020.105672

14. Hendy A. S., Zaky M. A., De Staelen R. H. A general framework for the numerical analysis of highorder finite difference solvers for nonlinear multi-term time-space fractional partial differential equations with time delay, Applied Numerical Mathematics, 2021, vol. 169, pp. 108-121. https://doi.org/10.1016/j.apnum.2021.06.010

15. Pimenov V. G. Raznostnye metody resheniya uravnenii v chastnykh proizvodnykh s nasledstvennost'yu (Difference methods for solving partial differential equations with heredity), Yekaterinburg: Ural State University, 2014.

16. Li C., Zeng F. Numerical methods for fractional calculus, Boca Raton: CRC Press, 2015.

17. Tadjeran C., Meerschaert M. M., Scheffler H.-P. A second-order accurate numerical approximation for the fractional diffusion equation, Journal of Computational Physics, 2006, vol. 213, no. 1, pp. 205-213. https://doi.org/10.1016/jjcp.2005.08.008

18. Kim A. V., Pimenov V. G. i-Gladkii analiz i chislennye metody resheniya funktsional'no-differentsi-al'nykh uravnenii (i-Smooth analysis and numerical methods for solving functional differential equations), Moscow-Izhevsk: Regular and Chaotic Dynamics, 2004.

19. Wilkinson J. H. The algebraic eigenvalue problem. Oxford, 1965.

Received 19.02.2022 Accepted 20.04.2022

Mohammad Ibrahim, Post-Graduate Student, Department of Computational Mathematics and Computer Science, Ural Federal University, pr. Lenina, 51, Yekaterinburg, 620000, Russia. ORCID: https://orcid.org/0000-0001-7991-497X E-mail: mokhammad.ibragim@urfu.ru

Vladimir Germanovich Pimenov, Doctor of Physics and Mathematics, Professor, Department of Computational Mathematics and Computer Science, Ural Federal University, pr. Lenina, 51, Yekaterinburg, 620000, Russia.

ORCID: https://orcid.org/0000-0002-4042-6079 E-mail: v.g.pimenov@urfu.ru

Citation: M. Ibrahim, V. G. Pimenov. Numerical method for system of space-fractional equations of superdiffusion type with delay and Neumann boundary conditions, Izvestiya Instituta Matematiki i Informatiki Udmurtskogo Gosudarstvennogo Universiteta, 2022, vol. 59, pp. 41-54.

М. Ибрагим, В. Г. Пименов

Численный метод для системы дробных по пространству уравнений супердиффузионного типа с запаздыванием и граничными условиями Неймана

Ключевые слова: супердиффузионные уравнения, условия Неймана, функциональное запаздывание, производные Рисса, аппроксимация Грюнвальда-Летникова, метод Кранка-Никольсон, порядок сходимости.

УДК: 519.63

DOI: 10.35634/2226-3594-2022-59-04

Рассматривается система двух дробных по пространству уравнений супердиффузии с функциональным запаздыванием общего вида и краевыми условиями Неймана. Для этой задачи конструируется аналог метода Кранка-Никольсон, основанный на сдвинутых формулах Грюнвальда-Летникова для аппроксимации дробных производных Рисса по пространственной переменной и применении кусочно-линейной интерполяции дискретной предыстории с экстраполяцией продолжением для учета эффекта запаздывания. С помощью теоремы Гершгорина доказана разрешимость разностной схемы и ее устойчивость. Получен порядок сходимости метода. Представлены результаты численных экспериментов.

Финансирование. Исследования второго автора выполнены при поддержке Российского Научного Фонда, проект № 22-21-00075.

СПИСОК ЛИТЕРАТУРЫ

1. Вольтерра В. Математическая теория борьбы за существование. М.: Наука, 1976.

2. Kayan §., Merdan H., Yafia R., Goktepe S. Bifurcation analysis of a modified tumor-immune system interaction model involving time delay // Mathematical Modelling of Natural Phenomena. 2017. Vol. 12. No. 5. P. 120-145. https://doi.org/10.1051/mmnp/201712508

3. Pindza E., Owolabi K. M. Fourier spectral method for higher order space fractional reaction-diffusion equations // Communications in Nonlinear Science and Numerical Simulation. 2016. Vol. 40. P. 112-128. https://doi.org/10.1016/j.cnsns.2016.04.020

4. Owolabi K. M. High-dimensional spatial patterns in fractional reaction-diffusion systems arising in biology // Chaos, Solitons and Fractals. 2020. Vol. 134. 109723. https://doi.org/10.1016/j.chaos.2020.109723

5. Tian W., Zhou H., Deng W. A class of second order difference approximations for solving space fractional diffusion equations, Mathematics of Computation, 2015, vol. 84, no. 294, pp. 1703-1727. https://doi.org/10.1090/S0025-5718-2015-02917-2

6. Jin X.-Q., Lin F.-R., Zhao Z. Preconditioned iterative methods for two-dimensional space-fractional diffusion equations // Communications in Computational Physics. 2015. Vol. 18. No. 2. P. 469-488. https://doi.org/10.4208/cicp.120314.230115a

7. Lin X.-L., Ng M. K., Sun H.-W. Stability and convergence analysis of finite difference schemes for time-dependent space-fractional diffusion equation with variable diffusion coefficients // Journal of Scientific Computing. 2018. Vol. 75. No. 2. P. 1102-1127. https://doi.org/10.1007/s10915-017-0581-x

8. Lin X.-L., Ng M. K. A fast solver multidimensional time-space fractional diffusion equation with variable coefficients // Computers and Mathematics with Applications. 2019. Vol. 78. No. 5. P. 1477-1489. https://doi.org/10.1016/j.camwa.2019.04.012

9. Hendy A. S., Macias-Diaz J. E. A conservative scheme with optimal error estimates for a multidimensional space-fractional Gross-Pitaevskii equation // International Journal of Applied Mathematics and Computer Science. 2019. Vol. 29. No. 4. P. 713-723. https://doi.org/10.2478/amcs-2019-0053

10. Yue X., Shu S., Xu X., Bu W., Pan K. Parallel-in-time multigrid for space-time finite element approximations of two-dimensional space-fractional diffusion equations // Computers and Mathematics with Applications. 2019. Vol. 78. No. 11. P. 3471-3484. https://doi.org/10.1016/j.camwa.2019.05.017

11. Du R., Alikhanov A. A., Sun Z.-Z. Temporal second order difference schemes for the multi-dimensional variable-order time fractional sub-diffUsion equations // Computers and Mathematics with Applications. 2020. Vol. 79. No. 10. P. 2952-2972. https://doi.org/10.1016/j.camwa.2020.01.003

12. Zaky M.A., Machado J. T. Multi-dimensional spectral tau methods for distributed-order fractional diffusion equations // Computers and Mathematics with Applications. 2020. Vol. 79. No. 2. P. 476-488. https://doi.org/10.1016/j.camwa.2019.07.008

13. Zaky M.A., Hendy A. S., Alikhanov A. A., Pimenov V. G. Numerical analysis of multi-term time-fractional nonlinear subdiffusion equations with time delay: What could possibly go wrong? // Communication in Nonlinear Science and Numerical Simulation. 2021. Vol. 96. 105672. https://doi.org/10.1016/j.cnsns.2020.105672

14. Hendy A. S., Zaky M. A., De Staelen R. H. A general framework for the numerical analysis of high-order finite difference solvers for nonlinear multi-term time-space fractional partial differential equations with time delay // Applied Numerical Mathematics. 2021. Vol. 169. P. 108-121. https://doi.org/10.1016/j.apnum.2021.06.010

15. Пименов В. Г. Разностные методы решения уравнений в частных производных с наследственностью. Екатеринбург: Изд-во Уральского ун-та, 2014.

16. Li C., Zeng F. Numerical methods for fractional calculus. Boca Raton: CRC Press, 2015.

17. Tadjeran C., Meerschaert M. M., Scheffler H.-P. A second-order accurate numerical approximation for the fractional diffusion equation // Journal of Computational Physics. 2006. Vol. 213. No. 1. P. 205-213. https://doi.org/10.1016/j.jcp.2005.08.008

18. Ким А. В., Пименов В. Г. i-Гладкий анализ и численные методы решения функционально-дифференциальных уравнений. М.-Ижевск: Регулярная и хаотическая динамика, 2004.

19. Уилкинсон Дж. Х. Алгебраическая проблема собственных значений. М.: Наука, 1978.

Поступила в редакцию 19.02.2022 Принята в печать 20.04.2022

Ибрагим Мохаммад, аспирант, кафедра вычислительной математики и компьютерных наук, Уральский федеральный университет, 620000, Россия, г. Екатеринбург, пр. Ленина, 51. ORCID: https://orcid.org/0000-0001-7991-497X E-mail: mokhammad.ibragim@urfu.ru

Пименов Владимир Германович, д. ф.-м. н., профессор, кафедра вычислительной математики и компьютерных наук, Уральский федеральный университет, 620000, Россия, г. Екатеринбург, пр. Ленина, 51.

ORCID: https://orcid.org/0000-0002-4042-6079 E-mail: v.g.pimenov@urfu.ru

Цитирование: М. Ибрагим, В. Г. Пименов. Численный метод для системы дробных по пространству уравнений супердиффузионного типа с запаздыванием и граничными условиями Неймана // Известия Института математики и информатики Удмуртского государственного университета. 2022. Т. 59. С. 41-54.