Meshless Local Weak form Method Based on a Combined Basis Function for Numerical Investigation of Brusselator Model and Spike Dynamics in the Gierer-Meinhardt System

Mohammad Ilatiand Mehdi Dehghan

1 Introduction

Reaction-diffusion models frequently arise in the study of chemical and biological systems and are naturally modeled by time dependent partial differential equations(PDEs).The nonlinear system of reaction-diffusion equations is composed of twodistinct terms:diffusion terms and reaction terms.Diffusion terms are involved with the random motion of each constituent due to turbulent nature of the flow field and reaction terms describe the interaction among the involved physical and biological species[Bhatt and Khaliq(2015);Zhao,Ovadia,Liu,Zhang,and Nie(2011)].A two-component reaction-diffusion system with general reaction termsLandNhas the following form:

subject to the no- flux boundary conditions

and initial conditions

whereu(x，t)andv(x，t)denote the concentration of two reactants,d1andd2(diffusion coefficients)are constants,u0andv0are known functions.Here,Ω represents the spatial domain of the equation with?Ω as the boundary and?2denotes the Laplacian associated with diffusion of the species.Based on many specific forms of the reaction terms,various models have been proposed for the pattern formation.Here,we investigate Brusselator model[Nicolis and Prigogine(1977)]and Gierer-Meinhardt(GM)model[Gierer and Meinhardt(1972)].

1.1 Brusselator model

One of the important reaction-diffusion equations is known as Brusselator system.The Brusselator system is used to describe mechanism of chemical reaction diffusion with non-linear oscillations.The importance of oscillations in biochemical systems has been emphasized by a number of authors.For instance,Turing[Turing(1952)]showed that when certain reactions are coupled with the process of diffusion,it is possible to obtain a stable spatial pattern(this laid the foundation of the theory of morphogenesis).The so-called Brussels school[Lefever and Nicolis(1971);Nicolis and Prigogine(1977);Prigogine and Lefever(1968);Tyson(1973)]developed and analysed the behaviour of a non-linear oscillator associated with the chemical system

whereδandρare input chemicals,DandEare output chemicals andUandVare intermediates.The chemical system(4),known as the“Brusselator”system,is important in that it admits limit-cycle oscillations and yet contains only two dependent variables(UandV)thus enabling the use of two-dimensional mathematical systems[Lefever and Nicolis(1971);Twizell,Gumel,and Cao(1999)].Letuandvbe the concentrations ofUandV,respectively,and assume that the concentrations of the input compoundsδandρare held constant during the reaction process.Then one obtains the following system of reaction-diffusion equations:

whereδ,ρ,and diffusion coefficientsα1andα2are positive constants.The parameterρis often chosen as a parameter for studying bifurcation.This model has been referred to as the trimolecular model or Brusselator.It represents a useful model for study of cooperative processes in chemical kinetics.Such a trimolecular reaction step arises in the formation of ozone by atomic oxygen via a triple collision.It arises also in enzymatic reactions,and in plasma and laser physics in multiple couplings between certain modes[Adomian(1995);Tyson(1973)].

It is well known that from[Twizell,Gumel,and Cao(1999)],for small values of the diffusion coefficientsα1andα2,the steady state solution of the Brusselator system(5)converges to equilibrium pointif 1-ρ+δ2＞0.

The analytical solution of the reaction-diffusion Brusselator system is not yet known and therefore it got interest from numerical point of view.In recent years,much attention has been paid in literature to the development of numerical schemes for the numerical solutions of reaction-diffusion Brusselator system such as secondorder finite difference scheme[Twizell,Gumel,and Cao(1999)],decomposition method of Adomian[Adomian(1995);Wazwaz(2000)],dual-reciprocity boundary element method[Ang(2003)],Runge-Kutta-Chebyshev method[Verwer,Hundsdorfer,and Sommeijer(1990)],collocation method using the radial basis functions[Siraj-ul-Islam and Haq(2010)],differential quadrature method[Mittal and Jiwari(2011a,b)],modified cubic B-spline differential quadrature method[Jiwari and Yuan(2014)],homotopy perturbation method[Kumar,Khanm,and Yildirim(2012)],alternating direction implicit(ADI)orthogonal spline collocation(OSC)method[Fernandes and Fairweather(2012)],local integral equation method[Shirzadi,Sladek,and Sladek(2013a)],moving finite element method[Hu,Qiao,and Tang(2012)],method of lines[Mohammadi,Mokhtari,and Schaback(2014)]and locally extrapolated exponential time differencing LOD scheme[Bhatt and Khaliq(2015)].

1.2 Gierer-Meinhardt model

