Transcript IMA, November 2004 Equation-free Modeling For Complex Systems or Enabling Microscopic Simulators to perform System Level Tasks or Solving Differential Equations Without the Equations Or A Systems Approach to Multiscale.

IMA, November 2004
Equation-free Modeling
For Complex Systems
or
Enabling Microscopic Simulators to perform
System Level Tasks
or
Solving Differential Equations
Without the Equations
Or
A Systems Approach to Multiscale Computations
I. G. Kevrekidis -C. W. Gear and several other good peopleDepartment of Chemical Engineering, PACM & Mathematics
Princeton University, Princeton, NJ 08544
Sell rate increases
bad news
arrival rate
SELL
Buy rate increases
good news
arrival rate
Bearish
Neutral
Bad News Arrives
No news – drift
to Center
jump left --
Sell
rate
-1
Bad News
arrival rate
(Poisson)
If moves to extreme left
agent sells and becomes
neutral
BUY
Bullish
Good News Arrives
jump right +
0
1
Buy
rate
Good News
arrival rate
(Poisson)
Influence of other investors
- or lemmings falling off a cliff
If moves to extreme right
agent buys and becomes
neutral
SIMULATION RESULTS
50000 agents, ε+=0.075, ε - = - 0.072, v0+=v0-=20
Open loop response, g = 42
SIMULATION RESULTS
50000 agents, ε+=0.075, ε - =-0.072, v0+ = v0- =20
Open loop response, g = 46.5
J. Ottino, NWU
“simple”
“complicated”
“complex”
Multiscale / Complex System Modeling
“Textbook” engineering modeling:
macroscopic behavior through
macroscopic models
(e.g. conservation equations augmented by
closures)
Alternative (and increasingly frequent) modeling
situation:
• Models
– at a FINE / ATOMISTIC / STOCHASTIC level
–
MD, KMC, BD, LB (also CPMD…)
• Desired Behavior
What I will tell you:
Solve the equations WITHOUT writing them down.
Write “software wrappers” around “fine level” microscopic codes
Top level:
all algorithms we know and love (e.g. AUTO)
Bottom level: MD, kMC, LB, BD, heterogeneous/ discrete media,
CPMD, hybrid
INTERFACE:
Trade Function Evaluation
for “on demand” experimentaton and estimation
Think of the microscopic simulator AS AN EXPERIMENT
That you can set up and run at will
“Equation Free” (motivated by “matrix free iterative linear algebra”)
Algorithms (coarse integration, patch dynamics, coarse RPM…)
Tasks (stability/ bifurcation, control, optimization, dynamic renormalization)
Examples (LB, KMC, BD, MD), and some nebulous thoughts
dx
 f (x)
