Solving Equations of Motion by Using Monte Carlo Metropolis: Novel Method Via Random Paths Sampling and the Maximum Caliber Principle

Abstract

A permanent challenge in physics and other disciplines is to solve Euler–Lagrange equations. Thereby, a beneficial investigation is to continue searching for new procedures to perform this task. A novel Monte Carlo Metropolis framework is presented for solving the equations of motion in Lagrangian systems. The implementation lies in sampling the path space with a probability functional obtained by using the maximum caliber principle. Free particle and harmonic oscillator problems are numerically implemented by sampling the path space for a given action by using Monte Carlo simulations. The average path converges to the solution of the equation of motion from classical mechanics, analogously as a canonical system is sampled for a given energy by computing the average state, finding the least energy state. Thus, this procedure can be general enough to solve other differential equations in physics and a useful tool to calculate the time-dependent properties of dynamical systems in order to understand the non-equilibrium behavior of statistical mechanical systems.

1. Introduction

The main objective of this work is to show a new computational framework based on Monte Carlo Metropolis, for the study of dynamical systems that are described by a Lagrangian, being a first approach for the understanding of non-equilibrium statistical mechanics (NESM) by using constraints, as performed in statistical mechanics.

Here, we propose a technique capable of simulating deterministic, dynamical systems through a stochastic formulation. This technique is based on sampling a statistical ensemble of paths defined as having the maximum path entropy (also known as the caliber) available under imposed time-dependent constraints. In Section 2, we introduce this approach, which is known as the maximum caliber principle (MaxCal) [1], a generalization of Jaynes’ principle of maximum entropy [2,3,4].

The principle of maximum entropy (MaxEnt) is a systematic method for constructing the simplest, most unbiased probability distribution function under given constraints, a conceptual generalization of Gibbs’ method of ensembles in statistical mechanics. The complete generality of the principle makes it widely used in several areas of science, such as astronomy, ecology, biology, quantitative finance, image processing, electronics, and physics, among others. According to Jaynes, choosing a candidate probability distribution by maximizing its entropy is a rule of inferential reasoning far beyond its original application in physics, which makes this rule a powerful tool for creating models in any context.

Monte Carlo Metropolis (MCM) is a rejection technique that can be used for sampling any probability distribution [5,6]. This probability distribution is usually constructed by using MaxEnt and optimized using the simulated annealing algorithm [7]. MCM and simulated annealing are used for the understanding of different systems in physics [8,9] such as spin models, material simulations, and thermodynamics, among others. In Section 3, an analogous methodology is presented, by assigning probabilities to paths instead of probabilities of states, allowing the minimization of functional quantities such as the classical action instead of the energy.

It is possible to directly use MaxEnt instead of MaxCal as an inference methodology for dynamical systems by constraining time-dependent functions, and this has been used to understand Newtonian dynamics [10] and the Schrödinger equation [11] in terms of traditional MaxEnt. A new approach has been explored recently for recovering frameworks for dynamical systems by using MaxCal and the path space [12,13,14,15]. In Section 4 and Section 5, we show a numerical implementation of this new approach, exposing the capability for solving complex problems in NESM. In summary, this constitutes a novel approach to the study of dynamical systems via Monte Carlo simulation.

2. Creating a Path Ensemble: The Maximum Caliber Principle

The maximum caliber principle allows the definition of a unique path ensemble given prior information and a number of dynamical constraints [1,14,16].

MaxCal is similar to the maximum entropy principle (MaxEnt) [13,15,17,18,19], but defined over the path space, allowing defining a probability functional $P\left[x\right]$ as follows. Consider an N-dimensional path x on the path space $\mathbb{X}$. In order to construct a probability functional for each path $P\left[x\right]$, given an initial probability (prior) $P_0\left[x\right]$ and an arbitrary constraint:

$$〈A\left[x\right]〉=a,$$

(1)

the caliber (or path entropy):

$$S[P_0\to P]=-{\int }_{\mathbb{X}}{D}_{x}P\left[x\right]ln\frac{P\left[x\right]}{P_0\left[x\right]},$$

(2)

