A Local Integral Equation Formulation Based on Moving Kriging Interpolation for Solving Coupled Nonlinear Reaction-Diffusion Equations

Themeshless local Pretrov-Galerkin method (MLPG) with the test function in view of the Heaviside step function is introduced to solve the system of coupled nonlinear reaction-diffusion equations in two-dimensional spaces subjected to Dirichlet and Neumann boundary conditions on a square domain. Two-field velocities are approximated by moving Kriging (MK) interpolation method for constructing nodal shape function which holds the Kronecker delta property, thereby enhancing the arrangement nodal shape construction accuracy, while the Crank-Nicolson method is chosen for temporal discretization. The nonlinear terms are treated iteratively within each time step. The developed formulation is verified in two numerical examples with investigating the convergence and the accuracy of numerical results.Thenumerical experiments revealing the solutions by the developed formulation are stable and more precise.


Introduction
Reaction-diffusion systems, which are proposed by Alan Turing [1], have an important application in mathematics, physics, chemistry, and biology.Turing or diffusion-driven instability is initiated by arbitrary random deviations of the stationary state and results in stationary spatially periodic variations in the chemical concentration, that is, chemical patterns.Intuitively turning instability can be understood by considering the long-range effects of chemicals, which are not equal due to the difference in the pace of diffusion and thus the instability arises [2].The Prigogine principle of minimum entropy production from [3] is not in general a necessary condition for the steady state and the most favorable state of the system cannot be determined based on the behavior in the vicinity of the steady state, but one must consider the global nonequilibrium dynamics.The reaction-diffusion equations are often solved by numerical methods and usually diffusion is thought to be stabilizing.The idea that diffusion could make a stable and uniform chemical state unstable was innovative.Shirzadi et al. [4] developed meshless local Petrov-Galerkin (MLPG) formulation for numerical solution of the nonlinear reaction-diffusion equations.The spatial variations are approximated by moving least squares and the nonlinear terms are treated iteratively within each time step.Constructing shape functions is one of the most important issues in the MLPG method.There are many methods for constructing shape functions such as the moving least squares (MLS) and the weighted least squares (WLS) method.The most popular method is MLS.Although the MLPG method and many other meshless methods have been gradually applied to different fields, there exists inconvenience because of the difficulty in implementing some essential boundary conditions; the shape function constructed by MLS approximation does not satisfy the Kronecker delta function property.In this research, the MLPG method based on moving Kriging approximation is developed to solve the system of two nonlinear partial differential equations of the parabolic type.The moving Kriging, which was proposed by Gu [5] for constructing shape function, has the Kronecker delta property that is a 2 Advances in Mathematical Physics good property for constructing the shape function.These systems are solved by local integral equation formulation and one-step time discretization method.The nonlinear terms are treated iteratively within each time step.The boundary and domain integrals are calculated using the Simpson and the Gauss-Legendre quadrature rules.Two numerical examples are considered in order to verify the proposed method with testing its accuracy and convergence.

Governing Equation
The numerical simulations of the coupled pair of nonlinear partial differential equations are as follows [4]: given initial and Dirichlet and/or Neumann boundary conditions in the two-dimensional region Ω, where  1 ,  2 ,  1 ,  2 , and  are given constants,  and  are functions of the field variables  and V, and  1 and  2 are assumed to be prescribed sources.In the case of two-component reaction system, (x, ) and V(x, ) stand for concentrations and  1 ,  2 stand for the diffusion coefficients of chemical species.

Moving Kriging Interpolation Method
The Kriging interpolation is a well-known geostatic technique for spatial interpolation in geology and mining [5].The formulation of the construction of meshless shape function by moving Kriging (MK) interpolation is introduced briefly in the following.Similar to the MLS approximation, consider the function (x) defined in the domain Ω discretized by a set of properly scattered nodes x  , ( = 1, 2, . . ., ), where  is the total number of nodes in the whole domain.It is assumed that only  nodes surrounding point x have the effect on (x).
The subdomain Ω x that encompasses these surrounding nodes is called the interpolation domain of point x.The MK interpolation  ℎ (x) at point x is defined as presented in [6].Therefore, the formulation of the meshless shape function using MK interpolation is given by where u = [(x 1 ) (x 2 ) ⋅ ⋅ ⋅ (x  )]  is a vector value of the function in the domain Ω.Φ(x) is a 1 ×  vector of shape functions, expressed as where matrices A and B are defined as in which I is a unit matrix of size  ×  and vector p(x) is In general, a linear basis in two-dimensional space is a quadratic basis is given as and a cubic basis is For matrix P with the size  × , values of the polynomial basis functions (5) at the given set of nodes are collected as follows: Matrices R and vector r(x) are defined by the following equations: where (x  , x  ) is the correlation function between any pair of nodes located at x  and x  , representing the covariance of the field value (x); that is, Similarly, the covariance [(x  ) (x  )] can be replaced by (x, x  ).It can be seen from the foregoing formulations that the values of matrices R and r play important roles in the computation.A simple and frequently used correlation function is a Gaussian function as where   = ‖x  − x  ‖ and  > 0 are the correlation parameters used to fit the model and are assumed to be given.The first-order partial derivatives of the shape function Φ(x) against the coordinates x  ,  = 1, 2, can be easily obtained from (3) as where (⋅) , denotes (⋅)/  .

