Models used in simulation may give accurate short-term trajectories but distort long-term (statistical) properties. In this work, we augment a given approximate model with a control law (a ‘thermostat’) that gently perturbs the dynamical system to target a thermodynamic state consistent with a set of prescribed (possibly evolving) observations. As proof of concept, we provide an example involving a point vortex fluid model on the sphere, for which we show convergence of equilibrium quantities (in the stationary case) and the ability of the thermostat to dynamically track a transient state.

1. Introduction

In applications of modern computational science, the underpinning physical laws (and equations of motion) are often well established, yet the detailed behaviour is unpredictable on long time scales due to the presence of deterministic chaos. Examples of this arise in molecular dynamics [1,2], in polymer simulation [3] and in the modelling of turbulent fluids in the atmosphere and ocean [4,5], where long simulations are used to extract statistical information (e.g. the statistics of rare transitions between basins in molecular dynamics, or slow relaxation processes in fluids). In this type of dynamical sampling, the common requirement is that simulated paths are sufficiently accurate to calculate measures of dynamical mixing such as 2-point temporal correlation functions [6].

In Hamiltonian systems such as molecular dynamics, it is common to run ensembles of microcanonical (i.e. constant energy) simulations. For such simulations, backward error analysis [7,8] suggests that a long simulation should not be viewed as the approximation of a particular trajectory of the system but rather as an exact solution of a perturbed continuum process described by modified equations. In the case of dynamical sampling of complex systems, the statistics of simulation data are biased because they sample an invariant measure of the modified equations (rather than those of the target system), i.e. the time-discretization error induces an effective statistical bias.

Statistical bias may also arise due to spatial discretization. For example, in the setting of geophysical fluid dynamics, a comparison of discretizations of the quasi-geostrophic equations reveals that the long-time mean potential vorticity field and pointwise fluctuation statistics are heavily dependent on discrete conservation laws such as energy, enstrophy and material conservation of vorticity [9–11]. The discretization bias may be controlled by reducing timestep or by incorporating a Metropolis condition [12,13], possibly increasing the computational cost. The use of a Metropolis condition can also increase the difficulty of computing accurate dynamical properties. The need for computations to be accurate with respect to the dynamical process as well as to the long-term statistics thus poses difficult challenges to the simulator.

In this paper, we consider a method of perturbing the dynamics of the system to correct for statistical bias. Our method is based on the concept of a thermostat, by which we mean a scheme based on a control law designed to facilitate sampling of a target probability density function (pdf). Although originally proposed in molecular dynamics (see, e.g. the Nosé-Hoover method [14,15]), thermostats can be extended to handle a wide variety of systems and have been used in a wide variety of applications. For example, in Dubinkina et al. [16], a thermostat was used as a model reduction technique for a vortex model of a fluid (suppressing the detailed interactions of a few strong vortices with a weak vortex field), an approach we draw on in the numerical experiments of this article. In another recent article [17], thermostats have been suggested as a means of sampling incompletely specified systems (with noisy gradients), with applications in machine learning for Bayesian inference in the ‘big data’ context.

In the general modelling scenario, thermostats are significantly hampered by the need to specify a functional form of the pdf a priori. In this article, we assume that, instead of the pdf, what is actually available is a partial set of expectations of observables with respect to the unknown probability measure, which may arise from experiment or other types of modelling. In this setting, information theory (in particular entropy maximization [18,19]) offers tools for constructing densities consistent with observations that are close to a given prior distribution (see [20–22] for examples in geophysical modelling, [23,24] for discussion in molecular modelling). Our algorithm for computing the parameters finds Lagrange multipliers (one for each observable) that modify the probability density, using iterative techniques like those in Agmon et al. [25], Haken [26], Davis & Gutiérrez [27], where each step of iteration involves ensemble simulation that reweights samples from the prior distribution appropriately. The density resulting from such an entropy maximization is then used as the basis for designing a thermostat that ensures the extended dynamical system samples the distribution and consequently exhibits the correct (long-term) statistical averages, while mildly perturbing the short-term dynamical behaviour of the system. The complexity and generality of our iterative procedure (sampling $\to$ parametrization $\to$ sampling…) means that convergence has to be verified through computer experiment.

Entropy maximization also forms the basis for closure approximations in the works of Ilg et al. [3] and Samaey et al. [28], where the result from entropy maximization is used in determining evolution equations for resolved macroscopic quantities. This approach is distinct in philosophy from the current work in which the result of entropy maximization is used to define a statistical ensemble and evolution equations for microscopic dynamics that are consistent with macroscopic observations.

While we motivate the method in the setting of equilibrium sampling of a fixed target distribution, we can easily apply our scheme to the case of a more general model with unknown or even undefined steady state. For this purpose, an ensemble of simulations are run in finite duration bursts, approximating the Lagrange multipliers via an iterative procedure based on the simulation history. This procedure allows entropy maximization to be performed with ensemble averages that evolve in time, providing a flexible method of nonlinear data assimilation.

The remainder of this article is organized as follows. In §2, we review the idea of a thermostat which is the key tool we employ to control the invariant distribution. Section 3 describes the combination of thermostats with the maximum entropy framework for correcting the density to reflect thermodynamic constraints. Section 4 describes the adaptive procedure. In §5, we apply and evaluate the method in the setting of a system of point vortices on the surface of a sphere, which represents a simple geophysical model with multiple statistically relevant first integrals.

2. Thermostats

We introduce the thermostat technique for generic Hamiltonian dynamical systems of the form

\frac{d y}{d t} = f (y) = B (y) \nabla H (y), y (t) \in D, B (y) = - B {(y)}^{T} and H (y) : D \to R,

2.1

possessing a divergence-free vector field¹ ∇⋅f≡0. Invariance of the Hamiltonian H along solutions of (2.1) follows from (d/dt)H(y(t))=∇H⋅dy/dt=∇H⋅B∇H=0, due to skew-symmetry of B(y). Additional first integrals may exist, denoted by I_ℓ(y):∇I_ℓ⋅f≡0, ℓ=1,…,L.

Thermostats were introduced in molecular dynamics to model the trajectories of molecules in a fluid at constant temperature. In the common statistical mechanical framework, it is assumed that the trajectories of a system of particles in thermal equilibrium with a reservoir at constant temperature sample the canonical distribution, with smooth invariant density proportional to exp(−βH) where β is a parameter reciprocal to temperature scaled by Boltzmann’s constant. The trajectories of the Hamiltonian system (2.1) are restricted to a level set of energy. Hence, to model a system at constant temperature, it is necessary to perturb the vector field to make trajectories ergodic with respect to the canonical distribution. The most common way of achieving this is by adding stochastic and dissipative terms satisfying a fluctuation–dissipation relation (Langevin dynamics), which is rigorously ergodic [29].

However, DelSole [30] points out that additive stochastic forcing of trajectories leads to inaccurate dynamical quantities as autocorrelation functions are strongly perturbed. For smooth deterministic Hamiltonian dynamics, normalized velocity autocorrelation functions are of the form 1−cτ², c>0 in the zero-lag limit $τ \to 0$ , whereas the autocorrelation of a variable that is directly forced by white noise must take the form $\exp (- κ τ)$ , κ>0 in the same limit. This implies that direct stochastic perturbation leads to autocorrelation functions that have non-zero slope and opposite curvature at zero lag [31].

