An exponential-free structure preserving integrator for stochastic Lie–Poisson systems

Erik Jansson

DAMTP, University of Cambridge and Wolfson College

 

 

Joint work with:

 

  • Annika Lang (Chalmers/University of Gothenburg)
  • Sagy Ephrati Imperial College London
  • Erwin Luesink (University of Amsterdam )

Made possible by:

 

  1. Klas Modin and Milo Viviani. Lie–Poisson methods for isospectral flows. Foundations of Computational Mathematics, 20:889–921, 2020.

  2. Grigori N. Milstein and Michael V. Tretyakov. Stochastic Numerics for Mathematical Physics. Springer Berlin Heidelberg, 2004.
  3. Charles-Edouard Bréhier, David Cohen, and Tobias Jahnke. Splitting integrators for stochastic Lie–
    Poisson systems. Mathematics of Computation, 92(343):2167–2216, 2023.

Talk based on:

asdasdasd

 

 

arxiv:2408.16701

Examples of systems

Stochastic rigid body

Stochastic (sine) Euler equations

Stochastic point vortex dynamics

Lie-Poisson systems

Lie group \(G\) with algebra \(\mathfrak{g}\)

Hamiltonian \(H\colon \mathfrak{g}^* \to \mathbb R\)

\,\dot X_t = \operatorname{ad}^*_{\nabla H(X_t)} X_t\,
\langle \operatorname{ad}_U^*X,V \rangle = \langle X,[U,V] \rangle

Two ingredients

Important geometric structure:

LP Systems evolve on coadjoint orbits

(symplectic manifolds that foliate space!)

+ Several preserved quantities, Hamiltonian, Casimirs etc.

\mathcal{O}_{X} = \{g^* X (g^*)^{-1}\colon g \in G\}

Coadjoint orbits

Dynamics on orbit

Lie-Poisson systems

\(G \subset \operatorname{GL}(n)\) is compact, simply connected and \(J\)-quadratic, i.e., \(A \in \mathfrak{g} \iff AJ + JA^* = 0 \)

Examples: \(\mathfrak{so}(N), \mathfrak{su}(N), \mathfrak{sp}(N)\).

In this setting, LP flow is isospectral:

\, \dot X_t = [\nabla H(X_t)^*, X_t] \,

Lie-Poisson systems

Example of LP systems

Classic example: rigid body (ODE)

  • Algebra: \(\mathfrak{so}(3) \cong \mathbb R^3\)
  • Bracket: \([v,w] = w \times v\)
  • Hamiltonian: \(H(v) = v \cdot Iv\)

Introducing stochasticisty

How to add noise to equations?

  • \(M\) independent scalar Brownian motions \(W^1, \ldots, W^M\)
  • \(M\) noise Hamiltonians \(H_k\colon\mathfrak{g}^*\to \R\)
~ \mathrm d X_t =[\nabla H_0(X_t)^*,X_t] \, \mathrm d t + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 
~ \mathrm d X_t =[\nabla H_0(X_t)^*,X_t] \, \mathrm d t + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

Why Strato?

Why these coefficients?

We need a chain rule

Same bracket-type coefficient: same geometric structure

Stochastic LP Systems evolve on coadjoint orbits and have Casimirs!

Introducing stochasticisty

Stochastic LP Systems evolve on coadjoint orbits and have Casimirs!

~ \mathrm d X_t =[\nabla H_0(X_t)^*,X_t] \, \mathrm d t + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

Coadjoint orbits

With transport noise

With non-transport noise

Introducing stochasticisty

Geometric structure is needed

~ \mathrm d X_t =[\nabla H_0(X_t)^*,X_t] \, \mathrm d t + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

What to assume to prove the existence of solutions?

E.g. Lipschitz is not realistic to assume

But! Some smoothness is sufficient for global existence

System remains on compact level sets of Casimirs
\(\rightarrow\) truncation argument to prove existence

Structure preservation: accurate numerics

Classic example: rigid body

  • Algebra: \(\mathfrak{so}(3) \cong \mathbb R^3\)
  • Bracket: \([v,w] = w \times v\)

Off the shelf integrator?

Non-physical behavior!

Structure preservation: accurate numerics

Structure-preserving integrators?

Structure preservation: accurate numerics

Wish list:

  • Preserve Casimirs a.s.
  • Preserve symplectic structure on coadjoint orbits a.s.

= Lie–Poisson integrator

Also: Avoid exponential map! (expensive, matrices in our applications are \(\sim2000\times2000\) )

Structure preservation: accurate numerics

Unreducing LP systems

LP systems on \(\mathfrak{g}^*\) are reductions of Hamiltonian systems on \(T^*G \cong G \times \mathfrak{g}^*\)

\(H\colon \mathfrak{g}^* \to \mathbb R\) lifts to left-invariant Hamiltonian!

\(G\) acts on \(T^*G\) by \(g.(Q,P) = (gQ,g^{-*}P)\)

Associated momentum map:

\mu(Q,P) = \frac{1}{2} Q^* P -\frac{1}{2c}J P^* Q J \in \mathfrak{g}^*
\tilde{H}(Q,P) = H(\mu(Q,P))