Among the reaction-diffusion models,Gierer-Meinhardt model is the most wellknown reaction-diffusion system of activator-inhibitor type[Qiao(2008)].It has been widely used to model localization processes in nature,such as cell differentiation and morphogenesis[Harrison and Holloway(1995);Meinhardt(1982)],and the formation of sea-shell patterns[Meinhardt(1995)].This model not only generates spatial patterns but also exhibits size regulation,a phenomenon that occurs in many developmental systems such as head development in the Hydra[Gierer and Meinhardt(1972);Ward,Mcinerney,and Houston(2002)].The dimensionless GM model can be written as[Fernandes and Fairweather(2012);Qiao(2008)]

whereu,v,ε,κ,λ(x),V(x)andμrepresent the scaled activator concentration,inhibitor concentration,activator diffusivity,inhibitor diffusivity,inhibitor decay rate,activator decay rate,and inhibitor time constant,respectively[Qiao(2008)].The exponents(p，q，r，s)are assumed to satisfy[Qiao(2008)]

Numerical studies by Meinhardt[Meinhardt(1982)]have revealed that whenεis small andκis finite,GM seems to have stable stationary states with the property that the activator is mainly concentrated inKpeaks which are each placed nearKdifferent points in Ω whose locations satisfy suitable conditions[Wei and Winter(2001)].Moreover,forε?1,many studies of GM model have shown the spike patterns become narrower and narrower whenε→0[Meinhardt(1982);Qiao(2008);Wei and Winter(2001)].In fact,their spatial extension is of the orderO(ε).And the spike patterns also have various dynamical behaviors,such as the drift of the center of the spikes,the oscillation of the height of the spikes,even the splitting of the spikes.So,in the mesh based methods,very fine meshes over the spatial extension of the spikes are needed to resolve this problem[Qiao(2008)].

For two-dimensional GM model,some analysis works can be found in[Chen and Kowalczyk(2001);Kolokolnikov and Ward(2004);Wei and Winter(1999,2001)].However,because of the extremely large computational cost,there are only a few works in numerical simulations for spike dynamics in 2D,see[Fernandes and Fairweather(2012);Harrison and Holloway(1995);Kolokolnikov,Sun,Ward,and Wei(2006);McCourt,Dovidio,and Gilbert(2008);Qiao(2008);Ward,Mcinerney,andHouston(2002)].Numerical difficulties in simulating two-dimensional GM model also lie in that there are different orders of errors:the error in spike height isthe error in spike location isO(ε),the critical threshold forκisand the time evolution for spikes is[Qiao(2008)].As mentioned in[Qiao(2008)],the traditional finite element method(FEM)can not resolve the spike dynamics for very smallε[Qiao(2008)].

1.3 The main aim of the current paper

In this paper we propose a meshless local weak form method based on new combined shape function for numerical solution of reaction-diffusion Brusselator and Gierer-Meinhardt systems.The meshless local weak form method is a truly meshless method,which requires no elements or background cells,for either the interpolation or the integration purposes.This concept was first proposed by Atluri and Shen[Atluri and Zhu(1998)],and later discussed in depth in[Atluri and Shen(2002a,b)].The most significant difference between this method and the finite element method or any other meshless method is that the local weak forms are generated on overlapping local sub-domains,instead of using the global weak form.Integration of the weak form is performed in local sub-domains with simple geometrical shapes,therefore no elements or background cells are necessary either for interpolation purposes or for integration purposes.For other investigations on the meshless local weak form methods we refer to[Atluri and Zhu(2000);Dehghan,Abbaszadeh,and Mohebbi(2014);Dehghan and Mirzaei(2008,2009);Dehghan and Salehi(2013,2014);Dong,Alotaibi,Mohiuddine,and Atluri(2014);Mirzaei and Dehghan(2010);Shirzadi,Sladek,and Sladek(2013b);Sladek,S-ladek,and Hon(2006);Sladek,Sladek,Zhang,and Schanz(2006);Sladek,Stanak,Han,Sladek,and Atluri(2013);Taleei and Dehghan(2014);Zhang,He,Dong,Li,Alotaibi,and Atluri(2014)].

It is well-known that,in the classical MLPG,the moving least squares approximation is used as a shape function.In the present paper,a new shape function is developed as a linear interpolating function of radial point interpolation(RPI),moving least squares(MLS)and moving Kriging(MK)shape functions.This new combined shape function inherits the properties of RPI,MLS and MK shape functions and is controlled by control parameters,which take different values in the domain[0，1].Therefore this shape function is used as a shape function in the proposed meshless local weak form method.

The main aim of this paper is to show that the new combined basis function can be used as a shape function in meshless local weak form methods and leads to flexibility of the method from the perspective of concurrent use of different shape functions.In this paper,the combined basis function is used with different valuesof control parameters and it is showed that better results can be obtained by using suitable values for control parameters.

