Finite Difference Implicit Schemes to Coupled Two-Dimension Reaction Diffusion System

In this research article, two finite difference implicit numerical schemes are described to approximate the numerical solution of the two-dimension modified reaction diffusion Fisher’s system which exists in coupled form. Finite difference implicit schemes show unconditionally stable and second-order accurate nature of computational algorithm also the validation and comparison of analytical solution, are done through the examples having known analytical solution. It is found that the numerical schemes are in excellent agreement with the analytical solution. We found, second-implicit scheme is much faster than the first with good rate of convergence also we used NVIDA devices to accelerate the computations and efficiency of the algorithm. Numerical results show our proposed schemes with use of HPC (High performance computing) are very efficient and reliable.


Introduction
Reaction diffusion (RD) equations rise up naturally in systems consisting of many interacting factors, such as chemical reactions and are widely used to identify pattern formation phenomena in diverseness of biological and physical systems [1].The primary ingredients of all these models are in the form of mathematical balance equation [1] which is in two dimensions where ( ) , , u u x y t = is a vector of concentration variables, R(u) describes a local ( ) ( ) ( )   , , with a nonlinear source term ( ) [4].A typical solution of the Equation ( 1), a propagating front, separating two non-equilibrium homogeneous states, one of which ( ) is stable and another one ( ) such fronts behavior is often said to be propagation into stable state and also other referred to as waves (or fronts) of transition from an unstable state [3] [4] [5].The interest in these fronts was stimulated in the early 1980s, such as these fronts can be found in various physical, chemical as well as biological systems [5] [6].u R [7] [8].
Researchers have studied these model problems such as the stability of symmetric traveling waves in the Cauchy problem for a more general case than Equation (2); also some researchers explained perturbation method and found an approximate solution by expanding the solution in terms of a power series and in terms of some small parameters [8].
In this paper, we suggest a reaction diffusion system, which agree to several physical phenomena, the most common is the change in space and time of the concentration of one or more chemical substances.Local chemical reactions in which the metamorphosed into each other, and diffusion which causes the substances to spread out over a surface in space.Reaction diffusion systems are naturally applied in chemistry.However, the system can also describe dynamical processes of non-chemical nature.Mathematically, reaction diffusion systems take the form of semilinear parabolic partial differential equations.The general form of our proposed reaction diffusion system in two dimension, is ( ) where β is diffusion coefficient and α is a reactive factor, and ( ) , u x t is concentration and ( ) , v x t is the velocity of the chemical reaction.The aim of this work is to look into the viability of finite difference schemes for the numerical solution of two-dimension coupled reaction diffusion system.The proposed finite difference schemes show good agreement with analytical solution along efficiency in time.Comparison of two finite difference (FD) schemes is also mentioned with CPU efficiency.
The outlook of the paper is in Section 2 analytical solution, Section 3 smoothness and uniqueness, Section 4 numerical methods, Section 5 numerical results and Section 6 discussion.

Analytical Solution
To derive the analytical solution of the given system in (3), (4), we assume the solution of the two dimension coupled reaction diffusion system, in the following form , , e , , e t x y v x y t

Smoothness and Uniqueness of the Reaction Diffusion System
In order to guarantee the smoothness and uniqueness of a positive solution and to obtain upper and lower bounds of the solution, it is necessary to impose some general assumptions on the various physical parameters and the reaction function [9] [10].Throughout this paper, we always assume that the diffusion coefficient β is positive in domain Ω also at 0 t = the initial values 0 0 , u v are non-negative, where the smoothness hypothesis is used only for the existence problem of the corresponding linear system, and the non-negative hypothesis on the data is to obtain non-negative solutions [9] [10].Let us consider the system, , motivated by the nonlinear reaction functions given by Equation ( 7), we make the following basic assumption on functions f and g [9] [10].

Assumption or Hypothesis (H)
f v ∂ ∂ exists and is bounded subsets of domain Ω and there exists a function with ( ) [9]   [10].This definition implies that the function f is monotone non-decreasing in v and is uniform bounded for 0 v ≥ [9] [10].Clearly this condition is satisfied by Journal of Applied Mathematics and Physics both functions f and g, thus this property leads to where above Equation ( 8) represents 1 2 , F F are quasi monotone non-increasing and quasi monotone non-decreasing functions in Ω respectively [9] [10].
According to classification of the reaction functions, 1 2 , F F are typed III functions [9] [10], which leads to the following definition of the solutions.

