# Structure Preserving Numerical Analysis of Reaction-Diffusion Models

Authors: Nauman Ahmed; Muhammad Aziz-ur Rehman; Waleed Adel; Fahd Jarad; Mubasher Ali; Muhammad Rafiq; Ali Akgül
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

In this paper, we examine two structure preserving numerical finite difference methods for solving the various reaction-diffusion models in one dimension, appearing in chemistry and biology. These are the finite difference methods in splitting environment, namely, operator splitting nonstandard finite difference (OS-NSFD) methods that effectively deal with nonlinearity in the models and computationally efficient. Positivity of both the proposed splitting methods is proved mathematically and verified with the simulations. A comparison is made between proposed OS-NSFD methods and well-known classical operator splitting finite difference (OS-FD) methods, which demonstrates the advantages of proposed methods. Furthermore, we applied proposed NSFD splitting methods on several numerical examples to validate all the attributes of the proposed numerical designs.


---

## Body

## 1. Introduction

A system of differential equations represents the dynamics of real-life phenomenon. These systems can be applied in the field of physics, chemistry, fluid dynamics, engineering, economics, biological sciences etc. Initially, only temporal systems were used to observe the dynamics of different situations. But, for studying the most general behaviors of such systems, we can include the additional quantities such as advection and diffusion in the systems of differential equations. These systems depict the generic situations of the temporal and spatial models [1–3]. Reaction-diffusion equations are considered as one of the most important equations that are used in the modeling of chemical and biological systems. Certain phenomena such as population densities, concentration of chemical substance, and pressure of fluid are described by the state variables of such systems. Therefore, the positivity of state variables is an important requirement for the discrete scheme and should not produce nonphysical oscillations and negative solutions. There is not much work in the literature for positivity preserving and chaos-free numerical methods for chemical reaction models. In this study, we have designed two operator splitting positivity preserving and chaos-free numerical schemes for various reaction diffusion models arising in chemistry and biology.In this paper, three different models of reaction diffusion equations in one-space dimension are considered for the study. The main aim of this work is to find the numerical solutions of these models with reliable numerical techniques. These proposed techniques preserve the structure of continuous systems in one dimension.First of all, among these important models, we consider the Brusselator system that is used in describing the simulations of nonlinear oscillations in chemical reaction-diffusion processes [4–7]. The importance of oscillation in a biochemical system is very significant and was first introduced by Turning et al. [8]. They showed that when some phenomena are considered with diffusion term, a stable spatial pattern is obtained. Also, this system appears in a wide variety of models, such as ozone formation during oxygen atoms collision and enzymatic reactions. The reaction diffusion Brusselator system contains a couple of variables which intervene with reactions, and this process produces chemicals whose concentrations are then controlled. This model has been revealed as the trimolecular model.Due to the importance of such models, many powerful and efficient techniques have been proposed for solving this system. Such models with nonlinear evolution terms are very hard to solve with the aid of analytical methods. Therefore, numerical methods are needed while dealing with such models. For example, Mittal et al. [9–11] investigated the solution of the one-dimensional Brusselator system using the differential quadrature technique. The behaviour of one-dimensional version of the Brusselator system is studied in [12] with the aid of collocation method based on the radial bases functions. A B-spline approach with a modification was introduced by Jiwari and Yuan in [13] for the Brusselator system with one and two dimensions. Also, a three dimensional form of this system was first solved in [14] using the modified B-spline differential quadrature method. Lin et al. [15] studied the inhomogeneous Brusselator model with cross diffusion process by using a finite volume element approximation and performed the stability analysis for the method and addressed the formation of turning patterns. A new exact solution for this system was driven using the exp-function method in [16]. For more details about the methods for solving the Brusselator model, see [17–20] and references therein.The Brusselator reaction-diffusion model in one dimension is given as follows:(1)∂ψ1∂t=εψ1∂2ψ1∂x2+ϑ2−ϑ1+1ψ1+ψ12ψ2,(2)∂ψ2∂t=εψ2∂2ψ2∂x2+ϑ1ψ1−ψ12ψ2,with initial conditions(3)ψ1x,0=αx,0≤x≤L,(4)ψ2x,0=βx,0≤x≤L,and homogeneous Neumann boundary conditions(5)∂ψ10,t∂x=∂ψ1L,t∂x=0,t>0,(6)∂ψ20,t∂x=∂ψ2L,t∂x=0,t>0.As discussed earlier thatψ1x,t and ψ2x,t are the concentrations of the chemical substances or species, ϑ1 and ϑ2 represent the concentrations’ constants that appear in the reaction process, and εψ1andεψ2 are the diffusion constants. The equilibrium point of systems (1) and (2) is ψ1∗,ψ2∗=ϑ2,ϑ1/ϑ2. The point ψ1∗,ψ2∗ is stable if 1−ϑ1+ϑ22≥0 and unstable if 1−ϑ1+ϑ22<0. The solution to the system represented by equations (1) and (2) describes the positivity criteria as ψ1andψ2 are the concentrations of the two species [21]. Therefore, negative values of the solution of this model are meaningless. The numerical technique applied to find the solution of system by the coupled equation in (1) and (2) must preserve the positivity.The finite-difference technique is an important tool to solve the nonlinear model involving differential equations. Because finding the analytical solution of such models is not an easy task. Therefore, various authors used finite-difference numerical schemes to solve several mathematical models involving ordinary differential equations or partial differential equations [22–28].The nonstandard finite difference method (NSFD) is a powerful technique for solving different type of nonlinear continuous models which was first presented by Mickens back in 1993 [29]. Over the years, this method proved that it can treat continuous dynamical models that should preserve positivity property. For example, Ahmed et al. [30] adapt the NSFD scheme for solving different types of such models including the SEIR reaction diffusion model. Detailed discussion regarding NSFD and positivity preserving techniques can be found in [31–37].The paper is organized as follows: in Section2, four numerical techniques are presented for solving the Brusselator model. In Section 3, the stability and accuracy of the presented methods are introduced. Section 4 is devoted for illustrating the positivity of the purposed schemes. The applications of various reaction-diffusion systems along with the comparison of the four presented methods are presented in Section 5. Lastly, Section 6 provides conclusion for the study.