Local Integral Equations
The local integral formulation of (1) can be written as where   is a Heaviside step used as the test function: and V are trial functions, and instead of the entire domain Ω   we have considered a subdomain Ω   located entirely at domain Ω, x = (, )  ∈ Ω ⊂  2 .The domain Ω is enclosed by Γ = Γ  ∪ Γ  , with boundary conditions where n = ( 1 ,  2 )  is an outward unit normal vector of the boundary and n⋅∇ ≡ /n.Condition ( 16) is often referred to as the Dirichlet boundary condition and ( 17) is referred to as the Neumann boundary condition.Let ũ(x, ) and Ṽ(x, ), which substitute (x, ) and V(x, ), respectively, be the trial solutions: For internal nodes, from local integral equation ( 14), and using the MK (2), we have the following nonlinear equations: The following abbreviations have been used for the integral term: The boundary and domain integrals are calculated using the Simpson and the Gauss-Legendre quadrature rules.We can rewrite (19) as (21) 4.1.Temporal Discretization.Equation ( 21) can be rewritten as Similarly, we have where The finite-difference approximation of the time derivatives of ( 22) and ( 23) in Crank-Nicolson method is given as follows: Because  and  are nonlinear functions of  and V, we solve them iteratively in each time step with replacing  +1 and  +1 by   and   , respectively, at zeroth iteration.Equation (25) converted into a set of nonlinear algebraic equations for unknowns Û+1 and V+1 .

Numerical Experiments
The analyzed domain is and û are the exact and computed values of  at point x  , respectively, and  is the number of nodes.
Example 1.The system of nonlinear PDEs in the region Ω = [0, 1] × [0, 1] is given as follows: where The initial and Dirichlet boundary conditions are chosen in such a way that the exact solution is In this example, the governing equations resemble those of ( Example 2 (application for Brusselator system).The developed formulation from this research can solve a real world application example.Let us consider the nonlinear reactiondiffusion Brusselator system in the two-dimensional region [7], with  = 0.002,  = 1/2,  = 1, initial conditions for which the exact solution is unknown.For small values of the diffusion coefficient , if 1 −  +  2 > 0 then the numerical solution of the Brusselator system converges to an equilibrium point (, /) (see [7]).The experimental results for maximum and minimum values of the exact solution are presented in Table 1.According to Table 1, it is found that the values of the exact solution tend to the steady state values of ( * , V * ) = (, /) = (1, 1/2).Figure 4 shows how the solution changes from initial condition to the steady state as  tends to the infinity.The experimental results are similar to those previously reported [7,8].

Conclusions
The developed formulation has been proposed for coupled nonlinear reaction-diffusion equation by using MK nodal shape function with temporal discretization by the Crank-Nicolson method.Cubic polynomial basis is the best for constructing the nodal shape function.The robust method works well in the sense of accuracy and satisfies the boundary condition.Moreover, the solutions by the developed formulation with temporal discretization by Crank-Nicolson with iterative methods are stable and more precise, so the increased Δ can be chosen for studying.The convergence testing is shown in the Brusselator model for which the exact solution was unknown, showing that the proposed method is competent at simulating coupled nonlinear reactiondiffusion problems.
with  = 0 and  1 =  2 = 1.The shown results have been obtained using  = 25, 81, 256, and 441 nodal points, respectively, and Δ = 0.1 at time instant  = 1.Figures1(a) and 1(b) show that the MRE of  and V increases as a function of the number of nodal points using  = 3; meanwhile, using  = 6 and 10, the MRE increases gradually as a function of the number of nodal points.Figures2(a) and 2(b) show that the RMSRE of  and V decreases as a function of the number of nodal points using  = 3, 6, and 10.Moreover, the errors of  and V using  = 6, 10 are less than the RMSRE of  and V when using  = 3.The profile of trial solution of  and V is similar to the profile of exact solution of  and V (see Figures3(a), 3(b), 3(c), and 3(d)).Figures3(e) and 3(f) reveal corresponding error profile of  and V.The errors of  and V satisfy the boundary conditions as well as the Kronecker delta property.

Table 1 :
The maximum and minimum values of the exact solution by using Crank-Nicolson method:  = 10, Δ = 0.1, and  = 25.