SLIDE 1 Computing ergodic limits for SDEs
M.V. Tretyakov School of Mathematical Sciences, University of Nottingham, UK
Talk at the workshop “Stochastic numerical algorithms, multiscale modelling and high-dimensional data analytics”, ICERM, 21st July 2016
SLIDE 2 Plan of the talk
ϕerg =
t→∞ Eϕ(X(t)) = lim t→∞
1 t
t
SLIDE 3 Plan of the talk
ϕerg =
t→∞ Eϕ(X(t)) = lim t→∞
1 t
t
Introduction Examples of Langevin-type equations and stochastic gradient systems Geometric integrators for Langevin equations and the gradient system ‘Non-Markovian’ scheme for stochastic gradient systems [Davidchack, Ouldridge&T. J Chem Phys 2015] [Leimkuhler, Matthews,T 2014]
SLIDE 4 Introduction
Hamiltonian H(r, p) canonical ensemble (NVT) ρ(r, p) ∝ exp(−βH(r, p)), where β = 1/(kBT) > 0 is an inverse temperature. ϕerg =
t→∞ Eϕ(R(t)) = lim t→∞
1 t
t
SLIDE 5 Rigid Body Dynamics
Following [Miller III et al J. Chem. Phys., 2002] H(r, p, q, π) = pTp 2m +
n
3
Vl(qj, πj) + U(r, q), (1) where r = (r 1
T
, . . . , r n
T
)T ∈ R3n, and p = (p1
T
, . . . , pn
T
)T ∈ R3n are the center-of-mass coordinates and momenta; q = (q1
T
, . . . , qn
T
)T ∈ R4n, qj = (qj
0, qj 1, qj 2, qj 3)T ∈ R4, |qj| = 1, are the rotational coordinates in the
quaternion representation, π = (π1
T
, . . . , πn
T
)T ∈ R4n, are the angular momenta;
3
Vl(q, π) = 1 8
3
1 Il
2 = 1 8πTS(q)DST(q)π, q, π ∈ R4, l = 1, 2, 3, Il – the principal moments of inertia and the constant 4-by-4 matrices Sl : S1q = (−q1, q0, q3, −q2)T, S2q = (−q2, −q3, q0, q1)T, S3q = (−q3, q2, −q1, q0)T, S0 = diag(1, 1, 1, 1)T. S(q) = [S0q, S1q, S2q, S3q], D = diag(0, 1/I1, 1/I2, 1/I3).
SLIDE 6 Langevin thermostat for Rigid Body Dynamics
dRj = Pj m dt, Rj(0) = r j, (2) dPj = f j(R, Q)dt − γPjdt +
β dw j(t), Pj(0) = pj, dQj = 1 4S(Qj)DST(Qj)Πjdt, Qj(0) = qj, |qj| = 1, (3) dΠj = 1 4
3
1 Il
SlΠjdt + F j(R, Q)dt − ΓJ(Qj)Πjdt +
β
3
SlQjdW j
l (t),
Πj(0) = πj, qj Tπj = 0, j = 1, . . . , n, where f j(r, q) = −∇r jU(r, q) ∈ R3, F j(r, q) = − ˜ ∇qjU(r, q) ∈ TqjS3, (wT, WT)T = (w 1 T, . . . , w n T, W 1 T, . . . , W n T)T is a (3n + 3n)-dimensional standard Wiener process; γ ≥ 0 and Γ ≥ 0 are the friction coefficients for the translational and rotational motions, and J(q) = MS(q)DST(q)/4, M = 4/
3
1/Il. Davidchack, Ouldridge&T. J Chem Phys 2015
SLIDE 7
Langevin-type equations
The Ito interpretation of the SDEs (2)–(3) coincides with its Stratonovich interpretation. The solution of (2)–(3) preserves the quaternion length |Qj(t)| = 1, j = 1, . . . , n , for all t ≥ 0. (4) The solution of (2)–(3) automatically preserves the constraint: Qj T(t)Πj(t) = 0 , j = 1, . . . , n , for t ≥ 0 (5) Assume that the solution X(t) = (RT(t), PT(t), QT(t), ΠT(t))T of (2)–(3) is an ergodic process on D = {x = (rT, pT, qT, πT)T ∈ R14n : |qj| = 1, qj Tπj = 0, j = 1, . . . , n}. Then it can be shown that the invariant measure of X(t) is Gibbsian with the density ρ(r, p, q, π) on D: ρ(r, p, q, π) ∝ exp(−βH(r, p, q, π)) (6) Davidchack, Ouldridge&T. J Chem Phys 2015
SLIDE 8 Langevin equations and quasi-symplectic integrators
dRj = Pj m dt, Rj(0) = r j, (9) dPj = f j(R, Q)dt − γPjdt +
β dw j(t), Pj(0) = pj, dQj = 1 4S(Qj)DST(Qj)Πjdt, Qj(0) = qj, |qj| = 1, (10) dΠj = 1 4
3
1 Il
SlΠjdt + F j(R, Q)dt − ΓJ(Qj)Πjdt +
β
3
SlQjdW j
l (t),
Πj(0) = πj, qj Tπj = 0, j = 1, . . . , n,
SLIDE 9 Langevin equations and quasi-symplectic integrators
dRj = Pj m dt, Rj(0) = r j, (9) dPj = f j(R, Q)dt − γPjdt +
β dw j(t), Pj(0) = pj, dQj = 1 4S(Qj)DST(Qj)Πjdt, Qj(0) = qj, |qj| = 1, (10) dΠj = 1 4
3
1 Il
SlΠjdt + F j(R, Q)dt − ΓJ(Qj)Πjdt +
β
3
SlQjdW j
l (t),
Πj(0) = πj, qj Tπj = 0, j = 1, . . . , n, Let D0 ∈ Rd, d = 14n, be a domain with finite volume. The transformation x = (r, p, q, π) → X(t) = X(t; x) = (R(t; x), P(t; x), Q(t; x), Π(t; x)) maps D0 into the domain Dt.
SLIDE 10 Langevin equations and quasi-symplectic integrators
Vt =
dX 1 . . . dX d (7) =
D(x1, . . . , xd)
The Jacobian J is equal to J = D(X 1, . . . , X d) D(x1, . . . , xd) = exp (−n(3γ + Γ) · t) . (8) [Milstein&T. IMA J. Numer. Anal. 2003 (also Springer 2004)], [Davidchack, Ouldridge&T. J Chem Phys 2015]
SLIDE 11 The stochastic gradient system
It is easy to verify that
exp(−βH(r, p, q, π))dpdπ (9) ∝ exp(−βU(r, q)) =: ˜ ρ(r, q), where (rT, qT)T ∈ D′ = {(rT, qT)T ∈ R7n : |qj| = 1} and the domain of conjugate momenta Dmom = {(pT, π
T)T ∈ R7n : qTπ = 0}.
SLIDE 12 The stochastic gradient system
It is easy to verify that
exp(−βH(r, p, q, π))dpdπ (9) ∝ exp(−βU(r, q)) =: ˜ ρ(r, q), where (rT, qT)T ∈ D′ = {(rT, qT)T ∈ R7n : |qj| = 1} and the domain of conjugate momenta Dmom = {(pT, π
T)T ∈ R7n : qTπ = 0}.
We introduce the gradient system in the form of Stratonovich SDEs: dR = υ mf(R, Q)dt +
mβ dw(t), R(0) = r, (10) dQj = Υ M F j(R, Q)dt +
Mβ
3
SlQj ⋆ dW j
l (t),
(11) Qj(0) = qj, |qj| = 1, j = 1, . . . , n, where “⋆” indicates the Stratonovich form of the SDEs, parameters υ > 0 and Υ > 0 control the speed of evolution of the gradient system (10)–(11), f = (f 1 T, . . . , f n T)T and the rest of the notation is as in (2)–(3). [Davidchack, Ouldridge&T. J Chem Phys 2015]
SLIDE 13
The gradient thermostat for rigid body dynamics
This new gradient thermostat possesses the following properties. As in the case of (2)–(3), the solution of (10)–(11) preserves the quaternion length (4). Assume that the solution X(t) = (RT(t), QT(t))T ∈ D′ of (10)–(11) is an ergodic process. Then, by the usual means of the stationary Fokker-Planck equation, one can show that its invariant measure is Gibbsian with the density ˜ ρ(r, q) from (9).
SLIDE 14 Langevin integrators
Davidchack, Ouldridge&T. J Chem Phys 2015 For simplicity we use a uniform time discretization of a time interval [0, T] with the step h = T/N. Goal: to construct integrators quasi-symplectic preserve | ¯ Qj(tk)| = 1, j = 1, . . . , n , for all t ≥ 0 automatically preserve ¯ Qj T(tk)¯ Πj(tk) = 0 , j = 1, . . . , n , for t ≥ 0 automatically
To this end: stochastic numerics+splitting techniques [see e.g. Milstein&T, Springer 2004] the deterministic symplectic integrator from [Miller III et al J. Chem. Phys., 2002]
SLIDE 15 ‘Langevin A’ integrator
The first integrator is based on splitting the Langevin system in dRj = Pj m dt, Rj(0) = r j, (12) dPj = f j(R, Q)dt +
β dw j(t), dQj = 1 4S(Qj)DST(Qj)Πjdt, (13) dΠj = 1 4
3
1 Il
SlΠjdt + F j(R, Q)dt +
β
3
SlQjdW j
l (t), j = 1, . . . , n,
SLIDE 16 ‘Langevin A’ integrator
The first integrator is based on splitting the Langevin system in dRj = Pj m dt, Rj(0) = r j, (12) dPj = f j(R, Q)dt +
β dw j(t), dQj = 1 4S(Qj)DST(Qj)Πjdt, (13) dΠj = 1 4
3
1 Il
SlΠjdt + F j(R, Q)dt +
β
3
SlQjdW j
l (t), j = 1, . . . , n,
and the deterministic system of linear differential equations ˙ p = −γp, ˙ πj = −ΓJ(qj)πj, j = 1, . . . , n . (14)
SLIDE 17 ‘Langevin B’ integrator
dPI = −γPI dt +
β dw(t), dΠj
I = −ΓJ(q)Πj Idt +
β
3
SlqdW j
l (t);
(15) dRII = PII m dt, dPII = f(RII, QII)dt, dQj
II = 1
4S(Qj
II)DST(Qj II)Πj IIdt , (16)
dΠj
II
= F j(RII, QII)dt + 1 4
3
1 Il
II)TSlQj II
IIdt , j = 1, . . . , n.
SLIDE 18 ‘Langevin B’ integrator
dPI = −γPI dt +
β dw(t), dΠj
I = −ΓJ(q)Πj Idt +
β
3
SlqdW j
l (t);
(15) dRII = PII m dt, dPII = f(RII, QII)dt, dQj
II = 1
4S(Qj
II)DST(Qj II)Πj IIdt , (16)
dΠj
II
= F j(RII, QII)dt + 1 4
3
1 Il
II)TSlQj II
IIdt , j = 1, . . . , n.
The SDEs (15) have the exact solution: PI(t) = PI(0) exp(−γt) +
β t exp(−γ(t − s))dw(s), (17) Πj
I(t)
= exp(−ΓJ(q)t)Πj
I(0) +
β
3
t exp(−ΓJ(q)(t − s))dW j
l (s).
To construct the method: half a step of (17) + one step of a symplectic method for (16) + half a step of (17).
SLIDE 19
‘Langevin C’ integrator
Based on the same spliting (15) and (16) as Langevin B, i.e., the determinisitic Hamiltonian system + OU. To construct the method: half a step of a symplectic method for (16) + step of (17) + half a step of a symplectic method for (16).
SLIDE 20
‘Langevin C’ integrator
Based on the same spliting (15) and (16) as Langevin B, i.e., the determinisitic Hamiltonian system + OU. To construct the method: half a step of a symplectic method for (16) + step of (17) + half a step of a symplectic method for (16). Various splittings are compared for a translational Langevin thermostat in [Leimkuhler&Matthews 2013]
SLIDE 21
Numerical experiments
Davidchack, Handel&T. J Chem Phys 2009 and Davidchack, Ouldridge&T. J Chem Phys 2015 Two objectives for the experiments: the dependence of Langevin dynamics properties on the choice of parameters γ and Γ errors of the numerical schemes. TIP4P rigid model of water (Jorgensen et. al J. Chem. Phys. 1983) Ah : = EA( ¯ X) = EA(X) + EAhp + O(hp+1) = A0 + EAhp + O(hp+1) p = 2 for Langevin integrators and p = 1 for the gradient thermostat integrator Talay&Tubaro Stoch.Anal.Appl. 1990
SLIDE 22
Accuracy of integrators
Langevin A (left), Langevin B (centre), and Langevin C (right) with γ = 5 ps−1 and Γ = 10 ps−1. Error bars denote estimated 95% confidence intervals.
SLIDE 23
Geometric integrator for the gradient thermostat
The main idea is to rewrite the components Qj of the solution to (10)–(11) in the form Qj(t) = exp(Y j(t))Qj(0) and then solve numerically the SDEs for the 4 × 4-matrices Y j(t)
SLIDE 24 Geometric integrator for the gradient thermostat
The main idea is to rewrite the components Qj of the solution to (10)–(11) in the form Qj(t) = exp(Y j(t))Qj(0) and then solve numerically the SDEs for the 4 × 4-matrices Y j(t) – Lie group methods [Hairer, Lubich, Wanner; Springer, 2002] To this end, we introduce the 4 × 4 skew-symmetric matrices: Fj(r, q) = F j(r, q)qj T − qj(F j(r, q))T, j = 1, . . . , n. Note that Fj(r, q)qj = F j(r, q) under |qj| = 1 and the equations (11) can be written as dQj = Υ M Fj(R, Q)Qjdt+
Mβ
3
SlQj⋆dW j
l (t), Qj(0) = qj, |qj| = 1.
(18) One can show that Y j(t + h) = h Υ M Fj(R(t), Q(t)) +
Mβ
3
l (t + h) − W j l (t)
+ terms of higher order.
SLIDE 25 Geometric integrator for the gradient thermostat
R0 = r, Q0 = q, |qj| = 1, j = 1, . . . , n, (19) Rk+1 = Rk + h υ mf(Rk, Qk) + √ h
mβ ξk, Y j
k = h Υ
M Fj(Rk, Qk) + √ h
Mβ
3
ηj,l
k Sl,
Qj
k+1 = exp(Y j k)Qj k,
j = 1, . . . , n, where ξk = (ξ1,k, . . . , ξ3n,k)T and ξi,k, i = 1, . . . , 3n, ηj,l
k , l = 1, 2, 3,
j = 1, . . . , n, are i.i.d. random variables with the same law P(θ = ±1) = 1/2. (20)
- Proposition. The numerical scheme (19)–(20) for (10)–(11) preserves
the length of quaternions, i.e., |Qj
k| = 1,
j = 1, . . . , n , for all k, and it is
Davidchack, Ouldridge&T. J Chem Phys 2015
SLIDE 26 Stochastic gradient system
dX = a(X)dt + σdw, X(0) = X0, (21) a(x) := −∇U(x), x ∈ Rd, (22) σ =
- 2/β, d = 3n, and w(t) is a standard d-dimensional Wiener
process.
SLIDE 27 Stochastic gradient system
dX = a(X)dt + σdw, X(0) = X0, (21) a(x) := −∇U(x), x ∈ Rd, (22) σ =
- 2/β, d = 3n, and w(t) is a standard d-dimensional Wiener
process. We use the following notation for the solution of (21): X(t) = Xt0,x(t) when X(t0) = x, t ≥ t0, and also we will write Xx(t) when t0 = 0.
SLIDE 28 Stochastic gradient system
dX = a(X)dt + σdw, X(0) = X0, (21) a(x) := −∇U(x), x ∈ Rd, (22) σ =
- 2/β, d = 3n, and w(t) is a standard d-dimensional Wiener
process. We use the following notation for the solution of (21): X(t) = Xt0,x(t) when X(t0) = x, t ≥ t0, and also we will write Xx(t) when t0 = 0. The process X(t) is exponentially ergodic [Hasminskii 1980] if for any x ∈ Rd and any function ϕ with a polynomial growth there are C(x) > 0 and λ > 0 such that |Eϕ(Xx(t)) − ϕerg| ≤ C(x)e−λt, t ≥ 0,
SLIDE 29 Stochastic gradient system
dX = a(X)dt + σdw, X(0) = X0, (21) a(x) := −∇U(x), x ∈ Rd, (22) σ =
- 2/β, d = 3n, and w(t) is a standard d-dimensional Wiener
process. We use the following notation for the solution of (21): X(t) = Xt0,x(t) when X(t0) = x, t ≥ t0, and also we will write Xx(t) when t0 = 0. The process X(t) is exponentially ergodic [Hasminskii 1980] if for any x ∈ Rd and any function ϕ with a polynomial growth there are C(x) > 0 and λ > 0 such that |Eϕ(Xx(t)) − ϕerg| ≤ C(x)e−λt, t ≥ 0, The solution X(t) of (21) is exponentially ergodic and ρ(x) ∝ exp
σ2 U(x)
- if there exist c0 ∈ R and c1 > 0 such that
(x, a(x)) ≤ c0 − c1|x|2. (23) [Hasminskii (1980), Mattingly, Stuart, Higham (2002)]
SLIDE 30 Stochastic gradient system
The Euler scheme: Xk+1 = Xk + ha(Xk) + σ √ hξk+1, (24) where ξk = (ξ1
k, . . . , ξd k)⊤ and ξi k, i = 1, . . . , d, k = 1, . . . , are i.i.d.
random variables with the law N(0, 1).
SLIDE 31 Stochastic gradient system
The Euler scheme: Xk+1 = Xk + ha(Xk) + σ √ hξk+1, (24) where ξk = (ξ1
k, . . . , ξd k)⊤ and ξi k, i = 1, . . . , d, k = 1, . . . , are i.i.d.
random variables with the law N(0, 1). Heun’s scheme: ˆ Xk+1 = Xk + ha(Xk) + σ √ hξk+1, Xk+1 = Xk + h 2
Xk+1) + a(Xk)
√ hξk+1. (25)
SLIDE 32 Stochastic gradient system
The Euler scheme: Xk+1 = Xk + ha(Xk) + σ √ hξk+1, (24) where ξk = (ξ1
k, . . . , ξd k)⊤ and ξi k, i = 1, . . . , d, k = 1, . . . , are i.i.d.
random variables with the law N(0, 1). Heun’s scheme: ˆ Xk+1 = Xk + ha(Xk) + σ √ hξk+1, Xk+1 = Xk + h 2
Xk+1) + a(Xk)
√ hξk+1. (25) Weak convergence: |Eϕ(X(T)) − Eϕ(XN)| ≤ Khp (26) N = T/h; Euler – p = 1, Heun – p = 2 [e.g. Milstein, T., Springer 2004]
SLIDE 33 Stochastic gradient system
The Euler scheme: Xk+1 = Xk + ha(Xk) + σ √ hξk+1, (24) where ξk = (ξ1
k, . . . , ξd k)⊤ and ξi k, i = 1, . . . , d, k = 1, . . . , are i.i.d.
random variables with the law N(0, 1). Heun’s scheme: ˆ Xk+1 = Xk + ha(Xk) + σ √ hξk+1, Xk+1 = Xk + h 2
Xk+1) + a(Xk)
√ hξk+1. (25) Weak convergence: |Eϕ(X(T)) − Eϕ(XN)| ≤ Khp (26) N = T/h; Euler – p = 1, Heun – p = 2 [e.g. Milstein, T., Springer 2004] Ergodic limits [Talay (1990); Talay, Tubaro (1990); Mattingly, Stuart, Higham (2002); Milstein, T (2007); Mattingly, Stuart, T. (2010)]: |ϕerg − Eϕ(XN)| ≤ Khp + Ce−λT (27)
SLIDE 34 Non-Markovian scheme
Xk+1 = Xk + ha(Xk) + σ √ h 2 (ξk + ξk+1), (28) where ξk = (ξ1
k, . . . , ξi k)⊤ defined on (Ω, P, F) and ξi k, i = 1, . . . , d,
k = 1, . . . , are i.i.d. random variables with the law N(0, 1) [Leimkuhler, Matthews (2013)]
SLIDE 35 Example
Let a(x) = −αx with α > 0, then X(t) from (21) is the Ornstein-Uhlenbeck process, which is Gaussian with EXx(t) = xe−αt, Cov(Xx(s), Xx(t)) = σ2 2α(e−α(t−s) − e−α(t+s)) for s ≤ t and Var(Xx(T)) = σ2
2α(1 − e−2αT).
SLIDE 36 Example
Let a(x) = −αx with α > 0, then X(t) from (21) is the Ornstein-Uhlenbeck process, which is Gaussian with EXx(t) = xe−αt, Cov(Xx(s), Xx(t)) = σ2 2α(e−α(t−s) − e−α(t+s)) for s ≤ t and Var(Xx(T)) = σ2
2α(1 − e−2αT).
For the Euler scheme (24): EXN = x0(1 − αh)N = x0e−αT(1 + O(h)), Var(XN) = σ2 2α 1 − (1 − αh)2N 1 + αh = σ2 2α(1 − e−2αT) − σ2 2 h + e−2αTO(h) + O(h2), αh < 1, where |O(hp)| ≤ Kh with K > 0 independent of T.
SLIDE 37 Example
Let a(x) = −αx with α > 0, then X(t) from (21) is the Ornstein-Uhlenbeck process, which is Gaussian with EXx(t) = xe−αt, Cov(Xx(s), Xx(t)) = σ2 2α(e−α(t−s) − e−α(t+s)) for s ≤ t and Var(Xx(T)) = σ2
2α(1 − e−2αT).
For the Euler scheme (24): EXN = x0(1 − αh)N = x0e−αT(1 + O(h)), Var(XN) = σ2 2α 1 − (1 − αh)2N 1 + αh = σ2 2α(1 − e−2αT) − σ2 2 h + e−2αTO(h) + O(h2), αh < 1, where |O(hp)| ≤ Kh with K > 0 independent of T. For the scheme (28): EXN = x0(1 − αh)N = x0e−αT(1 + O(h)), Var(XN) = σ2 2α
1 − αh
2α(1 − e−2αT) + e−2αTO(h).
SLIDE 38
Assumptions
Assumption 1 There exist c0 ∈ R and c1 > 0 such that (x, a(x)) ≤ c0 − c1|x|2. Assumption 2 The potential U(x) ∈ C 7(Rd), its first-order derivatives grow not faster than a linear function at infinity and higher derivatives are bounded. The function ϕ(x) ∈ C 6(Rd) and it and its derivatives grow not faster than a polynomial function at infinity. The most restrictive condition in Assumption 2 is the requirement for a(x) = −∇U to be globally Lipschitz: |a(x)|2 ≤ K(1 + |x|2), (29) where K > 0 is independent of x ∈ Rd, which can be relaxed via [Milstein, T. (2005) and (2007): the concept of rejecting exploding trajectories]
SLIDE 39 Introduce the operator L L := ∂ ∂t + L, where L L :=
d
ai(x) ∂ ∂xi + σ2 2
d
∂2 (∂xi)2 . (30) The function u(t, x) = Eϕ(Xt,x(T)) (31) satisfies the Cauchy problem for the backward Kolmogorov equation Lu = 0, (32) u(T, x) = ϕ(x).
SLIDE 40 Main theorem
Theorem (1) Let Assumptions 1-2 hold. Then the scheme (28) is first order weakly convergent and for all sufficiently small h > 0 its error has the form Eϕ(Xx(T)) − Eϕ(XN) = C0(T, x)h + C(T, x)h2, (33) C0(T, x) = E T B0(t, Xx(t))dt, (34) B0(t, x) = 1 2
d
aj(x)∂ai(x) ∂xj ∂u(t, x) ∂xi +σ2 2
d
∂ai(x) ∂xj ∂2u(t, x) ∂xi∂xj + σ2 2
d
∂2ai(x) (∂xj)2 ∂u(t, x) ∂xi , |C(T, x)| ≤ K(1 + |x|κe−λT), for some K > 0, κ ∈ N and λ > 0 independent of h and T.
SLIDE 41
Corollary
Theorem 1: Eϕ(Xx(T)) − Eϕ(XN) = C0(T, x)h + C(T, x)h2, C0(T, x) = E T B0(t, Xx(t))dt.
SLIDE 42
Corollary
Theorem 1: Eϕ(Xx(T)) − Eϕ(XN) = C0(T, x)h + C(T, x)h2, C0(T, x) = E T B0(t, Xx(t))dt. Theorem (2) Let Assumptions 1-2 hold. Then the coefficient C0(T, x) from (34) goes to zero as T → ∞ : |C0(T, x)| ≤ K(1 + |x|κ)e−λT (35) for some constants K > 0, κ ∈ N and λ > 0, i.e., over a long integration time the scheme (28) is of order two up to exponentially small correction.
SLIDE 43 Sketch of the proof
C0(T, x) = T EB0(t, Xx(t))dt = T
- Rd B0(t, y)p(t, x, y)dydt
(36) = T
+ T
- Rd B0(t, y)[p(t, x, y) − ρ(y)]dydt,
where p(t, x, y) is the transition density for (21) and ρ(y) is the invariant density.
SLIDE 44 Sketch of the proof
C0(T, x) = T EB0(t, Xx(t))dt = T
- Rd B0(t, y)p(t, x, y)dydt
(36) = T
+ T
- Rd B0(t, y)[p(t, x, y) − ρ(y)]dydt,
where p(t, x, y) is the transition density for (21) and ρ(y) is the invariant density.
σ2 U(y)
(37)
SLIDE 45 Discussion
1 We emphasize that the fact that the average of B0(t, x) with respect
to the invariant measure is equal to zero is the reason why the scheme (28) is second order accurate in approximating ergodic limits.
SLIDE 46 Discussion
1 We emphasize that the fact that the average of B0(t, x) with respect
to the invariant measure is equal to zero is the reason why the scheme (28) is second order accurate in approximating ergodic limits.
2 In the case of the Euler scheme (24) we get the same error
expansion as (33) for the scheme (28) but with a different B0(t, x) = BE
0 (t, x) (see [Milstein, T. (2004)]):
BE
0 (t, x) = 1
2
d
aj ∂u ∂xj ai ∂u ∂xi +σ2 2
d
∂2aj (∂xi)2 ∂u ∂xj + σ2 2
d
ai ∂3u ∂xi (∂xj)2 +σ2
d
∂aj ∂xi ∂2u ∂xj∂xi + σ4 6
d
∂4u (∂xi)2 (∂xj)2 . The average of BE
0 (t, x) with respect to the invariant measure is not
equal to zero and, consequently, the Euler scheme (24) approximates ergodic limits with order one – the same order as its weak convergence over a finite time interval.
SLIDE 47 Discussion
3 Let a one-step weak approximation ¯
Xt,x(t + h) of the solution Xt,x(t + h) of (21) generate a method of order p. The global error
R : = Eϕ(Xx(T)) − Eϕ(¯ Xx(T)) (38) = C0(T, x)hp + · · · + Cn(T, x)hp+n + O(hp+n+1) , where n ∈ N and the functions C0(T, x), . . . , Cn(T, x) are independent of h which can be presented in the form Ci(T, x) = T EBi(s, Xx(s))ds. One can deduce from the proof of Theorem 2 that if the averages of Bi(s, x) 0 ≤ i ≤ n, with respect to the invariant measure are equal to zero then in the limit of T → ∞ the scheme has p + n order of accuracy in h.
SLIDE 48 Discussion
3 Let a one-step weak approximation ¯
Xt,x(t + h) of the solution Xt,x(t + h) of (21) generate a method of order p. The global error
R : = Eϕ(Xx(T)) − Eϕ(¯ Xx(T)) (38) = C0(T, x)hp + · · · + Cn(T, x)hp+n + O(hp+n+1) , where n ∈ N and the functions C0(T, x), . . . , Cn(T, x) are independent of h which can be presented in the form Ci(T, x) = T EBi(s, Xx(s))ds. One can deduce from the proof of Theorem 2 that if the averages of Bi(s, x) 0 ≤ i ≤ n, with respect to the invariant measure are equal to zero then in the limit of T → ∞ the scheme has p + n order of accuracy in h. [Abdulle, Vilmart, Zygalakis 2014-15]
SLIDE 49 Numerical experiments
Anharmonic scalar model: the one-dimensional potential energy U(x) = cos(x) L2 error:
(ˆ ρi − ρi)2
0.2 0.4 0.8 Timestep L2 error
Euler-Maruyama Heun’s Method Candidate Scheme (1.7)
10-1 Error 10-3 10-4 10-2
Figure: The error in computed distributions is plotted for each scheme.
SLIDE 50 Error in finite time
L error
10-1 0.48
10-1
2
Euler-Maruyama Heun’s Method Candidate Scheme (1.7) Time Error using h = 0.16 10-2 10-3 10-4 1 2 3 4 5 6 7 8 9
10-4 10-2 10-3 0.16 0.24 0.32 0.48 0.16 0.24 0.32 0.48 0.16 0.24 0.32 0.48 0.16 0.24 0.32 0.48
Timestep Timestep Timestep Timestep Error at time t
10-4 10-1 10-2 10-3 0.16 0.24 0.32 0.48 0.16 0.24 0.32 0.48 0.16 0.24 0.32 0.48 0.16 0.24 0.32 0.48 0.16 0.24 0.32
Error at time t Timestep Timestep Timestep Timestep Timestep t=0.96 t=2.88 t=4.80 t=6.72 t=8.64 t=7.68 t=5.76 t=3.84 t=1.92
Figure: The lower plot shows the error in the distribution after time t, as computed using each scheme at h = 0.16. In the plots at the top, we compare the error growth with respect to stepsize h at multiples of t = 0.96.
SLIDE 51
Conclusions
Two new thermostats, one of Langevin type and one of gradient (Brownian) type, for rigid body dynamics are introduced. As in the deterministic case, it is important to preserve structural properties of stochastic systems for accurate long term simulations. Current work includes stochastic rigid body dynamics with hydrodynamic interactions. Development of more efficient methods for stochastic gradient systems.
SLIDE 52
Conclusions
Two new thermostats, one of Langevin type and one of gradient (Brownian) type, for rigid body dynamics are introduced. As in the deterministic case, it is important to preserve structural properties of stochastic systems for accurate long term simulations. Current work includes stochastic rigid body dynamics with hydrodynamic interactions. Development of more efficient methods for stochastic gradient systems. THANK YOU!