The organization of the rest of this paper is as follows:in Sections 2,the new shape function is introduced. A time stepping method is described in Section 3. In Section 4,full discritization based on meshless local weak forms is provided.In Section 5,results of numerical experiments are presented.Finally Section 6 completes the structure of this paper with a brief conclusion.

2 Development of a new combined shape function

In this section,a new shape function,which interpolates the radial point interpolation(RPI),moving least squares(MLS)and moving Kriging(MK)shape functions,is presented.Consider a local support domain with a set of arbitrarily scattered points xi，i=1，...，n,wherenis the number of nodes in the local support domain.The approximation of a functionu(x)in support domain can be defined by

whereψi(x)stands for the shape function and is defined as follows

A brief description of RPI,MLS and MK shape functions is as follows:

where

In the above relations,ri(x)is radial basis function(RBF),pj(x)is the monomial in the space coordinates x=[x，y]T,mis the number of monomial basis functions andR(xi，xj)is the correlation function between any pair of nodes located at xiand xj.In the present work,RBFs and correlation function are defined as follows:

wheredcis the average nodal spacing near the point of interest x,αcandqare two arbitrary real numbers of dimensionless shape parameters as suggested by Liu[Liu(2002)],andθis correlation parameter.As mentioned in[Liu and Gu(2005)],the recommended shape parameters for the MQ-RBF areq=1.03 and 1≤αc≤7.The quality of the shape function is heavily influenced by the correlation parameterθ.As mentioned in[Zheng and Dai(2011)],the optimalθis dependent on the number of nodes in the compact support,empirical formula is obtained

where ω is a constant,h is the average distance of the nodes in the support domain.It is a good choice to take ω∈[0.03，0.2].For more details about RPI,MLS and MK shape functions we refer to[Atluri and Zhu(1998);Bui,Nguyen,and Zhang(2011);Chen and Liew(2011);Dehghan,Abbaszadeh,and Mohebbi(2016);Dehghan and Ghesmati(2010);Lancaster and Salkauskas(1981);Liu and Gu(2005);Mirzaei and Dehghan(2010);Salehi and Dehghan(2013);Taleei and Dehghan(2015);Zheng and Dai(2011)].

The proposed shape function inherits the properties of RPI,MLS and MK shape functions.A key property of the combined shape function is that the proposed shape function can reproduce any function in the basis exactly.In particular,if a linear basis is employed to construct the shape functions,all constants and linear terms can then be reproduced exactly,i.e.

Eq.(25)is well known as the consistency property.

3 Discretizing the temporal space

For discretization of time variable,we need to some preliminary.We de fi ne

tk=kτ，k=0，1，2，...，M，

To deal with the time derivative,a time-stepping method based on famous Crank-Nicolson scheme[Dehghan(2006)]is employed.

where uk=u(x，kτ)and vk=v(x，kτ).

Since L1and L2are linear functions,they can be decomposed as follows

4 Full discritization based on meshless local weak form

The weak forms are used to derive a set of algebraic equations through a numerical integration process over the domain of the problem,globally or locally.The use of the global weak-form requires the system of equations in the global integral form to be satisfied over the entire problem domain,and hence,a set of background cells has to be used for the numerical integration.Therefore these methods are not truly meshless methods.To avoid the use of global background cells,the so-called local weak-form methods have been developed.When a local weak-form is used for a field node,the numerical integrations are carried out over a local quadrature domain defined for the node.Therefore,no global background mesh is required[Liu and Gu(2005)].

In the local weak form methods,around each xia small sub-domainis chosen such that integrations overare comparatively cheap.The local subdomains overlap each other,and cover the whole global domain Ω.The local subdomains could be of any geometric shape and size.For simplicity they are taken to be of circular shape.The local weak form of the approximate equations(28)for

whereθ?is a test function.Using

and the divergence theorem,Eqs.(29)yield the following expression

is chosen as the test function in each sub-domain,then the local weak forms(31)are transformed into the following simple local integral equations

For boundary point xi,respectively,whereis a part of the local boundary located on the global boundary andis the other part of the local boundary over which no boundary conditions are specified,Therefore the local weak form equationsfor boundary points are

By imposing the natural boundary conditions,we obtain

By substituting Eqs.(35)into(34),the local weak form for boundary points simplified as followsApplying the new combined approximation(8)for the unknown functions,the local integral equations(33)and(36)are transformed into a system of algebraic equations with unknown quantities.We suppose that,the unknownsuandvare approximated as follows:

Substituting the approximations(37)and(38)into the local integral equations(33)and(36)yields:

where