## 2. Numerical Techniques

In this section, we will introduce four techniques for solving systems (1)and(2) subjected to condition (3) and homogeneous Neumann boundary conditions. These numerical methods can be divided into two classical methods named as forward and backward operator splitting methods and two nonclassical methods named as explicit and implicit OS-NSFD methods. Each of these method has its own pros and cons. These methods are explained as follows.


### 2.1. Forward Euler Operator Splitting Method

In the current section, the first method is introduced which is a splitting method based on a finite difference scheme. The splitting techniques are computationally efficient and handle the complexity and nonlinearity of the differential equations. The main idea of these methods is to split the main equation into two-equation system. The nonlinear reaction steps are described as(7)12∂ψ1∂t=ϑ2−ϑ1+1ψ1+ψ12ψ2,12∂ψ2∂t=ϑ1ψ1−ψ12ψ2.Here,ϑ1 and ϑ2 represent the concentrations’s constants that appears in the reaction process. These last mentioned equations are then utilized for the solution of the linear part of the diffusion equation at the first half step for time as(8)12∂ψ1∂t=εψ1∂2ψ1∂x2,12∂ψ2∂t=εψ2∂2ψ2∂x2,that are used for the second equal time step.For the finite difference approximations, divide0,L×0,T into M×N with ω=L/M and κ=T/N.Grid points arexi=iω,i=0,1,2,…,M, tn=nκ,n=0,1,2,…,N, and ψ1in and ψ2in describe the difference approximations of ψ1iω,nκ and ψ2iω,nκ, respectively. Then, the proposed scheme is used to solve these equations in the form:(9)ψ1¯in+1/2=ψ1in+κϑ2−κϑ1+1ψ1in+κψ1in2u2in,ψ2¯in+1/2=ψ2in+κϑ1ψ1in−κψ1in2ψ2in,where ψ¯1in+1/2andψ¯2in+1/2 demonstrate the concentrations at the first equal step of time and ψ1tn+1/2=ψ1tn+1/2ω and ψ2tn+1/2=ψ2tn+1/2ω. At the next step of time, the form is(10)ψ1in+1=ψ¯1in+1/2+λ1ψ¯1i−1n+1/2−2ψ¯1in+1/2+ψ¯1i+1n+1/2,(11)ψ2in+1=ψ¯2in+1/2+λ2ψ¯2i−1n+1/2−2ψ¯2in+1/2+ψ¯2i+1n+1/2.The boundary condition (5) is incorporated in (10) and (11) by using central difference approximation as(12)ψ1i+1n−ψ1i−1n2ω=0,which implies thatψ11n=ψ1−1n,ψ1M+1n=ψ1M−1n,with initial conditions as(13)ψ1i0=αxi,ψ2i0=βxi.


### 2.2. Backward Euler Operator Splitting Method

Now, in this section, we will implement the backward Euler OS-FD technique to solve the Brusselator system. The procedure for the first half step of time is given as(14)ψ1¯in+1/2=ψ1in+κϑ2−κϑ1+1ψ1in+κψ1in2ψ2in,(15)ψ2¯in+1/2=ψ2in+κϑ1ψ1in−κψ1in2ψ2in,which is identical to the last scheme. At the next equal step of time, the final form for the above mentioned scheme is(16)−λ1ψ1i−1n+1+1+2λ1ψ1in+1−λ1ψ1i−1n+1=ψ¯1in+1/2,(17)−λ2ψ2i−1n+1+1+2λ2ψ2in+1−λ2ψ2i−1n+1=ψ¯2in+1/2.The boundary condition (5) is incorporated in (16) and (17) by using central difference approximation as(18)ψ1i+1n−ψ1i−1n2ω=0,which implies thatψ11n=ψ1−1n,ψ1M+1n=ψ1M−1n,with initial conditions as,(19)ψ1i0=αxi,ψ2i0=βxi.


### 2.3. Operator Splitting Nonstandard Explicit Finite Difference Method

Here, we turn our attention to construct a novel method based on a nonstandard finite difference and an operator splitting method named as (OS-NSFD) explicit scheme. For this, we apply the rules defined by Mickens [29] and the designed OS-NSFD technique at the first half time step is described as(20)ψ¯1in+1/2=ψ1in+κϑ2+κψ1in2ψ2in1+κϑ1+1,(21)ψ¯2in+1/2=ψ2in+κϑ1ψ1in1+κψ1in2.For the second equal time step, the strategy for aforementioned technique is(22)ψ1in+1=1−2λ1ψ¯1in+1/2+λ1ψ¯1i−1n+1/2+ψ¯1i+1n+1/2,(23)ψ2in+1=1−2λ2ψ¯2in+1/2+λ2ψ¯2i−1n+1/2+ψ¯2i+1n+1/2.The boundary condition (5) is incorporated in (22) and (23) by using central difference approximation as(24)ψ1i+1n−ψ1i−1n2ω=0,which implies thatψ11n=ψ1−1n,ψ1M+1n=ψ1M−1n,with initial conditions as(25)ψ1i0=αxi,ψ2i0=βxi.


### 2.4. Operator Splitting Nonstandard Implicit Finite Difference Method

