Abstract
The meshless 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. The numerical experiments revealing the solutions by the developed formulation are stable and more precise.
1. 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 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.
2. 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 , , , , and are given constants, and are functions of the field variables and , and and are assumed to be prescribed sources. In the case of two-component reaction system, and stand for concentrations and , stand for the diffusion coefficients of chemical species.
3. 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 defined in the domain discretized by a set of properly scattered nodes , where is the total number of nodes in the whole domain. It is assumed that only nodes surrounding point have the effect on .
The subdomain that encompasses these surrounding nodes is called the interpolation domain of point . The MK interpolation at point is defined as presented in [6]. Therefore, the formulation of the meshless shape function using MK interpolation is given by where is a vector value of the function in the domain . is a vector of shape functions, expressed as where matrices and are defined as in which is a unit matrix of size and vector is In general, a linear basis in two-dimensional space is a quadratic basis is given as and a cubic basis is For matrix with the size , values of the polynomial basis functions (5) at the given set of nodes are collected as follows: Matrices and vector are defined by the following equations: where is the correlation function between any pair of nodes located at and , representing the covariance of the field value ; that is, Similarly, the covariance can be replaced by . It can be seen from the foregoing formulations that the values of matrices and play important roles in the computation. A simple and frequently used correlation function is a Gaussian function as where and are the correlation parameters used to fit the model and are assumed to be given.
The first-order partial derivatives of the shape function against the coordinates , , can be easily obtained from (3) as where denotes .
4. Local Integral Equations
The local integral formulation of (1) can be written as where is a Heaviside step used as the test function: and are trial functions, and instead of the entire domain we have considered a subdomain located entirely at domain , . The domain is enclosed by , with boundary conditions where is an outward unit normal vector of the boundary and . Condition (16) is often referred to as the Dirichlet boundary condition and (17) is referred to as the Neumann boundary condition. Let and , which substitute and , 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
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 , we solve them iteratively in each time step with replacing and by and , respectively, at zeroth iteration. Equation (25) converted into a set of nonlinear algebraic equations for unknowns and .
5. Numerical Experiments
The analyzed domain is . The error of and , which are presented in the numerical results, is represented by maximum relative error (MRE) and root mean square of relative error (RMSRE) of and , respectively, where and are the exact and computed values of at point , respectively, and is the number of nodes.
Example 1. The system of nonlinear PDEs in the region 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 (1) with and . The shown results have been obtained using , 81, 256, and 441 nodal points, respectively, and at time instant . Figures 1(a) and 1(b) show that the MRE of and increases as a function of the number of nodal points using ; meanwhile, using and 10, the MRE increases gradually as a function of the number of nodal points. Figures 2(a) and 2(b) show that the RMSRE of and decreases as a function of the number of nodal points using , 6, and 10. Moreover, the errors of and using , 10 are less than the RMSRE of and when using . The profile of trial solution of and is similar to the profile of exact solution of and (see Figures 3(a), 3(b), 3(c), and 3(d)). Figures 3(e) and 3(f) reveal corresponding error profile of and . The errors of and satisfy the boundary conditions as well as the Kronecker delta property.
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
Example 2 (application for Brusselator system). The developed formulation from this research can solve a real world application example. Let us consider the nonlinear reaction-diffusion Brusselator system in the two-dimensional region [7], . Consider with , , , initial conditions and Neumann boundary conditions for which the exact solution is unknown. For small values of the diffusion coefficient , if 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 . 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].
(a)
(b)
6. 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 reaction-diffusion problems.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research is partially supported by Dhurakij Pundit University, Thailand.