Transcript PowerPoint Presentation - Variable Reduction in “Heat Bath

Dimension Reduction in
“Heat Bath” Models
Raz Kupferman
The Hebrew University
Part I
Convergence Results
Andrew Stuart
John Terry
Paul Tupper
R.K
Ergodicity results: Paul Tupper’s talk.
Set-up: Kac-Zwanzig Models
A (large) mechanical system consisting of a
“distinguished” particle interacting through springs
with a collection of “heat bath” particles.
The heat bath particles have random initial
data (Gibbsian distribution).
Goal: derive “reduced” dynamics for the
distinguished particle.
Motivation
•Represents a class of problems where dimension
reduction is sought. Rigorous analysis.
•Convenient test problem for recent dimension
reduction approaches/techniques:
•optimal prediction (Kast)
•stochastic modeling
•hidden Markov models (Huisinga-Stuart-Schuette)
•coarse grained time stepping (Warren-Stuart, Hald-K)
•time-series model identification (Stuart-Wiberg)
The governing equations
The Hamiltonian:
n
n
p2j
1 2
H  Pn  V(Qn )  
  k j (q j  Qn )2
2
2m j j1
j1
The equations of motion:
QÝn  Pn
PÝn  V '(Qn )   k j (q j  Qn )
j
qÝj  p j m j
pÝj  k j (q j  Qn )
(Pn,Qn): coordinates of
distinguished particle
(pj,qj): coordinates of jth heat bath particle
mj: mass of j-th particle
kj: stiffness of j-th spring
V(Q): external potential
Initial data
Heat bath particles have random initial data: Gibbs
distribution with temperature T:
 H( p,q,P,Q) 
f ( p,q)  Z exp



T
1
Initial data are independent Gaussian:

p j (0)  m j T  j
q j (0)  Q(0)  T k j  j
i.i.d N(0,1)
Generalized Langevin Equation
Solve the (p,q) equations and
substitute back into the (P,Q)
equation
Ý P
Q
n
n
PÝn  V '(Qn )   k j (q j  Qn )
j
qÝj  p j m j
Ýj  k j (q j  Qn )
p
t
Ý
Ý
Qn  V '(Qn )    K n (t  s)QÝn (s)ds  Fn (t)

0
Memory kernel:
(Ford-Kac 65, Zwanzig 73)
n (t)   k j cos(  j t),  j  k j m j
K
j
Random forcing:
Fn (t)  T  k1j 2  j cos(  j t)   j sin( j t)
j
“Fluctuation-dissipation”:
E Fn (t)Fn (s)  T K n (t  s)
Choice of parameters
Heat baths are characterized by broad and dense
spectra. Random set of frequencies:
 j ~ U[0,n a ], a  (0,1),   n a n  0
Ansatz for spring constants:

k j  f 2 ( j ) 
Assumption:

f2() is bounded
and decays faster than 1/.
Generalized Langevin Eq.
t
Ý
Ý
Qn  V '(Qn )    K n (t  s)QÝn (s)ds  Fn (t)
0
Memory kernel:

K n (t)   f 2 ( j )cos( j t) 
j



0
f 2 ( )cos(t) dt
(Monte-Carlo
approximation of the
Fourier transform of
f2)
Random forcing:
F (t)  T
 f ( j ) j cos( j t)   j sin( j t) 
n
j
 T


0
f ( )cos(t) dB1 ( )  sin(t) dB2 ( )
(Monte-Carlo approximation of a stochastic integral)
Lemma
1. For almost every choice of frequencies (-a.s.) Kn(t)
converges pointwise to K(t), the Fourier cosine transform
of f2().
2. KnK in L2(,L2[0,T])
Theorem:
(-a.s.) the sequence of random functions Fn(t) converges
weakly in C[0,T] to a stationary Gaussian process F(t)
with mean zero and auto-covariance K(t); (FnF).
Proof: CLT + tightness
can be extended to “long term” behavior: convergence of
empirical finite-dimensional distributions (Paul Tupper’s talk)
Example
If we choose
f 2 () 
2
1
 2  2
then

Kn (t)  K(t)  exp( t )
and, by the above theorem, Fn(t) converges weakly to
the Ornstein-Uhlenbeck (OU) process U(t) defined by
the Ito
 SDE:
dU  U dt  2T dB

Convergence of Qn(t)
t
Ý
Ý
Qn  V '(Qn )    K n (t  s)QÝn (s)ds  Fn (t)
0
Theorem:
 (-a.s.) Q (t) converges weakly in C2[0,T] to the
n
solution Q(t) of the limiting stochastic integro-diff,
equation:
t
Ý
Ý
Q  V '(Q)    K(t  s)QÝ(s)ds  F(t)
0
Proof: the mapping (Kn,Fn)Qn is continuous

Back to example
if
f 2 () 
2
1
 2  2