This section is devoted for the second proposed novel scheme named as OS-NSFD implicit scheme. The designed OS-NSFD implicit technique at the first half step of time is(26)ψ¯1in+1/2=ψ1in+κϑ2+κψ1in2ψ2in1+κϑ1+1,(27)ψ¯2in+1/2=ψ2in+κϑ1ψ1in1+κψ1in2.For the next step, the final form of current technique is(28)−λ1ψ1i−1n+1+1+2λ1ψ1in+1−λ1ψ1i−1n+1=ψ¯1in+1/2,(29)−λ2ψ2i−1n+1+1+2λ2ψ2in+1−λ2ψ2i−1n+1=ψ¯2in+1/2.The boundary condition (5) is incorporated in (28) and (29) by using central difference approximation as(30)ψ1i+1n−ψ1i−1n2ω=0,which implies thatψ11n=ψ1−1n,ψ1M+1n=ψ1M−1n,with initial conditions as(31)ψ1i0=αxi,ψ2i0=βxi.In all above numerical schemes, the value ofλ1=εψ1κ/ω2 and λ2=εψ2κ/ω2.


## 2.1. Forward Euler Operator Splitting Method

In the current section, the first method is introduced which is a splitting method based on a finite difference scheme. The splitting techniques are computationally efficient and handle the complexity and nonlinearity of the differential equations. The main idea of these methods is to split the main equation into two-equation system. The nonlinear reaction steps are described as(7)12∂ψ1∂t=ϑ2−ϑ1+1ψ1+ψ12ψ2,12∂ψ2∂t=ϑ1ψ1−ψ12ψ2.Here,ϑ1 and ϑ2 represent the concentrations’s constants that appears in the reaction process. These last mentioned equations are then utilized for the solution of the linear part of the diffusion equation at the first half step for time as(8)12∂ψ1∂t=εψ1∂2ψ1∂x2,12∂ψ2∂t=εψ2∂2ψ2∂x2,that are used for the second equal time step.For the finite difference approximations, divide0,L×0,T into M×N with ω=L/M and κ=T/N.Grid points arexi=iω,i=0,1,2,…,M, tn=nκ,n=0,1,2,…,N, and ψ1in and ψ2in describe the difference approximations of ψ1iω,nκ and ψ2iω,nκ, respectively. Then, the proposed scheme is used to solve these equations in the form:(9)ψ1¯in+1/2=ψ1in+κϑ2−κϑ1+1ψ1in+κψ1in2u2in,ψ2¯in+1/2=ψ2in+κϑ1ψ1in−κψ1in2ψ2in,where ψ¯1in+1/2andψ¯2in+1/2 demonstrate the concentrations at the first equal step of time and ψ1tn+1/2=ψ1tn+1/2ω and ψ2tn+1/2=ψ2tn+1/2ω. At the next step of time, the form is(10)ψ1in+1=ψ¯1in+1/2+λ1ψ¯1i−1n+1/2−2ψ¯1in+1/2+ψ¯1i+1n+1/2,(11)ψ2in+1=ψ¯2in+1/2+λ2ψ¯2i−1n+1/2−2ψ¯2in+1/2+ψ¯2i+1n+1/2.The boundary condition (5) is incorporated in (10) and (11) by using central difference approximation as(12)ψ1i+1n−ψ1i−1n2ω=0,which implies thatψ11n=ψ1−1n,ψ1M+1n=ψ1M−1n,with initial conditions as(13)ψ1i0=αxi,ψ2i0=βxi.


## 2.2. Backward Euler Operator Splitting Method

Now, in this section, we will implement the backward Euler OS-FD technique to solve the Brusselator system. The procedure for the first half step of time is given as(14)ψ1¯in+1/2=ψ1in+κϑ2−κϑ1+1ψ1in+κψ1in2ψ2in,(15)ψ2¯in+1/2=ψ2in+κϑ1ψ1in−κψ1in2ψ2in,which is identical to the last scheme. At the next equal step of time, the final form for the above mentioned scheme is(16)−λ1ψ1i−1n+1+1+2λ1ψ1in+1−λ1ψ1i−1n+1=ψ¯1in+1/2,(17)−λ2ψ2i−1n+1+1+2λ2ψ2in+1−λ2ψ2i−1n+1=ψ¯2in+1/2.The boundary condition (5) is incorporated in (16) and (17) by using central difference approximation as(18)ψ1i+1n−ψ1i−1n2ω=0,which implies thatψ11n=ψ1−1n,ψ1M+1n=ψ1M−1n,with initial conditions as,(19)ψ1i0=αxi,ψ2i0=βxi.


## 2.3. Operator Splitting Nonstandard Explicit Finite Difference Method

Here, we turn our attention to construct a novel method based on a nonstandard finite difference and an operator splitting method named as (OS-NSFD) explicit scheme. For this, we apply the rules defined by Mickens [29] and the designed OS-NSFD technique at the first half time step is described as(20)ψ¯1in+1/2=ψ1in+κϑ2+κψ1in2ψ2in1+κϑ1+1,(21)ψ¯2in+1/2=ψ2in+κϑ1ψ1in1+κψ1in2.For the second equal time step, the strategy for aforementioned technique is(22)ψ1in+1=1−2λ1ψ¯1in+1/2+λ1ψ¯1i−1n+1/2+ψ¯1i+1n+1/2,(23)ψ2in+1=1−2λ2ψ¯2in+1/2+λ2ψ¯2i−1n+1/2+ψ¯2i+1n+1/2.The boundary condition (5) is incorporated in (22) and (23) by using central difference approximation as(24)ψ1i+1n−ψ1i−1n2ω=0,which implies thatψ11n=ψ1−1n,ψ1M+1n=ψ1M−1n,with initial conditions as(25)ψ1i0=αxi,ψ2i0=βxi.


## 2.4. Operator Splitting Nonstandard Implicit Finite Difference Method

This section is devoted for the second proposed novel scheme named as OS-NSFD implicit scheme. The designed OS-NSFD implicit technique at the first half step of time is(26)ψ¯1in+1/2=ψ1in+κϑ2+κψ1in2ψ2in1+κϑ1+1,(27)ψ¯2in+1/2=ψ2in+κϑ1ψ1in1+κψ1in2.For the next step, the final form of current technique is(28)−λ1ψ1i−1n+1+1+2λ1ψ1in+1−λ1ψ1i−1n+1=ψ¯1in+1/2,(29)−λ2ψ2i−1n+1+1+2λ2ψ2in+1−λ2ψ2i−1n+1=ψ¯2in+1/2.The boundary condition (5) is incorporated in (28) and (29) by using central difference approximation as(30)ψ1i+1n−ψ1i−1n2ω=0,which implies thatψ11n=ψ1−1n,ψ1M+1n=ψ1M−1n,with initial conditions as(31)ψ1i0=αxi,ψ2i0=βxi.In all above numerical schemes, the value ofλ1=εψ1κ/ω2 and λ2=εψ2κ/ω2.