dt
Equations
Determinism (in a practical way)
Program
x1
x0
x
…..
x2
subroutine
f
t0
f(x)
Experimental Experimentalist
x1
x0
x0
t1
Estimate d/dt
t0
IN EFFECT, Do Forward Euler
not with the EQUATION, but with the (computational) EXPERIMENT
Look at the experiment & RESTART IT
t2
Legacy Code
Computational Experimentalist
x1
t0
t1
t1
RESTRICTION - a many-one mapping from a high-dimensional
description (such as a collection of particles in Monte Carlo
simulations) to a low-dimensional description - such as a finite
element approximation to a distribution of the particles.
LIFTING - a one-many mapping from low- to high-dimensional
descriptions.
We do the step-by-step simulation in the high-dimensional
description.
We do the macroscopic tasks in the low-dimensional description.
SIMULATION RESULTS
50000 agents, ε+=0.075, ε - = - 0.072, v0+=v0-=20
Open loop response, g = 42
Coarse projective integration: Accelerating things
Microscopic
Simulator
Microscopic
Simulator
0
1
x = x(ΙCDF)
x = x(ΙCDF)
x = x(ΙCDF)
x = x(ΙCDF)
0
1
0
0
1
1
Projection
[α1,…αns](k)
Simulation results at
g = 35, 200,000 agents
[α1,…αns](k+1)
[α1,…αns](k+2)
[α1,…αns](k+2+M)
Multiscale Modeling Challenges:
How to improve the spatial
technology?
Space
Proposal: detailed modeling in small spatial boxes with
interpolation between boxes - the “gap-tooth scheme”
Viscous Burgers equation:
kMC Realization
ut  uux  vuxx
Combined patch-dynamics
scheme
Solution
Time
Project
Integrate 2
Integrate 1
Restrict 2
Restrict 1
Initial Conditions 2
Initial Conditions 1
Lift
Boxes
Start from macroscopic description
Space
THE CONCEPT: What else can I do with an integration code ?
x  f (x)
Have equation
Do Newton
f (x )
Write Simulation
Compile
Do Newton on
x  (x)
x(t )
T
x(t   )
x  (x)  0
Estimate
Also
DΦ  ε
Φ(x ε)
xx x
( k 1 )
x 
xx
matrix-vector
product
Φ(x)
x
( k () k ) (k )
Matrix free
iterative linear algebra
The World
CG, GMRES
Newton-Krylov
“Coarse” Bifurcations (PNAS 2000, CChE 2002)
Bifurcation
Results
Parameter
Coarse Bifurcation Code
RPM-based
Local equilibrium assumption
e.g Maxwellian distribution
coarse IC
Averaging in time
And/or space and or nr.
of realizations and filtering
PDE-based
Timestepper
LIFT m
RESTRICT M
{
Microscopic IC’s
…
Microscopic
Timestepper
…
}
The Bifurcation Diagram
x( )
T
x(t   )
x  ( x, g )  0
Tracing the branch with arc-length continuation
H (x, g )  α(x x1 )   ( g  g1 )  S  0
400
m0 (x)
300
( x1  x 0 ) T
α
S
200
0.3

( g1  g 0 )
ΔS
100
0
-1
0.25
0
1
1 1.8
2
1.6
0.2
0.15
0.1
1.4
400
1.2
300
1
200
0.8
35
100
0
-1
0.05
30
0
35
40
1
40
45
50
g
45
g
50
The Bifurcation Diagram
50000 agents, g=35, ε+=0.075, ε-=-0.072, v0+=v0-=20
Open loop response. From unstable to stable markets
The Bifurcation Diagram
50000 agents, g=35, ε+=0.075, ε - = -0.072, v0+ = v0- = 20
Open loop response. Blow up
STABILIZING UNSTABLE M*****S
Feedback controller design
We consider the problem of stabilizing an equilibrium x*, p* of a dynamical system
of the form
x (t)  f ( x (t ), p (t )) , f : RN x R RN
where f and hence x* is characterized of uncertainty
To do this the dynamic feedback control law is implemented:
p  p *  K ( x D w )
w  x D w
Where w is a M-dimensional variable that satisfies
Choose matrices K, D such that the closed loop system is stable
At steady state:
w  0, x  0

p  p*
and the system is stabilized in it’s
“unknown” steady state
p  p*, x  x*
In the case under study the control variable is the exogenous arrival frequency of
“negative” information vex- and the controlled variables the coefficients of the orthogonal
polynomials used for the approximation of the ICDF
STABILIZING UNSTABLE MARKETS
Control variable: the exogenous arrival frequency of “negative” information v ex-
Coarse Self-Similar Solutions:
Dynamic Renormalization using Coarse Timesteppers.
An analogy: problems with traveling solutions (translational invariance)
Move along with the solution – it appears steady
u
 L(u)
Original Equation
(template-based reduction/reconstruction,
t
Rowley and Marsden, 2000)
v
v
Transformed equation
 c(v)
 L(v)
t
x
problems with focusing or collapsing solutions (scale invariance)
explode along with the solution - it appears steady
Original Equation
u
 L (u)
