Numerical Method for Solving the Nonlinear Superdiffusion Equation with Functional Delay

: For a space-fractional diffusion equation with a nonlinear superdiffusion coefﬁcient and with the presence of a delay effect, the grid numerical method is constructed. Interpolation and extrapolation procedures are used to account for the functional delay. At each time step, the algorithm reduces to solving a linear system with a main matrix that has diagonal dominance. The convergence of the method in the maximum norm is proved. The results of numerical experiments with constant and variable delays are presented


Introduction
Partial differential equations of both integer and fractional orders with various complicating effects are widely used in many mathematical models (gas dynamics, population dynamics, and others).These effects may include nonlinearities in the differentiation operator [1,2], delay effects [3,4], and the presence of space-fractional orders (superdiffusion) [5,6].
Linearization techniques initially introduced by Bellman [7] also provide iterative methods to overcome the nonlinear difficulty in differential equations.Later OHAM methods appeared and they have been extensively applied to several types of nonlinear differential equations [8,9].
Due to the complexity of the effects under study, the development of numerical algorithms for solving the problems posed comes to the fore.Analytical solutions can be obtained extremely rarely in such problems [10,11].For equations with fractional derivatives with respect to state, numerical methods are now being actively developed.We note the works [12,13], the results of which are used in this article.More accurate numerical methods for solving linear superdiffusion equations, including those with two or more spatial variables, were developed in [14][15][16][17][18][19][20].A numerical method for a space-fractional equation with a constant delay was developed in [21].
So, for partial differential equations with a functional delay effect, numerical methods were studied, in particular, in the articles [22][23][24][25].Algorithms for the numerical solution of linear space-fractional equations with a functional delay effect were studied in [26].
In this paper, we consider a quasilinear superdiffusion equation with a delay effect.In view of the special form of nonlinearity (quasilinearity), it is possible to construct an efficient algorithm for solving the considered equations.The idea of this algorithm was borrowed from [27]; for the diffusion equation with the delay effect, the idea was implemented in [28].In contrast to [28], where the algorithm reduces to solving a linear system with a tridiagonal matrix at each time step, for the superdiffusion equation with a two-sided Riemann-Liouville derivative, a linear system arises at each time step that does not have the tridiagonal property.However, this system has a diagonal predominance with positive diagonal elements, which makes it possible to solve it quite effectively.In addition, the properties of the system allow us to prove the stability of the method and, as a consequence, to obtain the convergence of the method in the maximum norm.
The structure of the work is as follows.In the second section, the problem is formulated and the main assumptions are made.Section 3 constructs a difference method that takes into account the effect of functional delay, the fractional nature of the derivative with respect to space, and the nonlinearity of the superdiffusion coefficient.In the next section, the local error of the method is studied.In the fifth section, the global error of the method is studied and the convergence theorem is proved.The results of numerical experiments on test examples are presented.In conclusion, the results are summarized and prospects for the development of the method are outlined.

Problem Statement and Basic Assumptions
Consider a nonlinear one-dimensional superdiffusion equation with a functional delay: is the history of the desired function up to the time t, τ is the delay value.K(u(x, t))-nonlinear superdiffusion coefficient, −1 q 1.The left-hand and right-hand fractional derivatives of order α, 1 < α < 2, are defined in the Riemann-Liouville sense Initial and boundary conditions are set Let us assume that there is a unique solution to the problem (1)-(3), while deriving error estimates, we assume that it is sufficiently smooth.
The set of functions q(s) that are piecewise continuous on [−τ, 0], with a finite number of discontinuity points of the first kind, and right-continuous at the discontinuity points is denoted by We will assume that the functional f x, t, u t (x, •) is Lipschitz with constant L f in the last argument, i.e., there is a constant L f , that for all Let for all x ∈ [0, X], t ∈ [t 0 − τ, t 0 ], the exact solution (1)-( 3) |u(x, t)| U is satisfied.We will assume that, in the domain |u| 2U, the following condition is satisfied: We will also assume that the function K(u) is Lipschitz in this domain, meaning that there exists a constant L K such that for any u and v from this domain, (5)