The second integrals in(42)and(43)are approximated using the combined approximation formula(8)and then the accuracy of these approximations becomes better by applying the predictor-corrector algorithm.

A brief description of predictor-corrector algorithm is as follows.

Predictor-corrector algorithm

It should be noted that,in the current section,the no- fl ux boundary conditions(2)are considered as boundary conditions.In the case of Dirichlet boundary condition,we can employ a penalty parameter to impose the Dirichlet boundary conditions[Atluri and Zhu(1998)]or impose it directly.

5 Numerical results

In numerical results,we use the quartic spline weight function in the constructing MLS shape function.

where andrsis the radius of the local support domain.The parameterrsshould be large enough to ensure the regularity of the moment matrixPWPTin MLS approximation[Dehghan and Mirzaei(2008)].Then the support sizersis a very important parameter in meshless methods.It is related to both accuracy of the solution,as well as the computational efficiency.

In this section,all of the numerical solutions are obtained by means of new combined shape function.Furthermore,the computational results are compared for different values of controlling parametersμ1,μ2andμ3.These comparisons show that,the new combined shape function allows to get better results.All of the computations are carried out with parametersq=1.03,ω=0.12,rq=0.75handrs=2.7h.

5.1 Brusselator model

5.1.1 Test problem 1

Consider the Brusselator system (5) together with the Dirichlet boundary conditions on the unit square domain Ω=[0，1]×[0，1].Let the parameter values beρ=1,δ=0 andα1=α2=0.25.The initial and boundary conditions are extracted from the exact solutions[Ang(2003)]

In this example,the numerical solutions are obtained withh=1/20,τ=0.001 andαc=7.The maximum error normL∞for the componentsuandvwith various values of the controlling parameters are shown in Tab.1.

5.1.2 Test problem 2

In this example we solve the Brusselator system(5)on the unit square domain Ω=[0，1]×[0，1]with no- flux boundary conditions for bothuandv.Initial conditions are taken as follows[Hu,Qiao,and Tang(2012);Verwer,Hundsdorfer,andSommeijer(1990)]

Table 1: L∞errors for different values of controlling parameters

Let the parameter values be given byρ=3.4,δ=1,α1=α2=0.002.In this example,the numerical solutions are obtained withh=1/30,τ=0.001,μ1=0.60,μ2=0.05,μ3=0.35 andαc=1.The time evolution of the concentration of activatoruandvat different times is shown in Fig.1,Fig.2 and Fig.3.It is found that the concentration pro files ofuandvare similar to[Hu,Qiao,and Tang(2012);Verwer,Hundsdorfer,and Sommeijer(1990)].

5.1.3 Test problem 3

Consider the Brusselator system(5)on the unit square domain Ω=[0，1]×[0，1]with no- flux boundary conditions for bothuandv.Initial conditions are considered as[Mohammadi,Mokhtari,and Schaback(2014)]

In this example,the numerical solutions are obtained withh=1/20,τ=0.001,μ1=0.6,μ2=0,μ3=0.4 andαc=1.The concentration pro files ofuandvatT=0 andT=5 with the parametersρ=1,δ=2,andα1=α2=0.002 are shown in Fig.4 and Fig.5.From Fig.5,it can be seen that the numerical values ofuandvat each collocation point approach to 2 and 0.5,respectively.These results show an agreement that for small values of the diffusion coefficientsα1andα2,astincreases,whenever 1-ρ+δ2＞0.The plots of the values ofuandvat the collocation point(0.3，0.3)versus time are shown in Fig.6.It canbe noted from Fig.6,that(u(0.3，0.3)，v(0.3，0.3))→(2，0.5)ast→∞.The results show an agreement with the results of[Mohammadi,Mokhtari,and Schaback(2014)].

Figure 1: Time evolution of the activator u at different times with parameter values ρ=3.4,δ=1,α1=α2=0.002.

Figure 2: Time evolution of the activator v at different times with parameter values ρ=3.4,δ=1,α1=α2=0.002.

Figure 3: Time evolution of the activator v at different times with parameter values ρ=3.4,δ=1,α1=α2=0.002.

Figure 4: Initial concentration pro files of u and v.

The algorithm is repeated withρ=3.4 andδ=1 up to timeT=40.The concentrations pro files ofuandvatT=40 are shown in Fig.7.The plots of the values ofuandvat the collocation point(0.3，0.3)versus time are shown in Fig.8.The computed results reveal that whenever the parametersρandδare chosen such that 1-ρ+δ2＞0,the concentration pro files ofuandvconverge to the fixed pointand for values ofρandδsuch that 1-ρ+δ2＜0,the numerical method is seen not to converge to any fixed concentration.The results show an agreement with the results of[Mohammadi,Mokhtari,and Schaback(2014)].