An alternative approach, due to Nosé [14,15] and Hoover [32], augments the phase space with an additional thermostat variable ξ(t) (driven by a differential equation) in order to ensure that the extended dynamics on R^d×R preserves an equilibrium density whose marginal on R^d is the target (e.g. Gibbs) density. The deterministic thermostats of Nosé and Hoover are not ergodic, but they can be combined with stochastic forcing of the thermostat variable ξ, leading to the so-called Nosé–Hoover–Langevin method [33]. For generic divergence-free dynamical systems, the generalized Bulgac–Kusnezov (GBK) [34,35] thermostat provides the augmented system:

d y = f (y) d t + ξ ε g (y) d t

2.2a

and

d ξ = ε h (y) d t - γ ξ d t + \sqrt{2 γ} d w,

2.2b

where ε>0 and γ>0 are parameters, w(t) is a scalar Wiener process, and g and h are vector fields which we discuss in more detail below. Given a target density

ρ (y) \propto \exp (- A (y))

A : D \to R

, denote the augmented product density by

\tilde{ρ} (y, ξ) = ρ (y) \cdot μ (ξ)

, with μ a univariate normal distribution with mean zero and standard deviation one. It is easily checked that

\tilde{ρ}

is stationary under the Fokker–Planck operator associated with (2.2) provided

h (y) = \nabla \cdot g - g \cdot \nabla A

2.3

and provided the ‘potential’ A is only a function of y through the first integrals of f. The (extended) target measure

\tilde{ρ}

is ergodic if the vector fields f and g satisfy a Hörmander condition [36]. Up to verification of the Hörmander condition, the choice of g is free. In the example of §5, we consider a specific form for g.

The parameter ε controls the strength of the thermostat relative to that of the vector field f. This affects the rate at which the invariant measure is sampled. It has been demonstrated in Bajars et al. [36,37] and Leimkuhler et al. [31] that GBK/NHL thermostating gives a weak perturbation of microcanonical dynamics: typical autocorrelation functions match the exact ones to leading order, i.e. have the form 1−cτ²+O(τ³), as $τ \to 0$ .

Designing a thermostat relies on knowing the functional form of the target distribution, which may not always be available. In this article, we relax this requirement by constructing the target distribution iteratively and adaptively, based on Jaynes’ principle of least-biased ensemble prediction.

3. Bias correction method

Suppose that we are given a simplified dynamical model for the evolution of some projection, i.e. a ‘coarse graining’ of the phase variables (coarse-grained variables y(t)∈R^d) which is in the form (2.1). Although the original system is complex and its details unknown, we assume that we can obtain in some way (e.g. through measurement) a collection of observations of mean values of functions of the reduced variables. That is, there are functions $C_{k} : R^{d} \to R$ , k=1,2…,K and given values c_k, k=1,2,…,K, such that

c_{k} = ⟨ C_{k} (y) ⟩, k = 1, \dots, K,

3.1

where 〈C_k(y)〉 represents averaging with respect to the true, empirical invariant measure of the dynamical system. The use of a GBK thermostat (2.2) limits our implementation to observables that are functions of the conserved quantities {H,I_ℓ,ℓ=1,…,L}, that is, C_k(y)=C_k(H(y),I₁(y),…,I_L(y)), k=1,…,K. This restriction is only due to the choice of thermostat, the entropy maximization has no such restriction. Our goal is to find a perturbed dynamical model for the reduced variables which (i) is compatible with the thermodynamic constraints (3.1), and (ii) mildly perturbs the dynamics compared to those of the native model.

(a) Entropy maximization

Empirical information theory generalizes the principle of insufficient reason, by proposing the least-biased probability density consistent with a set of observations. This idea was first proposed by Jaynes [18,19] and is treated in the monographs of Haken [26], Dewar et al. [38] and in the context of geophysical fluids by Majda & Wang [39]. The least-biased density is defined as the probability density ρ(y) that maximizes the information entropy functional $S [ρ] = - \int_{D} ρ (y) \log ρ (y) d y,$ subject to a set of constraints given by observations. When $D$ is a compact set and there are no observations, the minimizer is the uniform density $ρ \equiv | D |^{- 1}$ . The entropy $S$ is the unique measure of uncertainty that is positive valued, monotonically increasing as a function of uncertainty, and additive for independent random variables. With observable functions {C_k(y) | k=1,2,…K} let

E_{ρ} C_{k} = \int_{D} C_{k} (y) ρ (y) d y

3.2

denote expectation in the (as yet undetermined) density ρ. Defining Lagrange multipliers λ_k, k=1,…,K, associated with the observables C_k, we obtain the solution ρ from a stationarity condition for the Euler–Lagrange functional

W (ρ, λ_{1}, \dots λ_{K}) := [S [\hat{ρ}] - \sum_{k = 1}^{K} λ_{k} (E_{\hat{ρ}} C_{k} (y) - c_{k})] .

Conditions for a unique solution to exist are discussed in Balian [40]. When it exists, the maximum entropy solution satisfies

ρ (y) = λ_{0} \exp (- λ_{1} C_{1} (y) - \dots - λ_{K} C_{K} (y)),

where λ₀ is chosen to satisfy

\int_{D} ρ d y = 1

, and λ_k is chosen such that

E_{ρ} C_{k} (y) = c_{k}

In some cases, besides the observations, we may be given prior statistical information on the process y(t). The Kullback–Leibler divergence, or relative entropy,

S_{π} [ρ (y)] = \int ρ (y) \ln \frac{ρ (y)}{π (y)} d y

represents a (non-symmetric) distance between measures. Note the change of sign, which makes the relative entropy more like a distance measure; we stick to this convention. It quantifies the information lost in approximating ρ(y) by π(y).

Suppose y is a random variable with distribution (law) y∼ρ, where ρ is unknown. Suppose further, that we are given a prior distribution π, presumed to be close to ρ, and a set of K observations (3.1). Following Jaynes [18,19], the least-biased distribution ρ consistent with the observations c_k and prior π solves the constrained minimization problem

ρ = \underset{\hat{ρ}}{\arg \min} [S_{π} - λ_{0} (1 - \int \hat{ρ} (y) d y) - \sum_{k = 0}^{K} λ_{k} (c_{k} - \int C_{k} (y) \hat{ρ} (y) d y)],

where the λ_k are Lagrange multipliers to enforce the condition that the expectations (3.2) agree with the observations (3.1). The solution to the variational problem is

ρ (y) = λ_{0} \exp (- λ_{1} C_{1} (y) - \dots - λ_{K} C_{K} (y)) π (y),

3.3

where the Lagrange multipliers λ_k are chosen consistently with the observations (3.1) and λ₀ is a normalization constant so that ρ is a probability density function.

Calculation of the Lagrange multipliers is discussed in [25–27]. We use the following algorithm based on re-weighting. Assume we are given a sequence of samples yⁿ, n=1,…,N, distributed according to a known prior distribution π(y), i.e. yⁿ∼π. The expectation under π(y) of a function Φ(y) has the consistent and unbiased estimator

{\hat{Φ}}^{π} = \frac{1}{N} \sum_{n = 1}^{N} Φ (y^{n}) .

Given the posterior distribution ρ(y) of the form (3.3), compute the expectation $E_{ρ} Φ$ by re-weighting of the integral