## 3. Stability and Consistency of Underlying Techniques

In this section, the stability along with the consistency of the solution for all the above operator splitting methods is investigated. In all of these methods, the time derivative has an accuracy ofOκ when solved exactly for the reaction step. The diffusion step, in the same way, has an Oω2 accuracy, and the accumulative accuracy for all the techniques is first order and second order in time and space, respectively. The reaction steps are proven to be unconditionally stable for all the methods illustrated above when solving in an exact way [38,39]. In forward Euler OS-FD and explicit NSFD techniques represented in equations (22) and (23), the stability region is(32)λi≤12,i=1,2,3.In addition, for the Euler OS-FD and implicit NSFD methods, the stability for the diffusion process is found to be stable without any conditions.Next, we will study the positivity of the solution by considering all of the above schemes.


## 4. Positivity of the Proposed Schemes

This section is concerned with the validation of positivity of the designed NSFD explicit and implicit techniques. The following theorem verifies that the proposed explicit technique retains the positivity of the solution.Theorem 1.
The solution of the proposed explicit techniques in the given formulas (20), (21), (26), and (27) at the reaction step with the assumptions of nonnegative initial conditions, i.e.,(33)ψ1i,n≥0,ψ2in≥0⇒ψ1in+1/2≥0,ψ2in+1/2≥0.Theorem 1 verifies that the property of positivity of the solution is preserved by the proposed OS-NSFD schemes at the reaction step.Remark 1.
A positive solution is achieved by using the explicit NSFD method represented in equations (22) and (23) if(34)1−2λi≥0,i=1,2.
From the above expression, it is clear that(35)λi≤12,i=1,2,which is the same stability condition as (32) of the OS-NSFD explicit method. This proves that the OS-NSFD explicit method can preserve positivity within its stability region.We, then utilize the M-matrix theory [40] which helps in proving the positivity property for the proposed OS-NSFD implicit method (28)-(29).Theorem 2.
For anyω>0 and κ>0, systems (28) and (29) are positive, i.e., ψ1n>0 and ψ2n>0 for all n=0,1,2….Proof.
Systems (28) and (29) can be written as(36)Θψ1n+1=ψ1n,(37)Φψ2n+1=ψ2n.
In the above equations,Θ and Φ are the square matrices as follows:(38)Θ=Θ3Θ10⋯⋯⋯⋯0Θ2Θ3Θ2⋱⋮0Θ2Θ3Θ2⋱⋮⋮⋱⋱⋱⋱⋱⋮⋮⋱⋱⋱⋱⋱⋮⋮⋱Θ2Θ3Θ20⋮⋱Θ2Θ3Θ20⋯⋯⋯⋯0Θ1Θ3,(39)Φ=Φ3Φ10⋯⋯⋯⋯0Φ2Φ3Φ2⋱⋮0Φ2Φ3Φ2⋱⋮⋮⋱⋱⋱⋱⋱⋮⋮⋱⋱⋱⋱⋱⋮⋮⋱Φ2Φ3Φ20⋮⋱Φ2Φ3Φ20⋯⋯⋯⋯0Φ1Φ3.
The off-diagonal entries ofΘ are Θ1=−2λ1 and Θ2=−λ1, and diagonal entries are Θ3=1+2λ1. The off-diagonal entries of Φ are Φ1=−2λ1 and Φ2=−λ1, and diagonal entries are Φ3=1+2λ2. Thus, Θ, and Φ are M-matrices. Therefore, expressions (36) and (37) can be written as(40)ψ1n+1=Θ−1ψ1n,(41)ψ2n+1=Φ−1ψ2n.
If we consider thatψ1n>0 and ψ2n>0, then we get that ψ1n+1>0 and ψ2n+1>0, with the aid of the property of M-matrix and Theorem 1. So, the result is then proved by induction. From this theorem, we conclude that the implicit OS-NSFD preserves the property of positivity unconditionally. □


## 5. Application

To support our claims regarding the proposed OS-NSFD techniques, we will implement these schemes to several reaction diffusion models along with the famous Brusselator model.


### 5.1. Brusselator Model