5.1.4 Test problem 4

In this example,we consider the Brusselator system(5)on the unit square domain Ω=[0，1]×[0，1]with no- flux boundary conditions for bothuandv.Initial conditions are taken as follows[Ang(2003)]

Let the parameter values be given byρ=0.5,δ=1,α1=α2=0.002.In this example,computations are carried out with the parametersh=1/20,τ=0.001,μ1=1,μ2=μ3=0 andαc=1.The concentration pro files ofuandvatT=0

andT=10 are shown in Fig.9.Also the plots of the values ofuandvat the collocation point(0.5，0.5)versus time are shown in Fig.10.It is well known that for small values of the diffusion coefficients,if 1-ρ+δ2＞0 then the steady state solution of the Brusselator system converges to equilibrium pointThis fact is con firmed by Fig.9 and Fig.10.The results show an agreement with the results of[Ang(2003)].

Figure 5: Time evolution of the activators u and v at T=5 with ρ=1,δ=2,α1=α2=0.002.

Figure 6: Plots of u(0.3，0.3)and v(0.3，0.3)versus time with ρ=1,δ=2,α1=α2=0.002.

Figure 7: Plots of u and v at T=40 with ρ=3.4,δ=1,α1=α2=0.002.

Figure 8: Plots of u(0.3，0.3)and v(0.3，0.3)versus time with ρ=3.4,δ=1,α1=α2=0.002.

Figure 9: Plots ofuandvatT=0andT=10withρ=0.5,δ=1,α1=α2=0.002.

5.1.5 Test problem 5

Consider the diffusion-free Brusselator system corresponding toα1=α2=0 on the unit square domain Ω=[0，1]×[0，1]with no- flux boundary conditions for bothuandv.Computations are carried out with parametersh=1/20,τ=0.001,μ1=μ2=0,μ3=1 andαc=1.Parametersρandδtake different values and algorithm is tested up to timeT=60 with various values of the initial conditions taken from-8≤u0，v0≤8.Phase portraits forδ=0.5 withρ=2.5，2，1.2，0.2 are depicted in Fig.11,Fig.12,Fig.13 and Fig.14,respectively.From Fig.11 and Fig.12,it can be noted that the limit cycles do not exist astincreases.It can also be noted form Fig.13 and Fig.14 that limit cycle exists and the solution sequence converges to the fixed pointastincreases.The results show an agreement with the results of[Twizell,Gumel,and Cao(1999)]that limit cycle does not exist whenever 1-ρ+δ2＜0 and limit cycle exists if 1-ρ+δ2＞0.

Figure 10: Plots of u(0.5，0.5)and v(0.5，0.5)versus time with ρ=0.5,δ=1,α1=α2=0.002.

5.2 Gierer-Meinhardt model

For(p，q，r，s)=(2，1，2，0),m=1,V(x)=0 andλ(x)=1,we solve the system(6)on the Ω=[-1，1]×[-1，1]with no- flux boundary conditions applied to bothuandv.Let the parameter values be given byε=0.04 andμ=0.1.The initial conditions are chosen as[McCourt,Dovidio,and Gilbert(2008)]

Figure 11: Solution of diffusion-free system with ρ=2.5,δ=0.5.

Figure 12: Solution of diffusion-free system with ρ=2,δ=0.5.

Figure 13: Solution of diffusion-free system with ρ=1.2,δ=0.5.

Figure 14: Solution of diffusion-free system with ρ=0.2,δ=0.5.

Figure 15: Contour plots of time evolution of the activator u at different times with κ=0.0128.

Figure 16: Contour plots of time evolution of the activator u at different times with κ=0.0128.

The numerical solutions are obtained withh=2/40,τ=0.0001,μ1=0.5,μ2=0,μ3=0.5 andαc=3.

5.2.1 Test problem 6

In this example,we letκ=0.0128.Contour plots of time evolution of the activatoruat different times are shown in Fig.15 and Fig.16.From the numerical simulations,we found that the spike at the center will be splitting and the domain will be filled with spikes astincreases.In Fig.15 and Fig.16,we can find that the spike begins with a ring whent=20.The ring becomes bigger and bigger.Then the

ring becomes almost square and subsequently the four edges of the square become sunken in the middle and then the whole annulation splits to spikes which fi ll the domain with some symmetry.It is clear from our graphs that the dynamics ofufollow the same general pattern displayed in[Fernandes and Fairweather(2012);McCourt,Dovidio,and Gilbert(2008)].

Figure 17: Contour plots of time evolution of the activator u at different times with κ=0.0152.

Figure 18: Contour plots of time evolution of the activator u at different times with κ=0.0152.