Definition
A smooth pair of two vector functions ( ) are called upper and lower solutions respectively, if they satisfy the following inequalities , In the above definitions the smoothness of ( ) , u v is in the sense that these functions are continuously differentiable to the order appeared in Equations ( 7) and ( 8) respectively [9] [10].Hypothesis and above definition leads to the following theorem.

Theorem
Let f and g satisfy above hypothesis (H).If there exist upper and lower solutions ( ) u v converges monotonically from above and below, respectively, to a unique solution (u, v) of system (7

t x y u t x y u t x y t x y v t x y v t x y v t x y
The usefulness of the above theorem is that through suitable construction of upper and lower solutions, not only can the existence problem be ensured, but the stability and the asymptotic behavior of the time dependent solution can also be established from the behavior of the upper and lower solutions [9] [10].Thus, the definition of stability of a steady state solution is in the usual sense of Lyapunov function [9] [10].Unlike scalar systems or coupled systems with qua-si-monotone increasing function, the upper and lower solutions for the present system are interconnected and to be determined simultaneously from relations in Equation ( 9).This makes the determination of these functions more delicate, especially in relation to the stability property of non-homogeneous systems [9] [10].Nevertheless, for the global existence problem or the stability problem with homogeneous boundary conditions the construction of those functions is not very difficult [10].

Numerical Methods
We consider the numerical solution of the nonlinear system in ( 5), ( 6) and ( 7

Second Order Implicit Scheme
The Crank Nicolson scheme for the system in ( 3) and ( 4) can be displayed as follows: where [18].The scheme in Equation ( 11) is a nonlinear implicit Journal of Applied Mathematics and Physics scheme with block linear penta diagonal structure [17] [18].The Newton's iterative method is used to solve this linear system, such scheme is of second order accuracy in both directions, space and time respectively [17] [18].The scheme is unconditionally stable using the von Neumann stability analysis [18].

Computationally Efficient Implicit Scheme
In search of a time efficient alternate, we analyzed the naive version of the Crank Nicolson scheme for the two dimensional equation, and find out that that scheme is not time efficient such that to get time efficiency, the common name of Alternating Direction Implicit (ADI) method can be used [19] [20].In ADI scheme, the two steps are as follows: The trick used in constructing of ADI scheme, is to split time step into two sweeps and apply two different stencils in each half time step, therefore to increment time by one time step in grid point , we first compute both of these stencils, such that the resulting linear system is block tridiagonal [20] [21] [22].
The scheme in Equation ( 12) is a nonlinear implicit scheme of second order accuracy in space and time.According to Von Neumann stability analysis, such scheme is also unconditionally stable [22] [23] [24] [25].

Algorithm
The nonlinear system of Equation ( 12), can be written in the form: where ( ) , , , , , , , , , , , , , , n r r r r  are the nonlinear equations obtained from the system in Equation (12).The system of Equations in ( 12) is solved by Newton's iterative method using the following steps: 1) Specify ( ) 0 W as an initial approximation.
− Solve the linear system Jacobian matrix, which is computed analytically and is the correction vector [26] [27].In the iteration method solution at the previous time step is taken as the initial guess.Iteration at each time step is Journal of Applied Mathematics and Physics stopped when with Tol is a very small prescribed value.The linear system obtained from Newton's iterative method, is solved by Gauss elimination method with partial pivoting also convergence done with iterations along less CPU time [28] [29] [30].

Norms
The accuracy and consistency of the schemes is measured in terms of error norms specially 2 L and L ∞ which are defined as: ( ) Two more interesting error are listed below, ( ) Error

RMS
, where Total terms

Numerical Results
Numerical computations have been performed using the uniform grid, for the test problem, the approximated and analytical solutions such as ( ) , , u x y t and ( ) , , U x y t have been given in Table 1, Table 2 at different grids with some fixed parameters such as Nicolson scheme with different parameters while Figures 3-4 indicate ADI scheme with error profile in Figures 5-6.Also in Tables 3-8, we analyzed some more feathers using ADI scheme at different grids along the same parameters as we mentioned above.We discussed some interesting feathers regarding computer system using ADI scheme in Table 9, Table 10, along Self time which is the time spent in a function excluding the time spent in its child functions [29] [30] [31] [32].Self time, includes overhead resulting from the process of profiling, such as Child functions are involved in coding of such problems [29] [30] [31] [32].Thus later scheme works very fine for half spacing as we improve accuracy, which is lead by Richardon extrapolation method, see from Table 10.

Discussion
In this research article, Crank Nicolson scheme has been successfully applied to find the solutions of two-dimension nonlinear reaction diffusion system.The accuracy and stability of the scheme demonstrated by test problem with data tables and figures.According to T Lakoba [32], implicit CN scheme is not efficient in term of computational time, as we see, the linear algebraic system from Jacobean matrix which is Penta block diagonal [32] [33] [34], with two block diagonals adjacent and two are unpaired at the distance of L from main diagonal [33] [34] [35].During simulation, unpaired block diagonals increase   is that the implicit solver only requires a tridiagonal matrix algorithm to be solved, so that the difference between the true Crank Nicolson solution and ADI approximated solution has an order of accuracy of ( )

2
O k and hence can be ignored with a sufficiently small time-step [32]- [37].The ADI method is a predictor corrector scheme where part of the difference operator is implicit in the initial (prediction) step and another part is implicit in the final (correction) step, whereas in ADI scheme two diagonals are paired with main diagonal and no other diagonals as in Crank Nicolson, thus these schemes reduce the computational time and increase efficiency [32]- [37].The following diagrams indicate the Figure 6.Shows results at different grids, for simple error in the concentration of the diffusion reaction system, with step up in grid sizes, make significant change in error but incremental in time, increase error as we see from Figure 6. Figure 7. Shows simulations with capacity of the system along CPU usage and performance.Also processor calls and threads, by Crank Nicolson.Specification of the system is mentioned in computer applications header.performance of the CPU for two different schemes.The derivation of our ADI scheme for a nonlinear PDE system relies on a few key observations.Most importantly, using the solution at time levels previous to 1 n t t + = , the algorithm converts the nonlinear spatial operator into an implicit but linear operator with variable coefficients.The resulting approximately-factored equation is solved in sweeps along each of the Cartesian directions, including, as is common in ADI approaches, an intermediate 1 2 n t + step, so that all of the proposed algorithms are embodied in the two steps formula that every iteration updated the block tridiagonal linear algebraic system [32]- [39].

Computer Applications
Since last two decades, there is a challenging hasting among vendors and software development communities to bring improvement in performance for Journal of Applied Mathematics and Physics The reaction diffusion Equation (2) represents a model equation for the evolution of a neutron population in a nuclear reactor and also arises in chemical engineering applications, such equation allows for the effects of linear diffusion by means of xx yy + u u and nonlinear local multiplication or reaction through ( ) where the first step is to choose integers n and m to define step sizes[12].Partition the interval [a, b] into n equal parts of width h and the interval [c, d] into m equal parts of width k and place a grid on the rectangle R by drawing vertical and horizontal lines through the points with coordinates scheme.Graphically, Figures 1-2 explain Crank

Figure 4 .
Figure 4. Shows results using ADI scheme with 0.1 t = and 53 53 Grid = × with log

Figure 5 .
Figure 5. Shows results using CN scheme at 0.5 t = and 101 101 Grid = × with log

Figure 8 .
Figure 8. Shows simulations with ADI.CPU usage increase from 9% to 13% with efficiency.

Table 1 .
Error comparison by crank nicolson scheme at different grid sizes, at t = 1 and time step = k = 0.0001.

Table 2 .
Rate of convergence comparison by crank nicolson scheme at different grid sizes at t = 1 and time step = k = 0.0001.

Table 3 .
Solution by ADI scheme at different locations at t = 1, time step = k = 0.0001 and grid size = 25 × 25.

Table 4 .
Error comparison by ADI scheme at different grid sizes, at t = 1 and time step = k = 0.0001.

Table 5 .
Error comparison by ADI scheme at different time with time steps = k = 0.0001.

Table 6 .
Error comparison by ADI scheme at different time steps.

Table 7 .
Rate of Convergence Comparison by ADI scheme at different grid sizes.

Table 8 .
Error Comparison by ADI scheme at different space step sizes.

Table 9 .
Interesting feathers of ADI scheme at different grid sizes.

Table 10 .
Interesting feathers of ADI scheme at different grid sizes.