For the application of the proposed and classical methods under study, we chose the following numerical test. The graphical behaviour of the solution for each of the proposed four methods are demonstrated through the following figures for solving systems (1) and (2). As discussed above, systems (1) and (2) have a positive solution and converge toward the equilibrium point if the condition 1−ϑ1+ϑ22≥0 is satisfied and unstable if 1−ϑ1+ϑ22<0.First, Figures1 and 2 demonstrates the simulated results for the first method which is the forward Euler OS-FD technique at various step sizes. It can be noticed from these figures that this method gives nonconsistent data with systems (1) and (2). This proves that this method gives negative values of the concentrations which is meaningless. Also, the solution diverge from the equilibrium point which contradict with the physical behavior of the solution.Figure 1
Mesh graphs and plot graphs ofψ1andψ2 (concentration profile) using forward Euler OS-FD method at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Figure 2
Plot graphs ofψ1andψ2 implementing the forward Euler OS-FD technique at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.001,andεψ1=εψ1=10−4. (a) Plot graph of ψ1. (b) Plot graph of ψ2.
(a)(b)The second method is tested on the same example which is the backward Euler OS-FD method with the same value of parameters as in the first method for the sake of comparison. The behavior of the solution is illustrated in Figures3 and 4. These figures show that this method, like the first method, fails to preserve the positivity of the solution and also gives divergence.Figure 3
Solution graphs and plot graphs ofψ1andψ2 implementing backward Euler OS-FD technique at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Figure 4
Plot graphs ofψ1andψ2 implementing the backward Euler OS-FD method at κ=0.01,ϑ1=1,ϑ2=3.4,λ1=λ2=0.001,andεψ1=εψ1=10−4. (a) Plot graph of ψ1. (b) Plot graph of ψ2.
(a)(b)The graphical behaviour of the state variables show that the forward Euler and backward Euler OS-FD method are not the reliable techniques to solve the nonlinear auto-catalytic chemical reaction model. As they provide us with the negative solutions for the small step sizes.Due to the failure of the two abovementioned methods, the NSFD methods are presented to overcome this issue. Figures5 and 6 show the graphs of the solution using the explicit OS-NSFD method. For the sake of comparison with the other methods, we use the same values of the parameters which are used in the previous two methods and took the value of ϑ1 and ϑ2 such that 1−ϑ1+ϑ22≥0. From these figures, unlike the other methods, we observe that this method preserves the positivity and all the important properties of the glycolysis continuous model and converges toward the equilibrium point.Figure 5
Mesh graphs and plot graphs ofψ1andψ2 implementing NSFD explicit splitting method at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Figure 6
Mesh and plot graphs ofψ1andψ2 implementing NSFD explicit splitting technique at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.001,andεψ1=εψ1=10−4. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)The graphs in Figures7 and 8 demonstrate the concentrations by using the implicit OS-NSFD method. This behavior is also shown to be consistent with the explicit OS-NSFD method. This method like the previous method also preserves positivity and chaos-free properties and converges towards the equilibrium point which is stable under the condition 1−ϑ1+ϑ22≥0>0.Figure 7
Mesh graphs and plot graphs ofψ1andψ2 implementing NSFD implicit splitting technique at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Figure 8
Mesh and plot graphs ofψ1andψ2 implementing NSFD implicit splitting method at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.001,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Now, the simulations by using both designed OS-NSFD techniques are given by considering the values of parametersϑ1 and ϑ2 so that 1−ϑ1+ϑ22>0.As mentioned before, the stability of the Brusselator system is preserved with the condition1−ϑ1+ϑ22>0. Figures 9 and 10 prove this fact that the two new presented OS-NSFD methods are consistent with the continuous model when the stability condition holds and are inconsistent when 1−ϑ1+ϑ22<0. Figures 9 and 10 reflect the inconsistent behavior when the stability criterion does not hold, while the other while the other parameters and step sizes and kept fixed. This fact has been proven that stability criteria is very important in the stability analysis.Figure 9
Solution graphs ofψ1andψ2 implementing NSFD explicit splitting technique at κ=0.1,ϑ1=3.4,ϑ2=1,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Figure 10
Solution graphs ofψ1andψ2 implementing NSFD implicit splitting method at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Now, we considerT=1 and M=100. Figure 11 depicts log-log graphs of the error versus the computational time. We presented the four schemes, and the graphs correspond to the solutions when we estimate the concentration profiles ψ1 and ψ2 of the Brusselator system. The result demonstrates that the OS-NSFD explicit and implicit schemes are approximately as efficient as forward Euler OS-FD and backward Euler OS-FD schemes. Obviously, the former technique stands out in terms of its capability to preserve the structure of the relevant solutions of the Brusselator model under investigation. It is worth pointing out that the approximation to the exact solution was obtained using N=20000.Figure 11
Log-log graphs. The error versus the computational time when solving the system subjected to the initial data. Three schemes are employed: forward Euler OS-FD technique, backward Euler OS-FD technique, OS-NSFD explicit technique, and OS-NSFD implicit technique. The graphs correspond to the solutions when the functionsψ1andψ2 of the Brusselator model are estimated. The parameters employed are L=1, T=1, ϑ1=0.05,ϑ2=0.02, M=100, and εψ1=εψ1=10−4. The exact solution was approximated using N=20000. (a) Error time graph for ψ1. (b) Error time graph for ψ2.
(a)(b)Finally, four numerical methods are used to solve the underlying model, namely, the forward Euler operator splitting method, backward Euler operator splitting method, nonstandard finite difference explicit operator splitting method, and nonstandard finite difference implicit operator splitting method. The extensively used numerical Euler methods (backward and forward) are used to solve the underlying models. The obtained solutions are analyzed and compared with the newly developed methods. The failure of the classical methods motivated us to develop new numerical methods. The new techniques are structure-preserving and reliable numerical methods that give positive and bounded solutions. The computed solutions converge towards the exact steady-state. So, the numerical analysis demonstrates that the NSFD methods are the reliable tool to solve the nonlinear models.


### 5.2. Susceptible-Infected-Recovered Epidemic Model

