a r X i v :c o n d -m a t /0404583v 3 [c o n d -m a t .s t a t -m e c h ] 6 A u g 2004
Numerical schemes for continuum models of reaction-diffusion systems
subject to internal noi
Esteban Moro ∗
Grupo Interdisciplinar de Sistemas Complejos (GISC)and Departamento de Matem´a ticas,Universidad Carlos III de Madrid,Avda.de la Universidad 30,E-28911Legan´e s,Spain
(Dated:February 2,2008)We prent new numerical schemes to integrate stochastic partial differential equations which describe the spatio-temporal dynamics of reaction-diffusion (RD)problems under the effect of in-ternal fluctuations.The schemes conrve the nonnegativity of the solutions and incorporate the Poissonian nature of internal fluctuations at small densities,their performance being limited by the level of approximation of density fluctuations at small scales.We apply the new schemes to two different aspects of the Reggeon model namely,the study of its non-equilibrium pha transition and the dynamics of fluctuating pulled fronts.In the latter ca,our approach allows to reproduce quantitatively for the first time microscopic properties within the continuum model.
PACS numbers:05.40.-a,05.70.Ln,68.35.Ct
Continuum reprentations of the dynamics of spatially-extended systems subject to fluctuations is a very active area of rearch in statistical mechanics and nonlinear dynamics [1,2,3,4,5].This is becau they are frequently more tractable than discrete models,they can be put forward using simple symmetry arguments and applying conrvation laws,and therefore they pro-vide minimal reprentations of the obrved phenom-ena.Important instances are Langevin equations for the relaxational dynamics of equilibrium models [1],growth interface phenomena [4]or coar-grained descriptions of microscopic RD problems [5,6].Despite their apparent simplicity,most of the models can not be solved ana-lytically and one has to resort to approximate analytical techniques,or to numerical integration of the stochastic time-dependent t of equations using well established algorithms [7].In the important instance of RD sys-tems subject to internal fluctuations the configurations are given by a non-negative density field ρ(x,t )subject to fluctuations of typical strength
∂t
=D
∂2ρ
σρη(x,t ),(1)
where η(x,t )is a Gaussian white noi.The Reggeon model can be obtained under some approximations from the microscopic Master equations of RD microscopic models using well-known techniques [5,8].Heuristi-cally,Eq.(1)can also be considered as the simplest dy-namical equation for a coar-grained density field with鸡翅中怎么做好吃
σ=1/N ,N being the mean-field number of particles per site.The Reggeon model provides a minimal repre-ntation of the Directed Percolation (DP)universality class,which is currently regarded as paradigm of non-equilibrium systems with absorbing states [6]:if ¯ρ(t )is the mean density spatial average,there exists a critical value of σfor which (1)undergoes a transition between an active pha lim t →∞¯ρ(t )=0and an absorbing pha for which lim t →∞¯ρ(t )=0.
In addition,when σ=0Eq.(1)becomes the so called Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP)equa-tion [9],which displays pulled fronts in which the active pha invades the absorbing state [10,11,12,13].Simu-lations of microscopic particle models [11,14]have shown that the dynamics of pulled fronts are extremely nsitive to microscopic fluctuations at ρ≃1/N ,leading to strong correc
tions in the front properties when compared with tho of the FKPP equation.Since Eq.(1)is usually held as a continuum description of some particle models at the mesoscopic level (i.e.when ρ≫1/N )one might doubt that the Reggeon model describes correctly the behavior of pulled fronts subject to internal fluctuations.The ef-ficiency and accuracy of the numerical schemes propod here will allow us to show that Eq.(1)indeed incorpo-rates the ingredients to explain (even quantitatively)the phenomena obrved in particle models,thus providing also a minimal reprentation of pulled fronts subject to internal fluctuations.
To simplify the discussion,let us consider the simplest possible ca for the dynamics of a density subject to internal fluctuations:
dρ
σρη(t ).(2)Typical explicit or implicit methods bad on stochas-tic Taylor approximations of (2)immediately run into
problems,since they do not conrve the nonnegativity of ρ(t ).For example,the Euler approximation is [7]
ρt +∆t =ρt +aρt ∆t +
√王者荣耀法师
2
where ∆W t are random Gaussian numbers with zero mean and ∆t variance.Thus,there is a finite probability that ρt +∆t becomes negative,and the numerical integra-tion comes to a halt.In order to overcome this problem,Dickman propod an interesting solution bad on the
Euler
scheme (3)and the discretization of the possible values of ρt as
multiples
of ρ
∗=O (σ∆t )[16].Despite its success in reproducing the universality class exponents of DP using (1)and its application to other situations [16],Dickman’s algorithm is not really a numerical in-tegration of a continuum model.Moreover,no general study of its convergence and applicability for other situ-ations has been done yet.A more technical solution was propod by Schurz and coworker
外墙真石漆报价s [17]using Balanced Implicit Methods (BIM),in which implicit Euler meth-ods are ud to impo the nonnegativity of the solution.In the ca of Eq.(2)the BIM scheme reads [18]
ρt +∆t =
ρt +∆tρt +
√
1+
∆t )for approximations of individual trajec-tories and O (∆t )for moments of ρ(t )[7,17].
Another approach was taken by Pechenik and Levine [12]employing the exact conditional probability density (CDF)P (ρt +∆t |ρt )for the stochastic process satisfy-ing (2),which has been known for some time in econ-omy as the Cox-Ingersoll-Ross process [20].The CDF can be expresd in terms of modified Besl functions and,although it can be sampled numerically using re-jection or transformation methods [12],it is computa-tionally expensive.Here we propo a more efficient procedure,which is bad on the following:if we de-fine r d (t )= d i =1x 2i (t ),where x i (t )satisfies dx i /dt =
ax i /2+(σ/4)1/2
ηi (t )with ηi (t )independent white nois,then dr d /dt =dσ/4+ar d +(σr d )1/2η(t )which coin-cides with (2)in the limit d →0.Since the equation for each x i (t )is linear,r d (t )is the sum of squares of Gaussian random numbers with non-zero mean.Thus its probability distribution is related to the χ2distri-bution with d degrees of freedom [21].Specifically,we
find that ρt +∆t =r 0(t +∆t )=χ′2
0(λ)/(2k )where
k =2a/[σ(e a ∆t −1)],λ=ke a ∆t ρt and χ′2
0(λ)is a ran-dom number with a noncentral χ2distribution with zero degrees of freedom and noncentrality parameter λwho cumulative distribution function is given by [15,19,21]P [χ20(λ)≤x ]=
∞ j =1
(λ/2)j e −λ/22k
0if K =0, 2K
i =1z 2i
if K =0,
铁骨铮铮郭沫若(6)
where z i are independent Gaussian random numbers with
zero mean and unit variance.
Another interesting feature of the exact CDF for (2)is the fact that it converges asymptotically towards the Eu-ler approximation (3)when λ≃ρt /(σ∆t )≫1[19,21].However,for small λ,the Euler approximation under-estimates the large fluctuations prent in the exact so-lution of (2).This effect,which can be en in Fig.1,is related to the fact that the Gaussian approximation (3)of a Poisson random number (2)is only valid when the mean value is large enough [19].The failure of ap-proximations like (3)or (4)to reproduce large density fluctuations at small values of λintroduces an effective microscopic cutoffρ∗=O (σ∆t )in the numerical simu-lations below which the approximations break down.Although the scheme (4)can be easily generalized to integrate equations like (1),this is not the ca for the exact sampling of the CDF for (2).Thus,a splitting-step strategy for integrating Eq.(1)was propod in [12],where the time interval ∆t is split into two steps:(i)given ρt ,we u (6)to integrate (2)and get an intermediate value ˜ρt +∆t ;(ii)we take ˜ρt +∆t as the initial condition for
十八岁成人礼祝福语∂ρ/∂t =∂2ρ/∂x 2−ρ2
,producing ρt +∆t with the aid of any deterministic numerical algorithm.It can be proved that this splitting step method (SSM)converges towards the solutions of (1),its order of convergence being O (∆t )
3
1 / 2
σc (∆t )
t
∆t
ρFIG.2:Left:convergence analysis for the different algo-rithms.Points are the critical value of σas a fu
nction ∆t (below)and of ∆t 1/2(up)while dashed lines are linear fits to the data.System size is L =400.Right:time dependence of the mean density ¯ρ(t )at the critical point σ=σc (∆t )
with ∆t =10−2
obtained using the different algorithms (lines are shifted vertically for clarity).Thin line is the power law ¯ρ(t )∼t −δwith δ=0.1595.System size is L =1000.
both for realizations and for moments of ρ(t )[19].This means that the splitting-step method provides better ap-proximations than tho bad on Euler methods (like the Dickman and BIM algorithms)for any realization of the noi.This has significant conquences when charac-terizing the critical point,as will be shown below.In the following we apply the two methods propod here [BIM and the SSM using (6)]and the Dickman algorithm to the two problems for which Eq.(1)is archetypal [22].Study of the DP pha transition.To test the propod algorithms,we study the well known non-equilibrium pha transition that Reggeon model displays for mod-erate values of σ[16].At the critical point,the mean average density ¯ρ(t )≡1
∆t ),while the SSM has O (∆t )order of convergence.The improvement in the order of convergence comes with a price:the computer time needed for our numerical sim-ulations at the critical point (e
table I)indicate that methods bad on Euler approximations,despite having an effective microscopic cutoffat ρ∗=O (σ∆t ),are faster than the SSM,and thus could provide better strategies for integrating numerically equations for RD models clo to the critical point,where only accurate approximations of large length and time scales are needed.
Dynamics of fluctuating pulled fronts.When σ=0,equation (1)displays a wave-like solution (front)which travels with velocity v 0=2
√Dickman 2.551BIM
1.61 1.2Splitting-Step 1.457.6
initial conditions are given)[9,10].The dynamics of this
pulled front is verely affected when microscopic fluctu-ations clo to the absorbing state ρ=0are considered.Specifically,it has been obrved in particle models who mean field limit is given by the FKPP equation,that the front speed is universally modified as [10,11,14]
v N ≡lim
t →∞
x f (t )
ln 2N
,(7)
where x f (t )is the instantaneous position of the front,C is a positive constant and N is the number of particles per site [14].Moreover,the pulled front diffus with diffusion constant
D f,N ≡lim
t →∞
(x f (t )−v N t )2
ln 3N
,(8)
where C ′is a positive constant.Whereas the velocity correction can be easily understood becau microscopic fluctuations at ρ≃N −1provide an effective cutoffin the dynamics [11],the diffusion coeffic
ient ems to depend on the existence of relatively large fluctuations in the density at ρ≃N −1and on their slow relaxation by the pulled front dynamics [14].
As mentioned in the introduction,one might doubt that the large microscopic fluctuations at ρ=N −1ob-rved in particle models are correctly reproduced by such type of equation like (1).Note,however,that the relationship between particle models and the Reggeon field model is deeper than at the coar-grained level.Specifically,in [13]it was shown that there is an exact duality transformation between the A ↔A +A micro-scopic particle model and the so called stochastic FKPP equation,which is similar to the Reggeon model but with a σρ(1−ρ)≃
√
静待黎明
西兰花蒸多久熟
4
log 10 N
10
101010
10
10
10D f , N
FIG.3:Front diffusion coefficient as a function of σfor the
different algorithms and different ∆t ,compared with hybrid MC simulations of the microscopic model A ↔A +A [14].Solid (dashed)line is the log −3N (log −6N )power law.
scaling which,interestingly,can be obtained through standard perturbation techniques bad on Gaussian ap-proximations for the fluctuations of the front position [10].The reason for this difference among the various schemes is related to the fact that both the Dickman and BIM algorithms are bad on Gaussian approximations for the density fluctuations which are underestimated for ρ<ρ∗=O (σ∆t ),while only the SSM reproduces exactly the large density Poissonian fluctuations (also obrved in particle models)when ρis small.This does not mean that the Dickman and BIM algorithms do not converge in this ca:specifically,if we take ∆t →0we obrve that the value of the diffusion coefficient approaches that of the hybrid MC simulations for the A ↔A +A (e Fig.3).Thus t
he applicability of the Dickman and BIM algo-rithms is limited in this ca since they fail to reproduce fluctuations at small density and time scales .
In summary,we have prented new strategies for in-tegrating stochastic (partial)differential equations for models of RD subject to internal fluctuations.While all of them prerve the nonnegativity of the solution,algorithms bad on Gaussian approximations introduce a microscopic cutoffbelow which density fluctuations are not correctly accounted for.This is not important when the system properties are dominated by the dynamics of large length and time scales (as in critical behavior),and thus,schemes bad on Euler approximation suffice to integrate numerically equations like (1).However,when the obrved phenomena are nsitive to microscopic fluc-tuations,only algorithms which take into account the exact sampling of density fluctuations at small scales are computationally efficient.Moreover,our results vali-date continuum models like (1)to study the dynamics of fluctuating pulled fronts and corroborate the importance of Poissonian large fluctuations of the density at small scales.We hope that our results will be ud in future
for the analytical understanding of pulled front dynamics [10,14].
Finally,we mention that the methods prented here (the BIM and SSM)can be easily extended to ot
her situations in which the relevant degrees of freedom are non-negative [15,19],like the study of density fluctu-tions in more general RD problems [5],the understand-ing of critical phenomena of systems subject to exter-nal/multiplicative noi (e.g.with ρη(x,t )nois)[2,3],or the nonlinear modelling of the behavior of interest rates in economy [17,20].
We are grateful to E.Brunet,R.Cuerno,C.Doering,H.Schurz,and P.Smereka for comments and discussions.Financial support is acknowledged from the Ministerio de Ciencia y Tecnolog´ıa (Spain).
∗
Electronic address:emoro@math.uc3m.es
[1]M.C.Cross and P.C.Hohenberg,Rev.Mod.Phys.65,851(1993).[2]J.Garc ´ıa-Ojalvo and J.Sancho,Noi in spatially Ex-tended Systems (Springer-Verlag,New York,1999).[3]M.A.Mu˜n oz,in Advances in Condend Matter and Sta-tistical Mechanics ,E.Korutcheva and R.Cuerno Eds.(Nova Sience,New York,2004).[4]A.-L.Barab´a si and H.E.Stanley,Fractal concepts in sur-face growth (Cambridge Univ.Press,Cambridge,1995);J.Krug,Adv.Phys.46,129(1997).
[5]C.W.Gardiner,Handbook of Stochastic Methods ,(Springer,Berlin,1996).
马由页
[6]H.Hinrichn,Adv.Phys.49,815(2000).
[7]P.E.Kloeden,E.Platen,Numerical Solution of Stochas-tic Differential Equations (Springer-Verlag,1992).
[8]M.Doi,J.Phys.A 9,1479(1976);L.Peliti,J.Phys.(Paris)46,1469(1985);D.C.Mattis and M.L.Glasr,Rev.Mod.Phys.70,979(1998).
[9]R.A.Fisher,Ann.Eugenics VII ,355(1936);A.Kol-mogorov,I.Petrosvky,and N.Piscounov,Moscow Univ.Bull.Math.A 1,1(1937);
[10]W.van Saarloos,Phys.Rep.38629(2003);D.Panja,ibid ,393,87(2004).
[11]E.Brunet and B.Derrida,Phys.Rev.E 56,2597(1997);J.Stat.Phys.103,269(2001).
[12]L.Pechenik and H.Levine,Phys.Rev.E 59,3893(1999).[13]C.R Doering,C.Mueller,and P.Smereka,Physica A 325,243(2003).
[14]E.Moro,Phys.Rev.Lett.87238303(2001);Phys.Rev.E 69,060101(2004)
[15]
While writing this paper we became aware of a re-cent preprint,I.Dornic,H.Chat´e ,M. A.Mu˜n oz,cond-mat/0404105,in which a splitting-step like the one in [12]and a similar procedure to compute ρt +∆t using (5)is given.
[16]R.Dickman,Phys.Rev.E 50,4409(1994);C.L´o pez and M.A.Mu˜n oz,ibid 56,4864(1997).
[17]G.N.Milshtein,E.Platen,and H.Schurz,SIAM J.Nu-mer.Anal.35,1010(1998);H.Schurz,Dyn.Syst.Appl.5,323(1996).
[18]
Convergence of scheme (4)requires a cut-offat ρ≃ε,where εis chon small enough (e [17]).
5
[19]E.Moro,in preparation.
[20]J.Cox,E.Ingersoll,and S.A.Ross,Econometrics53,
385(1985).
[21]N.L.Johnson,S.Kotz and N.Balakrishnan,Continuous
univariate distributions(Vol.II)(John Wiley&Sons,
New York,1994);A.F.Siegel,Biometrika66,381(1979).
[22]In our simulations we have ud D=1and a spatial
discretization with∆x=1in a one dimensional lattice with L nodes.
