A Second-Order Stochastic Leap-Frog Algorithm for Multiplicative Noise Brownian Motion

8 Pages·0.22 MB·English

by  Ji Qiang

A Second-Order Stochastic Leap-Frog Algorithm for Multiplicative Noise Brownian
Motion
Ji Qiang1,⋆ and Salman Habib2,†
1LANSCE-1, MS H817, Los Alamos National Laboratory, Los Alamos, NM 87545
2T-8, Theoretical Division, MS B285, Los Alamos National Laboratory, Los Alamos, NM 87545
(February 2, 2008)
0
Astochasticleap-frogalgorithmforthenumericalintegrationofBrownianmotionstochasticdiffer-
0
entialequationswithmultiplicativenoiseisproposedandtested. Thealgorithm hasasecond-order
0
convergenceofmomentsinafinitetimeintervalandrequiresthesamplingofonlyoneuniformlydis-
2
tributedrandomvariablepertimestep. Thenoisemaybewhiteorcolored. Weapplythealgorithm

n to a study of the approach towards equilibrium of an oscillator coupled nonlinearly to a heat bath
a and investigate the effect of the multiplicative noise (arising from the nonlinear coupling) on the
J
relaxation time. This allows us to test the regime of validity of the energy-envelopeapproximation
4
method.

]
PACS Numbers: 05.10.-a, 05.40.-a, 02.60.Cb, 02.50.Ey LAUR99-5263
h
p
-
p
m I. INTRODUCTION the behavior of low-order moments, but the systematics
of these approximations remains to be elucidated. Com-
o
c Stochastic differential equations with multiplicative pared to the Fokker-Planck equation, stochastic differ-
s. noise have not only found many applications in physics ential equations are not difficult to solve, and with the
c but also have interesting mathematical properties. Con- advent of modern supercomputers, it is possible to run
si sequently they have attracted substantial attention over very large numbers of realizations in order to compute
y the years [1–11]. The key point lies in the fundamen- low-ordermomentsaccurately. (Wemaymentionthatin
h taldifference betweenadditive andmultiplicative noises: applicationstofieldtheoriesitisessentiallyimpossibleto
p
Additive noise does not couple directly to the system solvethecorrespondingFokker-Planckequationsincethe
[
variablesanddisappearsfromthenoise-averagedformof probability distribution is now a functional.) However,

2 the dynamical equations. However, in the case of multi- the extractionofthe probabilitydistribution functionit-
v plicativenoise,thesystemvariablesdocoupledirectlyto selfisverydifficultdue to the samplingnoiseinherentin
5 the noise (alternatively, we may say that the noise am- a particle representation of a smooth distribution.
5
plitude depends on the system variables). This fact can Numerical algorithms to solve stochastic differential
0
2 leadtodramaticchangesofsystembehaviorthatcannot equations have been discussed extensively in the litera-
1 occurinthepresenceofadditivenoisealone. Twoclassic ture[14–19]. The simplest,fastest,andstillwidely-used,
9 illustrationsaretheKubooscillator[12]andtheexistence is Euler’s method which yields first-orderconvergenceof
9 of long-time tails in transport theory [13]. In this paper moments for a finite time interval. Depending on the
s/ we willinvestigateanotherexample,that ofanoscillator controloverstatisticalerrorsarisingfromthe necessarily
c nonlinearly coupled to a heat bath, in which the effects finite number of realizations, in the extraction of statis-
si of multiplicative noise can significantlyalter the qualita- tical information it may or may not pay to use a higher
y tive nature, as well as the rate [2], of the equilibration orderalgorithmespeciallyifitiscomputationallyexpen-
h process (relative to that of an oscillator subjected only sive. Because of this fact, it is rare to find high-order
p to additive noise). schemes being put to practical use for the solution of
:
v The dynamical behavior of systems subjected to noise stochastic differential equations, and second-order con-
Xi canbestudiedintwodifferentways: wemayeithersolve vergence is usually considered a good compromise be-
stochasticdifferentialequationsandaverageoverrealiza- tweenefficiency andaccuracy. Apopularalgorithmwith
ar tionstoobtainstatisticalinformation,orwemaydirectly second-order convergence of moments for additive noise
solve the Fokker-Planck equation which describes the butwithonlyfirst-orderconvergenceofmomentsformul-
evolution of the corresponding probability distribution tiplicative noise is Heun’s algorithm(also calledstochas-
function. Bothapproacheshavetheirshareofadvantages tic RK2 by some authors) [14,17,20]. A stochastic leap-
and disadvantages. Fokker-Planck equations are partial frog algorithm which has the same order convergence of
differential equations and their mathematical properties momentsasHeun’s method wassuggestedinRef.[21]to
are still not fully understood. Moreover, they are very study particle motion in a stochastic potential without
expensive to solve numerically even for dynamical sys- damping. Severalotheralgorithmsforparticlemotionin
tems possessing only a very modest number of degrees a quasi-conservative stochastic system were proposed in
of freedom. Truncation schemes or closures (such as cu- Ref.[16]andinthe bookby Allen andTildesley [22]. At
mulanttruncations)havehadsome successinextracting every time step, these methods all require sampling two
1
t
Gaussian random variables which adds to the computa-
tional cost. A modified algorithm suggested in Ref. [19] xi(t)=xi(0)+ dsFi(x1(s),···,xn(s))
Z0
requires only one Gaussian random variable but applies t
only to white Gaussian noise. In the following sections, + dsσij(x1(s), ,xn(s))ξj(s) (4)
···
we present a new stochastic leap-frog algorithm for mul- Z0
tiplicativeGaussianwhite noiseandOrnstein-Uhlenbeck
where x (0) is a given sharp initial condition at t = 0.
i
colored noise which not only has second-order conver-
The infinitesimal update form of this equation may be
gence of moments but also requires the sampling of only
derived by replacing t with an infinitesimal time step h:
one random uniform variable per time step.
The organization of this paper is as follows: General h t′
numerical integration of a system of stochastic differ- x (h)=x (0)+ dt′ F x (0)+ dsF (x(s))
i i i k k
ential equations with Gaussian white noise is discussed Z0 " Z0
in Section II. The stochastic leap-frog algorithms for t′
Brownianmotionwith Gaussianwhite noise and colored + dsσkl(x(s))ξl(s)
Ornstein-Uhlenbeck noise are given in Section III. Nu- Z0 #
mericaltestsofthesealgorithmsusingaone-dimensional h t′
harmonic oscillator are presented in Section IV. A phys- + dt′ σij xk(0)+ dsFk(x(s))
ical application of the algorithm to the multiplicative- Z0 " Z0
noiseBrownianoscillatorisgiveninSectionV.SectionVI t′
containsthefinalconclusionsandandashortdiscussion. + dsσkl(x(s))ξl(s) ξj(t′) (5)
Z0 #
SinceF andσ aresmoothfunctionsofthex ,theymay
II. NUMERICAL INTEGRATION OF i ij i
be expanded about their values at t = 0, in which case
STOCHASTIC DIFFERENTIAL EQUATIONS
we can write the exact solution for x (h) as
i
A general system of continuous-time stochastic differ- x (h)=D (h)+S (h) (6)
i i i
ential equations (Langevin equations) can be written as
where D (h) and S (h) denote the deterministic and
i i
x˙i =Fi(x1, ,xn)+σij(x1, ,xn)ξj(t) (1) stochastic contributions respectively. The deterministic
··· ···
contribution D (h) is
i
where i = 1, ,n and ξ (t) is a Gaussian white noise
j
···
with 1
D (h)=x (0)+hF + h2F F +O(h3) (7)
i i i i,k k
2
ξ (t) =0 (2)
j
ξ (t)hξ (t′)i=δ(t t′) (3) whereFi,k =∂Fi/∂xk,thesummationconventionforthe
j j
h i − repeated indices having being employed. The stochastic
contribution S (h) is
and the symbol represents an average over realiza- i
h···i
tions of the inscribed variable (ensemble average). The
S (h)=σ W (h)+σ σ C (h)+F σ Z (h)
noise is said to be additive when σ is not a function i ij j ij,k kl lj i,k kl l
ij
1
of the x , otherwise it is said to be multiplicative. In
i +σ F (hW (h) Z (h))+ σ σ σ H (h)
ij,k k j j ij,kl km ln mnj
thecaseofmultiplicativenoises,amathematicalsubtlety − 2
arises in interpreting stochastic integrals, the so-called 1 1
+ F σ σ G (h)+ F σ σ K (h)
i,kl ks lt st k ij,kl lm mj
Ito-Stratonovich ambiguity [23]. It should be stressed 2 2
that this is a point of mathematics and not of physics. 1 1
+ F σ σ K (h)+ σ σ σ σ I
Once it is clear how a particular Langevin equation has 2 l ij,kl km mj 6 ij,klm kn lo mp nopj
been derived and what it is supposed to represent, it +O(h5/2) (8)
should either be free of this ambiguity (as in the case of
the example we study later) or it should be clear that The quantities W , C , H , Z , G , K , and I
i ij ijk i ij ij ijkl
there must exist two different stochastic equations, one are random variables which can be written as stochastic
written in the Ito form, the other in Stratonovich, both integrals over the Gaussian white noise ξ(t):
representingthesamephysicalprocessandhenceyielding
identical answers for the variables of interest. (Another h
W (h)= dtξ (t) O(h1/2) (9)
wayto state this is that there should be only one unique i i ∼
Fokker-Planckequation.) Itisimportanttonotethatthe Z0
h
vast majority of numerical update schemes for Langevin C (h)= dtW (t)ξ (t) O(h) (10)
ij i j
equations use the Ito form of the equation. Z0 ∼
The integral representation of the set of equations (1) h
H (h)= dtW (t)W (t)ξ (t) O(h3/2) (11)
is ijk i j k
∼
0
Z
2
h
Z (h)= dtW (t) O(h3/2) (12) III. STOCHASTIC LEAP-FROG ALGORITHM
i i ∼ FOR BROWNIAN MOTION
0
Z
h
Gij(h)= dtWi(t)Wj(t) O(h2) (13) The approach to modeling Brownian motion that we
∼
Z0 considerhereisthatofaparticlecoupledtothe environ-
h
K (h)= tdtW (t)ξ (t) O(h2) (14) ment through its position variable [1]. When this is the
ij i j
∼ case, noise terms enter only in the dynamical equations
Z0
h fortheparticlemomenta. Inthecaseofthreedimensions,
Iijkl(h)= dtWi(t)Wj(t)Wk(t)ξl(t) O(h2) (15) the dynamical equations take the general form:
∼
0
Z
Ito integration has been employed in the derivation of x˙1 =F1(x1,x2,x3,x4,x5,x6)+σ11(x2,x4,x6)ξ1(t)
the above equations. x˙2 =F2(x1)
The nth moment of the x is
i x˙3 =F3(x1,x2,x3,x4,x5,x6)+σ33(x2,x4,x6)ξ3(t)
x (h)n = (D (h)+S (h))n
h i i h i i i x˙4 =F4(x3)
=D (h)n+nD (h)n−1 S (h)
i i h i i x˙5 =F5(x1,x2,x3,x4,x5,x6)+σ55(x2,x4,x6)ξ5(t)
+C2D (h)n−2 (S (h))2 + (16)
n i h i i ··· x˙6 =F6(x5) (27)
where
The convention used here is that the odd indices corre-
i n!
Ci = = (17) spond to momenta, and the even indices to the spatial
n n i!(n i)! coordinate. Inthedynamicalequationsforthemomenta,
(cid:18) (cid:19) −
and the first termonthe righthand side is a systematic drift
termwhichincludestheeffectsduetoexternalforcesand
1
hSi(h)i= 4F,iklσksσlsh2+O(h3) (18) damping. The second term is stochastic in nature and
describes a noise force which, in general, is a function of
1
S (h)S (h) =σilσjlh+ σimσklσjmσplh2 position. The noise ξ(t) is first assumed to be Gaussian
h i j i 2 ,k ,p
and white as defined by Eqns. (2)-(3). The stochastic
1 1
+ σilFjσklh2+ σjlFi σklh2 leap-frog algorithm for the Eqns. (27) is written as
2 ,k 2 ,k
+1σilσjlFkh2+ 1σjlσilFkh2 x¯i(h)=D¯i(h)+S¯i(h) (28)
2 ,k 2 ,k
1 ThedeterministiccontributionD¯ (h)canbeobtainedus-
+ σipσjpσkmσlmh2 i
4 ,kl ingthedeterministicleap-frogalgorithm. Thestochastic
1 contributionS¯(h)canbeobtainedbyapplyingEq.(8)on
+ σjpσipσkmσlmh2+O(h3) (19) i
4 ,kl Eq. (27). The stochastic integration defined by Eqs. (9)
S (h)S (h)S (h) =O(h3) (20) to (15) can be approximated so that the moment rela-
i j k
h i tionshipsdefinedbyEqs.(18)to (22)aresatisfied. After
Si(h)4 =3(σii)4+O(h3) (21) some calculation, the deterministic contribution D¯ (h)
h i i
h(Si(h))5i=O(h3) (22) andthe stochasticcontributionS¯i(h)ofthe aboverecur-
sion formula for one-step integration are found to be
Supposethattheresultsfromanumericalalgorithmwere
written as D¯ (h)=x¯ (0)+hF (x¯∗,x¯∗,x¯∗,x¯∗,x¯∗,x¯∗);
i i i 1 2 3 4 5 6
x¯i(h)=D¯i(h)+S¯i(h) (23) i=1,3,5
{ }
where the x¯ are approximations to the exact solutions D¯ (h)=x¯∗
i i i
xi. Thx¯e(nht)hnm=om(eD¯nt(hof)x+¯iS¯is(h))n +21hFi[xi−1+hFi−1(x¯∗1,x¯∗2,x¯∗3,x¯∗4,x¯∗5,x¯∗6)];
i i i
h i=Dh¯i(h)n+nD¯i(h)ni−1 S¯i(h) {i=2,4,6}
h i 1
+C2D¯ (h)n−2 (S¯(h))2 + (24) S¯(h)=σ √hW (h)+ F σ h3/2W˜ (h)
n i h i i ··· i ii i 2 i,k kk i
ComparingEqns.(16)and(24),weseethatifD (h)and 1
D¯i(h), and Si(h) and S¯i(h) coincide up to h2,iwe will +2σii,jFjh3/2W˜i(h)
have 1
+ F σ σ h2W˜ (h)W˜ (h);
x (h) x¯(h)=O(h3) (25) 4 i,kl kk ll i i
i i
− i=1,3,5; j =2,4,6; k,l=1,3,5
and for a finite time interval { }
1
xi(t)n x¯i(t))n =O(h2) (26) S¯i(h)= √3Fi,jσjjh3/2W˜j(h)
h i−h i
3
1 1
+4Fi,jjσj2jh2W˜j(h)W˜j(h) S¯i(h)= √3σii(x¯2,x¯4,x¯6)kih3/2W˜i(h);
i=2,4,6; j =1,3,5 i=1,3,5
{ } { }
x¯∗i =x¯i(0)+ 12hFi(x¯1,x¯2,x¯3,x¯4,x¯5,x¯6) S¯i(h)=0;
i=2,4,6
i=1,2,3,4,5,6 (29) { }
{ } 1
S¯ =k √hW˜ (h) k2h3/2W˜ (h);
where W˜ (h)is aseries ofrandomnumbers withthe mo- ξi i i − 2 i i
i
ments i=1,3,5 (36)
{ }
where
W˜ (h) = (W˜ (h))3 = (W˜ (h))5 =0 (30)
(W˜hi(hi))2i=h1, i (W˜i(ih))4h =i3 i (31) x¯∗i =x¯i(0)+ 12h[Fi(x¯1,x¯2,x¯3,x¯4,x¯5,x¯6)
h i h i
+σii(x¯2,x¯4,x¯6)ξi];
This can not only be achieved by choosing true Gaus-
i=1,3,5
sian random numbers, but also by using the sequence of
{ }
random numbers following: x¯∗i =x¯i(0)+ 21hFi(x¯1,x¯2,x¯3,x¯4,x¯5,x¯6);
√3, R<1/6 i=2,4,6
W˜i(h)= −√03,, 1/65≤/6R<R5/6 (32) ξi∗ ={ξi(0)exp(−}12kih);
≤
i=1,3,5 (37)
where R is a uniformly distributed random number on { }
the interval (0,1). This trick significantly reduces the
computational cost in generating random numbers.
2.24
Next we consider the case that the noise in Eqs. (27)
2.22
is a colored Ornstein-Uhlenbeck process which obeys
2.2
ξi(t) =0 (33) 2.18
h i
ξ (t)ξ (t′) = ki exp( k t t′ ) (34) <x^2>2.16
i i i
h i 2 − | − | 2.14
2.12
where the correlation factor k is the reciprocal of the
i
2.1
correlation time. In the limit of k , the Ornstein-
i
Uhlenbeck process reducesto Gaussi→an∞white noise. The 2.080 0.1 0.2 0.3 0.4 0.5 0.6 0.7
h
aboveprocesscanbegeneratedbyusingawhiteGaussian
noise from a stochastic differential equation 2.095
2.09
ξ˙(t)= k ξ (t)+k ζ (t) (35)
i i i i i
− 2.085
where ζi(t) is a standard Gaussian white noise. The ini- 2.08
tialvalueξ (0)ischosentobeaGaussianrandomnumber <x^2>
i 2.075
with ξ (0) =0 and ξ (0)2 =k /2.
i i i
Forhthesitochasticphrocessiwithcolorednoise,theleap- 2.07
frog algorithmfor Eqns.(27) is of the same form as that 2.065
for white noise (Cf. Eqn. (29)), but with
2.06
0 0.05 0.1 0.15 0.2 0.25 0.3
h
D¯ (h)=x¯ (0)+hF (x¯∗,x¯∗,x¯∗,x¯∗,x¯∗,x¯∗)
i i i 1 2 3 4 5 6 FIG. 1. Zero damping convergence test. Top: hx2(t)i at
+hσ (x¯∗,x¯∗,x¯∗)ξ∗;
ii 2 4 6 i t = 6 as a function of step size with white Gaussian noise.
i=1,3,5 Bottom: hx2(t)iatt=6asafunctionofstepsizewithcolored
{ }
D¯ (h)=x¯∗ Ornstein-Uhlenbecknoise. Solidlinesrepresentquadraticfits
i i
to thedata points (diamonds).
1
+ hFi[x¯i−1+hFi−1(x¯∗1,x¯2∗,x¯∗3,x¯4∗,x¯∗5,x¯6∗)
2
+hσi−1i−1(x¯∗2,x¯∗4,x¯∗6)ξi∗−1 ; IV. NUMERICAL TESTS
i=2,4,6
{ } (cid:3)
D¯ (h)=ξ (0)exp( k h); Theabovealgorithmsweretestedonaone-dimensional
ξi i − i stochastic harmonic oscillator with a simple form of the
i=1,3,5
{ } multiplicative noise. The equations of motion were
4
p˙ =F1(p,x)+σ(x)ξ(t) gorithm, we computed x2 as a function of t using
h i
x˙ =p (38) 100,000 numerical realizations for a particle starting
from(0.0,1.5)inthe(x,p)phasespace. Theresultsalong
where F1(p,x)= γp η2x and σ(x)= αx. with the analyticalsolutionand a numericalsolutionus-
As a firsttest, w−e co−mputed x2 asa f−unction of time ing Heun’s algorithm are given in Fig. 3. Parameters
h i
step size. To begin, we took the case of zero damping used were h = 0.1, η = 1.0, and α = 0.1. The ad-
constant (γ = 0), where x2 can be determined analyt- vantageinaccuracyofthe stochasticleap-frogalgorithm
ically. The top curve in hFigi. 1 shows x2 at t = 6.0 as overHeun’s algorithmis clearlydisplayed,bothinterms
h i
a function of time step size with white Gaussian noise. of error amplitude and lack of a systematic drift.
Here, the parameters η and α are set to 1.0 and 0.1. WenotethatwhileingeneralHeun’salgorithmisonly
The ensemble averages were taken over 106 independent linearformultiplicativenoiseapplications,forthepartic-
simulations. The analytically determined value of x2 ular problem at hand it turns out to be quadratic. This
h i
at t = 6.0 is 2.095222 (The derivation of the analytical is due to a coincidence: the stochastic term of x does
results is given in the Appendix). The quadratic con- not contain W(h) but does posses a higher order term
vergence of the stochastic leap-frog algorithm is clearly hW(h). However, this higher order term has a larger
seen in the numerical results. We then considered the coefficient compared with our stochastic leap-frog algo-
caseofcoloredOrnstein-Uhlenbecknoiseasafunctionof rithm,andthis accountsforthe largererrorsobservedin
time step size using the same parametersas in the white Fig. 3.
Gaussian noise case and with the correlation parameter
k = 0.16. The result is shown as the bottom curve in
Fig. 1 and the quadratic convergence is again apparent.
12
0.54
0.53 Exact
8
0.52
0.51 <x^2>(t)
<x^2>0.5
4
0.49
Error: Heun
0.48
0.47
0
0.46
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
h Error: Leapfrog
-2
0 100 200 300 400 500
0.454
t
0.452
FIG. 3. Comparing stochastic leap-frog and the Heun al-
0.45 gorithm: hx2(t)i as a function of t. Errors are given relative
0.448
to theexact solution.
<x^2>0.446
0.444
0.442 V. A PHYSICAL APPLICATION: THE
MECHANICAL OSCILLATOR
0.44
0.438
0 0.05 0.1 0.15 0.2 0.25 0.3
h Inthissection,weapplyouralgorithmtostudyingthe
approachto thermal equilibrium of an oscillatorcoupled
FIG. 2. Finite damping (γ = 0.1) convergence test. Top:
hx2(t)iatt=12asafunctionofstepsizewithwhiteGaussian nonlinearly to a heat bath modeled by a set of noninter-
noise. Bottom: hx2(t)i at t = 12 as a function of step size acting harmonic oscillators [1]. The nonlinear coupling
leads to the introduction of multiplicative noise into the
with colored Ornstein-Uhlenbeck noise. Solid lines represent
quadratic fitsto thedata points (diamonds). systemdynamics. LindenbergandSeshadrihavepointed
out that, at weak coupling, multiplicative noise may sig-
We verified that the quadratic convergence is present nificantly enhance the equilibration rate relative to the
for nonzero damping (γ = 0.1). At t = 12.0, and with rateforweaklinearcoupling(additivenoise)[2]. Wewill
all other parameters as above, the convergence of x2 choose the same form of the coordinate couplings as in
as a function of time step is shown by the top and bhoti- Ref. [2], in which case the additive noise equations are
tom curves in Fig. 2 (white Gaussian noise and colored
Ornstein-Uhlenbeck noise, respectively). p˙ = ω02x λ0p+ 2D0ξ0(t)
− −
As a comparison against the conventional Heun’s al- x˙ =p (39)
p
5
and for the system with multiplicative noise only: λ0
1; additive noise (45)
ω0 ≪
p˙= ω02x λ2x2p 2D0xξ2(t) kTλ2
− − − 1; multiplicative noise (46)
x˙ =p (40) ω3 ≪
p 0
where the diffusion coefficients D =λ kT, i=0,2,λ is Asafirstcheck,weperformedsimulationswithω0 =1.0,
i i i
thecouplingconstant,kisBoltzmann’sconstant,T isthe λ0 = λ2 = 0.01, and kT = 4.5, in which case both
the above conditions are satisfied. Moreover, with these
heat bath temperature, and ω0 is the oscillator angular
choices of parameter values, and within the energy en-
frequency without damping. The approach to thermal
velope approximation, the relaxation time predicted for
equilibrium is guaranteed for both sorts of noises by the
multiplicative noise is substantially smaller than for the
fluctuation-dissipation relation
case of additive noise. At the same time we also ran a
ξ (t)ξ (s) =δ δ(t s) (41) simulation at kT = 200 to see how the energy envelope
i j ij
h i − approximation for multiplicative noise breaks down at
written here for the general case when both noises are high temperatures.
simultaneously present. While in all cases, it is clear
that the final distribution is identical and has to be the
1.2
thermal distribution, the precise nature of the approach
to equilibrium can certainly be different. We wish to
I
explore this issue in more detail. An important point to 1
keepinmindisthatinthisparticularsystemofequations
II
thereisnonoise-induceddriftintheFokker-Planckequa- 0.8
tionobtainedfromtheStratonovichformoftheLangevin III
<E(t)>/kT
equation, i.e., there is no Ito-Stratonovichambiguity.
0.6
It is a simple matter to solve the Langevin equations
given above applying the algorithm from Eqs. (29). As
0.4
our primary diagnostic,we computed the noise-averaged
energy E(t) of the oscillator as a function of time t,
h i
where 0.2
1 1
E(t)= 2p2+ 2ω02x2. (42) 00 50 100 150 200 250 300 350 400
t
Intheweakcouplinglimitandemployingorbit-averaging
FIG. 4. Temporal evolution of the scaled average energy
(validpresumablywhenthedynamicaltimescaleismuch
hE(t)i/kT with additive noise and multiplicative noise. The
smaller than the relaxation time scale), one finds [2]
dashed lines I and II are the predictions from Eqn. (44) for
E(t) =kT (kT E0)e−λ0t (43) kT = 200 and kT = 4.5 respectively. The dashed line III is
h i − − the theoretical prediction for additive noise with kT = 4.5.
in the case of additive noise (a result which can also be Aspredicted,therelaxationproceedsmuchfasterwithmulti-
directlyobtainedasalimiting casefromthe knownform plicative noise: The solid lines are numerical results for mul-
tiplicativenoiseatkT =200andkT =4.5. Itisclearthatat
of the exact solutiongiven, e.g., in Ref. [24]). The corre-
higher temperatures, the theory grossly underestimates the
sponding formofthe approximatesolutioninthe caseof
relaxation time.
multiplicative noise is
E0kT In Fig. 4, we display the time evolution of the aver-
E(t) = . (44)
h i E0+(kT E0)exp( λ2kTt/ω02) age energy (scaled by kT for convenience) with additive
− − and multiplicative noise both from the simulations and
While in the case of additive noise, the exponential na- the approximate analytical calculations. In the case of
tureoftherelaxationisalreadyclearfromtheformofthe weak coupling to the environment (small λ0, λ2), the
exact solution (cf. Ref. [24]), the situation in the case of rate atwhich the averageenergy approachesequilibrium
multiplicativenoiseisnotobviouslyapparentasnoexact issignificantlygreaterforthe caseofmultiplicativenoise
solutionisknowntoexist. Thepredictionofarelaxation relative to the case of additive noise more or less as ex-
processcontrolledbyasingleexponentialasfoundin(44) pected. In addition, the analytic approximation result-
is a consequence of the assumption x2(t) kT/ω2 at ing from the application of the energy-envelope method
0
h i ≃
“late”times,thisimplyingaconstantdampingcoefficient (44) is seen to be in reasonable agreement with the nu-
in the Langevin equation (40). merical simulations for kT = 4.5. The slightly higher
The timescale separations necessary for the energy- equlibrationratefromtheanalyticalcalculationisdueto
envelope methodto be applicableare encodedin the fol- the truncationinthe energyenvelope equationusing the
lowing inequalities [2]: E2(t) 2 E(t) 2 relationwhichyieldsanupper bound
h i≈ h i
6
on the rate of equilibration of the average energy [2]. oscillator-heat-bathsysteminordertoinvestigatetheef-
Note that in the case of high temperature (kT = 200) fectofmultiplicativenoiseonthenatureoftherelaxation
the relaxaton time computed from the energy envelope process.
method is much smaller than the numerical result, con-
sistent with the violation of the condition (46).
While the results shown in Fig. 4 do show that the VII. ACKNOWLEDGMENTS
energy envelope approximation is qualitatively correct
within its putative domain of validity, it is clear that We acknowledge helpful discussions with Grant Lythe
the actual relaxation process is not of the precise form and Robert Ryne. Partial support for this work came
(44). In Fig. 5 we illustrate this point by plotting from the DOE Grand Challenge in Computational Ac-
celeratorPhysics. Numericalsimulationswereperformed
E0(kT −hE(t)i) =exp( λ2kTt/ω02) (47) on the SGI Origin2000 systems at the Advanced Com-
E(t) (kT E0) − puting Laboratory (ACL) at Los Alamos National Lab-
h i −
oratory, and on the Cray T3E at the National En-
[equivalent to (44)] against time on a log scale: the re-
ergy Research Scientific Computing Center (NERSC) at
laxation is clearly nonexponential. The reason for the
Lawrence Berkeley National Laboratory.
failure of the approximationis that despite the fact that
equipartition of energy does take place on a relatively
shorttimescale,itisnottruethat x2(t) canbetreated
h i
as a constant even at relatively late times.
1 ⋆
Electronic address: [email protected]
† Electronic address: [email protected]
Additive Noise [1] R. Zwanzig, J. Stat. Phys. 9, 215 (1973).
[2] K. Lindenberg and V. Seshadri, Physica 109 A, 483
0.1 (1981).
[3] A.Careta and F.Sagues, Phys.Rev.A44, 2284 (1991).
[4] S.HabibandH.Kandrup,Phys.Rev.D46,5303(1992).
[5] S. Habib, Ann.N.Y. Acad. Sci. 706, 111 (1993).
Multiplicative Noise
[6] G.Efremov,L.Mourokh,andA.Smirnov,Phys.Lett.A
0.01 175, 89 (1993).
[7] A. Becker and L. Kramer, Phys. Rev. Lett. 73, 955
(1994).
[8] H. Leung, Physica A 221, 340 (1995).
[9] J. Bao, Y. Zhuo, and X. Wu, Phys. Lett. A 217, 241
0.001
0 50 100 150 200 250 300 (1996).
t [10] S.Mangioni, R.Deza, H.Wio, andR.Toral, Phys.Rev.
Lett. 79, 2389 (1997).
FIG.5. TheLHSof(47)asafunctionoftime(straightline)
[11] W. Genovese, M. Munoz, and J. Sancho, Phys. Rev. E
compared with numerical results for kT = 4.5. Also shown
57, R2495 (1998).
isanumericalresult for thecaseof additivenoise whichisin
[12] R. Kubo,J. Math. Phys 4, 174 (1963).
excellentagreementwiththepredictedexponentialrelaxation
[13] R.W. Zwanzig, in Statistical mechanics; new concepts,
with therelaxation timescale =1/λ0.
new problems, new applications edited by S.A. Rice,
K.F.Freed,andJ.C.Light(UniversityofChicago Press,
Chicago, 1972).
VI. CONCLUSIONS
[14] A.Greiner, W.Strittmatter,andJ.Honerkamp,J.Stat.
Phys. 51, 94 (1988).
We have presented a stochastic leap-frog algorithm [15] R. Mannella, and V. Palleschi, Phys. Rev. A 40, 3381
for single particle Brownian motion with multiplicative (1989).
noise. This method has the advantages of retaining the [16] R. Mannella, in Noise in Nonlinear Dynamical Systems,
symplecticpropertyinthedeterministiclimit,easeofim- vol.3,F.MossandP.V.E.McClintock,Eds.(Cambridge
plementation, and second-orderconvergenceof moments University Press, Cambridge, 1989).
formultiplicativenoise. Samplingauniformdistribution [17] R.L. Honeycutt,Phys. Rev.A 45, 600 (1992).
instead of a Gaussian distribution helps to significantly [18] P.E. Kloeden and E. Platen, Numerical Solution of
reduce the computational cost. A comparison with the Stochastic Differential Equations (Springer, New York,
conventionalHeun’salgorithmhighlightsthegaininacu- 1992).
racyduetothenewmethod. Finally,wehaveappliedthe [19] R.Mannella,inSupercomputation inNonlinear andDis-
stochastic leap-frog algorithm to a nonlinearly coupled
7
ordered Systems, L. Vazuez, F. Tirado, and I. Marun, 2α2 4η2x x3 =0 (A6)
− −
Eds., p. 101 (World Scientific, 1996).
[20] S. Habib, H.E. Kandrup, and M.E. Mahon, Phys. Rev. which gives
E 53, 5473 (1996).
[21] M. Seesselberg, H.P. Breuer, H. Mais, F. Petruccione, r1 = 64/27η6+α4+α2 1/3
and J. Honerkamp,Z. Phys. C 62, 63 (1994).
[22] M.P. Allen, and D.J. Tildesley, Computer Simulation of (cid:16)p (cid:17) 1/3
64/27η6+α4 α2
Liquids(Clarendon Press, Oxford,1987). − −
[23] C.W. Gardiner, Handbook of Stochastic Methods for 1 (cid:16)p (cid:17) 1/3
Physics, Chemistry, and the Natural Sciences (Springer, r2 = (1+√3i) 64/27η6+α4 α2
2 −
NewYork, 1983).
1 (cid:16)p (cid:17) 1/3
[24] H.Risken,The Fokker-Planck Equation: Methods of So- (1 √3i) 64/27η6+α4+α2
lution and Applications (Springer,New York,1989). −2 −
r3 =r2∗ (cid:16)p (cid:17) (A7)
where the superscript represents complex conjugation.
APPENDIX A: Thepositiverealrootr∗1 impliesthat x2(t) willhavean
h i
exponential growth in time.
The analytic solution of Eqns. (38) for x2(t) (with
h i
whiteGaussiannoise)asafunctionoftimeinthespecial
case of zero damping, i.e. γ = 0, can be obtained by
solvingtheequivalentFokker-Planckequation[24]forthe
probability density f(x,p,t):
∂
f(x,p,t)=
∂t
p ∂ ∂F1(p,x) + 1σ2(x) ∂2 f(x,p,t) (A1)
− ∂x − ∂p 2 ∂p2
(cid:20) (cid:21)
The expectation value of any function M(x,p;t) can be
written as
+∞
M(x,p) = dxdpM(x,p)f(x,p,t) (A2)
h i
Z−∞
Equations(A1)and(A2)canbeusedtoyieldaBBGKY-
like heirarchy for the evolution of phase space moments.
Since the systemweareconsideringis linear,this heirar-
chy truncates exactly and yields a group of coupled lin-
ear ordinarydifferential equations for the moments x2 ,
xp ,and p2 . Theseequationscanbewrittenasashinglie
hthirid-ordehr tiime evolution equation for x2 :
h i
d3 x2 d x2
h i = 4η2 h i +2α2 x2 (A3)
dt3 − dt h i
subject to the initial conditions
x2(0) =x2(0)
h i
x˙2(0) =2x(0)p(0)
h i
x¨2(0) =2p2(0) 2η2x2(0) (A4)
h i −
This equation has an analytical solution written as
x2(t) =c1exp(r1t)+c2exp(r2t)+c3exp(r3t) (A5)
h i
where c1, c2, and c3 are constants depending on initial
conditions,andr1,r2andr3aretherootsofathirdorder
alegbraic equation
8