A novel time efficient structure-preserving splitting method for the solution of two-dimensional reaction-diffusion systems

In this article, the first part is concerned with the important questions related to the existence and uniqueness of solutions for nonlinear reaction-diffusion systems. Secondly, an efficient positivity-preserving operator splitting nonstandard finite difference scheme (NSFD) is designed for such a class of systems. The presented formulation is unconditionally stable as well as implicit in nature and even time efficient. The proposed NSFD operator splitting technique also preserves all the important properties possessed by continuous systems like positivity, convergence to the fixed points of the system, and boundedness. The proposed algorithm is implicit in nature but more efficient in time than the extensively used Euler method.

Many numerical methods have gained considerable attention, as exact solution of many chemical and biological models involving partial differential equations is an uphill task. Therefore many authors proposed different numerical techniques to find the numerical solution of ordinary, partial, and fractional differential equations arising in the models of chemistry, physics, biology, engineering, and many other fields of different sciences [24][25][26][27][28][29][30][31][32][33][34][35][36]. The proposed NSFD method introduced in this work is concerned with time efficiency as well as preserving of the properties possessed by the continuous reactiondiffusion models. The proposed method is an operator splitting NSFD method. The operator splitting FD technique introduced by Yanenko and his collaborators is a remarkable achievement for the numerical solution to the problems in theoretical mechanics [37]. An explicit forward Euler FD method is compared with the proposed method. The proposed operator splitting NSFD method maintains the same numerical accuracy and reduces the computation time efficiently. Reflecting numerically, more accuracy requires less CPU time; furthermore, easy implementation enhances the applicability of the proposed NSFD scheme.
The major motivation for the current research is to construct the solutions under special physical traits that lead to the computational difficulties. The developed algorithm not only preserves the necessary properties of the solution, but also reduces computational cost. Furthermore, the existence of solutions of the underlying problem is assured.

Existence of solution
The goal of the paper is to solve the following system in a square domain: subject to the initial conditions w i (x, y, 0) = g i (x, y), Ω = a ≤ x, y ≤ b (2) and the homogeneous Neumann boundary conditions where d w i are the diffusion coefficients and are constants, while g i (x, y) are continuous in Ω. The conventional integration leads to the inversion of the partial differential operator ∂ t by simple integration in the following new analytic form: In the operator form we can write where F = F 2 w i + f (x, y, t, w i ) .
Clearly the solutions w i of system (1) are twice continuously differentiable in the domain Ω and continuously differentiable in the time domain [0, T], T > 0. Hence F in (6) is continuous and bounded, i.e., |F | < M, M is a positive real number. For the existence of the solution of system (1)- (2), Hence T i are bounded, and further if we consider a closed and convex subset B r (Θ) (centered at zero element) in the space of continuous functions with radius r, thus which gives restriction on the bounds κ i of the initial values. Now consider the family of images T j i of pre-images w j i , then from (5) we get the following difference: This implies that T j i is equicontinuous, so by the Arzela-Ascoli theorem there exists a subsequence T j l i of T j i which is uniformly convergent, and hence the operator T in (5) turns out to be relatively compact, so the Schauder fixed point theorem is applicable. Summing up the above arguments, the following theorems have been proved. Suppose the right-hand side of F(w i ) of (1) satisfies the following condition: which is a necessary condition to show the unique existence under condition (7). Further, The operator T will be a contraction provided Finally, we establish the following important theorem for the unique existence of the solution of (1)-(2), i.e., we have the following.

Theorem 3
Suppose that the right-hand side of (1) satisfies the Lipschitz condition of the form (8), then problem (1)-(2) is uniquely solvable provided conditions (7) and (9) are satisfied. The notation w i n l,m is used for the finite difference approximations of w i (lh, mh, nτ ). We present two different FD techniques, the forward Euler FD technique and the operator splitting NSFD technique. We propose the positivity preserving NSFD operator splitting technique for the solution of two-dimensional reaction-diffusion systems. In order to verify all the attributes of the proposed NSFD technique, we also use the forward Euler FD technique for the purpose of comparison. The forward Euler FD technique is explicit in nature, time efficient, and easy to implement. But the forward Euler FD technique has the conditional stability. The proposed operator splitting NSFD scheme is implicit in nature and has unconditional stability. The splitting techniques are very efficient numerical techniques for the solution of differential equations. The splitting techniques are widely used numerical techniques for the solution of different nonlinear physical problems [38][39][40][41][42][43].

