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

Authors: Kanittha Yimnak; Anirut Luadsong
Journal: Advances in Mathematical Physics
 (2014)

---

## 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.


---

## Body

## 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]:
(1)∂u∂t=D1∇2u+α1u+A(u,v)+f1(x,t),∂v∂t=D2∇2v+α2v+βu+B(u,v)+f2(x,t)
given initial and Dirichlet and/or Neumann boundary conditions in the two-dimensional region Ω, where D1, D2, α1, α2, and β are given constants, A and B are functions of the field variables u and v, and f1 and f2 are assumed to be prescribed sources. In the case of two-component reaction system, u(x,t) and v(x,t) stand for concentrations and D1, D2 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 functionu(x) defined in the domain Ω discretized by a set of properly scattered nodes xi,(i=1,2,…,n), where n is the total number of nodes in the whole domain. It is assumed that only N nodes surrounding point x have the effect on u(x).The subdomainΩx that encompasses these surrounding nodes is called the interpolation domain of point x. The MK interpolation uh(x) at point x is defined as presented in [6]. Therefore, the formulation of the meshless shape function using MK interpolation is given by
(2)uh(x)=∑I=1NϕI(x)uI=Φ(x)u,x∈Ωx,
where u=[u(x1)u(x2)⋯u(xN)]T is a vector value of the function in the domain Ω. Φ(x) is a 1×N vector of shape functions, expressed as
(3)Φ(x)=pT(x)A+rT(x)B,
where matrices A and B are defined as
(4)A=(PTR-1P)-1PTR-1,B=R-1(I-PA),
in which  I is a unit matrix of size N×N and vector p(x) is
(5)pT(x)=[p1(x1)⋯pm(xN)].
In general, a linear basis in two-dimensional space is
(6)pT(x)=(1,x,y),m=3,
a quadratic basis is given as
(7)pT(x)=(1,x,y,x2,xy,y2),m=6,
and a cubic basis is
(8)pT(x)=(1,x,y,x2,xy,y2,x3,x2y,xy2,y3),m=10.
For matrix P with the size N×m, values of the polynomial basis functions (5) at the given set of nodes are collected as follows:
(9)P=[p1(x1)⋯pm(x1)⋮⋮⋮p1(xN)⋯pm(xN)].
Matrices R and vector r(x) are defined by the following equations:
(10)R=[γ(x1,x1)⋯γ(x1,xN)⋮⋮⋮γ(xN,x1)⋯γ(xN,xN)],rT(x)=[γ(x,x1)⋯γ(x,xN)],
where γ(xi,xj) is the correlation function between any pair of nodes located at xi and xj, representing the covariance of the field value u(x); that is,
(11)γ(xi,xj)=E[u(xi)u(xj)].
Similarly, the covariance E[u(xi)u(xj)] can be replaced by γ(x,xj). 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
(12)γ(xi,xj)=e-θrij2,
where rij=∥xi-xj∥ 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 xi, i=1,2, can be easily obtained from (3) as
(13)Φ,i(x)=p,iT(x)A+r,iT(x)B,
where (·),i denotes ∂(·)/∂xi.


## 4. Local Integral Equations