t
v
v
 G(v)(v  
)  L (v )
Transformed equation t

Templates: Aronson, Betelu and Kevrekidis, 2001
Rowley, Kevrekidis, Marsden and Lust, 2003)
Discrete time, lift –run – restrict:
Coarse renormalization flow integration / bifurcation analysis
Formulation
ut  Dx(u); Dx ( Bf ( x ))  Aa Bb Dy ( f ( y)), y  x
A
A
If u ( x, s)  B( s ) w( x , ( s)) 

x
If u ( x, t )  s U (  ) 
A( s )
s
  1  a  b
Dy ( w)  w  B w  A ywy (3)
(1)
B
A
Dy (U )   U   yUy (2)
(4)
 s  Aa Bb 1
Compare Eqns (2) (3):
 lim
 
 lim
Eqns (1) (5) determines  and 
A / A
 B / B

(5)
Time evolution of an initial
rectangular density profile
by the 1d diffusion equation
Regular
Dynamically renormalized
Schematic view of the dynamically
renormalized coarse timestepper
1D glass compaction model
Simulation Method
• 100000 particles at given density placed in
1D simulation box with periodic boundary
condition.
• Particles interact through hard-core
potential.
• Monte Carlo random walk is performed in
each step.
• Once a gap of unit size opened up
between two adjacent particles, an
additional block will deposit.
Glassy dynamics
Dynamics in linear time
Dynamics in logarithmic time

~ ln(t ) Eq. 7 Stinchcombe
1 
PRL 88, 125701 (2002)
Self-Similar Distribution
Dynamic Rescaling Fixed Point
Calculation at =0.9
GMRES Fixed Point Calculation
Dynamics With Limited
Distribution and Adjusted t0

1 
 k ln(t  t 0)  b
Reverse Projective Integration –
a sequence of outer integration steps backward;
based on forward steps + estimation
1
2
3
Reverse Integration:
and then
a little forward,
a lot backward !
We are studying the accuracy and stability of these methods
Forward and Backward
Integration
Forward Integration
Backward Integration
NLDbyMDofH2OinCNT
Water Confinement in Carbon Nanotube Molecular Channels
Model system for studying H20 transport in channel pores of membrane proteins
Classical Molecular Dynamics (MD)
- Flexible CNT (L = 13.5 Å; R = 8.1 Å)
- Graphite parameters
- 1034 TIP3P water molecules
- AMBER 6.0
- Particle-mesh Ewald electrostatics
- Gerhard Hummer
G. Hummer, J. C. Rasaiah, and J. P. Noworyta, Nature 414, 188-190 (2001)
Filling-Emptying Transition of H2O in CNT
Water
CNT
G. Hummer, J. C. Rasaiah, and J. P. Noworyta, Nature 414, 188-190 (2001)
Filling-Emptying Transition of H2O in CNT
“Equation Free” Approach
Analyse effect of CNT-Water Interactions on the H20 Occupancy
G. Hummer, J. C. Rasaiah, and J. P. Noworyta, Nature 414, 188-190 (2001)
The Heart of the Method – The “Timestepper”
1
 N ,   
M
i 1
i

xk
Choose
Parameters
 N  ; N , t  0;    N  ; N0 , t  0;  
M
0
d N
D N

dt
k BT
Xk+1
 G N 
 N
d
var N   2 D N
dt
“Lifting”

Constrained MD
“Coarse”
Initial Condition
var N
“FREE”
.
.
.
MD
.
.
.
N0
 = 1 ps
50 Copies
(Coarse “fine”)
• Constrain for  = 15 ps and sample configurations:  > 10 ps
• Initialize replicas: velocities from Maxwell-Boltzmann distribution
N