For the second application, we consider nonlinear reaction-diffusion system of infectious disease dynamics. The system is known as susceptible-infected-recovered (SIR) epidemic [41],(42)∂ψ1∂t=εψ1∂2ψ1∂x2+ϑ1N−ϑ1ψ1−ϑ3ψ1ψ2,(43)∂ψ2∂t=εψ2∂2ψ2∂x2−ϑ1+ϑ2ψ2+ϑ3ψ1ψ2,with initial conditions(44)ψ1x,0=325000x,0≤x≤12,3250001−x,12≤x≤1,ψ2x,0=7500x,0≤x≤12,75001−x,12≤x≤1,and homogeneous Neumann boundary conditions. In systems (42) and (43), ψ1=ψ1x,t and ψ2=ψ2x,t are susceptible and infected population densities. ϑ1 is the parameter which represents natural birth and death rate. ϑ2 indicates the rate of recovery from infected to recovered class. The transmission coefficient from susceptibility to disease is denoted by ϑ3.The epidemic systems (42) and (43) have two stable fixed points, disease-free point (DFP) and endemic point (EP). DFP is the point when disease eradicates from the population. The DFP of systems (42)and(43) ϵ0=N,0. The EP is the point when disease persists in the population. The EP of systems (42)and(43) is ϵ∗=ψ1∗,ψ2∗, where ψ1∗=ϑ1+ϑ2/ϑ3 and ϑ2∗=ϑ1N/ϑ1+ϑ2−ϑ1/ϑ3.The basic reproductive number of the epidemic system (42)and(43) is ℬ0=Nϑ3/ϑ1+ϑ2 when dψi=0,i=1,2. ℬ0 is very important quantity which decides whether disease is eradicated or persisted. If ℬ0<1, then disease will wipe out, and the disease is present if ℬ0>1.In the above experiment, we use the following numerical values of parameters involved in the above system.N=3×10−05, ϑ1=0.04, ϑ2=24, and εψ1=εψ2=0.01. First, we depict the simulations of OS-NSFD explicit and implicit schemes at DFP. For the DFP, we take the value ϑ3=9.2×10−07 so that ℬ0<1.The simulations in the Figures12 and 13 describe that the proposed splitting methods show positive solution as well as preserves the stability of fixed point DFP ϵ0=N,0=3×10−05 under the condition ℬ0<1. This shows that the designed NSFD splitting techniques do not show the contrived chaos.Figure 12
Solution graphs ofψ1andψ2 (susceptible and infected population densities) using the NSFD explicit splitting method at κ=0.1,λ1=λ2=0.05. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Figure 13
Solution graphs ofψ1andψ2 (susceptible and infected population densities) using NSFD implicit splitting method at κ=0.002,λ1=λ2=125. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Also, the graphs for susceptible and infected individuals reflect the positive and bounded solutions. It is notable that those graphs converge towards the exact fixed points for the step sizes, and other parameter values are mentioned above.Next, we present the simulations of OS-NSFD explicit and implicit schemes at EP. For the EP, we use the valueϑ3=1.4×10−03 so that ℬ0>1.As we have taken the values of parameters so thatℬ0>1. This implicates that SIR epidemic model converges toward EP ϵ∗=ψ1∗,ψ2∗. The simulation executed in Figures 14 and 15 illustrate the positive behavior and convergence to the EP ϵ∗=ψ1∗,ψ2∗ of NSFD explicit and implicit splitting schemes. In the light of the above discussion, it can be concluded that OS-NSFD explicit and implicit methods are reliable numerical methods for the solution of reaction diffusion models. Because these methods confine all the important features of reaction diffusion systems like as positivity, boundedness, and stability at equilibrium points.Figure 14
Solution graphs ofψ1andψ2 (susceptible and infected population densities) using the NSFD explicit splitting method at κ=0.1,λ1=λ2=0.05. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Figure 15
Solution graphs ofψ1andψ2 (susceptible and infected population densities) using the NSFD implicit splitting method at κ=0.002,λ1=λ2=125. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)


## 5.1. Brusselator Model