The local integral formulation of (1) can be written as
(14)∫Ωsi∂u∂twidΩ=D1∫Ωsi(∇2u)widΩ+α1∫ΩsiuwidΩ+∫ΩsiA(u,v)widΩ+∫Ωsif1(x,t)widΩ,∫Ωsi∂v∂twidΩ=D2∫Ωsi(∇2v)widΩ+α2∫ΩsivwidΩ+β∫ΩsiuwidΩ+∫ΩsiB(u,v)widΩ+∫Ωsif2(x,t)widΩ,
where wi is a Heaviside step used as the test function:
(15)wi(x)={1,x∈Ωsi0,x∉Ωsi.u and v are trial functions, and instead of the entire domain Ωsi we have considered a subdomain Ωsi located entirely at domain Ω, x=(x,y)T∈Ω⊂R2. The domain Ω is enclosed by Γ=ΓD∪ΓN, with boundary conditions
(16)u=u-,v=v-on⁡ΓD,(17)n·∇u=q-u,n·∇v=q-von⁡ΓN,
where n=(n1,n2)T is an outward unit normal vector of the boundary and n·∇u≡∂u/∂n. Condition (16) is often referred to as the Dirichlet boundary condition and (17) is referred to as the Neumann boundary condition. Let u~(x,t) and v~(x,t), which substitute u(x,t) and v(x,t), respectively, be the trial solutions:
(18)u~(x,t)=∑j=1Nϕj(x)u^j(t),v~j(x,t)=∑j=1Nϕj(x)v^j(t),A~=A(u~(x,t),v~(x,t)),B~=B(u~(x,t),v~(x,t)).
For internal nodes, from local integral equation (14), and using the MK (2), we have the following nonlinear equations:
(19)∑j=1N[∫Ωsiϕj(x)dΩ]∂u^j∂t=∑j=1N[D1(∫Lsiϕj,n(x)dΓ+∫ΓsDiϕj,n(x)dΓ)+(α1∫Ωsiϕj(x)dΩ)]u^j+D1∫LsNiq-udΓ+∫ΩsiA~dΩ+∫Ωsif1(x,t)dΩ,∑j=1N[∫Ωsiϕj(x)dΩ]∂v^j∂t=∑j=1N[D2(∫Lsiϕj,n(x)dΓ+∫ΓsDiϕj,n(x)dΓ)+(α2∫Ωsiϕj(x)dΩ)]v^j+∑j=1N[β∫Ωsiϕj(x)dΩ]u^j+D2∫LsNiq-vdΓ+∫ΩsiB~dΩ+∫Ωsif2(x,t)dΩ.
The following abbreviations have been used for the integral term:
(20)kij=∫Ωsiϕj(x)dΩ,hij1=D1(∫Lsi(∂ϕj(x)∂n)dΓ+∫LsDi(∂ϕj(x)∂n)dΓ)+α1∫Ωsiϕj(x)dΩ,hij2=D2(∫Lsi(∂ϕj(x)∂n)dΓ+∫LsDi(∂ϕj(x)∂n)dΓ)+α2∫Ωsiϕj(x)dΩ,lij=β∫Ωsiϕj(x)dΩ,f1i=D1∫LsNiq-udΓ+∫ΩsiA~dΩ+∫Ωsif1(x,t)dΩ,f2i=D2∫LsNiq-vdΓ+∫ΩsiB~dΩ+∫Ωsif2(x,t)dΩ.
The boundary and domain integrals are calculated using the Simpson and the Gauss-Legendre quadrature rules.We can rewrite (19) as
(21)∑j=1Nkij∂u^j∂t=∑j=1nhij1u^j+f1i,∑j=1Nkij∂v^j∂t=∑j=1n(hij2v^j+liju^j)+f2i.


### 4.1. Temporal Discretization

Equation (21) can be rewritten as
(22)K∂U^∂t=H1U^+F1.
Similarly, we have
(23)K∂V^∂t=H2V^+LU^+F2,
where
(24)K=[kij]N×N,H1=[hij1]N×N,H2=[hij2]N×N,L=[lij]N×N,F1=[f1i]N×1,F2=[f2i]N×1,U^=[u^1u^2⋯u^N]T,V^=[v^1v^2⋯v^N]T.The finite-difference approximation of the time derivatives of (22) and (23) in Crank-Nicolson method is given as follows:
(25)(2K-ΔtH1)U^k+1=(2K+ΔtH1)U^k+Δt(F1k+1+F1k),(2K-ΔtH2)V^k+1=(2K+ΔtH2)V^k(2K-ΔtH1)U^k+1+Δt(LU^k+1+LU^k)(2K-ΔtH1)U^k+1+Δt(F2k+1+F2k).
Because A and B are nonlinear functions of u and v, we solve them iteratively in each time step with replacing Ak+1 and Bk+1 by Ak and Bk, respectively, at zeroth iteration. Equation (25) converted into a set of nonlinear algebraic equations for unknowns U^k+1 and V^k+1.


## 4.1. Temporal Discretization

Equation (21) can be rewritten as
(22)K∂U^∂t=H1U^+F1.
Similarly, we have
(23)K∂V^∂t=H2V^+LU^+F2,
where
(24)K=[kij]N×N,H1=[hij1]N×N,H2=[hij2]N×N,L=[lij]N×N,F1=[f1i]N×1,F2=[f2i]N×1,U^=[u^1u^2⋯u^N]T,V^=[v^1v^2⋯v^N]T.The finite-difference approximation of the time derivatives of (22) and (23) in Crank-Nicolson method is given as follows:
(25)(2K-ΔtH1)U^k+1=(2K+ΔtH1)U^k+Δt(F1k+1+F1k),(2K-ΔtH2)V^k+1=(2K+ΔtH2)V^k(2K-ΔtH1)U^k+1+Δt(LU^k+1+LU^k)(2K-ΔtH1)U^k+1+Δt(F2k+1+F2k).
Because A and B are nonlinear functions of u and v, we solve them iteratively in each time step with replacing Ak+1 and Bk+1 by Ak and Bk, respectively, at zeroth iteration. Equation (25) converted into a set of nonlinear algebraic equations for unknowns U^k+1 and V^k+1.


## 5. Numerical Experiments

The analyzed domain isΩ=[0,1]×[0,1]. The error of u and v, which are presented in the numerical results, is represented by maximum relative error (MRE) and root mean square of relative error (RMSRE) of u and v, respectively, where
(26)MREu=max⁡|ui-u^iui|,i=1,2,…,N,RMSREu=1N∑i=1N(ui-u^iui)2,i=1,2,…,N.ui and u^i are the exact and computed values of u at point xi, respectively, and N is the number of nodes.Example 1.
The system of nonlinear PDEs in the regionΩ=[0,1]×[0,1] is given as follows:
(27)∂u∂t=(∂2u∂x2+∂2u∂y2)+u-u2v2+g1(x,y,t),∂v∂t=(∂2v∂x2+∂2v∂y2)+v-u2+uv+g2(x,y,t),
where
(28)g1(x,y,t)=e4x-6t-5e-2t+x+y,g2(x,y,t)=-4e-t+x-y+e-4t+2x+2y-e2x-3t.
The initial and Dirichlet boundary conditions are chosen in such a way that the exact solution is
(29)u(x,y,t)=e-2t+x+y,v(x,y,t)=e-t+x-y.
In this example, the governing equations resemble those of (1) with β=0 and D1=D2=1. The shown results have been obtained using N=25, 81, 256, and 441 nodal points, respectively, and Δt=0.1 at time instant t=1. Figures 1(a) and 1(b) show that the MRE of u and v increases as a function of the number of nodal points using m=3; meanwhile, using m=6 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 u and v decreases as a function of the number of nodal points using m=3, 6, and 10. Moreover, the errors of u and v using m=6, 10 are less than the RMSRE of u and v when using m=3. The profile of trial solution of u and v is similar to the profile of exact solution of u and v (see Figures 3(a), 3(b), 3(c), and 3(d)). Figures 3(e) and 3(f) reveal corresponding error profile of u and v. The errors of u and v satisfy the boundary conditions as well as the Kronecker delta property.Both errors ofu and v against the number of nodal points: m=10, t=1, and Δt=0.1. (a) MRE of u; (b) MRE of v.
(a)
(b)Both errors ofu and v against the number of nodal points: m=10, t=1, and Δt=0.1. (a) RMSRE of u; (b) RMSRE of v.
(a)
(b)Exact and trial solutions ofu and v using N=441, m=10, t=1, and Δt=0.1. (a) exact solution of u; (b) exact solution of v; (c) trial solution of u; (d) trial solution of v; (e) corresponding error profile of u; (f) corresponding error profile of v.
(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], Ω=[0,1]×[0,1]. Consider
(30)∂u∂t=D(∂2u∂x2+∂2u∂y2)-(A+1)u+u2v+B,∂v∂t=D(∂2v∂x2+∂2v∂y2)+Au-u2v,
with D=0.002, A=1/2, B=1, initial conditions
(31)u(x,y,0)=12x2-13x3,v(x,y,0)=12y2-13y3,
and Neumann boundary conditions
(32)∂u(0,y,t)∂x=∂u(1,y,t)∂x=∂u(x,0,t)∂y=∂u(x,1,t)∂y=0,∂v(0,y,t)∂x=∂v(1,y,t)∂x=∂v(x,0,t)∂y=∂v(x,1,t)∂y=0,
for which the exact solution is unknown. For small values of the diffusion coefficient D, if 1-A+B2>0 then the numerical solution of the Brusselator system converges to an equilibrium point (B,A/B) (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 (u*,v*)=(B,A/B)=(1,1/2). Figure 4 shows how the solution changes from initial condition to the steady state as t tends to the infinity. The experimental results are similar to those previously reported [7, 8].Table 1
The maximum and minimum values of the exact solution by using Crank-Nicolson method:m=10,Δt=0.1, and N=25.
t
u
v
u
max
⁡
u
min
⁡
v
max
⁡
v
min
⁡
1
0.5871
0.5258
0.3173
0.1384
2
0.7555
0.6944
0.4658
0.3458
3
0.8701
0.8060
0.5328
0.4794
4
0.9499
0.8985
0.5395
0.5306
5
0.9918
0.9631
0.5309
0.5212
6
1.0053
0.9953
0.5163
0.5072
7
1.0060
1.0048
0.5052
0.5007
8
1.0046
1.0031
0.5003
0.4990
9
1.0024
1.0011
0.4991
0.4990
10
1.0008
1.0001
0.4995
0.4993
↓
↓
↓
↓
↓
∞
1
1
0.5
0.5The approximate solutions profiles ofu and v for the Brusselator model at t=10 obtained by using Crank-Nicolson method: m=10, N=25, and Δt=0.1. (a) u at t=10; (b) v at t=10.
(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Δt 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.



---

*Source: 196041-2014-06-04.xml*