Implicit-Explicit Difference Method
Let us divide the segments [0, X], [t 0 , θ] into parts with steps h = X/N and ∆ = (θ − t 0 )/M, respectively, and introduce the points x i = ih, i = 0, N, t k = t 0 + k∆, k = 0, M. Without loss of generality, we assume that the value τ/∆ = m is an integer.
If τ/∆ is a non-integer, which may happen when (θ − t 0 )/τ is a non-integer, then the step can be introduced as follows.Let us define the step ∆ = τ/m, where m is an integer.Then the number of steps M can be determined through the integer part of the relation: The approximation of the exact solution u(x i , t k ) at the nodes of the grid (x i , t k ) will be denoted by u i k .For each fixed i = 0, N, we introduce a discrete history up to the moment Definition 1.The operator of interpolation-extrapolation of discrete prehistory {u i l } k is the mapping I : In what follows, we will use the piecewise constant interpolation with extrapolation by continuation This method of interpolation with extrapolation has the first order of ∆, i.e., if the exact solution u(x, t) is continuously differentiable with respect to t on [t 0 − τ, θ], then there are constants C 1 = 1 and C 2 , so that for all i, k and t ∈ [t k − τ, t k+1 ], the following inequality holds Note also that the operator of piecewise constant interpolation with extrapolation by continuation is Lipschitz with the Lipschitz constant L I = 1 in the following sense: if w i k (t) and v i k (t) are the results of piecewise constant interpolation with extrapolation by continuation of two discrete prehistories, respectively, {w i j } k and {v i j } k , then for all t ∈ [t k − τ, t k+1 ] the following inequality holds To approximate the left-hand fractional derivative on the k + 1-th time layer, we will use the right-shifted Grunwald-Letnikov formula [5] where the normalized Grunwald weights are defined by the relations , j = 0, 1, . . . .
In particular, g α,0 = 1, g α,1 = −α, Similarly, to approximate the right-hand fractional derivative on the k + 1-th time layer, we will use the left-shifted Grunwald-Letnikov formula [12] We introduce the difference operator Lemma 1.If the condition ( 5) is satisfied, the left-hand and right-hand derivatives of order α + 1 of the exact solution, as well as their Fourier transforms, are continuous.Thus, Proof.By virtue of the assumptions about the corresponding smoothness of the solution u(x, t) of the problem ( 1)-( 3) and assumptions about the function K(u), we obtain and, as follows from [5] The same assumptions also imply that are bounded, whence follows the conclusion of the lemma.
For k = 0, M − 1, consider the implicit-explicit difference scheme with initial and boundary conditions Without interpolation and extrapolation procedures, difference methods would arise in an infinite-dimensional space.The procedure of interpolation and extrapolation with given properties allows us to make the numerical method implicit only in a finite-dimensional space.The use of interpolation and extrapolation makes it possible to explicitly calculate the functional f (x i , t k+1 , v i k (•)); therefore, we call the method implicit-explicit.
The idea of taking the value of the non-linear superdiffusion coefficient from the previous time layer makes the implicit-explicit method linear.
The values of the function K(u) in the system ( 8) are calculated at the points of the time layer t k .Thus, the scheme ( 8)-( 10) is a linear system with respect to the values u i k+1 on the time layer t k+1 .
Let us rewrite the system (8) as Let us write out the matrix A of the coefficients of the unknowns of the system (11), elements of matrix A of dimension N − 1 × N − 1 have the form where Lemma 2. The coefficients of the matrix A of the system (11) have strict diagonal dominance with positive diagonal elements; hence, the system is solvable and has a unique solution.
Proof.Note the properties of the coefficients g Also, we performed the following calculation: Then the diagonal elements are positive; moreover, In addition, all off-diagonal elements of the matrix are negative.Let us show strict diagonal dominance.Let, for example, i = 1, then This follows from the fact that The system (11) is solved using LU factorization, which is computed by Gaussian elimination with partial pivoting.

