Computing the quasipotential for nongradient SDEs in 3D Samuel F. - - PowerPoint PPT Presentation

Computing the quasipotential for nongradient SDEs in 3D

Samuel F. Potter

Joint work with Maria Cameron and Shuo Yang

University of Maryland, College Park, MD

Mid-Atlantic Numerical Analysis Day, 2019

What are rare events in SDEs? dxt = b(xt) dt + √ǫ dwt

C1 vector field small parameter Brownian motion stable equilibria stable limit cycles more complex attractors

What do we want to estimate?

expected escape times from basins of attraction maximum likelihood escape paths quasi-invariant probability densities in the neighborhood of attractors

Gradient SDEs dxt = −∇V (xt) dt +

• 2β−1 dwt

Min #1 Min #2 Saddle

Gradient SDEs

Arrhenius equation (1884) Gibbs measure Kramer’s/Langer’s formula

(1940) (1969)

reaction rate ∝ exp

• Vsaddle−Vmin

kBT

• f(x) = Z−1 exp
• −V (x)

kBT

• reaction rate ≅ |λ|

2π

• det Hmin

| det Hsaddle|e−β(Vsaddle −Vmin)

Large deviation theory

Friedlin-Wentzell action functional Quasipotential Expected escape time ST(φ) = 1

2

T

0 ˙

φ − b(φ)2dt UA(x) = infφ,T {ST(φ)|φ(0) ∈ A, φ(T) = x} limβ→∞

1 β log E[τexit] = 1 2 infx∈∂B(A) UA(x) basin of attraction of A

attractor A basin of attraction B(A) transition state

Computing the quasipotential

Geometric action Hamilton-Jacobi equation Minimum action paths S(ψ) = L

• ψ′b(ψ) − ψ′ · b(ψ)
• ds

∇UA2 + 2b(x) · ∇UA = 0, UA(x) = 0, x ∈ A

dψ∗ ds ∇U(ψ∗) + b(ψ∗)

Orthogonal decomposition of b

b(x) = −∇U(x)/2 + ℓ(x) ℓ(x) = ∇U(x)/2 + b(x) ℓ(x) ⊥ ∇U(x)

(Heymann and Vanden-Eijnden, 2008) Nonlinear PDE for the quasipotential

Computing the quasipotential on a mesh: goal

Develop numerical methods for computing: ◮ Quasipotential: UA(x) ◮ Minimum action paths: ψ∗ What can we do with UA? ◮ More complete information about the dynamics in a region of interest ◮ A useful visualization tool

Computing the quasipotential on a mesh: key ideas

Motivated by Sethian’s fast marching method (1996) for solving the eikonal equation: F(x)∇U = 1

The ordered upwind method for solving HJ PDEs

F

• x,

∇U ∇U

• ∇U(x) = 1

(Sethian and Vladimirsky, 2001, 2003)

where 0 < Fmin ≤ F ≤ Fmax < ∞

Key idea and motivation

Dynamic programming principle U(x, y) x y

Computing the quasipotential on a mesh: previous methods

OUM for the quasipotential

∇U2 + 2b(x) · ∇U = 0 = ⇒

1 −2b(x)· ∇U

∇U ∇U = 1

(Cameron, 2012)

unbounded

Ordered line integral methods

(Dahiya and Cameron, 2017)

Solve the minimization problem directly: S(ψ) = L

• ψ′b(ψ) − ψ′ · b(ψ)
• ds

UA(x) = infψ

• S(ψ)
• ψ(0) ∈ A, ψ(L) = x
• Apply a time saving hierarchical update strategy

faster convergence and 100–1000 times more accurate than OUM 1.5–4 times faster than OUM

play video

Ordered line integral methods in 3D

Developed a 3D solver that: ◮ is 1–2 orders of magnitude faster than the OUM extended to 3D, ◮ improves the 2D hierarchical update strategy, ◮ skips higher-dimensional updates using the KKT conditions for constrained optimization, ◮ and searches for updates intelligently, ◮ without compromising accuracy.