E_{ρ} Φ = \int Φ (y) ρ (y) d y = λ_{0} \int Φ (y) e^{- \sum_{i = 1}^{K} λ_{i} C_{i} (y)} π (y) d y = λ_{0} E_{π} {Φ (y) λ_{0} e^{- \sum_{i = 1}^{K} λ_{i} C_{i} (y)}},

yielding an unbiased estimator for

E_{ρ} Φ

given by

{\hat{Φ}}^{ρ} = \frac{λ_{0}}{N} \sum_{n = 1}^{N} Φ (y^{n}) e^{- \sum_{i = 1}^{K} λ_{i} C_{i} (y^{n})} .

We wish to ensure that the estimators ${\hat{C}}_{k}^{ρ}$ for observables C_k match their empirical values c_k, i.e.

c_{k} = {\hat{C}}_{k}^{ρ} = \frac{λ_{0}}{N} \sum_{n = 1}^{N} C_{k} (y^{n}) e^{- \sum_{i = 1}^{K} λ_{i} C_{i} (y^{n})}, k = 1, \dots, K .

We can use this fact to define a Newton–Raphson iteration to determine the Lagrange multipliers λ_k. Define the residual r with components

r_{k} (λ) = c_{k} - \frac{λ_{0}}{N} \sum_{n = 1}^{N} C_{k} (y^{n}) e^{- \sum_{i = 1}^{K} λ_{i} C_{i} (y^{n})}, k = 1, \dots, K,

with λ=(λ₁,…,λ_K) and r=(r₁(λ),…,r_K(λ)). Note that λ₀ can be viewed as a function of λ₁,λ₂,…,λ_K chosen from the normalization condition, i.e.

λ_{0} = {[\sum_{n = 1}^{N} e^{- \sum_{j = 1}^{K} λ_{j} C_{j} (y^{n})}]}^{- 1} .

The Jacobian matrix J=(J_kj) of the vector function r is determined as

J_{k j} (λ) := \frac{\partial r_{k}}{\partial λ_{j}} = \frac{λ_{0}}{N} \sum_{n = 1}^{N} C_{k} (y^{n}) C_{j} (y^{n}) e^{- \sum_{i = 1}^{K} λ_{i} C_{i} (y^{n})} j, k = 1, \dots, K .

The iteration then proceeds as

λ^{α + 1} \leftarrow λ^{α} - J^{- 1} (λ^{α}) r (λ^{α})

4. Adaptive determination of Lagrange multipliers

When the observations deviate far from the prior, the Lagrange multipliers take larger values and the weights of the samples diverge, i.e. some weights become very large and others become very small. This results in a larger variance on the estimates and consequently in poorer results for the Lagrange multipliers. Moreover, the approach of the previous section excludes applications where (i) the statistical knowledge is expected to improve as the simulation progresses, (ii) the average observables are known to vary slowly with time or (iii) it is unfeasible to construct a large enough ensemble distributed in the prior. For these cases we consider using the simulation data of a small ensemble (propagated in short bursts of M timesteps) for updating the Lagrange multipliers for mean observation data. This results in an adaptive algorithm for obtaining the Lagrange multipliers on-the-fly during simulation. In this situation, the sample set used is drawn from a distribution closer to the target, resulting in better estimations for the residuals and Jacobian.

Let us assume that we have a (moderately sized) ensemble of initial conditions according to some distribution

ρ_{0} (y) = π (y) \exp (- \sum_{k = 1}^{K} λ_{k}^{0} C_{k} (y)) .

4.1

In our numerical experiments, we assume

λ_{k}^{0} = 0

, ∀ k and hence ρ₀=π, but we present the general form here nevertheless. Rather than first computing the Lagrange multipliers corresponding to observations and then running the simulation over the time interval of interest, we perform a short burst of the thermostatted simulation for P ensemble members, using the Lagrange multipliers

λ_{k}^{0}

in the thermostat. The system states at the end of this burst, denoted by y^p(MΔt) for ensemble member p at time MΔt, are used for an estimator of the observables in the distribution ρ₀:

{\hat{C}}_{k}^{(1)} = \frac{1}{P} \sum_{p = 1}^{P} C_{k} (y^{p} (M Δ t)) λ_{0}^{1} \exp (\sum_{i = 1}^{K} (λ_{i}^{0} - λ_{i}^{1}) C_{i} (y^{p} (Δ t))), k = 1, \dots, K,

4.2

where the superscript (1) indicates that it is an estimator for a distribution with Lagrange multipliers

λ_{k}^{1}

. Note that the samples being distributed in ρ₀ and not in π leads to the appearance of the Lagrange multipliers

λ_{k}^{0}

, this is a simple importance sampling procedure. A Newton–Raphson procedure similar to that in §3 is then used to find

λ_{k}^{1}

using the estimator (4.2). The resulting multipliers are used in the second burst of thermostatted simulations, sampling

ρ_{1} (y) = π (y) \exp (- \sum_{k = 1}^{K} λ_{k}^{1} C_{k} (y))

. Now, crucially, this estimator from the ensemble after two simulation bursts is combined with that at the end of the first burst. As both estimators are unbiased, we may use any weighted sum of the two (with the sum of the weights equal to one) to find another estimator. For simplicity, we weight both estimators equally. We write the combination of estimators for the general case after L bursts of the simulation, using the system states at the end of each previous burst, that is y^p(ℓMΔt) for ℓ=1,2,…,L, as an estimator for the densities ρ_ℓ−1 with the Lagrange multipliers

λ_{k}^{ℓ - 1}

. This estimator takes the form

{\hat{C}}_{k}^{(L)} = \frac{1}{L P} \sum_{p = 1}^{P} \sum_{ℓ = 1}^{L} C_{k} (y^{p} (ℓ M Δ t)) λ_{0}^{ℓ - 1} \exp (\sum_{i = 1}^{K} (λ_{i}^{ℓ - 1} - λ_{i}^{L}) C_{j} (y^{p} (ℓ M Δ t))), k = 1, \dots, K .

4.3

As with the non-adaptive scheme, a Newton–Raphson procedure is used to find Lagrange multipliers such that

{\hat{C}}_{k}^{(L)} = c_{k}

with the residuals

r_{k}^{L} = {\hat{C}}_{k}^{(L)} - c_{k} k = 1, \dots, K,

4.4

and the gradients

\frac{\partial r_{k}^{m}}{\partial λ_{j}^{m}} = \frac{1}{m P} \sum_{p = 1}^{P} \sum_{ℓ = 0}^{m - 1} C_{k} (y_{ℓ M}^{p}) C_{j} (y_{ℓ M}^{p}) λ_{0}^{ℓ} \exp (\sum_{i = 1}^{K} (λ_{i}^{l} - λ_{i}^{m}) C_{i} (y_{ℓ M}^{p})), j, k = 1, \dots, K .

4.5

By performing this iteration after each burst, and using the result for the subsequent burst, the Lagrange multipliers are computed ‘on-the-fly’ without the need for constructing a large ensemble in the prior.

Let us consider two limiting cases to justify the approach. (i) By increasing the number of ensemble members P, one recovers the non-adaptive scheme, where the Lagrange multipliers are found based on the initial ensemble after a single burst. On the other hand, (ii) if we have a modest ensemble P, it is necessary that the burst is long enough for the thermostatted dynamics with the adaptive parameters to equilibrate on this time scale. The error due to the transient of the thermostat is removed in the limit of many timesteps per burst, i.e. as $M \to \infty$ . We present numerical verification of the practical use of the scheme in a specific application setting in §6c.

