Numerical Simulation of Reaction-Diffusion Systems of Turing Pattern Formation

Differential method and homotopy analysis method are used for solving the two-dimensional reaction-diffusion model. And the structure of the solutions is analyzed. Finally, the homotopy series solutions are simulated with the mathematical software Matlab, so the Turing patterns will be produced. Overall analysis and experimental simulation of the model show that the different parameters lead to different Turing pattern structures. As time goes on, the structure of Turing patterns changes, and the final solutions tend to stationary state.


Introduction
In time or space, patterns have nonuniform macroscopic structure with regularity.From the thermodynamic point of view, the nature of the pattern formation can be divided into two categories.One is presenting in the thermodynamic equilibrium conditions, such as the crystal structure of inorganic chemistry, the self-organic pattern formation of organic polymers and so on.The other is the station far from the thermodynamic equilibrium conditions, such as the ripples of the sea, the surface patterns of the animal, the strip clouds in the sky and so on.Reaction-diffusion system is one of the fundamental equations which describe the motion of the nature.It not only has a wide practical background, but also is used in many fields, for example, predator-prey model, spread of infectious diseases, migration of population and spread of forest fires.Its mathematical model is a special kind of parabolic partial differential equations.As for reaction-diffusion systems, the coupling of nonlinear dynamical and linear diffusion leads to spontaneously producing a variety of ordered or disordered pattern of the system.This is the pattern dynamics of the reaction-diffusion systems [1].According to the pattern dynamics, the ordered patterns can be divided into two categories: stationary state (such as Turing patterns) and traveling wave (such as spiral pattern).As for the Turing patterns of reaction-diffusion systems, in 1952, Turing [2] for the first time showed that the homogeneous system will be stable with little disturbance for the absence of proliferation, while it becomes unstable with the spatial disturbance for joining the steady-state diffusion.This is the Turing instability of reaction-diffusion equations.
The classical method to study the Turing patterns of reaction-diffusion is the analysis of linear stability method [3]- [8].Firstly, the nonlinear reaction-diffusion model is turned into linear model.Then, perturbation theory is utilized to study its solution which can be seen as the little perturbation of equilibrium state (i.e.uniform and stable state).And we can obtain the conditions of generating Turing pattern through studying the linear equations with perturbation by the analysis of stability method.Finally, through the numerical simulation, the reaction-diffusion model can obtain different Turing pattern structures.This method shows the relationship between the parameters of generating Turing pattern.But for the strong nonlinear problems, this method may not be applicable.
The range of parameters of the Turing pattern can be obtained by the analysis of linear stability method.Based on the parameters which limit to the range, this article solves the reaction-diffusion model by the use of the combination of differential method and homotopy analysis method.Then the changes of the mechanism of Turing patterns can be under control by the simulation of experimental data and analysis of utilizing this method.Finally, the study of the concrete case shows the feasibility and effectiveness of this method.

Homotopy Analysis Solution of Reaction-Diffusion Model
The general mathematical representation of two-dimensional reaction-diffusion model is as follows, ( Their definition are as follows: u and v are different reactant concentrations vector; µ is a control parameter; f and g represent the nonlinear dynamics function of the system; u D and v D describe the diffusion coefficient; . The discretization of the model takes the following form, And, , , , , , According to the thought of homotopy analysis method, we take the value and impose the initial guess solutions , ,0 i j u and , ,0 .Auxiliary adjusting parameters meet the conditions 0 u ≠  and 0 v ≠  , while auxiliary control functions meet the conditions ( ) ( ) , zero-order deformation equation is as follows: , The problem (8) shows that two cases.One is ,0  2), while embedded variable q is continuous changing from 0 to 1.

( )
, , i j t q ϕ and ( ) , , i j t q φ express the form of Taylor expansion as follows, ( !
The two sides of zero-order Equations ( 8) are solved m-order derivative about q and divided !m at the same time.Take 0 q = can obtain the m-order deformation equation as follows, ( , 0 0 The solutions of problem (10) are, ( , , , , 1 , , 1 0 , , , , 1 , , 1 0 e d e d In order to make Taylor series are convergent at 1 q = , we can choose the value of , . Obtain the solutions of problem (2) as follows, ,1 In summary, the m-order homotopy approximate solutions of problem ( 1) are , ,

Homotopy Analysis of Brusselator Reaction Diffusion Model
The typical Brusselator reaction diffusion model is as follows, ( ) where , a b are control parameters for the system; u D and v D represent the coefficient of diffusion and 0 ,0 x l y l ≤ ≤ ≤ ≤ ; h describes the space step (i.e. the step of x and y are the same length h) and l M h = .The uniform steady-state solution to the reaction-diffusion system is ( ) . Combined the significance of reaction-diffusion model and the steady-state solution of the homogeneous perturbation analysis, the guess of initial solution can be as follows, , ,0 1 2 0 cos cos 0 cos cos where δ is a perturbation parameter and 1 k and 2 k represent the wave number of x-direction and y-direc- tion, respectively.Base on the method above, three order approximate solution for the problem of two dimensional reaction diffusion model is as follows, ( , , 0 0 ( ) . The Figures 1-4 show the effective ranges of  in different conditions.The effective ranges in Table 1 show that we will get different effective ranges under different condition of different parameters.

Simulation of Turing Patterns about Reaction-Diffusion
Parameter  have the ability to regulate and control the convergent region of series solutions and the speed of convergence of the series solutions.The analysis of the basis function in the homotopy analysis method shows that when u h and v h are limited value, the homotopy series solutions are convergent.Choose the step 1 h = and other parameters

Conclusions
In this paper, a new method based on differential method and homotopy analysis method is used to solve the typical two-dimensional reaction diffusion model.Different shapes of Turing pattern can be obtained through Matlab mathematical software on experimental data simulation of the structure of the solution.And it is proved that the proposed method which to solving nonlinear reaction diffusion problems is feasible and effective.The new method not only reduces the dimension about differential in space, but also gets the analytical expression with physical parameters through the homotopy in time.It will facilitate the analysis of the influence of parameters variation on Turing pattern structure.The method which is differential in space and homotopy in time domain has enormous potential for solving definite solution of nonlinear partial differential, and it has great promotional value.
The innovation point of this article is making use of the new method to solve the two-dimensional reaction diffusion model.The new method is combining the differential method and homotopy analysis method.

2 ∇
is space of the research for the model.Discrete the reaction-diffusion model in discrete nodes , zero-order approximate solution of problem(2).Introduce the embedded variable

2 , 2 .
and n d are the functions which depend on the parameters 1 , , , a b k k δ .Brusselator model is a classical dissipative structure model, which has been studied by many researchers.The conditions for Turing bifurcation of Brusselator model is as follows, The Effective Range of  Parameters  have the ability to regulate and control the convergent region of series solution and the speed of convergence of the series solution.We can obtain an appropriate  through the effective range about curve  of physical quantity

Figure 1 .
Figure 1.The effective area about  of u, when 3, 9.5 a b = =.

Figure 2 .
Figure 2. The effective area about  of v, when 3, 9.5 a b = = .

Figure 3 .
Figure 3.The effective area about  of u, when 4, 10 a b = =.

Figure 4 .Table 1 .
Figure 4.The effective area about  of v, when 4, 10 a b = =.Table1.The effective area about  under different parameters.ParametersThe effective area about  of uThe effective area about  of v 3, 9.5 a b = = this case, as time goes by, the changes of Turing patterns structures show as Figures 5-12.

Figure 20 .
Figure 20.The patterns in space of v(t), when t = 10.