must be maximized under the constraint in Equation (1) and the requirement that probability be normalized,

$${\int }_{\mathbb{X}}{D}_{x}P\left[x\right]=1.$$

The operation ${\int }_{\mathbb{X}}{D}_x$ is a path integral, where $D_x$ is a differential in the path space $\mathbb{X}$. Then, the probability functional obtained is,

$$P\left[x\right|\beta ]=\frac{1}{Z\left(\beta \right)}{P}_0\left[x\right]exp(-\beta A\left[x\right]),$$

(3)

where $Z\left(\beta \right)$ is the partition function and $\beta$ is the Lagrange multiplier, given by the constraint equation:

$$-\frac{\partial }{\partial \beta }lnZ\left(\beta \right)=a.$$

(4)

In this formalism, analogous to Gibbs’ method, $\beta$ is identified as the inverse of the temperature ($\beta =1/k_{B}T$) similar to the canonical ensemble of statistical mechanics. Here, the expected value of an arbitrary functional $F\left[x\right]$ is given by:

$${〈F\left[x\right]〉}_{\beta }=\frac{1}{Z\left(\beta \right)}{\int }_{\mathbb{X}}{D}_{x}exp(-\beta A\left[x\right])F\left[x\right],$$

(5)

but, perhaps more importantly, the expectation value of a function over time can be defined similarly as:

$${〈f(\dot{x},x,t)〉}_{\beta ,t}=\frac{1}{Z\left(\beta \right)}{\int }_{\mathbb{X}}{D}_{x}exp(-\beta A\left[x\right])f(\dot{x},x,t).$$

(6)

This last relation shows that MaxCal can be used to understand the macroscopic properties of time-dependent systems, which are the main elements in NESM.

3. Least Action Principle and Most Probable Path

From classical mechanics, it is well known that the path followed by a particle under a potential $V(x;t)$ under the boundary conditions $x\left(0\right)=0$ and $x\left(T\right)=x_f$ is the one given by the least action principle [20,21], which in practice leads to an equation of motion describing the evolution of the particle from zero to T.

The classical action is a functional defined as $A\left[x\right]={\int }_0^{T}dtL(\dot{x},x;t)$, where $L(\dot{x},x;t)$ is called the Lagrangian of the system, which for classical systems is:

$$L(\dot{x},x;t)=\frac{m{\dot{x}}^2}{2}-V(x;t).$$

For a MaxCal framework where the classical action is constrained, following an analogous treatment to the constraint in Equation (1), the probability functional is of the form given in Equation (3). The most probable path can be obtained by finding the extrema of the functional $P\left[x\right]$, by solving the equation $\frac{\delta P\left[x\right]}{\delta x\left(t^{\prime }\right)}=0$. This is because the exponential function in Equation (3) is convex and monotonically increasing, and so, maximum probability is equivalent to imposing that x should be an extremum for the argument of the exponential,

$$-\beta \frac{\delta A\left[x\right]}{\delta x\left(t^{\prime }\right)}=0.$$

(7)

If the Lagrange multiplier is positive ($\beta >0$), the requirement is that the action is actually a minimum. This equation is precisely the Euler–Lagrange equation of motion for the Lagrangian [20,22]. In summary, the most probable path and the least action path in a MaxCal framework are the same. According to this, sampling trajectories from the probability distribution in Equation (3) using Monte Carlo methods and computing the averages of quantities should converge to a description of the dynamical properties of a classical system evolving in time.

Another important consequence of the use of MaxCal and the form of the probability is that the most probable path in general coincides with the average path according to the central limit theorem.

4. Computational Method

In order to define the elements on the path space $\mathbb{X}$, for an N-dimensional path x, it is always possible to write it in an orthonormal basis ${\varphi }_i$ of the form:

$$x\left(t\right)=\sum _i^N{a}_i{\varphi }_i\left(t\right)=x(t;\mathit{a}).$$

(8)