(a) Adaptive algorithm

There are two important modifications to the algorithm described above that are included in the numerical implementation of this method:

— If the Lagrange multipliers change rapidly, the thermostat may require a long time to equilibrate. This requires a larger value for M, increasing the simulation time required before including new samples. The effect is notable at the beginning of a simulation, due to two factors: (i) the small sample size leads to inaccurate expectations for the observables, and (ii) the initial values for Lagrange multipliers may be far from their equilibrium. By limiting the rate of change of the Lagrange multipliers, these problems are circumvented.
— In equations (4.4) and (4.5), all previous values $λ_{k}^{ℓ}$ are included. In a long simulation, this leads to a growing computational demand. By taking only a fixed number (q) of recent steps the cost is reduced. In the case that the initial samples cannot accurately be drawn from the prior, this has the further advantage that these inaccuracies are eventually forgotten.

The resulting method is summarized in algorithm 1.

Inline Graphic

5. Application to reduced modelling of point vortices

Here we apply the least-biased correction method to a system of point vortices on the sphere, which has a Poisson structure and multiple conserved quantities, including total energy and angular momentum and a set of Casimirs.

(a) Point vortex system

Point vortices represent singular solutions to two-dimensional, incompressible fluid flow in which each vortex, positive or negative, is centred at a specified position. Point vortices have been studied extensively in [41–43]. We work on a spherical domain [44–48], viewing the vortex position as a vector in R³ with unit norm. The equations of motion have a Lie–Poisson structure

Γ_{i} {\dot{x}}_{i} = x_{i} \times \nabla_{x_{i}} H i = 1, 2, \dots, M,

5.1

where the Hamiltonian is given by

H = - \sum_{i = 1}^{M} \sum_{j = 1}^{i - 1} \frac{Γ_{i} Γ_{j}}{4 π} \ln (2 - 2 x_{i} \cdot x_{j})

and each Γ_i represents the circulation of a single point vortex.

By introducing $y = {(x_{1}^{T}, x_{2}^{T}, \dots, x_{M}^{T})}^{T}$ , equation (5.1) has the form (2.1) with the block-diagonal structure matrix

B (y) = [\begin{matrix} Γ_{1}^{- 1} {\hat{x}}_{1} \\ Γ_{2}^{- 1} {\hat{x}}_{2} \\ ⋱ \\ Γ_{M}^{- 1} {\hat{x}}_{M} \end{matrix}],

where

{\hat{x}}_{i}

is the 3×3 skew-matrix satisfying

{\hat{x}}_{i} a := x_{i} \times a

, for all a∈R³. The Poisson bracket for the system is given equivalently by

{F, G} = \sum_{i = 1}^{M} \frac{1}{Γ_{i}} \nabla_{x_{i}} F \cdot (x_{i} \times \nabla_{x_{i}} G) or {F, G} = \nabla F {(y)}^{T} B (y) \nabla G (y) .

This Poisson structure generalizes the rigid body Poisson structure and also occurs in ferromagnetic spin lattices [49–51] and elastic rods (e.g. [52]).

The vortex positions are defined in Cartesian coordinates, but initial positions x_i(0) are chosen on the sphere. Because each |x_i| is a Casimir of the Poisson bracket, the vortices remain on the sphere. Furthermore, the rotational symmetry of the sphere gives rise to three Noether momenta, which are expressed by the angular momentum vector

J = \sum_{i = 1}^{M} Γ_{i} x_{i} .

5.2

The GBK thermostat (2.2) is only applicable to divergence free systems ∇⋅f≡0. It is straightforward to check that this condition holds for the spherical point vortex model.

A numerical integrator for the point vortex system is constructed as in Myerscough & Frank [48] based on splitting the differential equations into integrable subproblems (see related ideas in [53,54]). The backward error analysis of symplectic integrators has an analogous development for Poisson systems, implying approximate conservation of the Hamiltonian [8]. The Casimirs (|x_i|≡1) and the momentum $J = \sum Γ_{i} x_{i}$ are exactly preserved by the integrator.

(b) Thermostat perturbation vector

The choice of the perturbation vector field g in (2.2) is flexible, but the vector fields f and g should satisfy a Hörmander condition [36]. It is vital that the Casimirs should be respected by the perturbation vector field. Double-bracket dissipation [55] preserves Casimirs of the original system and is thus an ideal candidate for g:

{\tilde{g}}_{i} (x_{i}) = \sum_{j \neq i} x_{i} \times x_{i} \times \frac{Γ_{j}}{4 π} \frac{x_{j}}{1 - x_{i} \cdot x_{j}} .

5.3

This ‘double-bracket thermostat’ also has the advantage that it can be split into pairwise interactions along with the original dynamics f, simplifying computation. The denominator in (5.3) leads to a stiff differential equation when like-signed vortices approach one another, restricting the step size of an explicit splitting method. For this reason, we use a slightly modified scheme (but still Casimir-preserving) defined by

g_{i} (x_{i}) = \sum_{j \neq i} x_{i} \times x_{i} \times \frac{Γ_{j}}{4 π} x_{j} .

5.4

Details of the numerical integration of these dynamics may be found in appendix A.

The thermostat (2.2) is designed to sample a target density ρ(y)∝e^−A(y) on the phase space of y. The thermostat variable ξ is normally distributed, yielding the extended distribution $ρ \propto e^{- A (y) - \frac{1}{2} ξ^{2}}$ . The perturbation vector field g must additionally ensure that the thermostatted system is ergodic in the target density. Because the target measure is positive for all open sets on the phase space, hypoellipticity of the Fokker–Planck equation associated with (2.2) is sufficient to prove uniqueness of the invariant measure [56]. Hypoellipticity follows from Hörmander’s controllability condition [57]. The condition has been tailored to GBK thermostats in Bajars et al. [36], but it is difficult to verify in practice. Here we verify using simulation that single trajectories have statistics that agree with the target distribution.

(c) Maximum entropy model

To apply the methods proposed in §3 in the setting of a reduced model for point vortices, we use point vortices distributed evenly over the surface of the sphere as the prior π. The use of a thermostat for the bias correction limits us to using observables that are function of the first integrals. Furthermore, we assume the angular momentum has no directional preference, a choice we motivate in the numerical comparisons of §6. In that case, only the magnitude |J| of the momentum vector needs to be considered.

To accurately reproduce statistics from a full model with a moderate number of point vortices, it is necessary to modify the canonical density with a term quadratic in the Hamiltonian, that is, a density of the form $ρ (y) \propto \exp (- β H (y) - γ H {(y)}^{2})$ [16]. Motivated by Dubinkina et al. [16], we choose observations that include linear and quadratic functions of H and |J|².

We consider the following set of observables:

C_{1} = H, C_{2} = | J |^{2}, C_{3} = H^{2}, C_{4} = | J |^{4} and C_{5} = H | J |^{2}

5.5

and denote the corresponding Lagrange multipliers by β_H, β_J, γ_H, γ_J and γ_HJ.

The least-biased density consistent with observations of the $E C_{k}$ is given by

\tilde{ρ} (H) = e^{- β_{H} H - β_{J} | J |^{2} - γ_{H} H^{2} - γ_{J} | J |^{4} - γ_{H J} H | J |^{2}},