Residual of the Difference Method
Definition 2. Residual without interpolation of the method ( 8) is called Lemma 3. Let the conditions of Lemma 1 be satisfied, and, moreover, the exact solution be twice continuously differentiable with respect to t.Then the residual without interpolation of the method (8) has the order h + ∆, i.e., there exists a constant C 7 , that Proof.According to the numerical differentiation formula, we have Then, from the statement of Lemma 1, taking into account the fact that u(x i , t k+1 ) is the exact solution of the equation ( 1), we obtain the assertion of the lemma.

Definition 3.
Residual with piecewise constant interpolation and extrapolation by continuation of the method (8) is called where v(x i , t) for t ∈ [t k − τ, t k+1 ] is the result of piecewise constant interpolation and extrapolation by continuation (6) of the discrete prehistory of the exact solution at the nodes {u(x i , t l )} k .
Lemma 4.Under the condition of the previous lemma the residual with piecewise constant interpolation and extrapolation by continuation of the method (8) has the order h + ∆, i.e., there exists a constant C 9 such that Proof.Residual with piecewise constant interpolation and extrapolation by continuation of the method ( 8) is related to residual without interpolation of the method (8) by Using the Lipschitz conditions (4) and the fact that piecewise constant interpolation with extrapolation by continuation is of the first order of ∆ (7), we obtain From here, from (12) and the assertion of Lemma 3, the required assertion follows.

Error Analysis
Let us determine the error of the method ( 8) We say that the method converges with order h p + ∆ q , if there exists a constant C independent of h and ∆, that |ε i j | C(h p + ∆ q ) for all i = 0, 1, . . ., N and j = 0, 1, . . ., M. Let us define for each time layer with the number j = 0, 1, . . ., M the layer-by-layer error by the vector ) with the norm ε j = max In addition, we determine the accumulated prehistory of the layer-by-layer error by the time t k , k = 0, 1, . . ., M : {ε j } k = {ε j , 0 j k} with norm {ε j } k = max 0 j k Proof.Let us write the method (8) as Let us rewrite the definition of the residual with interpolation in the form Then the equation for the error has the form The modulus of the left side of this relation for the index i 0 is rewritten in the form Due to the property of the coefficient g α,1 = −α, properties of the quantities Let us estimate the right-hand side of the inequality (15) from below using the properties of the coefficients g α,j > 0, j = 0, 2, 3, . . ., ) and the definition of the number i 0 From this inequality and the inequality (15), we obtain the estimate Let us estimate the modulus of the right side of the relation ( 14) for the index i = i 0 .From the definition of the operator Λ, the fact that the function K is Lipschitz, and the boundedness of δ α,x [u(x i , t k+1 )] and δ α,−x [u(x i , t k+1 )], it follows that Using also Lemma 4, we obtain an estimate for the right-hand side of the relation ( 14) for the index i = i 0 : 14), ( 16) and ( 17), the assertion of the lemma follows.
The proved statement means the stability of the difference scheme.In the following statement, the accumulated prehistory of the layer-by-layer error by the time t k+1 is estimated in terms of the accumulated prehistory of the layer-by-layer error by the time t k .Lemma 6.Under the conditions of the previous lemma, we have the estimate 13) we obtain the relation Due to the Lipschitz property of the function f with respect to the last argument and the properties of piecewise constant interpolation with extrapolation by continuation, we obtain From ( 18) and ( 19), the assertion of the lemma follows.
Theorem 1. Letting the exact solution u(x, t) of the problem ( 1)-( 3) satisfy the conditions of Lemma 1 and Lemma 3, then the method ( 8) converges with the order h + ∆.
Proof.From Lemma 6, we have where Using the geometric progression formula, we obtain for all time layers with the number n M, Let us substitute the expressions for A and B into this estimate, and also use the relation ∆M = θ − t 0 : From this, we obtain the estimate uniform over all n = 1, 2 . . ., M. This estimate means the convergence of the method with order h + ∆.