Do for all Hamiltonians!

\begin{split} \mathrm d Q_t &= Q_t\nabla H_0(\mu(Q_t,P_t)) \, \mathrm d t + \sum_{i=1}^M Q_t \nabla H_i (\mu(Q_t,P_t))\circ \mathrm d W_t^i, \\ \mathrm d P_t &= -P_t\nabla H_0(\mu(Q_t,P_t))^* \, \mathrm d t - \sum_{i=1}^M P_t \nabla H_i (\mu(Q_t,P_t))^* \circ \mathrm d W_t^i. \end{split}

Note: No a priori guarantee that this system remains on \(T^*G\)!

Embedded into \((R^{n \times n})^2\)

However! Remains on \(T^*G\)

Unreducing LP systems

\begin{split} \mathrm d Q_t &= Q_t\nabla H_0(\mu(Q_t,P_t)) \, \mathrm d t + \sum_{i=1}^M Q_t \nabla H_i (\mu(Q_t,P_t))\circ \mathrm d W_t^i, \\ \mathrm d P_t &= -P_t\nabla H_0(\mu(Q_t,P_t))^* \, \mathrm d t - \sum_{i=1}^M P_t \nabla H_i (\mu(Q_t,P_t))^* \circ \mathrm d W_t^i. \end{split}

\(X_t = \mu(Q_t,P_t) \in \mathfrak{g}^*\) satisfies

~ \mathrm d X_t =[\nabla H_0(X_t)^*,X_t] \, \mathrm d t + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

Unreducing LP systems

\(\Phi_h \colon T^*G \to T^* G\): \(G\)-equivariant symplectic integrator \(\rightarrow\)
Reduces to a Lie–Poisson integrator by \(\mu\)

Exploit the unreduction!

\begin{split} \mathrm d Q_t &= Q_t\nabla H_0(\mu(Q_t,P_t)) \, \mathrm d t + \sum_{i=1}^M Q_t \nabla H_i (\mu(Q_t,P_t))\circ \mathrm d W_t^i, \\ \mathrm d P_t &= -P_t\nabla H_0(\mu(Q_t,P_t))^* \, \mathrm d t - \sum_{i=1}^M P_t \nabla H_i (\mu(Q_t,P_t))^* \circ \mathrm d W_t^i. \end{split}
~ \mathrm d X_t =[\nabla H_0(X_t)^*,X_t] \, \mathrm d t + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

Deriving the integrator

Do we know such an integrator?

Implicit midpoint!

\begin{split} \Phi_h^{(1)} &\colon \begin{cases} Q_n & = \tilde{Q}- \frac{1}{2}\left( f_0(\tilde{Q}, \tilde{P}) h + \sum_{i=1}^M f_i(\tilde{Q}, \tilde{P})(\zeta_{i})_n \sqrt{h}\right), \\ P_n & = \tilde{P}-\frac{1}{2}\left(k_0(\tilde{Q}, \tilde{P}) h + \sum_{i=1}^M k_i(\tilde{Q}, \tilde{P})(\zeta_{i})_n\sqrt{h} \right), \end{cases} \\ \Phi_h^{(2)} &\colon \begin{cases} Q_{n+1} & = \tilde{Q} + \frac{1}{2}\left( f_0(\tilde{Q}, \tilde{P}) h + \sum_{i=1}^M f_i(\tilde{Q}, \tilde{P})(\zeta_{i})_n \sqrt{h}\right), \\ P_{n+1} & = \tilde{P}+\frac{1}{2}\left(k_0(\tilde{Q}, \tilde{P}) h + \sum_{i=1}^M k_i(\tilde{Q}, \tilde{P})(\zeta_{i})_n\sqrt{h} \right), \end{cases} \end{split}\\ \Phi_h = \Phi_h^{(2)} \circ \Phi_h^{(1)}
\zeta_h = \begin{cases} & A_h \quad \text{if } ~ \xi>A_h, \\ &-A_h \quad \text{if } ~ \xi<-A_h,\\ & \xi \quad \text{if} ~ |\xi| \leq A_h, \end{cases}
  • Symplectic
  • Equivariant
  • Well studied

Deriving the integrator

Stochastic isospectral midpoint:

Take \(Q_n,P_n\) such that \(\mu(Q_n,P_n) = X_n\)

\begin{split} \tilde{\Psi}_{h, n}(\tilde{X}) &= \nabla H_0(\tilde{X})^*h + \sum_{i=1}^M\nabla H_i(\tilde{X})^*(\zeta_i)_n\sqrt{h}, \\ X_n &= \left(I- \frac{1}{2}\tilde{\Psi}_{h, n}(\tilde{X})\right) \tilde{X} \left(I+\frac{1}{2}\tilde{\Psi}_{h, n}(\tilde{X})\right), \\ X_{n+1} &= \left(I+\frac{1}{2}\tilde{\Psi}_{h, n}(\tilde{X})\right) \tilde{X} \left(I- \frac{1}{2}\tilde{\Psi}_{h, n}(\tilde{X})\right). \end{split}

Take \(X_{n+1} = \mu(\Phi_h(Q_n,P_n))\)