5.2.2 Test problem 7

In this example,we letκ=0.0152.Contour plots of time evolution of the activatoruat different times are shown in Fig.17 and Fig.18.From Fig.17 and Fig.18,we can find that,the spike splits into two spikes spreading in thexdirection and becomes symmetric(t=140).Then,each of the spikes splits,spreading in theydirection,and maintains symmetry(t=220).Next,each of the four spikes splits into two along thexdirection,and the eight spikes arrange themselves symmetrically about the center(t=520).Finally,att=570,the inner spikes split and the 12 spikes arrange themselves symmetrically about the center.It is clear from our graphs that the dynamics ofufollow the same general pattern displayed in[McCourt,Dovidio,and Gilbert(2008)].

6 Conclusion

In this paper,we have proposed a new combined shape function.Based on this shape function,the meshless local weak form method has been applied for solving two-dimensional time-dependent non-linear Brusselator and Gierer-Meinhardt systems.The new combined shape function is developed as a linear interpolating function of radial point interpolation(RPI),moving least squares(MLS)and moving kriging(MK)shape functions.This new shape function inherits the properties of RPI,MLS and MK shape functions and is controlled by control parameters,which take different values in the domain[0，1].The results show that good accuracy can be obtained with using different values for controlling parameters.Finally,we believe that the new combined shape function provides the opportunity of using different shape functions,simultaneously and this leads to fl exibility of the method.

Acknowledgement:The authors thank the reviewers for their useful comments and suggestions that improved the paper.

Adomian,G.(1995):The diffusion-Brusselator equation.Comput.Math.Appl.,vol.29,pp.1–3.

Ang,W.T.(2003):The two-dimensional reaction-diffusion Brusselator system:a dual-reciprocity boundary element solution.Eng.Anal.Bound.Elem.,vol.27,pp.897–903.

Atluri,S.N.;Shen,S.P.(2002):The meshless local Petrov-Galerkin(MLPG)method.Tech.Science Press.

Atluri,S.N.;Shen,S.P.(2002):The meshless local Petrov-Galerkin(MLPG)method:A simple and less-costly alternative to the finite element methods.CMES:Computer Modeling in Engineering&Sciences,vol.3,no.1,pp.11–51.

Atluri,S.N.;Zhu,T.(1998):A new meshless local Petrov-Galerkin(MLPG)approach in computational mechanics.Comput.Mech.,vol.22,pp.117–127.

Atluri,S.N.;Zhu,T.(2000):The meshless local Petrov-Galerkin(MLPG)approach for solving problems in elasto-statics.Comput.Mech.,vol.25,pp.169–179.

Bhatt,H.P.;Khaliq,A.Q.M.(2015):The locally extrapolated exponential time differencing LOD scheme for multi-dimensional reaction-diffusion systems.J.Comput.Appl.Math.,vol.285,pp.256–278.

Bui,T.Q.;Nguyen,M.N.;Zhang,C.(2011):A moving Kriging interpolationbased element-free Galerkin method for structural dynamic analysis.Comput.Methods Appl.Mech.Eng.,vol.200,pp.1354–1366.

Chen,L.;Liew,K.M.(2011):A local Petrov-Galerkin approach with moving Kriging interpolation for solving transient heat conduction problems.Comput.Mech.,vol.47,pp.455–467.

Chen,X.;Kowalczyk,M.(2001):Dynamics of an interior spike in the Gierer-Meinhardt system.SIAM J.Math.Anal.,vol.33,pp.172–193.

Dehghan,M.(2006):Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices.Math.Comput.Simul.,vol.71,pp.16–30.

Dehghan,M.;Abbaszadeh,M.;Mohebbi,A.(2014):Numerical solution of system of N-coupled nonlinear Schrodinger equations via two variants of the meshless local Petrov-Galerkin(MLPG)method.CMES:Computer Modeling in Engineering&Sciences,vol.100,no.6,pp.399–444.

Dehghan,M.;Abbaszadeh,M.;Mohebbi,A.(2016):The use of element free Galerkin method based on moving Kriging and radial point interpolation techniques for solving some types of Turing models.Eng.Anal.Bound.Elem.,vol.62,pp.93–111.

Dehghan,M.;Ghesmati,A.(2010):Numerical simulation of two-dimensional sine-Gordon solitons via a local weak meshless technique based on the radial point interpolation method(RPIM).Comput.Phys.Commun.,vol.181,pp.772–786.

Dehghan,M.;Mirzaei,D.(2008):The meshless local Petrov-Galerkin(MLPG)method for the generalized two-dimensional non-linear Schrodinger equation.Eng.Anal.Bound.Elem.,vol.32,pp.747–756.