Then, by changing the parameters $a_i$, it is possible to map the entire space of paths x. In other words, there is a one-to-one correspondence between an arbitrary path x and its parameter vector $\mathit{a}$, so the action becomes a function of $\mathit{a}$, namely $\mathcal{A}\left(\mathit{a}\right):=A\left[x\right(·;\mathit{a}\left)\right]$, and the problem of path sampling reduces to ordinary sampling [23] of N-dimensional states $\mathit{a}$,

$$P\left(\mathit{a}\right|\lambda )=\frac{1}{Z\left(\lambda \right)}exp(-\lambda \mathcal{A}\left(\mathit{a}\right)).$$

(9)

The choice of the basis functions ${\varphi }_i$ is, in principle, arbitrary. However, a convenient choice can be made related to the particular conditions of the problem to be solved. In this case, the target is the study of classical dynamical systems by using the MCM where most of the problems have well-defined boundary conditions; therefore, it is important to find a basis set in which one can easily generate paths in the desired path space holding the required, fixed boundary conditions. For this reason, we considered the Bézier curves as a basis set.

Bézier curves are defined by “control points” $c_i$, where the first $c_0$ and the last $c_n$ control point determine the boundary conditions of the curve, allowing the mapping of the path space with well-defined boundary conditions.

For an N-dimensional path x with boundary conditions $x\left(t_0\right)=c_0$ and $x\left(t_f\right)=c_n$, a Bézier curve is defined of the form:

$$x\left(t\right)=\sum _{i=0}^N{c}_i{B}_i(t;n),$$

(10)

where the basis functions $B_i(t;n)$ are the Bernstein polynomials,

$$B_i(t;n)=\left(\genfrac{}{}{0pt}{}{n}{i}\right)\frac{{(t-t_0)}^i{(t_f-t)}^{n-i}}{{(t_f-t_0)}^n}.$$

(11)

Following these definitions, it is clear that all Bézier curves automatically follow the specified boundary conditions at $t=t_0$ and $t=t_f$.

5. Implementation of the Monte Carlo Metropolis for the Sampling Path Space

The Monte Carlo Metropolis (MCM) implementation is usually employed for sampling a multidimensional space governed by a probability distribution [5,6]. The following MCM implementation is used for sampling the path space $\mathbb{X}$ governed by a probability functional obtained via MaxCal. This procedure makes it possible to find the minimum action path.

For a given action $A\left[x\right]$ and boundary conditions $x\left(t_0\right)=c_0$ and $x\left(t_f\right)=c_n$, the MCM evolution for sampling the path space is implemented as shown in Figure 1.

By performing this process in an iterative way, the path space is sampled, allowing the calculation of the properties for the system which is determined by the classical action used. As shown in Section 3, the probability distribution obtained when the classical action is constrained allows the sampling of the path space where the most probable path and the least action path coincide. Finally, $\lambda$ can be understood as an analogue of the inverse of temperature in the usual MCM, due to the fact that, as $\lambda \to 0$, the sampled paths will be randomly more spread over the space, while taking $\lambda \to \infty$ constrains the sampled paths closer to the least action path. The value for $\lambda$ in an MCM is related to the change from x to $x^{\prime }$, and empirically, this change must be adjusted to have approximately a $80\%$ acceptance rate. We used the value $\lambda =100$ for all the MCM simulations performed in this work.

6. Results and Discussion

6.1. Free Particle Action

The equation of motion for a free particle is obtained by minimizing the classical action:

$$A\left[x\right]={\int }_{t_0}^{t_f}dt\frac{m\dot{x}{\left(t\right)}^2}{2}$$

according to the least action principle. For a free particle with mass $m=1$ and boundary conditions $x\left(0\right)=0$ and $x\left(1\right)=1$, the analytical solution for the least action path is the straight line $x\left(t\right)=t$.

By using MCM simulation as shown in Figure 2, we sampled the path space and calculated the simple average position $\overline{x}\left(t\right)$ at each time. We see that the simulation readily converges to the correct least action path.

For this simulation, five control points were used and were sufficient to obtain the exact result ($R^2=0.99$) in less than 10,000 Monte Carlo steps.

6.2. Harmonic Oscillator Action

In the case of the harmonic oscillator, the action is of the form:

$$A\left[x\right]={\int }_{t_0}^{t_f}dt\left(\frac{m\dot{x}{\left(t\right)}^2}{2}-\frac{kx{\left(t\right)}^2}{2}\right).$$

(12)

Without loss of generality, for numerical simulations, $m=1$ and $k=1$ were used. The solution for this problem was divided into two parts. The first solution found was the least action path for a short time interval. More precisely, we simulated a particle with boundary conditions $x\left(0\right)=0$ and $x\left(t_f\right)=asin\left(\omega t_f\right)$, with $t_f$ less than the half period $\frac{T}{2}$.

In this case, the least action path also converges to the analytical solution, correctly solving the equation of motion as shown in Figure 3. For this simulation, five control points were used and were sufficient to obtain the exact result ($R^2=0.99$) in less than 10,000 Monte Carlo steps.

The second case corresponds to the harmonic oscillator with boundary conditions $x\left(0\right)=0$ and $x\left(t_f\right)=asin\left(\omega t_f\right)$, but where $t_f$ is larger than the half period $\frac{T}{2}$. In other words, the end condition is past the first node of the harmonic oscillator. Under these boundary conditions, an unexpected result is obtained: the Monte Carlo sampling procedure diverges. This result can be understood as due to the fact that, for this case, the action extremum is not a global minimum [24]. More precisely, the second functional derivative for the action of the harmonic oscillator action shows that the extremum is a saddle point in the case where the total time is longer than half a period, and a “true” minimum only for paths with total time less than half a period. We solved this convergence problem by considering an additional constraint to the action solved, suggested by the work of Gray and Taylor [24] in classical mechanics. The constraint involves the so-called kinetic foci, defined by the condition:

$$\frac{\partial x}{\partial v_0}=0.$$

(13)

where $v_0$ is the initial velocity for the path x. As it turns out, the action extremum is guaranteed to be a minimum if the paths used pass close enough to a kinetic focus $x_i$. This can be implemented in the Monte Carlo simulation by including a quadratic constraint in the probability functional, leading to:

$$P\left[x\right|\lambda ,\beta ]=\frac{1}{\eta (\lambda ,\beta )}{P}_0\left[x\right]exp(-\lambda A\left[x\right]-\beta \sum _i{(x\left(t_i\right)-x_i)}^2),$$

(14)

where $(t_i,x_i)$ are the set of kinetic foci (solutions of Equation (13)) and $\beta >>\lambda$, in order to stop the system from drifting away from the action extremum.

As an example, we solved the harmonic oscillator for a total time close to one and a half periods $\frac{3T}{2}$, sampling all the paths that cross the first two kinetic foci of the harmonic oscillator, namely the points $\{(0,\frac{T}{2}),(0,T)\}$. Now, the Monte Carlo procedure does converge to the expected solution as shown in Figure 4.

For this simulation, eight control points were used and were sufficient to reach the exact result ($R^2=0.95$) in less than 20,000 Monte Carlo steps. An important remark here is to note the number of Bézier basis elements, or control points, because this is closely related to the computing time. For this reason, the main goal in an efficient simulation is to find the minimum number of control points to use without sacrificing precision, needed to map any solution of the differential equation of interest.

7. Concluding Remarks

In summary, we described a technique for implementing Monte Carlo sampling of dynamical trajectories in classical Lagrangian systems under the maximum caliber formalism. We demonstrated its usefulness by applying this technique to the case of the classical free particle and harmonic oscillator, recovering in both cases a statistical distribution of paths centered on the classical solution of the Euler–Lagrange equation. For the case of the harmonic oscillator, we noted the need for fixing additional points known as the kinetic foci of the system in order for the simulation to converge properly. Our proof-of-concept implementation could be the starting point for a complete computational scheme of the simulation of classical systems under uncertainty. It remains to be seen how this method scales to multidimensional and many-particle interacting systems. This framework could be used for the study of dissipative systems with a known Lagrangian [25], being a promising line of work for the study of non-equilibrium physical problems. Finally, one of the main proposed uses of this framework is to obtain the instantaneous probability density of positions at each time, which would allow obtaining the instantaneous macroscopic properties of classical systems under uncertainty [26].