5.6

where the Lagrange multipliers are to be found via the procedure detailed in §3a or that of §4.

6. Numerical comparison

Here we apply algorithm 1 to a reduced model of point vortices similar to the configuration used in Bühler [58] and Dubinkina et al. [16]. We distinguish between three models. The full model consists of a system (5.1) of M_full=288 point vortices, of which eight strong vortices of circulation Γ_j=±1 and 280 weak vortices of circulation $Γ_{j} = \pm \frac{1}{5}$ . Both strong and weak classes comprised equal numbers of positively and negatively oriented point vortices. The reduced model consists of (5.1) with just M=8 strong vortices. Finally, the corrected model consists of a thermostatted system (2.2) with unperturbed vector field f given by (5.1) for M=8 strong vortices, perturbation vector field g given by (5.4), and equilibrium measure defined by the least-biased density (5.6). Additionally, we compare with Metropolis–Hastings samples from the least-biased density (5.6) to help distinguish between errors incurred due to the maximum-entropy model and those due to sampling bias of the (discretized) thermostatted dynamics.

We run five long simulations of the full model with total energies chosen from the set H_full∈{−2,−1,0,1,2}. The total angular momentum vector is fixed at J_full=0 such that there is no directional preference for the momentum of the reduced model embedded in the full model. For each run, we determine the time averages of the observables (5.5) for the subset of strong vortices. When computing the Hamiltonian H of this subsystem, we include only the internal coupling between strong vortices. The time averages are tabulated in table 1.

Table 1.Full model observations and (in parentheses) corrected values of first integrals.
Collapse
	〈H〉	〈\|J\|²〉	〈H²〉	〈\|J\|⁴〉	〈H\|J\|²〉
H_full=−2	−0.33 (−0.38)	4.59 (4.45)	0.22 (0.23)	−0.63 (−0.98)	34.58 (31.45)
H_full=−1	−0.11 (−0.18)	4.78 (4.68)	0.10 (0.12)	0.38 (−0.01)	37.55 (35.88)
H_full=0	0.02 (−0.04)	4.63 (4.56)	0.08 (0.08)	0.90 (0.60)	35.26 (34.30)
H_full=1	0.17 (0.15)	4.74 (4.75)	0.13 (0.12)	1.73 (1.61)	37.76 (37.44)
H_full=2	0.31 (0.28)	4.87 (5.00)	0.22 (0.21)	2.49 (2.46)	39.26 (41.74)

Using these averages, we compute the Lagrange multipliers using the algorithm described in §4 with prior distribution π taken to be uniform on the sphere. The Lagrange multipliers are also recorded in table 2. The magnitude of γ_{H,J,HJ} indicates that including these observables impacts the resulting density. In other words, all observables add information to the least-biased density, with the possible exception of the lowest energy case with H_full=−2.

Table 2.Lagrange multipliers for each energy level.
Collapse
	β_H	β_J	γ_H	γ_J	γ_HJ
H_full=−2	5.98	−0.20	0.69	0.41×10⁻³	−0.04
H_full=−1	2.89	−0.03	2.67	9.77×10⁻³	−0.33
H_full=0	−0.76	0.20	3.38	9.97×10⁻³	−0.37
H_full=1	−3.54	0.37	4.29	15.31×10⁻³	−0.54
H_full=2	−6.42	0.53	4.45	14.05×10⁻³	−0.51

We then run simulations of the corrected model using the computed parameters. Table 1 also records expectations from the thermostat-corrected model.

By analogy with canonical statistical mechanics, we may think of the weak vortices that are ignored in the reduced model as forming a reservoir with which our reduced model exchanges energy and angular momentum. Experience with canonical statistical mechanics of point vortices in the plane [16,58] suggests that for small reservoir sizes, the canonical distribution must be modified with higher order terms to agree with the full system statistics. Table 3 gives the multipliers for different numbers M_full of vortices in the full system, confirming that those corresponding to γ_H, γ_J and γ_HJ are more significant for smaller M_full; as the number of weak vortices grows, the strong vortex system approaches the canonical density.

Table 3.Lagrange multipliers as a function of M_full, all for H_full=0.
Collapse
	β_H	β_J	γ_H	γ_J	γ_HJ
M_full=36	−1.51	1.35	27.75	117.15×10⁻³	−3.07
M_full=72	−4.27	0.82	8.71	37.79×10⁻³	−1.12
M_full=144	−0.97	0.32	6.70	25.80×10⁻³	−0.82
M_full=288	−0.76	0.20	3.38	9.97×10⁻³	−0.37
M_full=576	−1.09	0.13	0.87	3.08×10⁻³	−0.10

Remark 6.1

The form of the distribution (5.6) implies there is a non-zero probability for any configuration of the strong vortex system with finite energy and momentum. When considering the strong vortex system in contact with a reservoir of weak vortices, its energy is bounded by the energy in the total system plus or minus how much energy can be stored in the reservoir and in the interactions between strong and weak vortices. For the momentum, a similar condition holds, but here there are no interaction terms. If it is not possible for the reservoir to add or subtract sufficient energy from the reservoir, it will not be possible for the subsystem of strong vortices to attain all possible configurations with non-zero probability.

For the energy, it is readily verified that as long as the reservoir contains three vortices not all of equal sign, the reservoir may supply or remove an arbitrary amount of energy, as close of approaches of two like-signed (respectively, opposite-signed) vortices produces arbitrarily large-positive (respectively, negative) energies.

For the momentum, the result is more interesting. The system of M strong vortices, all with strength ±Γ_strong, may exist in configurations with an angular momentum $sup | J_{red .} | = M Γ_{strong}$ . Note the supremum, as the configuration with exactly $sup | J_{red .} | = M Γ_{strong}$ will have infinite energy. For the reservoir, it holds that $sup | J_{weak} | \leq (N - M) Γ_{weak}$ . For the the reservoir to supply sufficient angular momentum, it is necessary that

M Γ_{strong} \leq (M_{full} - M) Γ_{weak} \Leftrightarrow \frac{M_{full} - M}{M} \geq \frac{Γ_{strong}}{Γ_{weak}} .

In the thermal bath simulations discussed in this section M=8 and Γ_strong/Γ_weak=5, this means M_full should satisfy M_full≥48. The smallest system considered (M_full=36) does not, explaining its eccentric parameter values in table 3.

(a) Equilibrium results

In this section, we compare statistical properties of the corrected model with those of the full and reduced models. In figures 1–3, we show histograms of a number of solution features for the eight vortex model: the distributions of H and |J|, as well as typical distances between like- and opposite-signed vortices, a metric also used by Bühler [58]. In each histogram, the statistics corresponding to the strong vortices in the full model, the reduced model, thermostat-corrected reduced model and Metropolis–Hastings samples are displayed. Figures 1–3 correspond to approximate total energies H_full≈−2, 0 and 2, respectively.

Figure 1. Histograms for H_full≈−2. (a,b) Compare strong vortex energy and angular momentum, respectively. Panels (c,d) Compares the distance between like (resp. opposite) signed strong vortices. The parameters are specified in the text. (Online version in colour.)

Figure 2. Histograms for H_full≈0. Panel layouts same as figure 1. (Online version in colour.)