Forward Euler FD method
In the forward Euler FD method, time derivative is approximated by the forward difference, and space derivatives are approximated by central differences. The forward Euler FD scheme for system (1) is given as follows: where The stability region of the forward Euler FD scheme after linearizing the nonlinear reaction-diffusion is approximately d w i (τ /h 2 ) ≤ 1/4.

NSFD method
The NSFD schemes give numerical approximations to differential equations by constructing discrete models. These schemes preserve the important physical properties of the continuous model and consequently give reliable results. The following are the rules defined by Mickens in [1] for designing the NSFD scheme: Rule 1. The discrete derivative should be of the same order as the order of continuous derivative appearing in the model. After applying these rules, equations (11)-(13) become an NSFD operator splitting scheme.

Stability and accuracy of the proposed NSFD method
The stability and consistency of operator splitting schemes depend on the split solutions [38,44]. The time derivative has O(τ ) accuracy, and the reaction step is exactly solved so it is unconditionally stable. In the similar way, the space derivative has O(h 2 ) accuracy, and the diffusion step is also unconditionally stable.
In the next sections we consider two different chemical reaction models for the application of the positivity-preserving NSFD operator splitting method and the forward Euler method.

Brusselator reaction-diffusion model
The Brusselator model is an auto-catalytic chemical reaction developed by Prigogine [45]. This Brusselator reaction-diffusion system plays a significant role in the application of cooperative processes of chemical kinetics. This system arises in a large number of physical problems. The Brusselator reaction-diffusion system appears in the formation of ozone by atomic oxygen through a triple collision as well as in enzymatic reactions, and in plasma and laser physics. The system of two-dimensional Brusselator model is with the initial condition and the homogeneous Neumann boundary conditions. Here, w 1 = w 1 (x, y, t) and w 2 = w 2 (x, y, t) are the concentrations of the chemical substances. C and D are constant concentrations. The equilibrium point of system (14)- (15) is (w * 1 , w * 2 ) = (D, C/D). Twizell et al. [49] concluded that system (14)-(15) has a stable equilibrium point (w * 1 , w * 2 ) under the condition 1 -C + D 2 ≥ 0 and unstable if 1 -C + D 2 < 0. Now we apply the forward Euler FD scheme (10) on equations (14)-(15), we have Now we apply the NSFD operator splitting technique (11)-(13) on equation (14). At the first time step, the reaction term is solved: After simplification, we have At the second time step, the space derivative w.r.t. x is approximated In a similar way, the NSFD operator splitting scheme for equation (15) at all steps is

Matrix representation
where B 1 and B 3 are real matrices of size (M + 1) × (M + 1) Again let w n+1 where (·) t represents the transposition of a vector. Along with the homogeneous Neumann boundary conditions, equations (22) and (25) are equivalent to the matrix form where B 2 and B 4 are real matrices of size (M + 1) × (M + 1)  27), (30), and (31) can be recovered using the inversion of the matrices B i , i = 1, 2, 3, 4. Since these are M-matrices, by Lemma 1 these, in general, are nonsingular. Therefore these matrices are invertible, and the inverse can be obtained by computational software such as MATLAB.

Positivity of the NSFD operator splitting method
In order to discuss the positivity of the NSFD operator splitting scheme, we use M-matrix theory. A matrix which is real and strictly diagonally dominant is called M-matrix if it has nonpositive off diagonal entries and positive diagonal entries. If a matrix is M-matrix, then it is nonsingular and its inverse matrix has entries of positive numbers [46]. (20) and (23) guarantee the solutions to be positive under the assumptions of nonnegative initial functions, i.e., (20) and (23) as g 1 (x, y) ≥ 0, g 2 (x, y) ≥ 0, and all the terms involved in the right-hand sides of (20) and (23)

Numerical example
Simulations of the example are performed by selecting different values of parameters involved in designing the two-dimensional Brusselator chemical reaction model. Different numerical parametric values are chosen for the better illustration and understanding of the claimed discrepancies about the widely used existing numerical scheme and salient features of the proposed numerical scheme. The following example is considered for system (14)- (15) in two dimensions with initial conditions [47,48]: and homogeneous Neumann boundary conditions.