Dehghan,M.;Mirzaei,D.(2009):Meshless local Petrov-Galerkin(MLPG)method for the unsteady magneto hydrodynamic(MHD) flow through pipe with arbitrary wall conductivity.Appl.Numer.Math.,vol.59,pp.1043–1058.

Dehghan,M.;Salehi,R.(2013):A mesh free weak-strong(MWS)form method for the unsteady magneto hydrodynamic(MHD) flow in pipe with arbitrary wall conductivity.Comput.Mech.,vol.52,pp.1445–1462.

Dehghan,M.;Salehi,R.(2014):A meshless local Petrov-Galerkin method for the time-dependent Maxwell equations.J.Comput.Appl.Math.,vol.268,pp.93–110.

Dong,L.;Alotaibi,A.;Mohiuddine,S.A.;Atluri,S.N.(2014):Computational methods in engineering:a variety of primal and mixed methods,with global and local interpolations,for well-posed or ill-Posed BCs.CMES:Computer Modeling in Engineering&Sciences,vol.99,no.1,pp.1–85.

Fernandes,R.I.;Fairweather,G.(2012):An ADI extrapolated Crank-Nicolson orthogonal spline collocation method for nonlinear reaction-diffusion systems.J.Comput.Phys.,vol.231,pp.6248–6267.

Gierer,A.;Meinhardt,H.(1972):A theory of biological pattern formation.Kybernetik,vol.12,pp.30–39.

Harrison,L.;Holloway,D.(1995):Order and localization in reaction-diffusion pattern.Physica A,vol.222,pp.210–233.

Hu,G.;Qiao,Z.;Tang,T.(2012):Moving fi nite element simulations for reaction-diffusion systems.Advances in Applied Mathematics and Mechanics,vol.4,pp.365–381.

Jiwari,R.;Yuan,J.(2014):A computational modeling of two-dimensional reaction-diffusion Brusselator system arising in chemical processes.J.Math.Chem.,vol.52,pp.1535–1551.

Kolokolnikov,T.;Sun,W.;Ward,M.J.;Wei,J.(2006):The stability of a stripe for the Gierer-Meinhardt model and the effect of saturation.SIAM J.Appl.Dyn.Syst.,vol.5,pp.313–363.

Kolokolnikov,T.;Ward,M.J.(2004):Bifurcation of spike equilibria in the nearshadow Gierer-Meinhardt model.Dis.Cont.Dyn.Sys.B,vol.4,pp.1033–1064.

Kumar,S.;Khanm,Y.;Yildirim,A.(2012):A mathematical modeling arising in the chemical systems and its approximate numerical solution.Asia Pac.J.Chem.Eng.,vol.7,pp.835–840.

Lancaster,P.;Salkauskas,K.(1981):Surfaces generated by moving least squares methods.Math.Comp.,vol.37,pp.141–158.

Lefever,R.;Nicolis,G.(1971):Chemical instabilities and sustained oscillations.J.Theor.Biol.,vol.30,pp.267–284.

Liu,G.R.(2002):Meshfree Methods:Moving Beyond the Finite Element Method.CRC Press,Boca Raton,USA.

Liu,G.R.;Gu,Y.T.(2005):An Introduction to Meshfree Methods and Their Programming.Berlin:Springer.

McCourt,M.;Dovidio,N.;Gilbert,M.(2008):Spectral methods for resolving spike dynamics in the Gierer-Meinhardt model.Commun.Comput.Phys.,vol.3,pp.659–678.

Meinhardt,H.(1982):Models of Biological Pattern Formation.Academic Press,London.

Meinhardt,H.(1995):The Algorithmic Beauty of Sea Shells.Springer-Verlag,Berlin.

Mirzaei,D.;Dehghan,M.(2010):Meshless local Petrov-Galerkin(MLPG)approximation to the two dimensional sine-Gordon equation.J.Comput.Appl.Math.,vol.233,pp.2737–2754.

Mittal,R.C.;Jiwari,R.(2011):Numerical solution of two-dimensional reactiondiffusion Brusselator system.Appl.Math.Comput.,vol.217,pp.5404–5415.

Mittal,R.C.;Jiwari,R.(2011):Numerical study of two-dimensional reactiondiffusion Brusselator system by differential quadrature method.Int.J.Comput.Methods Eng.Sci.Mech.,vol.12,pp.14–25.

Mohammadi,M.;Mokhtari,R.;Schaback,R.(2014):A meshless method for solving the 2D Brusselator reaction-diffusion system.CMES:Computer Modeling in Engineering&Sciences,vol.101,no.2,pp.113–138.

Nicolis,G.;Prigogine,I.(1977):Self-Organization in Non-Equilibrium Systems.Wiley/Interscience,New York.