Numerical Experiments
Example 1.Let us consider the following test equation with constant concentrated delay with respect to the variable t: Coefficient β is taken equal to 1. Initial and boundary conditions are set as Let us denote the maximum error in nodes as The method (8) was tested when the spatial step h was changed from the value 1   10   to the value 1 80 with a fixed time step ∆ = 1 4000 .The convergence order of the space is characterized in the experiment by ). Table 1 contains the values E(∆, h), order h and CPU time (sec.)for different parameters α, q and step h.
The simulation of the test equations was performed on the computer Intel Core i5-2467M, 4 cores, CPU 1.6 GHz, 8 Gb RAM.
Table 1.Dependence of the values E(∆, h), order h and CPU time from the spatial step and parameters.To study the dependence of the error on time, the time step ∆ varied from 1 10 to 1 80 with a fixed spatial step h = 1 1000 .The convergence order of the time is characterized in the experiment by

h E(∆, h) order h CPU Time E(∆, h) order h CPU Time
). Table 2 contains the values E(∆, h), order ∆ and CPU time (sec.)for different parameters α, q and step ∆.
Coefficient β is taken equal to 1. Initial and boundary conditions are set as The exact solution is the function u Table 3 for fixed time step ∆ = 1 1000 contains the values E(∆, h), order h and CPU time (sec.)for different parameters α, q and step h.

Conclusions and Directions for Further Research
In this paper, for the first time, a numerical method for solving the superdiffusion differential equation with a nonlinear superdiffusion coefficient is presented.This equation contains several effects that complicate the solution: the fractional nature of the space derivative, the presence of a functional delay, and, most importantly, the nonlinearity of the diffusion coefficient (in this case, superdiffusion).
The difference numerical algorithm is based on three methods.Accounting for lefthanded and right-handed fractional derivatives is carried out using the shifted Grunwald-Letnikov formulas [5].The delay effect is taken into account using interpolation, and implicitness becomes finite-dimensional after additional extrapolation.The nonlinearity of the superdiffusion coefficient is overcome by using the value of this coefficient in the previous time layer, this technique is described in [27].
As a result, the algorithm is reduced to solving a system of linear equations of a special form at each time layer.Following [12], the main matrix of this system is written out.It is shown that it has the diagonal dominance, which implies that the system is uniquely solvable.
The local error (residual) of the algorithm is written out without taking into account interpolation and taking into account the piecewise constant interpolation with extrapolation by continuation.It is shown that the residual values are of the first order of smallness with respect to the partitioning steps in time and space.The grid value of the method error, the layer-by-layer error vector, and the vector of the accumulated error history are introduced.The main result of the paper is the proof of Lemma 5, which plays the role of a statement about the stability of the algorithm in the maximum norm.From this lemma, we derive an estimate of the norm of the accumulated error history at the next time layer in terms of the norm of the accumulated error history at the previous time layer.This also implies the theorem on convergence of the algorithm with the first order of smallness with respect to the steps of partitioning in time and space.It should be noted that the proof of convergence does not rely on the embedding of the algorithm in a general difference scheme with heredity [23,25,26]; however, the ideas of the proof are similar.
The last section of the work presents the results of numerical experiments on test examples with exact solutions.Note that the selection of such examples is also a difficult task, since all three effects must be taken into account: the fractional nature of the derivatives, the nonlinearity of the superdiffusion coefficient, and the presence of delays of various types.The tests confirmed the theoretical conclusions about the convergence of the algorithms.
However, a small order of convergence requires a large number of operations to achieve high accuracy, which, in turn, leads to an increase in the computational error.Therefore, the main issue for further research in solving this problem is the development of algorithms of a higher order of convergence.
The research methodology proposed in this paper can also be applied to other types of equations, primarily to equations with two and three spatial variables.For linear superdiffusion equations, different authors have obtained efficient numerical algorithms, such as ADI method.We hope that this methodology can be modified for non-linear multidimensional space-fractional equations with functional delay.

Table 2 .
Dependence of the values E(∆, h), order ∆ and CPU time from the time step and parameters.

Table 3 .
Dependence of the values E(∆, h), order h and CPU time from the spatial step and parameters.

Table 4
for fixed spatial step h = 1 2000 contains the values E(∆, h), order ∆ and CPU time (sec.)for different parameters α, q and step ∆.

Table 4 .
Dependence of the values E(∆, h), order ∆ and CPU time from the time step and parameters.