Deriving the integrator

Error analysis: essentially without crying

Exploit the unreduction!

\begin{split} \mathrm d Q_t &= Q_t\nabla H_0(\mu(Q_t,P_t)) \, \mathrm d t + \sum_{i=1}^M Q_t \nabla H_i (\mu(Q_t,P_t))\circ \mathrm d W_t^i, \\ \mathrm d P_t &= -P_t\nabla H_0(\mu(Q_t,P_t))^* \, \mathrm d t - \sum_{i=1}^M P_t \nabla H_i (\mu(Q_t,P_t))^* \circ \mathrm d W_t^i. \end{split}
~ \mathrm d X_t =[\nabla H_0(X_t)^*,X_t] \, \mathrm d t + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

Check the literature for error analysis for implicit midpoint

Error analysis: actually with some crying

 

1) Show that the implicit midpoint method does not travel too far:

\(\sup_{h \geq 0} \sup_{n \geq 0} \|Q_n,P_n\| \leq R(Q_0,P_0). \)

2) Truncate the CHS system on \(T^*G\)

Upstairs

Downstairs

3) Apply error analysis from literature

\begin{split} \mathrm d Q_t &= Q_t\nabla H_0(\mu(Q_t,P_t)) \, \mathrm d t + \sum_{i=1}^M Q_t \nabla H_i (\mu(Q_t,P_t))\circ \mathrm d W_t^i, \\ \mathrm d P_t &= -P_t\nabla H_0(\mu(Q_t,P_t))^* \, \mathrm d t - \sum_{i=1}^M P_t \nabla H_i (\mu(Q_t,P_t))^* \circ \mathrm d W_t^i. \end{split}

Reduction by \(\mu\)

~ \mathrm d X_t =[\nabla H_0(X_t)^*,X_t] \, \mathrm d t + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

4) Convergence of method for truncated LP system

  1. Strong convergence: \(\mu\) is Lipschitz
  2. Weak convergence: \(\phi \circ \mu\) is valid test func

5) Truncated LP system = LP system

Recipe:

  • Assume: \(H_0 \in C^4\) and \(H_1,H_2, \ldots, H_M \in C^5\).
  • \(T>0\): be a fixed final time,
  • \(N \in \mathbb{N}\): number of steps
  • \(X_0 \in \mathfrak{g}^*\): Deterministic initial condition
  • \(\phi \in C^5\): Test function

 

~\sup_{0 \leq n \leq N} \mathbb E\left[ \|X(hn)-X_n\|^2\right]^{1/2} \leq \kappa h^{1/2}\\ \sup_{0 \leq n \leq N} \left|\mathbb E[\phi(X(hn))-\phi(X_n)]\right| \leq \kappa h~

Error analysis: actually with some crying

~\sup_{0 \leq n \leq N} \mathbb E\left[ \|X(hn)-X_n\|^2\right]^{1/2} \leq \kappa h^{1/2}\\ \sup_{0 \leq n \leq N} \left|\mathbb E[\phi(X(hn))-\phi(X_n)]\right| \leq \kappa h~

Geometric structure of equations and its preservation is central. Used to prove

  1. Existence and uniquness
  2. Convergence

Without structure preservation, no convergence guarantees

Error analysis: actually with some crying

Mandatory error plots

Strong error

Weak error

Rigid body

500 realizations

10^7 realizations

\(\phi(X) = \sin(2\pi x_1) + \sin(2\pi x_2) + \sin(2\pi x_3)\).

What's all this for?

e_1,e_2,e_3 \in C^\infty(\mathbb{S}^2)
\Delta = \sum_{i=1}^3 \{e_k,\{e_k,\cdot \}\}
\dot \omega = \{\Delta_{\mathbb S^2}^{-1} \omega, \omega \}, \omega \in C^\infty_0(\mathbb S^2)
\dot W = [\Delta_{N}^{-1} W, W ], W \in \mathfrak{su}(N)
\Delta_N = \sum_{i=1}^3 [X_i,[X_i,\cdot]]
X_1,X_2,X_3 \in \mathfrak{u}(N)
\dot X = [\Delta_{N}^{-1} X, X ], X \in \mathfrak{su}(N)
\Delta_N = \sum_{i=1}^3 \{X_k,\{X_k,\cdot \}\}
X_1,X_2,X_3 \in \mathfrak{u}(N)
~ \mathrm d X_t = [\nabla H_0(X_t)^*,X_t] \mathrm dt + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

Zeitlin is a convenient computational method for 2D ideal spherical hydrodynamics. 

What's all this for?

\dot X = [\Delta_{N}^{-1} X, X ], X \in \mathfrak{su}(N)
\Delta_N = \sum_{i=1}^3 \{X_k,\{X_k,\cdot \}\}
X_1,X_2,X_3 \in \mathfrak{u}(N)
~ \mathrm d X_t = [\nabla H_0(X_t)^*,X_t] \mathrm dt + \sum_{k=1}^M [\nabla H_k(X_t)^*,X_t] \circ \mathrm d W_t^k ~ 

Matrices are very large, makes exponential map/algebra to group mappings too expensive.

What's all this for?

What's all this for?