For the application of the proposed and classical methods under study, we chose the following numerical test. The graphical behaviour of the solution for each of the proposed four methods are demonstrated through the following figures for solving systems (1) and (2). As discussed above, systems (1) and (2) have a positive solution and converge toward the equilibrium point if the condition 1−ϑ1+ϑ22≥0 is satisfied and unstable if 1−ϑ1+ϑ22<0.First, Figures1 and 2 demonstrates the simulated results for the first method which is the forward Euler OS-FD technique at various step sizes. It can be noticed from these figures that this method gives nonconsistent data with systems (1) and (2). This proves that this method gives negative values of the concentrations which is meaningless. Also, the solution diverge from the equilibrium point which contradict with the physical behavior of the solution.Figure 1
Mesh graphs and plot graphs ofψ1andψ2 (concentration profile) using forward Euler OS-FD method at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Figure 2
Plot graphs ofψ1andψ2 implementing the forward Euler OS-FD technique at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.001,andεψ1=εψ1=10−4. (a) Plot graph of ψ1. (b) Plot graph of ψ2.
(a)(b)The second method is tested on the same example which is the backward Euler OS-FD method with the same value of parameters as in the first method for the sake of comparison. The behavior of the solution is illustrated in Figures3 and 4. These figures show that this method, like the first method, fails to preserve the positivity of the solution and also gives divergence.Figure 3
Solution graphs and plot graphs ofψ1andψ2 implementing backward Euler OS-FD technique at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Figure 4
Plot graphs ofψ1andψ2 implementing the backward Euler OS-FD method at κ=0.01,ϑ1=1,ϑ2=3.4,λ1=λ2=0.001,andεψ1=εψ1=10−4. (a) Plot graph of ψ1. (b) Plot graph of ψ2.
(a)(b)The graphical behaviour of the state variables show that the forward Euler and backward Euler OS-FD method are not the reliable techniques to solve the nonlinear auto-catalytic chemical reaction model. As they provide us with the negative solutions for the small step sizes.Due to the failure of the two abovementioned methods, the NSFD methods are presented to overcome this issue. Figures5 and 6 show the graphs of the solution using the explicit OS-NSFD method. For the sake of comparison with the other methods, we use the same values of the parameters which are used in the previous two methods and took the value of ϑ1 and ϑ2 such that 1−ϑ1+ϑ22≥0. From these figures, unlike the other methods, we observe that this method preserves the positivity and all the important properties of the glycolysis continuous model and converges toward the equilibrium point.Figure 5
Mesh graphs and plot graphs ofψ1andψ2 implementing NSFD explicit splitting method at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Figure 6
Mesh and plot graphs ofψ1andψ2 implementing NSFD explicit splitting technique at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.001,andεψ1=εψ1=10−4. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)The graphs in Figures7 and 8 demonstrate the concentrations by using the implicit OS-NSFD method. This behavior is also shown to be consistent with the explicit OS-NSFD method. This method like the previous method also preserves positivity and chaos-free properties and converges towards the equilibrium point which is stable under the condition 1−ϑ1+ϑ22≥0>0.Figure 7
Mesh graphs and plot graphs ofψ1andψ2 implementing NSFD implicit splitting technique at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Figure 8
Mesh and plot graphs ofψ1andψ2 implementing NSFD implicit splitting method at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.001,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2. (c) Plot graph of ψ1. (d) Plot graph of ψ2.
(a)(b)(c)(d)Now, the simulations by using both designed OS-NSFD techniques are given by considering the values of parametersϑ1 and ϑ2 so that 1−ϑ1+ϑ22>0.As mentioned before, the stability of the Brusselator system is preserved with the condition1−ϑ1+ϑ22>0. Figures 9 and 10 prove this fact that the two new presented OS-NSFD methods are consistent with the continuous model when the stability condition holds and are inconsistent when 1−ϑ1+ϑ22<0. Figures 9 and 10 reflect the inconsistent behavior when the stability criterion does not hold, while the other while the other parameters and step sizes and kept fixed. This fact has been proven that stability criteria is very important in the stability analysis.Figure 9
Solution graphs ofψ1andψ2 implementing NSFD explicit splitting technique at κ=0.1,ϑ1=3.4,ϑ2=1,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Figure 10
Solution graphs ofψ1andψ2 implementing NSFD implicit splitting method at κ=0.1,ϑ1=1,ϑ2=3.4,λ1=λ2=0.0006,andεψ1=εψ1=10−4. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Now, we considerT=1 and M=100. Figure 11 depicts log-log graphs of the error versus the computational time. We presented the four schemes, and the graphs correspond to the solutions when we estimate the concentration profiles ψ1 and ψ2 of the Brusselator system. The result demonstrates that the OS-NSFD explicit and implicit schemes are approximately as efficient as forward Euler OS-FD and backward Euler OS-FD schemes. Obviously, the former technique stands out in terms of its capability to preserve the structure of the relevant solutions of the Brusselator model under investigation. It is worth pointing out that the approximation to the exact solution was obtained using N=20000.Figure 11
Log-log graphs. The error versus the computational time when solving the system subjected to the initial data. Three schemes are employed: forward Euler OS-FD technique, backward Euler OS-FD technique, OS-NSFD explicit technique, and OS-NSFD implicit technique. The graphs correspond to the solutions when the functionsψ1andψ2 of the Brusselator model are estimated. The parameters employed are L=1, T=1, ϑ1=0.05,ϑ2=0.02, M=100, and εψ1=εψ1=10−4. The exact solution was approximated using N=20000. (a) Error time graph for ψ1. (b) Error time graph for ψ2.
(a)(b)Finally, four numerical methods are used to solve the underlying model, namely, the forward Euler operator splitting method, backward Euler operator splitting method, nonstandard finite difference explicit operator splitting method, and nonstandard finite difference implicit operator splitting method. The extensively used numerical Euler methods (backward and forward) are used to solve the underlying models. The obtained solutions are analyzed and compared with the newly developed methods. The failure of the classical methods motivated us to develop new numerical methods. The new techniques are structure-preserving and reliable numerical methods that give positive and bounded solutions. The computed solutions converge towards the exact steady-state. So, the numerical analysis demonstrates that the NSFD methods are the reliable tool to solve the nonlinear models.


## 5.2. Susceptible-Infected-Recovered Epidemic Model

For the second application, we consider nonlinear reaction-diffusion system of infectious disease dynamics. The system is known as susceptible-infected-recovered (SIR) epidemic [41],(42)∂ψ1∂t=εψ1∂2ψ1∂x2+ϑ1N−ϑ1ψ1−ϑ3ψ1ψ2,(43)∂ψ2∂t=εψ2∂2ψ2∂x2−ϑ1+ϑ2ψ2+ϑ3ψ1ψ2,with initial conditions(44)ψ1x,0=325000x,0≤x≤12,3250001−x,12≤x≤1,ψ2x,0=7500x,0≤x≤12,75001−x,12≤x≤1,and homogeneous Neumann boundary conditions. In systems (42) and (43), ψ1=ψ1x,t and ψ2=ψ2x,t are susceptible and infected population densities. ϑ1 is the parameter which represents natural birth and death rate. ϑ2 indicates the rate of recovery from infected to recovered class. The transmission coefficient from susceptibility to disease is denoted by ϑ3.The epidemic systems (42) and (43) have two stable fixed points, disease-free point (DFP) and endemic point (EP). DFP is the point when disease eradicates from the population. The DFP of systems (42)and(43) ϵ0=N,0. The EP is the point when disease persists in the population. The EP of systems (42)and(43) is ϵ∗=ψ1∗,ψ2∗, where ψ1∗=ϑ1+ϑ2/ϑ3 and ϑ2∗=ϑ1N/ϑ1+ϑ2−ϑ1/ϑ3.The basic reproductive number of the epidemic system (42)and(43) is ℬ0=Nϑ3/ϑ1+ϑ2 when dψi=0,i=1,2. ℬ0 is very important quantity which decides whether disease is eradicated or persisted. If ℬ0<1, then disease will wipe out, and the disease is present if ℬ0>1.In the above experiment, we use the following numerical values of parameters involved in the above system.N=3×10−05, ϑ1=0.04, ϑ2=24, and εψ1=εψ2=0.01. First, we depict the simulations of OS-NSFD explicit and implicit schemes at DFP. For the DFP, we take the value ϑ3=9.2×10−07 so that ℬ0<1.The simulations in the Figures12 and 13 describe that the proposed splitting methods show positive solution as well as preserves the stability of fixed point DFP ϵ0=N,0=3×10−05 under the condition ℬ0<1. This shows that the designed NSFD splitting techniques do not show the contrived chaos.Figure 12
Solution graphs ofψ1andψ2 (susceptible and infected population densities) using the NSFD explicit splitting method at κ=0.1,λ1=λ2=0.05. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Figure 13
Solution graphs ofψ1andψ2 (susceptible and infected population densities) using NSFD implicit splitting method at κ=0.002,λ1=λ2=125. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Also, the graphs for susceptible and infected individuals reflect the positive and bounded solutions. It is notable that those graphs converge towards the exact fixed points for the step sizes, and other parameter values are mentioned above.Next, we present the simulations of OS-NSFD explicit and implicit schemes at EP. For the EP, we use the valueϑ3=1.4×10−03 so that ℬ0>1.As we have taken the values of parameters so thatℬ0>1. This implicates that SIR epidemic model converges toward EP ϵ∗=ψ1∗,ψ2∗. The simulation executed in Figures 14 and 15 illustrate the positive behavior and convergence to the EP ϵ∗=ψ1∗,ψ2∗ of NSFD explicit and implicit splitting schemes. In the light of the above discussion, it can be concluded that OS-NSFD explicit and implicit methods are reliable numerical methods for the solution of reaction diffusion models. Because these methods confine all the important features of reaction diffusion systems like as positivity, boundedness, and stability at equilibrium points.Figure 14
Solution graphs ofψ1andψ2 (susceptible and infected population densities) using the NSFD explicit splitting method at κ=0.1,λ1=λ2=0.05. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)Figure 15
Solution graphs ofψ1andψ2 (susceptible and infected population densities) using the NSFD implicit splitting method at κ=0.002,λ1=λ2=125. (a) Solution graph of ψ1. (b) Solution graph of ψ2.
(a)(b)