Figure 3. Histograms for H_full≈2. Panel layouts same as figure 1. (Online version in colour.)

The full and reduced model simulations are performed with a timestep of 5×10⁻³ and run up to T=5×10⁴, taking 10⁵ samples spaced evenly in time. For the Metropolis–Hastings method, we use 10⁶ samples. The same figures also show results from the thermostatted system (dashed-dotted lines), run with a timestep of 10⁻³ up to T=10⁶, taking 10⁶ samples. The parameters in (2.2) were set to be ϵ=10 and γ=0.1. These results confirm that the thermostatted system samples the least-biased density closely. The complete simulation trajectories are available in the electronic supplementary material.

The reduced model is Hamiltonian and the Poisson integrator ensures that the energy is conserved with a standard deviation of order 10⁻³ and the angular momentum constant to machine precision. Both cases correspond to approximate delta-distributions in the upper histograms in figures 1–3. Note that due to the high skewness of the distribution for |J|, the observed mean differs significantly from the median and mode, implying some ambiguity in choosing the angular momentum for an appropriate initial condition for the reduced model.

A simple Hamiltonian-reduced model is incapable of sampling the energy and angular momentum spectra, as these are first integrals, hence the reduced model shows significant bias in statistics such as vortex separation. On the other hand, the thermostat-corrected model faithfully samples the least-biased probability density, as indicated by the good agreement in the histograms of the corrected model and Metropolis–Hastings samples. The least-biased density approximates the strong-vortex statistics well, particularly in the negative to moderate total energy regime. At large positive total energies, the strong vortex energy and angular momentum distributions are still well represented by the least-biased PDF, but some bias in the vortex separations can be observed. The closeness of the thermostat results to those from the Metropolis–Hastings sampling indicate the error lies in the choice of least-biased density, not in the thermostat sampling.

(b) Dynamic consistency

The results in the previous section confirm that the thermostatted simulations lead to equilibrium distributions of observables H and |J| similar to those of the full system. In this section, we address the issue of the degree to which our equilibrium correction mechanism disturbs dynamics, as encoded in autocorrelation functions and diffusivity. Diffusivity was considered by Chavanis [59] for a system of identical point vortices and by Cotter & Pavliotis [60] for a wide array of problems with scale separation. We emphasize that the values of the thermostat parameters ε and γ have no impact on the equilibrium statistics presented in the previous section, and only affect the rate at which the least-biased PDF is sampled. Parameters with a larger deviation from the unperturbed dynamics lead to faster equilibration of the distribution, giving rise to a modelling choice.

(i) Autocorrelation functions

Given a sequence of L equally spaced observation times t_i∈[0,T] for i∈[0,L], and the values of the relevant observable (in our case vortex position) u_i=u(t_i) at those times, the discrete autocorrelation function is defined by

ν_{i}^{u} = \frac{1}{L - i} \sum_{j = i}^{L} u (t_{j}) u (t_{j - i}) .

A normalized autocorrelation function

{\hat{ν}}^{u}

is given by dividing each

ν_{i}^{u}

ν_{0}^{u}

, i.e.

{\hat{ν}}_{i}^{u} = ν_{i}^{u} / ν_{0}^{u}

We average the autocorrelation function over all 3M (strong) vortex coordinates {x_m,y_m,z_m}_m=1…M. Three symmetries in the problem justify this averaging: the vortex numbering is arbitrary; the choice of reference frame is arbitrary and the sign of the vortices appears in the dynamics as a reversal of time, to which the autocorrelation is insensitive. Additionally, the observables H and |J| are isotropic.

Furthermore, we ensure that the phase space is well sampled by averaging the autocorrelation functions over an ensemble of P solutions. We then find the average autocorrelation function

ν_{i} = \frac{1}{3 M P (L - i)} \sum_{p = 1}^{P} \sum_{m = 1}^{M} \sum_{j = i}^{L} x_{m}^{p} (t_{j}) x_{m}^{p} (t_{j - i}) + y_{m}^{p} (t_{j}) y_{m}^{p} (t_{j - i}) + z_{m}^{p} (t_{j}) z_{m}^{p} (t_{j - i}),

6.1

where a superscript p represents the solution from ensemble member p. The normalized average autocorrelation function is given by

{\hat{ν}}_{i} = \frac{1}{M P (L - i)} \sum_{p = 1}^{P} \sum_{m = 1}^{M} \sum_{j = i}^{L} x_{m}^{p} (t_{j}) \cdot x_{m}^{p} (t_{j - i}),

6.2

where we have used that the Casimirs C_i=x_i(t)⋅x_i(t)=1 ∀ i,t are conserved exactly, also in the discretized equations.

In figure 4, we compare autocorrelation functions for the strong vortices in the full and reduced models as well as for the thermostat-corrected model over a range of parameters ε and γ. The thick solid black line represents the result for an (unthermostatted) system in contact with 280 weak ( $Γ_{B} = \pm \frac{1}{5}$ ) vortices, with a total energy H_full=0. The results present the average over an ensemble of 1000 runs. For each simulation, the initial placement of each strong vortex was taken uniformly over the sphere and the weak vortices were placed such that the full system satisfied H_full=0 and J_full=0. The thick dashed black line represents the results for an ensemble of simulations of the isolated system, with everything else unchanged. The other lines represent results for thermostatted simulations using the Lagrange multipliers as given in table 2 for the case of H=0. As expected, when the thermostat is weak the autocorrelation functions are similar to those of the isolated system. The stronger thermostat parameters presented here match the autocorrelation well for a short-time period, but all result in excessive decorrelation for longer lag times. In the following section on the diffusivity constant, we discuss the performance of different parameter values in more detail.

(ii) Diffusivity

For general multiscale dynamical systems with a separation of slow and fast dynamics, it is often desirable to model fast forces by a diffusion process, resulting in stochastic differential equation of the form [61,62]

d X = f (X) d t + K (X) d W,

where f represents the slow dynamics, W is a Wiener process and K is the diffusivity. The value of the diffusivity K can be estimated by sampling solutions to the original, multiscale, problem and applying Kubo’s formula

K_{Δ τ} = \frac{⟨ Δ X Δ X ⟩}{2 Δ τ},

where ΔX represents displacement during the sampling interval Δτ. Choosing the correct sampling interval is a notorious problem; for a comparison, see [61].

If we take the average diffusivity for each vortex coordinate, we find

K_{Δ τ} = \frac{1}{6 M Δ τ} \sum_{m = 1}^{M} ⟨ Δ x_{m} \cdot Δ x_{m} ⟩ .

We assume the observations are given at the same times t_i as before and that the sampling time is an integer multiple of the observation interval, i.e. Δτ=k(T/L). With an ensemble of P simulations the diffusivity estimator would then be

\begin{aligned} K_{Δ τ} & = \frac{1}{6 M P Δ τ} \sum_{p = 1}^{P} \sum_{m = 1}^{M} (x_{m}^{p} (Δ τ) - x_{m}^{p} (0)) \cdot (x_{m}^{p} (Δ τ) - x_{m}^{p} (0)) \\ = \frac{1}{6 M P Δ τ} \sum_{p = 1}^{P} \sum_{m = 1}^{M} 2 - 2 x_{m}^{p} (0) \cdot x_{m}^{p} (Δ τ) \\ = \frac{1}{3 Δ τ} - \frac{1}{3 M P Δ τ} \sum_{p = 1}^{P} \sum_{m = 1}^{M} x_{m}^{p} (0) \cdot x_{m}^{p} (Δ τ), \end{aligned}