then Qn(t) converges to the solution of
t  (ts)
Ý
Ý
Q  V '(Q)    e
QÝ(s)ds  U(t)

0
which is equivalent to the (memoryless!!) SDE
dQ  P dt

dP  V '(Q)  R dt
dR  (R  P) dt  2T dB

Numerical validation
Empirical distribution
of Qn(t) for n=5000
and various choices of
V(Q) compared with
the invariant measure
of the limiting SDE
double well
single well
triple well
extremely long correlation time
•“Unresolved” component of the solution are
modeled by an auxiliary, memoryless, stochastic
variable.
•Bottom line: instead of solving a large, stiff system
in 2(n+1) variables, solve a Markovian system of 3
SDEs!
•Similar results can be obtained for nonlinear
interactions. (Stuart-K ‘03)
Part II
Fractional Diffusion
Fractional (or anomalous) diffusion:
E Q(t) ~ t 
2
 1
Found in a variety of systems and models:
(e.g., Brownian particles in polymeric fluids,
continuous-time random walk)
In all known cases, fractional diffusion
reflects the divergence of relaxation
times; extreme non-Markovian behaviour.
Question: can we construct a heat bath
models that generated anomalous diffusion?
Reminder
t
Ý
Ý
Qn  V '(Qn )    K n (t  s)QÝn (s)ds  Fn (t)
0
Memory kernel:
Parameters:
K n (t)   k j cos(  j t),  j  k j m j
 j ~ U[0,n a ]
k j  f 2 ( j ) 
j
Random forcing:
Fn (t)  T  k1j 2  j cos(  j t)   j sin( j t)
j
If we take
f () 
2
C
1


C
then K n (t)  K(t)  
t
power law decay of
memory kernel
The limiting GLE
Theorem:
(-a.s.) Qn(t) converges weakly in C1[0,T] to the
solution Q(t) of the limiting stochastic integro-diff,
equation:
t
Ý
Ý
Q  V '(Q)    K(t  s)QÝ(s)ds  F(t)
0
C
K(t)  
t


F(t) is a Gaussian process with
covariance K(t); derivative of
fractional Brownian motion (1/fnoise)
(Interpreted in distributional sense)
Solving the limiting GLE
Ý(s)
Q
Ý
Ý V '(Q)   
Q
ds F(t)
0 (t  s)
t
For a free particle, V’(Q)=0, and a particle in a quadratic
potential well, V’(Q)=Q, the SIDE can be solved using the
Laplace
 transform.
Free particle: Gaussian profile, variance given by subdiffusive (Mittag-Leffler) function of time, var(Q)~t.
Quadratic potential: sub-exponential approach to the
Boltzmann distribution.
Numerical results
Variance of an ensemble of 3000 systems, V(Q)=0
(compared to exact solution of the GLE)
Quadratic well: evolving distribution of 10,000
systems (dashed line: Boltzmann distribution)
What about dimensional reduction?
Even a system with power-law memory can be
well approximated by a Markovian system with a
few (less than 10) auxiliary variables.
How? Consider the following Markovian SDE:
Ý
Ý V '(Q)  GT u
Q
uÝ GQÝ Au CWÝ
u(t) : vector of size m
A: mxm constant matrix
C: mxm constant matrix
G: constant m-vector


Ý
Ý V '(Q)  GT u
Q
uÝ GQÝ Au CWÝ
Solve for u(t) and substitute into Q(t) equation:

Ý
Ý V '(Q)    K (t  s)QÝ(s)ds F (t)
Q
0
t ~
~
where
~
K (t)  GT exp(At) G
~
F (t)  G exp(At) u(0) 
T
T
G
 0 exp(A(t  s))C dWs
t
Goal: find G,A,C so that fluc.-diss. is satisfied and
the kernel approximates power-law decay.
It is easier to approximate in Laplace domain:
Kˆ (s)  GT (A  sI)1 G
(and the Laplace transform of a power is a power).
The RHS is a rational function of degree (m-1) over m.
Pade approximation of the Laplace transform of the
memory kernel (classical methods in linear sys.
theory).
Even nicer if kernel has continued-fraction
representation
1

2
A  2



2
2
2
2
5
6
6
6
9
10
6 10 14


6 
10 



4m  3
2
1
 
1
G  1
 
 
 
1
Laplace transform of memory kernel (solid line) compared
with continued-fraction approximation for 2,4,8,16 modes
(dashed lines).
Variance of GLE (solid line) compared with Markovian
approximations with 2,4,8 modes.
Fractional diffusion scaling observed over long time.
Comment:
Approximation by Markovian system is not only a
computational tools. Also an analytical approach to study
the statistics of the solution (e.g. calculate stopping times).
Controlled approximation (unlike the use of a “Fractional
Fokker-Planck equation”).
Bottom line:
Even with long range memory system can be reduced (with
high accuracy) into a Markovian system of less than 10
variables (it is “intermediate asymptotics but that what we
care about in real life).