Simulations of the forward Euler FD scheme
In this section we present numerical simulations of Brusselator reaction-diffusion system  Figure 1 (a) and (b) shows the mesh graph of concentration profile w 1 and w 2 at the grid point (0.9, :) along z-direction vs time t using the forward Euler FD scheme. Figure 1 (c) and (d) depicts the two-dimensional plots of profile w 1 and w 2 at x = y = 1 vs time t. It is clear from Fig. 1 (b) and (d) that w 2 demonstrates the negative solution of concentration, which is useless.

Simulations of the proposed NSFD operator splitting method
This section is devoted to the numerical simulation of the Brusselator model with the same values as given in the above section using the proposed NSFD operator splitting method. Figure 3 (a) and (b) shows the mesh graph of concentration profile w 1 and w 2 at the grid point (0.9, :) along z-direction vs time t using the proposed NSFD operator splitting method. Figure 3 (c) and (d) depicts the plot of concentration profile w 1 and w 2 at x = y = 1 vs time t. These graphs verify that the proposed NSFD operator splitting technique preserves the positive solution of the continuous system and also exhibits the convergence towards the point of equilibrium of continuous system (14)- (15).
In Fig. 4, all parts of the figure clearly describe that the proposed NSFD operator splitting method shows the convergence of the system towards the equilibrium point (w * 1 , w * 2 ) compared with the forward Euler FD scheme.
In addition, we chose C = 3.4 and D = 1.0 and see the behavior of the solution using the proposed NSFD operator splitting method.
It is clearly observed from Figs. 5(a), 5(b) and 6(a), 6(b) that the solutions found using the proposed NSFD operator splitting method are oscillatory and did not converge to the equilibrium point (w * 1 , w * 2 ) with the combination of C = 3.4 and D = 1.0 so that 1 -C + D 2 < 0. This result agrees with the conclusion made by Twizzel et al. in [49]. The comparison of the proposed method and the forward Euler method is presented in Table 1. In this table we show that the exec. (execution) time of NSFD operator splitting technique is better than that of the forward Euler technique.

Glycolysis reaction-diffusion model
A classic and representative system in the biochemical reaction is glycolysis model. It is a basic biochemical reaction appearing in living cells. The fundamental model was presented by Sel'kov [50] and is represented by two coupled first-order differential equations [51,52]. We consider the following two-dimensional glycolysis reaction-diffusion system: with the initial condition Strogatz [52] and Mickens [51] discussed that the equilibrium point has the following stability properties: ξ < 0 : stable.

Example
In this section, we consider a suitable example of two-dimensional glycolysis model of auto-catalytic chemical reactions. This example elucidates all the attributes and essential features about the proposed scheme and system. The example is given as follows: and nonnegative initial conditions with homogeneous Neumann boundary conditions.

Simulations of the forward Euler FD scheme
This section is devoted to the numerical simulations of glycolysis reaction-diffusion system (36)-(37) with initial conditions (51) Fig. 7 (b) and (d) that the concentration profile w 2 gives the negative solution of concentration profile, which is meaningless.

Simulations of the proposed NSFD operator splitting method
This section is devoted to the numerical simulation of glycolysis model with the same values as given in the above section using the proposed NSFD operator splitting method. Figure 9 (a) and (b) shows the mesh graph of the concentration profile w 1 and w 2 at the grid point (0.6, :) along z-direction vs time t using the proposed NSFD operator splitting method. Figure 9 (c) and (d) depicts the plot of concentration profile w 1 and w 2 at x = y = 1 vs time t. These graphs show that the proposed NSFD splitting scheme preserves the positivity property of the solution of the continuous system and also converges to the equilibrium points of glycolysis continuous system (36)- (37).
In Fig. 10, all parts of the figure clearly describe that the proposed NSFD operator splitting method shows the convergence of system towards the equilibrium point (w * 1 , w * 2 ) comparative to the forward Euler FD scheme.
In addition, we chose C = 0.008 and D = 0.6 and see the behavior of the solution using the proposed NSFD operator splitting method. It is clearly observed from Figs. 11(a), 11(b) and 12(a), 12(b) that the solutions found using the proposed NSFD operator splitting method were oscillatory and did not to converge to the equilibrium point (w * 1 , w * 2 ) with the combination of C = 0.008 and D = 0.6 so that ξ > 0. Table 2 demonstrates the efficiency of the proposed technique in various aspects for the glycolysis model.