where again a superscript p denotes the solution from ensemble member p. Averaging over all time series data yields the estimator:

K_{Δ τ} = \frac{1}{3 Δ τ} - \frac{1}{3 M P (L - k) Δ τ} \sum_{p = 1}^{P} \sum_{j = 1}^{L - k} \sum_{m = 1}^{M} x_{m}^{p} (t_{j}) \cdot x_{m}^{p} (t_{j} + Δ τ) = \frac{1 - {\hat{ν}}_{k}}{3 Δ τ},

with Δτ=k(T/L)=t_k. This shift-averaged estimator incorporates all possible time intervals of length Δτ available from the data; it is shown by Cotter & Pavliotis [60] to improve the quality of the estimator.

Figure 5 presents the diffusivity constant as a function of the interval time Δτ for different parameter values for ε and γ. The results are taken from the same simulations used for the autocorrelation functions, as described in the section.

For ε small, the thermostat perturbation is weak, and both autocorrelation functions and diffusivity approach those of the reduced model with constant H, J. Also, the autocorrelations are insensitive to the parameter γ in this regime. For values of ε<1, the autocorrelations and diffusivities become indistinguishable from those of the isolated model. Hence even though the dynamics samples the least-biased density on long time scales, its short-time dynamics is similar to an unperturbed model. For moderate ε, dependence on γ becomes more pronounced, and a diffusivity closer to that of the full model can be achieved. For even larger values of ε, the diffusivity becomes much more sensitive to the value of γ, as indicated in figure 5c.

Figure 5d has been included to illustrate two important properties by plotting the diffusion parameter against the sampling interval on a log–log scale. First, for large sampling interval the estimator shows an inverse linear tendency. This corresponds simply to the decorrelation of the vortex dynamics. Second, as the sampling interval goes to zero, the estimator of the diffusivity constant shows polynomial behaviour. This is in agreement with known results for the GBK thermostat [63] and is an improvement on Langevin thermostats, which tend to a constant value for short sampling intervals, in disagreement with a deterministic reference [30]. See also the discussion in §2.

(c) Adaptive determination of multipliers

In this section, we apply the scheme of §4 to the point vortex system. Consider the same reduced point vortex model of eight vortices with Γ=±1. We perform a thermostatted simulation of this system, which we feed with observation data taken from a simulation of the 288-vortex system including the weak vortices that form the thermal bath. The observation from simulation with the thermal bath may be replaced by observation from another source–if available–such as empirical observations or a prediction of the trend for the observables. In particular, we use the observed data in three different ways for computing the Lagrange multipliers that determine the equilibrium distribution used by the thermostat:

(i) The first approach mirrors the non-adaptive method of §3. We use the long time averages from equilibrated thermal bath simulation. The equilibrium distribution varies during the thermostatted simulation, but during the whole simulation it converges to the same equilibrium density as before (like those denoted in table 2).
(ii) The second approach does not assume long time averages are available a priori. Instead, the running mean from the thermal bath simulation is used. This means the iterative procedure (4.4) drives the system towards a varying target during simulation. Eventually, the target approaches the long time average used above.
(iii) In the final implementation, we relax the notion of statistical equilibrium for the strong system and consider the statistical state to be slowly varying. This is implemented by using a time-local averaged observations from the thermal bath simulations. As a result, there is no stationary equilibrium density.

For testing the adaptive scheme, we only consider the thermal bath simulations where H_full=0. Also, we consider just two observables C₁=H and C₂=|J|², associated with the Lagrange multipliers β_H and β_J. Including more observables is possible, but the Newton–Raphson iteration is sensitive to noise in the estimators, requiring a larger ensemble. We initialize with an ensemble of P=100 initial conditions drawn from the uniform prior, applying algorithm §4a The timestep in the thermostatted simulation is chosen as 1×10⁻² and the method described in §4 for updating the Lagrange multipliers is applied every time unit, i.e. M=100. Between subsequent updates of the multipliers, the maximum difference is limited by |Δλ_k|≤0.1. When using equilibrium statistics, this limit only affects the beginning of the simulation, when the small sample size used leads to a large variance in the estimators.

We present the results for the three different uses of the observed data. The complete simulation trajectories are available in the electronic supplementary material.

(i) In figure 6, the long time mean is taken and used throughout.
(ii) In figure 7, the running mean is used. This reflects the situation where we have no a priori knowledge of the observations, and are continuously feeding new real-time data into the simulation.
(iii) In figure 8, a time-localized average of the observable is used. The averaging has a time-scale of a 100 time units.

Figure 6. Results when using *long-time mean* observations as a target while adaptively determining the Lagrange multipliers. Target observations for Hamiltonian (a) and momentum magnitude (b) are overlaid with the instantaneous ensemble mean (black dotted line)and the running ensemble mean (red solid line) from simulation. (c) Lagrange multipliers, the red dots indicate time steps at which their rate of change was limited. (Online version in colour.)

Figure 7. Results when using *running mean* observations as a target while adaptively determining the Lagrange multipliers. (Online version in colour.)

Figure 8. Results when using *time-local averaged* observations as a target while adaptively determining the Lagrange multipliers. (Online version in colour.)

In each of the figures, panel (a) shows the target mean for the energy solid cyan (light), the instantaneous ensemble mean in black dots and the running mean of the ensemble in dashed red (dark). Panel (b) shows the same means for the momentum magnitude and (c) shows the Lagrange multipliers for the energy (cyan, light) and momentum (black). Superimposed on the Lagrange multiplier for the energy is an indicator (red, darker dots) showing when the rate of change in the estimate for the Lagrange multiplier has been limited.

When using either a long time mean observation or a running mean observation, the simulation results tend towards the correct long-time averages. When using time-local averages the simulation averages appear to tend towards a similar value. In all three cases, the instantaneous ensemble mean remains close to the (moving) target for both energy and momentum. This is especially notable for the third case, where the target varies over time, but the simulation ensemble mean follows closely, with only a little lag.

The inaccuracies during the first approximately 100 time units indicate that the prior does not match the observed state well. This results in the (negative) growth of β_H being limited briefly at the beginning of each simulation. Subsequently, both Lagrange multipliers appear to oscillate irregularly about some mean value for the first two cases. In the case of a shifting target, the Lagrange multipliers vary in time more erratically, resulting in the limiter being active for a few brief periods of the simulation.

7. Conclusion

We have proposed a thermostat-based method for perturbing trajectories of numerical simulations to correct for stationary or time-dependent thermodynamic observations. We have also described an adaptive procedure that uses the thermostatted simulation data in finding the least-biased density, removing the need to compute the density a priori.

Our numerical experiments with the vortex model confirm that the distributions of the observed energy H and angular momentum magnitude |J| can be well approximated using the thermostat technique. Other equilibrium metrics such as the distribution of distances between like- and opposite-signed vortices are also in agreement across a range of total energy values of the full system, although some discrepancies occur at large positive energies.

We also investigated the degree to which correction of trajectories for expectations may affect dynamical information in the form of autocorrelation functions and diffusivity. By decreasing the perturbation parameter ε of the thermostat, the autocorrelation functions of the unperturbed, reduced system may be precisely recovered. As ε is increased, one may increase the diffusivity to values that agree with the full system. This is consistent with results reported in Frank & Gottwald [63], in the context of molecular dynamics where it was shown that the GBK thermostat used here approaches Langevin dynamics in the limit of large stochastic forcing.