Bifurcation Diagram: N vs 
“Vestibule”
7
Fully Filled States
“Vestibule”
6
5
N
*
4
3
2
1
0
Fully Empty States
0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
w/ Gerhard Hummer, NIDDK / J.Chem.Phys. 03
Alanine Dipeptide
In 700 tip3p waters
The waters
The dipeptide
and the Ramachandran plot
Parametrizing nonlinear mainfolds
Given: (noisy) data in space, unordered, in high
dimension. Want: discover meaningful parameters.
Euclidean distances between any two points usually
not meaningful.
Euclidean distances between very close points is
meaningful. Use local Euclidean distances and glue
them to find paths between any two points.
Diffusion distance: Compute ALL paths between A
and B and take a weighted average (long paths count
less, but there are more of them). This is stable
under noise and behaves well with bottlenecks. It is
a diffusion process.
0.5
0
-0.5
B
-1
A
-1.5
-0.5
0
0.5
1
1.5
-0.6
Slowly communicating states:
2
2.5
B
-0.8
-1
A
-1.2
-1.4
-1.6
-1.8
geod.dist(A,B)≈geod.dist(B,C), diff.dist(A,B)>>diff.dist(B,C)
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
Graphs, Laplacian Eigefunctions, Embeddings, Heat Parametrization
We want to compute diffusion distances, and then deduce parametrization from those.
Data
Vertices of graph
Local distances
Conductivity of edges
B
Want to solve an heat/electrical network
equation: look at Laplacian on graph,
A
compute eigenfunctions
-4
x 10
-0.6
2
-0.8
-1
-1.2
1
-1.4
The eigenfunctions (long time behaviour) can be
used to compute diffusion distances:
and to parametrize, mapping the set into R n
-1.6
1.2
1.4
1.6
1.8
2
2.2
2.4
0
Color proportional to diffusion distance,
related to how much heat flows between A
and B in a certain time
A
B
C
** Belkin, Nyogi, Lafon
Computations:
- Order n, nlog(n) via eigenfunctions
- Order nlog(n) via diffusion wavelets (full
multiscale organization)
Slowly communicating states are mapped by the
All depend on decay of eigenvalues, i.e. whether
eigefunctions into far away points.
there is separation of time scales
Alanine Dipeptide Data
12-atom dipeptide fragments in water.
Long time coarse molecular dynamics by projective
integration.
Data consists of the coordinates of the atoms in the
molecule.
The eigenfunctions as parameters
They do a good job at parametrizing the set of
configurations of the molecule during the
simulation. In here we can see the image of the
trajectories in the coordinates given by the top
3 eigenfunctions, each point is a configuration,
red points are from the first (comprehensive)
simulation, blue points from second (mostly
stuck in the wells) simulation.
Eigenvalues of Laplacian on the set of structures
0.07
0.06
0.05
0.04
0.03
0.02
Big gap
0.01
0
-0.01
0
2
4
6
8
10
12
14
16
18
20
4 almost disconnected
components corresponding to
potential wells.
Potential wells
Transitions
Just from the gap in the spectrum,
one can see the set of
configurations embeds very well
in 4 dimensions.
The parametrization is dynamically
significant if the equation of motion
(differential/stochastic) are local and
smooth
(modulo
small
random
perturbations).
Rare Events: Escaping from stability
Open loop response, g = 45.8
Histogram of escape times
( x(t=0)=0 )
80
180 runs
70
60
50
40
30
20
10
0
0
50
100
150
200
250
300
Recursive Projection Method (RPM)
Reconstruct solution:
u = p+q = PN(p,q)+QF
Newton
iterations
Initial state un
F.P.I.
Timestepping
Legacy Code
F(un)
Subspace
Q =I-P
Subspace P of few NO
slow eigenmodes
• Treats timstepping routine, as a
“black-box”
– Timestepper evaluates un+1= F(un)
•
Recursively identifies subspace of
slow eigenmodes, P
•
Substitutes pure Picard iteration with
–Newton method in P
–Picard iteration in Q = I-P
Picard
iteration
Convergence?
YES
Final state uf
•
Reconstructs solution u from sum of
the projectors P and Q onto subspace
P and its orthogonal complement Q,
respectively:
–u = PN(p,q) + QF
Rapid Pressure Swing Adsorption
1-Bed 2-Step Periodic Adsorption Process
t=0 to T/2
•Isothermal operation
Ci(z=0)=PfYf/(RTf)
•Modeling Equations (Nilchan & Pantelides)
P(z=0)=Pf
z=L
z=0
Mass balance in ads. bed
Ci
qi
 (vCi )
 2Ci