Ordered line integral methods in 3D: simplex updates

Three types of updates needed: line, triangle, and tetrahedron updates

x xm

2

xm xm

1

x0 x1 x2 xλ xm

λ

x xm

1

xm

λ

xm x0 x1

Ordered line integral methods in 3D: local minimization problem

E.g., the minimization problem solved to compute U(x): minimize Uλ + Q(xλ, x) subject to λ1 ≥ 0, λ2 ≥ 0, λ1 + λ2 ≤ 1 Where:

Uλ = U(x0) + λ1(U(x1) − U(x0)) + λ2(U(x2 − x0)) xλ = x0 + λ1(x1 − x0) + λ2(x2 − x0) Q(xλ, x) = bm

λ x − xλ − bm λ · (x − xλ)

bm

λ = bm 0 + λ1(bm 1 − bm 0 ) + λ2(bm 2 − bm 0 )

bm

i = b

x + xi 2

• ,

i = 0, 1, 2.

Ordered line integral methods in 3D: simplex search

Ordered line integral methods in 3D: hierarchical update strategy

The idea: ◮ Line updates are incident in adjacent triangle updates ◮ Triangle updates incident in adjacent tetrahedron updates ◮ For h small enough, strictly convex minimization problems ◮ So: do line updates first, check KKT conditions for adjacent triangle updates, ditto adjacent tetrahedron updates

Ordered line integral methods in 3D: complexity

For N × N × N grid with N 3 grid points: ◮ Heap sort: O(N 3 log N) ◮ Line updates: O(f(N)N 3) ◮ Triangle updates: O(N 3) ◮ Tetrahedron updates: O(N 3) The function f(N) is not O(1)—more likely, f(N) = O(N) or O(N 2) with a very small constant

Numerical results: overview

Five examples to test accuracy ◮ Three examples using linear SDEs ◮ Two nonlinear examples with analytic UA Two examples by Tao (2018)—two point attractors with a hyperbolic periodic orbit as a transition state A genetic switch model with asymptotically stable equilibria

Numerical results: errors and runtimes

Extensive numerical tests in our JCP paper, brief synposis here

Example #1 #2 #3 T vs K (N = 513) [s] 31.4K1.79 8.57K1.80 93.8K1.78 E∞ vs h (K opt.) 0.181h0.750 0.0808h0.658 12.3h1.49 T vs N (K opt.) [ns] 27.4N4.05 61.0N4.11 131.0N3.98 Example #4 #5 T vs K (N = 513) [s] 66.3K1.80 1.94K1.80 E∞ vs h (K opt.) 0.635h0.989 19.4h1.54 T vs N (K opt.) [ns] 49.2N4.10 92.6N4.00

Numerical results: example E

Numerical results: a system with hyperbolic periodic orbits

Numerical results: a system with hyperbolic periodic orbits

Numerical results: genetic switch model

Numerical results: genetic switch model

Conclusion

◮ Extended the 2D OLIM to 3D, preserving the improved accuracy compared with the OUM ◮ Developed a hierarchical update strategy which reduces runtime in 3D on large meshes from days to hours ◮ Presented guidelines for choosing the neighborhood truncation parameter K

Future work

◮ Compute the quasipotential on (or around) the lower dimensional manifold occupied by the dynamics ◮ Optimal time complexity: O(N) ◮ Can the solver be made higher-order accurate?

References

Shuo Yang, Samuel F. Potter, and Maria K. Cameron. “Computing the quasipotential for nongradient SDEs in 3D.” Journal of Computational Physics 379 (2019): 325–350. Samuel F. Potter and Maria K. Cameron. “Ordered line integral methods for solving the eikonal equation.” Journal

• f Scientific Computing (2019): 1–41.

Cameron, Maria, and Shuo Yang. ”Computing the quasipotential for highly dissipative and chaotic SDEs an application to stochastic Lorenz’63.” Communications in Applied Mathematics and Computational Science 14.2 (2019): 207-246.