Prigogine,I.;Lefever,R.(1968):Symmetries breaking instabilities in dissipative systems II.J.Phys.Chem.,vol.48,pp.1695–1700.

Qiao,Z.(2008):Numerical investigations of the dynamical behaviors and instabilities for the Gierer-Meinhardt system.Commun.Comput.Phys.,vol.3,pp.406–426.

Salehi,R.;Dehghan,M.(2013):A moving least square reproducing polynomial meshless method.Appl.Numer.Math.,vol.69,pp.34–58.

Shirzadi,A.;Sladek,V.;Sladek,J.(2013):A local integral equation formulation to solve coupled nonlinear reaction-diffusion equations by using moving least square approximation.Eng.Anal.Bound.Elem.,vol.37,pp.8–14.

Shirzadi,A.;Sladek,V.;Sladek,J.(2013):A meshless simulations for 2D nonlinear reaction-diffusion Brusselator system.CMES:Computer Modeling in Engineering&Sciences,vol.95,no.4,pp.259–282.

Siraj-ul-Islam,A.A.;Haq,S.(2010):A computational modeling of the behavior of the two-dimensional reaction-diffusion Brusselator system.Appl.Math.Model.,vol.34,pp.3896–3909.

Sladek,J.;Sladek,V.;Hon,Y.C.(2006):Inverse heat conduction problems by meshless local Petrov-Galerkin method.Eng.Anal.Bound.Elem.,vol.30,pp.650–661.

Sladek,J.;Sladek,V.;Zhang,C.;Schanz,M.(2006):Meshless local Petrov-Galerkin method for continuously nonhomogeneous linear viscoelastic solids.Comput.Mech.,vol.37,pp.279–289.

Sladek,J.;Stanak,P.;Han,Z.D.;Sladek,V.;Atluri,S.N.(2013):Applications of the MLPG method in engineering and sciences:a review.CMES:Computer Modeling in Engineering&Sciences,vol.92,no.5,pp.423–475.

Taleei,A.;Dehghan,M.(2014):Direct meshless local Petrov-Galerkin method for elliptic interface problems with applications in electrostatic and elastostatic.Comput.Methods in Appl.Mech.Eng.,vol.278,pp.479–498.

Taleei,A.;Dehghan,M.(2015):An efficient meshfree point collocation moving least squares method to solve the interface problems with nonhomogeneous jump conditions.Numer.Methods Partial Differential Eq.,vol.31,pp.1031–1053.

Turing,A.M.(1952):The chemical basis of morphogenesis.Philos.Trans.Roy.Soc.London B,vol.237,pp.37–72.

Twizell,E.H.;Gumel,A.B.;Cao,Q.(1999):A second-order scheme for the Brusselator reaction-diffusion system.J.Math.Chem.,vol.26,pp.297–316.

Tyson,J.(1973):Some further studies of nonlinear oscillations in chemical systems.J.Chem.Phys.,vol.58,pp.3919–3930.

Verwer,J.G.;Hundsdorfer,W.H.;Sommeijer,B.P.(1990):Convergence properties of the Runge-Kutta-Chebyshev method.Numer.Math.,vol.57,pp.157–178.

Ward,M.J.;Mcinerney,D.;Houston,P.(2002):The dynamics and pinning of a spike for are action-diffusion system.SIAMJ.Appl.Math.,vol.62,pp.1297–1328.

Wazwaz,A.M.(2000):The decomposition method applied to systems of partial differential equations and to the reaction-diffusion Brusselator model.Appl.Math.Comput.,vol.110,pp.251–264.

Wei,J.;Winter,M.(1999):On the two-dimensional Gierer-Meinhardt system with strong couping.SIAM J.Math.Anal.,vol.30,pp.1241–1263.

Wei,J.;Winter,M.(2001):Spikes for the two-dimensional Gierer-Meinhardt system:the weak coupling case.J.Nonlinear Sci.,vol.11,pp.415–458.

Zhang,T.;He,Y.;Dong,L.;Li,S.;Alotaibi,A.;Atluri,S.N.(2014):Meshless local Petrov-Galerkin mixed collocation method for solving cauchy inverse problems of steady-state heat transfer.CMES:Computer Modeling in Engineering&Sciences,vol.97,no.6,pp.509–553.

Zhao,S.;Ovadia,J.;Liu,X.;Zhang,Y.-T.;Nie,Q.(2011):Operator splitting implicit integration factor method for stiff reaction-diffusion-advection systems.J.Comput.Phys.,vol.230,pp.5996–6009.

Zheng,B.;Dai,B.(2011):A meshless local moving Kriging method for twodimensional solids.Appl.Math.Comput.,vol.218,pp.563–573.