## 6. Schnakenberg Model

Schnakenberg model is autocatalytic in nature given in 1979 [42]. The Schnakenberg model is a coupled reaction diffusion system [43] given as(45)∂ψ1∂t=∂2ψ1∂x2+ϑ3ϑ1−ψ1+ψ12ψ2,(46)∂ψ2∂t=εψ2∂2ψ2∂x2+ϑ3ϑ2−ψ12ψ2,with initial conditions(47)ψ1x,0=e−1+0.001x,(48)ψ2x,0=19,and no flux boundary conditions.Here,ψ1 and ψ2 demonstrate the chemical concentrations of the two species. ϑ1 and ϑ2 are chemical kinetic positive constants, and ϑ3 is the positive constant for the model which is dimensionless. The equilibrium point of systems (45)and(46) is ψ1∗,ψ2∗, where ψ1∗=ϑ1+ϑ2 and ψ2∗=ϑ2/ϑ1+ϑ2. The equilibrium point is stable under the condition ϑ2−ϑ1<ϑ1+ϑ23 [43]. If this condition violates, then the point ψ1∗,ψ2∗ is unstable.Now, we present the simulations of the above experiment. First, we take the values of parametersϑ1=2, ϑ2=3, ϑ3=0.5, and εψ2=10−04 such that the condition ϑ2−ϑ1<ϑ1+ϑ23 satisfies. For these numerical values, the equilibrium point is ψ1∗,ψ2∗=5,0.12.The graphs in Figures16 and 17 show the consistent behavior of both OS-NSFD explicit and implicit schemes. Both schemes show that the graphs converge towards the point ψ1∗,ψ2∗=5,0.12. We are not presenting the graphical solution with forward Euler and backward Euler splitting schemes, but it is confirmed that the schemes will present the inconsistent behavior for different values of parameters as shown earlier for the Brusselator model.Figure 16
Mesh graphs ofψ1andψ2 (chemical concentration) using NSFD explicit splitting method at κ=0.1,λ1=0.1,λ2=0.00001. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2.
(a)(b)Figure 17
Mesh graphs ofψ1andψ2 (chemical concentration) using NSFD implicit splitting method at κ=0.01,λ1=1000,λ2=0.1. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2.
(a)(b)Now, we choose the values of parametersϑ1=0.1, ϑ2=0.4, ϑ3=0.5, and εψ2=10−04 such that the condition ϑ2−ϑ1<ϑ1+ϑ23 violates.As we discussed earlier that systems (45) and (46) show the unstable behavior if the condition ϑ2−ϑ1<ϑ1+ϑ23 does not satisfy. The graphs in Figures 18 and 19 clarify that OS-NSFD explicit and implicit schemes also reveal the unstable behavior which is possessed by the continuous system.Figure 18
Mesh graphs ofψ1andψ2 (chemical concentration) using the NSFD explicit splitting method at κ=0.1,λ1=0.1,λ2=0.00001. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2.
(a)(b)Figure 19
Mesh graphs ofψ1andψ2 (chemical concentration) using the NSFD implicit splitting method at κ=0.01,λ1=1000,λ2=0.1. (a) Mesh graph of ψ1. (b) Mesh graph of ψ2.
(a)(b)It is important to note that all the parameters and step sizes are kept same during the simulations. It is the stability condition that also plays an important role in describing the stability of the system.


## 7. Conclusion

In this article, four numerical schemes based on finite difference approximation are presented for solving different reaction-diffusion models in one-space dimension. It has been observed that two classical methods fail to provide accurate solution. Also, these schemes do not hold the positivity condition of the unknown variables in the continuous system. To overcome this issue, two new positivity preserving techniques have been proposed, based on OS-NSFD schemes. Our proposed methods not only provide the positive solution but also retain the essential physical attributes of the state variables. The designed schemes are applied on the Brusselator model, Schnakenberg model, and SIR epidemic model. The simulations are carried out to obtain the graphical solutions. The numerical results ascertained that the newly designed NSFD splitting schemes have some prominent features of the solution such as positivity of the solutions and stability at equilibrium points of the continuous system. One of the main significance is that the designed schemes do not generate the contrived chaos.The OS-NSFD methods grant the positive solution irrespective of the step sizes. So, these schemes are unconditionally positivity preserving. Also, these schemes are time efficient and consistent.The successful implementation of NSFD splitting schemes on three different problems show that our proposed schemes are futuristic and these techniques can be applied on various physical reaction diffusion problem in the fields of physical science, engineering, fluid mechanics, economics, and many more.



---

*Source: 5128343-2022-03-20.xml*