Data accessibility

The code for thermostated simulation of the point vortex system can be found at https://www.dropbox.com/sh/ekhj4ohtetj35ri/AABNMpHki63cn16hpnaQ.OYa?dl=0.

Authors' contributions

The authors have contributed equally to this work.

Competing interests

We declare we have no competing interests.

Funding

The work of K.W.M. was supported by a grant financed by the Netherlands Organisation for Scientific Research (NWO), as part of research programme no. 613.001.009. B.L. was supported by grant no. EP/K039512/1, grant no. SI2-CHE: ExTASY (Extensible Tools for Advanced Sampling and analYsis) from the Engineering and Physical Sciences Research Council (UK) and a faculty fellowship from the Alan Turing Institute (ATI).

Footnotes

1 The latter condition is automatic for systems (2.1) with constant B. Strictly speaking, the approach described here is applicable to any system with divergence-free vector field ∇⋅f≡0 possessing one or more first integrals.

Appendix A. Time integration of the thermostatted system

In this appendix, we detail the method of integration for the perturbation dynamics and thermostat variable.

The unperturbed dynamics are formed by the combination of pairwise interactions:

f (y) = B (y) \nabla_{y} H (y) = - B (y) \sum_{i = 1}^{M} \sum_{j = 1}^{i - 1} \nabla_{y} \frac{Γ_{i} Γ_{j}}{4 π} \ln (2 - 2 x_{i} \cdot x_{j}) = \sum_{i < j} f_{i j} .

A 1

This fact is exploited in our earlier work [48] to construct a splitting method for integration of the point vortex system with favourable (geometric) properties. The modified double bracket thermostat is also a combination of pairwise vortex interactions. Hence it is natural to apply the same splitting to g as was applied to f. In fact, the pairwise thermostat interaction can be executed along with the unperturbed dynamics. Let us denote by

ϕ_{Δ t}^{i, j}

the time Δt flow map associated with f_ij. Similarly, we define by

γ_{Δ t}^{i, j}

the time Δt flow map associated with the perturbation dynamics of a vortex pair (i,j). The flow map of the dynamics of the thermostat variable ξ is represented by χ_Δt. A symmetric composition of these flows is given by

Φ_{Δ t} = \prod_{(i, j) \in C} (ϕ_{Δ t / 2}^{i, j} \circ γ_{Δ t / 2}^{i, j}) \circ χ_{Δ t} \circ \prod_{(i, j) \in C^{*}} (γ_{Δ t / 2}^{i, j} \circ ϕ_{Δ t / 2}^{i, j}),

where C is an ordering of pairs and C* is the reverse ordering. While advancing ξ by χ_Δt, it is assumed that y is fixed and, therefore, h(y) is a constant. This means the dynamics of ξ(t) is just an Ornstein–Uhlenbeck process with mean h/γ and unit variance. The exact solution is given by

χ_{Δ t} ξ_{0} = ξ_{0} e^{- γ Δ t} + \frac{h}{γ} (1 - e^{- γ Δ t}) + e^{- γ Δ t} W (e^{2 γ Δ t} - 1) .

The integration of the perturbed dynamics for a vortex pair (i,j) means integrating the system of ordinary differential equations given by

\begin{aligned} {\dot{x}}_{i} & = \frac{Γ_{j}}{4 π} x_{i} \times x_{i} \times x_{j} = \frac{Γ_{i} + Γ_{j}}{4 π} x_{i} \times x_{i} \times x_{j} + \frac{Γ_{i} - Γ_{j}}{4 π} x_{i} \times x_{i} \times x_{j}, \\ {\dot{x}}_{j} & = \frac{Γ_{i}}{4 π} x_{j} \times x_{j} \times x_{i} = \frac{Γ_{i} + Γ_{j}}{4 π} x_{j} \times x_{j} \times x_{i} - \frac{Γ_{i} - Γ_{j}}{4 π} x_{j} \times x_{j} \times x_{i} . \end{aligned}

The equations have been split into a symmetric and an antisymmetric part. A symmetric composition may once more be applied to integrate the two parts. In our case, however, the thermostatted system consists only of vortices with strength ±Γ_strong and thus each vortex pair interaction is either fully symmetric or fully antisymmetric.

By introducing α_S=(Γ_i+Γ_j)/4π, the symmetric part may be written as

{\dot{x}}_{i} = α_{S} x_{i} \times x_{i} \times x_{j} = α_{S} (x_{j} \times x_{i}) \times x_{i}

A 2

and

{\dot{x}}_{j} = α_{S} x_{j} \times x_{j} \times x_{i} = - α_{S} (x_{j} \times x_{i}) \times x_{j} .

A 3

These dynamics are symmetric with respect to the plane equidistant between the two vortices. This is because the dynamics are rotations in opposite direction about the same vector x_j×x_i and this vector must lie in said plane. Furthermore, the dynamics of (7.2)–(7.3) does not allow the vortex pair to pass through the position where they are exactly opposed, as in this case their cross product is zero. All in all, this means the final position of the two vortices is uniquely determined by their chord distance.

The change in the distance between the two vortices can be represented by the change of their inner product

\begin{aligned} \frac{\partial}{\partial t} (x_{i} \cdot x_{j}) & = {\dot{x}}_{i} \cdot x_{j} + x_{i} \cdot {\dot{x}}_{j} \\ = 2 α_{S} ({(x_{i} \cdot x_{j})}^{2} - 1) . \end{aligned}

The solution to this differential equation is given by

x_{i} \cdot x_{j} |_{t = Δ t} = \tanh (2 α Δ t - artanh (x_{i} \cdot x_{j}) |_{t = 0}) .

In the antisymmetric part, we introduce α_A=(Γ_i−Γ_j)/4π to write

\begin{aligned} {\dot{x}}_{i} & = α_{A} x_{i} \times x_{i} \times x_{j} \\ and {\dot{x}}_{j} & = - α_{A} x_{j} \times x_{j} \times x_{i} . \end{aligned}

By rearranging the order of the cross product, we achieve

{\dot{x}}_{i} = α_{A} (x_{j} \times x_{i}) \times x_{i} = {\hat{a}}_{A} x_{i}

A 4

and

{\dot{x}}_{j} = α_{A} (x_{j} \times x_{i}) \times x_{j} = {\hat{a}}_{A} x_{j} .

A 5

The vector a_A implicity defined by (7.4) and (7.5) can simply be shown to be constant under the antisymmetric flow. This means the dynamics of the antisymmetric part may be integrated by using Rodigues’s formula.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3659462.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

VIEW FULL TEXT DOWNLOAD PDF

This Issue

Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences cover image

January 2017

Volume 473Issue 2197

Article Information

DOI:https://doi.org/10.1098/rspa.2016.0730
PubMed:28265197
Published by:Royal Society
Online ISSN:1471-2946

History:

Manuscript received25/09/2016
Manuscript accepted09/12/2016
Published online01/01/2017
Published in print31/01/2017

License:

Citations and impact

Keywords

Subjects

Your password must have 8 characters or more and contain 3 of the following:

Observation-based correction of dynamical models using thermostats

Abstract