Susceptible-infected-recovered (SIR) epidemic model with saturated incident rate
In this section, we present an SIR epidemic reaction-diffusion model in two space dimensions, and both techniques are applied on this model. Epidemic models depict the spread of communicable diseases in population through mathematical modeling and give future ideas to control the infection. In these models, an important term is incidence rate which has been considered to ensure that the system of differential equations demonstrates understandable qualitative explanation of the transmission dynamics of diseases [53][54][55][56].
A commonly used incidence rate is the bilinear incidence rate defined as the rate at which susceptible population gains infection. These days various modifications of the standard bilinear incidence rate have been used; for instance, saturation incidence rate. In this section, the SIR epidemic model with saturated incidence rate proposed by Suryato [57] is extended to the following reaction-diffusion system: Here, w 1 , w 2 , and w 3 represent the susceptible, infected, and recovered individuals. The infection force is measured by βw 2 and the inhibition effect is measured by 1 1+σ w 2 . The death rate and the birth rate are demonstrated by the parameter μ, the recovery rate from infection is denoted by γ , and the inhibition parameter is assumed by σ . All the parameters involved in system (53)-(55) are positive. Since w 3 is not present in the first two equations, the above system (53)-(55) can be written as with the initial condition and homogeneous Neumann boundary conditions. System (56)-(57) possesses two steady states, disease-free steady state (DFSS) and endemic steady state (ESS). DFSS of system (56)-(57) is (1, 0) and ESS is (w * 1 , w * 2 ), where w * 1 = μσ +γ +μ μσ +R 0 (γ +μ) and w * 2 = μ(R 0 -1) μσ +R 0 (γ +μ) . The quantity R 0 is the reproductive value that concludes that the disease is erased if R 0 is less than one and disease is presented in population if R 0 is greater than one.

Example
In this section, we consider system (56)-(57) with the following initial conditions: The values of the parameters involved in this experiment are μ = 0.04, γ = 24, d w 1 = d w 2 = 0.001, and σ = 1.

Simulations of the forward Euler FD scheme
In this section, we present the behavior of forward Euler method graphically. First we take β = 20 so that the value R 0 = 0.83195, which is less than one, and the system is stable at DFSS. Figure 13 shows the graphs of forward Euler method at DFSS. In this figure, the graph of infected population describes the negative solution, which is meaningless. For the next figure (Fig. 14), we consider the value β = 40 such that the reproductive value R 0 = 1.66389, which is greater than one. Therefore the system is stable at ESS. But graphs clearly depict that the forward Euler method cannot retain the stability of EES and diverges.

Simulations of NSFD operator splitting scheme
This section is devoted to the presentation of simulations for the proposed method. Initially we consider β = 20 so that the value R 0 = 0.83195, which is less than one, and the system is stable at DFSS. Figure 15 shows the graphs of the proposed splitting method at DFSS. This figure clearly demonstrates that the proposed technique preserves the positive solution of the system under study. For the next figure (Fig. 16), we take the value β = 40 such that the reproductive value R 0 = 1.66389, which is greater than one. Therefore the system is stable at ESS. The NSFD operator splitting method sustains the stability of ESS and converges to the ESS (w * 1 , w * 2 ) = (0.6013, 0.000663). Table 3 also validates the efficacy of the NSFD operator splitting technique in many ways for this model.

Conclusion
In this work, we design a novel numerical method. This method preserves all the essential conditions demonstrated by the continuous reaction-diffusion systems. The proposed method is the operator splitting NSFD scheme that is unconditionally convergent. The proposed NSFD scheme is implicit in nature but efficient in computation time as compared to the forward Euler FD explicit scheme. For the applications, we take three reaction-diffusion models in two space dimensions and observe that the proposed method shows good agreement with the positivity property and convergence to the stable equilibrium points of reaction-diffusion systems. On the other hand, the forward Euler numerical method fails to retain the property of positivity as well as convergence towards equilibrium points of the given reaction-diffusion system. The figures and tables are presented in this work to verify all the attributes of the proposed NSFD operator splitting method.