t
 b

 Di
t
t
z
z 2
n
P
  Ci
RT i 1
P
180mv (1   b ) 2

z
d p2
 b3
qi
 ki (mi pi  qi )
t
Step 1 :
Pressurisation
t= T/2 to T
Ci
( z  0)  0
z
P(z=0)=Pw
Darcy’s law
Rate of ads
.
Step 2:
Depressurisation
RPSA simulation results
C1 (mole/kg-1)
70
65
60
COARSE INTEGRATION
55
50
50
45
45
40
34
35
40
32
30
0
30
35
0.2
28
0.4
0.6
Axial (m)
0.8
1
2
0
4
6
8
Time (sec)
10
12
30
26
25
24
22
20
20
160
180
200
220
240
260
280
15
10
5
0
500
1000
1500
2000
2500
3000
3500
Coming full circle
No equations ?
Isn’t that a little medieval ? Equations = “Understanding”, right ?
AGAIN matrix free iterative linear algebra
Ax=b
PRECONDITIONING,
BAx =Bb
B approximate inverse of A
Use “the best equation you have”
to precondition equation-free computations.
With enough initialization authority:
equation free laboratory experiments
Computer-Aided Analysis
of Nonlinear Problems in Transport Phenomena
Robert A. Brown, L. E. Scriven and William J. Silliman
in HOLMES, P.J., New Approaches to Nonlinear Problems in
Dynamics, 1980
ABSTRACT
The nonlinear partial differential equations of mass, momentum, energy,
Species and charge transport…. can be solved in terms of functions of limited differentiability,
no more than the physics warrants, rather than the analytic functions of classical analysis…
….. basis sets consisting of low-order polynomials. …. systematically generating and
analyzing solutions by fast computers employing modern matrix techniques.
….. nonlinear algebraic equations by the Newton-Raphson method. … The Newton-Raphson
technique is greatly preferred because the Jacobian of the solution is a treasure trove, not only
for continuation, but also for analysing stability of solutions, for detecting bifurcations of
solution families, and for computing asymptotic estimates of the effects, on any solution, of
small changes in parameters, boundary conditions, and boundary shape……
In what we do, not only the analysis, but the equations themselves are obtained on the
computer, from short experiments with an alternative, microscopic description.
EXISTING WORK
Coarse kMC of catalytic surface reactions
Coarse instabilities in gas-liquid multiphase flows (LB)
Coarse rheology of nematic liquid crystals (BD)
Homogenization / Effective Equations for random media
Coarse Optimal Paths / Coarse Optimization and Control (LB, kMC)
Coarse Molecular Dynamics (w/NIH, MD, CPMD)
Dislocation motion w/ diffusing impurities (kMC)
Legacy Simulators-Pressure Swing Adsorption, (gPROMS)
Coarse Computational Epidemiology
Influenza A virus evolution, (kMC)
Coarse granular flow ( MD)
Structural Transitions in Metal Crystals (MS, MD)
Coupled Oscillator Fields / Computational Neuroscience (w/ NIH; kMC)
Equilibrium Complex Fluids-Micelle Formation, (MC)
Bifurcations and Coarse Bifurcations in Renormalization Flows
Agent-based modeling in environmental and social sciences (kMC)
Coarse modeling of transport in model turbulent flows (SDEs)
Growth experiments – KPZ type equations (SDEs, kMC)