An Improved Adaptive Minimum Action Method for the Calculation of Transition Path in Non-gradient Systems

TL;DR

This paper introduces an enhanced adaptive minimum action method that effectively addresses path-tangling and non-smoothness issues in calculating transition paths in non-gradient stochastic systems, improving accuracy and robustness.

Contribution

The paper proposes a new monitor function and a generalized Euler-Lagrange scheme, combined with WENO interpolation, to improve the adaptive minimum action method for non-gradient systems.

Findings

Reduces path-tangling near transition states.

Handles non-smooth tangent directions at transition points.

Demonstrates improved accuracy in numerical examples.

Abstract

The minimum action method (MAM) is to calculate the most probable transition path in randomly perturbed stochastic dynamics, based on the idea of action minimization in the path space. The accuracy of the numerical path between different metastable states usually suffers from the "clustering problem" near fixed points. The adaptive minimum action method (aMAM) solves this problem by relocating image points equally along arc-length with the help of moving mesh strategy. However, when the time interval is large, the images on the path may still be locally trapped around the transition state in a tangle, due to the singularity of the relationship between arc-length and time at the transition state. Additionally, in most non-gradient dynamics, the tangent direction of the path is not continuous at the transition state so that a geometric corner forms, which brings extra challenges for the aMAM. In this note, we improve the aMAM by proposing a better monitor function that does not contain the numerical approximation of derivatives, and taking use of a generalized scheme of the Euler-Lagrange equation to solve the minimization problem, so that both the path-tangling problem and the non-smoothness in parametrizing the curve do not exist. To further improve the accuracy, we apply the Weighted Essentially non-oscillatory (WENO) method for the interpolation to achieve better performance. Numerical examples are presented to demonstrate the advantages of our new method.

Tables

Table 1. Table 1: L 2 subscript 𝐿 2 L_{2} error and L ∞ subscript 𝐿 L_{\infty} error of our numerical solution when T = 200 𝑇 200 T=200 . The parameters are Δ ​ τ = 0.1 Δ 𝜏 0.1 \Delta\tau=0.1 , r = 0.5 𝑟 0.5 r=0.5 . δ k subscript 𝛿 𝑘 \delta_{k} is the difference between two iterations when the algorithm stops.

$N$	${\delta }_k$	max-iter	$L_2$ Error	$L_{\infty }$ Error
50	1.4245e-05	1000	2.9000e-03	6.6000e-03
100	3.6864e-06	1000	4.9028e-04	1.4000e-03
200	1.7743e-05	1000	1.1779e-04	2.6127e-04

Table 2. Table 2: L 2 subscript 𝐿 2 L_{2} error of the numerical solution when T = 200 𝑇 200 T=200 . The number of points: N = 50 , 100 , 200 𝑁 50 100 200 N=50,100,200 . Number of iterations: max-it = 1000 absent 1000 =1000 . Step size: Δ ​ τ = 0.1 Δ 𝜏 0.1 \Delta\tau=0.1 .

$N$	$r=0.5$		$r=1$
	$L_2$ Error	$L_{\infty }$ Error	$L_2$ Error	$L_{\infty }$ Error
50	2.9000e-03	6.6000e-03	2.9000e-03	4.4000e-03
100	4.9028e-04	1.4000e-03	7.3990e-04	3.3000e-03
200	1.1779e-04	2.6127e-04	5.6424e-04	2.4000e-03

Table 3. Table 3: L 2 subscript 𝐿 2 L_{2} error and L ∞ subscript 𝐿 L_{\infty} error of our numerical solution when T = 100 𝑇 100 T=100 . The number of points N = 50 , 100 , 200 𝑁 50 100 200 N=50,100,200 . The parameters are selected as follows: Δ ​ τ = 0.05 Δ 𝜏 0.05 \Delta\tau=0.05 , r = 0.5 𝑟 0.5 r=0.5 .

$N$	${\delta }_k$	max-iter	$L_2$ Error	$L_{\infty }$ Error
50	3.6647e-05	1000	7.4775e-04	9.0732e-04
100	3.5343e-06	1000	1.7256e-04	2.0253e-04
200	6.6296e-08	1000	4.7831e-05	5.6775e-05

Table 4. Table 4: L 2 subscript 𝐿 2 L_{2} error of the numerical solution when T = 100 𝑇 100 T=100 . The number of points: N = 50 , 100 , 200 𝑁 50 100 200 N=50,100,200 . Number of iterations: max-it = 1000 absent 1000 =1000 . Step size: Δ ​ τ = 0.05 Δ 𝜏 0.05 \Delta\tau=0.05 .

$N$	$r=0.5$		$r=1$
	$L_2$ Error	$L_{\infty }$ Error	$L_2$ Error	$L_{\infty }$ Error
50	7.4775e-04	9.0732e-04	6.8663e-04	1.1000e-03
100	1.7256e-04	2.0253e-04	1.6042e-04	6.9453e-04
200	4.7831e-05	5.6775e-05	4.1854e-05	3.2541e-04

Equations

$$\dot{\phi }(t)=b(\phi (t))+\sqrt{ϵ}\sigma (\phi (t))\xi (t),$$

$$\lim_{\epsilon\downarrow 0}\epsilon\log P(X_{T}\in A)=-\inf_{\mbox{ \tiny$\begin{array}[]{c}\varphi(0)=\phi_{0},\\ \varphi(T)\in A\end{array}$ }}S_{T}(\varphi),$$

$$S_T(\phi )=\frac{1}{2}{\int }_0^T|\dot{\phi }-b(\phi (t)){|}_a^{2}dt,$$

$$V({\varphi }_0,{\varphi }_1)$$

$$k_{0\to 1}\asymp \exp (-{ϵ}^{-1}V({\varphi }_0,{\varphi }_1)),k_{1\to 0}\asymp \exp (-{ϵ}^{-1}V({\varphi }_1,{\varphi }_0)).$$

$$\underset{\phi \in {AC}_{{\varphi }_0}^{{\varphi }_1}(0,T)}{\inf }{S}_T(\phi ).$$

$$\{\begin{array}{l}{\phi }_{tt}-(\nabla b(\phi )-(\nabla b(\phi ){)}^T){\phi }_t-(\nabla b(\phi ){)}^{T}b(\phi )=0,t\in [0,T]\\ \phi (0)={\varphi }_0,\phi (T)={\varphi }_1.\end{array}$$

$$\min_{\mbox{ \tiny$\begin{array}[]{c}\varphi_{1}=\phi_{0},\\ \varphi_{N}=\phi_{1}\end{array}$ }}\frac{1}{2}\sum_{i=1}^{N}|\frac{\varphi_{i}-\varphi_{i-1}}{\Delta t_{i}}-b(\frac{\varphi_{i}+\varphi_{i-1}}{2})|^{2}\Delta t_{i},$$

$$\{\begin{array}{l}\frac{d}{dt}\left(\frac{1}{\omega (t)}\frac{d\alpha }{dt}\right)=0,\\ \alpha (0)=0,\alpha (T)=1,\end{array}$$

$$\frac{{\int }_{t_i}^{t_{i+1}}\omega (t)dt}{\hat{K}}=\Delta \alpha ,$$

$$\omega (t)=\sqrt{1+\tilde{K}|{\phi }_t{|}^2},$$

$$\omega =\frac{|b(\phi ){|}^r}{C},$$

$$\frac{d\alpha }{dt}=\omega (t)=\frac{|b(\phi (t)){|}^r}{C}.$$

$$\alpha (t)=\frac{{\int }_0^t|b(\phi ){|}^{r}d{t}^{\prime }}{C},$$

$$\{\begin{array}{l}{\phi }_{tt}-(\nabla b(\phi )-(\nabla b(\phi ){)}^T){\phi }_t-(\nabla b(\phi ){)}^{T}b(\phi )=0,t\in [0,T]\\ \phi (0)={\varphi }_0,\phi (T)={\varphi }_1.\end{array}$$

$${\phi }_t=\omega {\phi }^{\prime },$$

$${\phi }_{tt}={\omega }^2{\phi }^{\prime \prime }+\omega {\omega }_{\alpha }{\phi }^{\prime }.$$

$$\{\begin{array}{l}0={\omega }^2{\phi }^{\prime \prime }-\omega (\nabla b(\phi )-(\nabla b(\phi ){)}^T-{\omega }_{\alpha }){\phi }^{\prime }-((\nabla b(\phi ){)}^T)b(\phi ),\\ \phi (0)={\varphi }_0,\varphi (1)={\varphi }_1.\end{array}$$

$$S_T(\phi )=\frac{1}{2}{\int }_0^T|{\phi }_t-b(\phi ){|}^{2}dt=\frac{1}{2}{\int }_0^1\frac{1}{\omega }|\omega {\phi }^{\prime }-b(\phi ){|}^{2}d\alpha .$$

$$0=H({\phi }^{\ast },{\phi }_t^{\ast }-b({\phi }^{\ast }))=⟨b({\phi }^{\ast }),{\phi }_t^{\ast }-b({\phi }^{\ast })⟩+\frac{1}{2}⟨{\phi }_t^{\ast }-b({\phi }^{\ast }),{\phi }_t^{\ast }-b({\phi }^{\ast })⟩=\frac{1}{2}(|{\phi }_t^{\ast }{|}^2-|b({\phi }^{\ast }){|}^2),$$

$$\left|{\phi }_{\alpha }^{\ast }\right|=\left|{\phi }_t^{\ast }\right|\frac{dt}{d\alpha }=C\cdot {\left|b({\phi }^{\ast })\right|}^{1-r},$$

$${\varphi }^1(s)-{\varphi }^1(s_{\ast })\propto {\left|s-s_{\ast }\right|}^{\gamma },\gamma \in (0,1],$$

$$\frac{d{\varphi }^1}{ds}(s)\propto \text{mbox}sign(s-s_{\ast }){\left|s-s_{\ast }\right|}^{\gamma -1},$$

$$d\alpha /dt\propto {\left|b(\phi )\right|}^r,ds/dt\propto \left|b(\varphi )\right|.$$

$$\begin{array}{rl}\frac{d{\phi }^1}{d\alpha }(\alpha )& =ds/d\alpha \times \frac{d{\varphi }^1}{ds}(s)\\ & \propto {\left|b(\phi )\right|}^{1-r}{\left|s-s_{\ast }\right|}^{\gamma -1}\text{mbox}sign(s-s_{\ast })\\ & \propto {\left|s-s_{\ast }\right|}^{m-mr+\gamma -1}\text{mbox}sign(s-s_{\ast })\end{array}$$

$$m-mr+\gamma -1>0,i.e.,r<1-(1-\gamma )/m.$$

$$r<\gamma .$$

$$r<\underset{i}{\min }\left\{{\gamma }_i\right\}.$$

$${\phi }^{\ast }(\alpha )={\phi }_0{L}_{n,0}(\alpha )+{\phi }_1{L}_{n,1}(\alpha )+...+{\phi }_n{L}_{n,n}(\alpha ),$$

$$L_{n,k}(\alpha )=\frac{(\alpha -{\alpha }_0)(\alpha -{\alpha }_1)...(\alpha -{\alpha }_{k-1})(\alpha -{\alpha }_{k+1})...(\alpha -{\alpha }_{n-1})(\alpha -{\alpha }_n)}{({\alpha }_k-{\alpha }_0)({\alpha }_k-{\alpha }_1)...({\alpha }_k-{\alpha }_{k-1})({\alpha }_k-{\alpha }_{k+1})...({\alpha }_k-{\alpha }_{n-1})({\alpha }_k-{\alpha }_n)}.$$

Full text

\ams

82C26, 60H30, 34F05

1 Introduction

The calculation of quasi-potential and the most probably transition path between stable equilibria in metastable systems is of interest to researchers in the study of dynamics of complex and stochastic systems in long time scales [1]. The minimum action method (MAM) [2] was introduced to find such optimal paths by directly minimizing the action functional, which is the rate function in the large deviation theory [1]. The action minimizer within a certain rare-event subset in the path space, i.e., the minimum action path, carries the dominant contribution to the probability for the corresponding rare events. So, conditioned on the occurrence of rare events, the minimum action path is the most probable path under the influence of small noise in the long run. In many applications, the path represents the progressive physical process of important rare events such as phase transformation, chemical reaction, etc. Therefore, the numerical study of how to calculate the path efficiently is of great importance.

There have been quite many developments of numerical methods for the minimum action path. Firstly, when the system is of gradient type, i.e., the dynamics is gradient flow driven by a potential energy, the variational problem for the minimum action path gets simplified and it turns out that essentially the optimal transition path is simply the time-reversed trajectory of the gradient flow. The joint location of the “uphill” path from one well and the “downhill” path toward the other well in the phase space is an index-1 saddle point which serves as transition state. The min-mode eigen-direction of the saddle point collapses with the tangential direction of the path from the both sides. This result reveals a second important feature for the gradient system: the tangential direction of the path is always continuous (belonging to $C^{1}$ curve), even when it crosses the separatrix via the saddle point. In practice, the path-finding algorithms, such as the string method [3] never searches the “uphill” and “downhill” paths separately since the saddle point is unknown a priori, but these methods search the whole path once for all between local minima and then locate the saddle point from the numerical path and split the path into “uphill” and “downhill” segments for interpretations.

But in the non-gradient systems, the above features of the path no longer holds, due to the lack of detailed balance for the stationary probability distribution. The “uphill” path and the “downhill” path are distinctively different in nature. More importantly, they meet at the saddle point from different directions in two sides of the separatrix: a sharp corner is usually formed where the path crosses the separatrix. The “uphill” path escapes the characteristic boundary by choosing a direction different from the eigen-direction of the saddle point. The path then exhibits non-smoothness, meaning that if the path $\varphi$ is written in terms of arc length parameter $s\in[0,1]$, then the tangent vector $\varphi^{\prime}(s)$ is not continuous at $s_{*}$, even though $\left|\varphi^{\prime}(s)\right|\equiv const$, where $\varphi(s_{*})$ is the location of the saddle point. This is a generic phenomenon for transition path in non-gradient system and it is the origin of the non-Gaussian skewed distribution of the exist point on boundary. It also affects the prefactor estimation of mean exist time for non-gradient systems [4, 5]. Refer to the work of [6, 7, 9] for theoretical analysis of the connection to effect of focusing and caustics. To visualize this non-smooth feature of the paths, we shall present two examples including the Maier-Stein model in [8] in the figures shown later.

To locate the saddle points, on the other hand, the direct search of saddle point, such as the dimer method [10] or the gentlest ascent dynamics[11, 12], sounds an alternative strategy. Actually, in many existing practical applications [13, 14, 15], after applying the standard MAM, one can usually obtain a reasonable result of the path and then by checking the force along the path, one can further guess the locations of critical points along the curve to serve as initial guess for a further independent run of saddle-point search algorithms (from the most straightforward Newton method to the new method of GAD [11]). After successfully splitting the path into “uphill” and “downhill” pieces , one can then run the MAM again separately for each piece to further refine the path and the quasi-potential. This approach only works when the numerical path in the first round of calculation is not too far away from the true solution and also critically relies on the success of finding saddle points. So it should be regarded as a post-process step rather than a true path-finding approach.

It seems there has no much work of a natural and accurate numerical method to directly address the issue of non-smoothness of the path, without extra interference or post-processes mentioned above. The discontinuous tangential direction across the boundary arising in non-gradient systems imposes very subtle challenges for the numerical calculations of the transition path. In the Freidlin-Wentzel theory, the action functional is defined in the absolute continuous function space $AC[0,T]$ for a fixed $T>0$, while the optimal $T$ that further minimizes the action is infinitely. Actually, at the saddle point, say $x_{*}$, the uphill path, denoted by $\varphi_{+}(t)$, goes to $x_{*}$ as $t\to+\infty$ while the downhill path, denoted by $\varphi_{-}(t)$, emits from $x_{*}$ as $t\to-\infty$. The time derivatives $\dot{\varphi}_{+}$ and $\dot{\varphi}_{-}$ both vanish to zero as $t$ tends to $+\infty$ and $-\infty$ respectively. In practice, the two jointly infinite-time intervals are approximated by a single finite time interval $[0,T]$ for a very large $T$. When $T$ is large, the problem actually becomes ill-posed: the grid points are extremely dense in geometry near fixed points. Numerous computational results of this clustering effect show the critical importance of calculating the path by adopting the arc-length parametrization. To capture the geometric shape of the path efficiently, it is very important to represent the path by the grid points with equal geometric distance along the curve. The adaptive MAM, “aMAM”, [16] and the geometric MAM, “gMAM”, [17] were developed for this purpose by resolving the arc-length parametrization either numerically or analytically. The numerical approach in aMAM is based on the moving mesh strategy by choosing a monitor function equivalent to the derivative in time variable. The analytical approach in gMAM is based on the reformulation of the Freidlin-Wentzell action functional by using the Maupertuis’ principle. The main difference between these two methods is that aMAM works for arbitrary but finite $T$, but gMAM works for the optimal case $T=\infty$. The subsequent research work to improve the numerical efficiency of the MAM includes [18, 15, 29], which covers the topics of high order time discretization and spatial discretization for applications to Navier-Stokes equation.

As we said, the dynamics is simply zero at the saddle points, where the non-smooth corners of the path lie. So, one can find that the coefficient for the second spatial derivative in the associated Euler-Lagrange equation written in form of arc-length parametrization has a degeneracy at these zero-dynamics locations. This is the underlying reason why the a second-order equation will give arise to the non-smooth solution. But numerically, this actually might not be a severe challenge because we empirically find that the simple finite-difference scheme in the aMAM [16] can also guarantee the second order convergence rate, in terms of the grid size, for the $L_{2}$ error of the path. The real challenges in practice for the version of the aMAM proposed in [16] are: (1) for a large time interval, i.e., a large domain of $t$, the tangling phenomena usually come up and the effect for the moving mesh breaks down. (2) the cubic spline interpolation used in the reparametrization step lowers its accuracy around the corner. In short, the scheme used in [16] is not efficient to address both issues when a very large domain $[0,T]$ is used.

The tangling phenomena, illustrated in Figure 1 and 2 below, mean a lot of points (images) on the path are in a tangle near the fixed point, which usually arise for a large $T$ and does not appear for the small $T$ or a small number of discrete points representing the whole path. This implies that the moving mesh strategy failed to take effect near the saddle points and the reason is related to the monitor function used in [16]: $\omega(t)=\left|\varphi_{t}\right|$, for which the derivative ${\varphi_{t}}\approx(\varphi_{n+1}-\varphi_{n-1})/2\Delta t$ has to be calculated numerically. Near the saddle point, the condition number in the optimization problem is very large, the images $\{\varphi_{n}\}$ there behave nearly randomly so that any numerical derivatives based on these “tangled” images generates a huge error in approximating the true value, failing to correct the random-like distribution of the images back into order. The larger domain for $t$, the higher possibility for this “chaotic” tangling to occur.

Here we want to propose a new form of the monitor function $\omega$ for the aMAM to use, so that the $L_{\infty}$ norm of the error in numerical path can be controlled with good accuracy. In principle, one could introduce the idea like ENO or WENO (Weighted Essentially non-oscillatory) to approximate the derivative used by $\omega(t)$ in a better way. However, our idea is even simpler and motivated by the heuristic and the deep connection to the underlying Hamiltonian dynamics governing the path (which is a key ingredient in gMAM). Our new monitor function only depends on the dynamics at each images — no any derivative is involved, and has a scaling parameter $0<r\leq 1$ to further control the parametrizaton. It can also be shown that at the limit $T\to\infty$, the parametrization corresponding to the new monitor function with $r=1$ is still in arc-length parameter. In this way, the tangling problem occurred in numerical solution of aMAM is avoided and the moving mesh strategy becomes more robust. The second important consequence of our new monitor function is that with the proper choice of $r<1$ in our method, the path associated with the corresponding parametrization, $\varphi(\alpha)$ for $\alpha\in[0,1]$ , actually does lie in $C^{1}([0,1])$: the singularity is removed into the relation between the arc-length and our new parametrization, not in the path-finding algorithm (see Figure 4 in Section 3 later for this point). To further improve the interpolation accuracy, we also make use of the well-known WENO method to do reparameterization, which has been widely applied in capturing the discontinuous solutions such as shock waves and imaging sciences [19, 20, 21, 22, 23, 24, 25]. Interestingly, this idea of using the WENO scheme for path calculation has already been applied in the string method for gradient systems [26] to construct the higher order scheme to improve the discretizing accuracy in handling the Euler-Lagrange equation. However, as we already emphasized, only the path in non-gradient systems has the non-smooth structure and calls for the application of the WENO purely for the purpose of reconstruction of more accurately reparametrized paths. The traditional cubic spline interpolation for the gradient system used in [26] is not desired for our purpose here. We prefer the use of the high order WENO interpolation since it works very well regardless of the regularity of the path (as a function of $\alpha$).

The rest of the paper is organized as follows. In Section 2, we review the minimum action method and adaptive minimum action method. And in Section 3, we give a derivation of our improved scheme of adaptive minimum action method, and we review the WENO scheme for estimating function values in reparameterization step, then we present our improved numerical scheme. In Section 4, we show numerical results of the tunnel-diode model and the Maier-Stein model. The last is the conclusion session.

2 Theoretical background and the Adaptive Minimum Action Method

In this section we describe the minimum action method [2] and adaptive minimum action method [16] for the calculation of minimum action path.

2.1 Theoretical background

Consider the system modeled by the following It$\hat{o}$ stochastic differential equation,

[TABLE]

where $b$ is the deterministic drift field, $\sigma$ is the noise amplitude tensor, and $\xi$ is a Gaussian white noise with zero mean and covariance $\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$. In the framework of the large deviation (cf. [1]), for an SDE like (1) in the vanishing noise limit $\epsilon\to 0$, the probability of a path starting from $\phi_{0}$ at time [math] and end in $A$ at time $T$ is given by

[TABLE]

where the action functional

[TABLE]

and $a(x):=\sigma(x)\sigma(x)^{T}$ is the diffusion tensor, $\langle u,v\rangle_{a}=\langle u,a^{-1}v\rangle$ and $|u|_{a}=\sqrt{\langle u,u\rangle_{a}}$. Also, for generalization of time $T$ goes to infinity, the quasi-potential from point $\phi_{0}$ to point $\phi_{1}$ is defined by

[TABLE]

where ${AC}_{x}^{y}(0,T)$ is the space of absolutely continuous functions $f:[0,T]\rightarrow R^{n}$ with $f(0)=x$ and $f(T)=y$. All the starting points $\phi_{0}$ and the ending points $\phi_{1}$ in concern throughout this article are two stable fixed points of the field $b(\varphi)$. And the transition rates between two states $\phi_{0}$, $\phi_{1}$ are described by the following equations

[TABLE]

To find the minimum of action $S_{T}(\varphi)$ is essential in calculating probability of the event $(X_{T}\in A)$ in equation (2) and transition rate in (5). The Minimum Action Method calculates the following minimization problem:

[TABLE]

And the minimizer $\varphi^{*}$ is the minimum action path (MAP), which characterize the most probable path from $\phi_{0}$ to $\phi_{1}$ under the influence of small noise. We can obtain the minimum action path by optimization algorithms to solve the problem (6). To solve the minimization problem (6) in practice, one usually uses the standard optimization method such as evolving the gradient flow $-\delta S/\delta\varphi$, whose steady state equation is the Euler-Lagrange equation of the following boundary value problem

[TABLE]

For problem (6), in many examples, the minimum action path passes through the saddle points with two neighbouring stable states as the starting and ending points. For gradient system, this statement is rigorous [30]. For non-gradient systems, there is no rigorous statement, since the existence of saddle point is still a question. But for many known examples where some saddle point lies on the separatrix, the optimal path does pass the saddle points [31]. As $\varphi$ crosses a fixed point, it is known that the time derivative $\varphi_{t}\rightarrow 0$ [31, 17] . If the original MAM in[2], where the time interval is discretized by equal time step size, is applied, this property of slowdown in dynamics near saddle point brings the“clustering problem” to the calculation of the MAP i.e., too many points are around the three critical states while too few points are allocated for shaping the path. The clustering problem increases errors and significantly reduces the accuracy of the algorithm. Another issue for the MAM in [2] is the existence of error for the quasi-potential problem (6) due to the truncation of the time interval. Numerical results show, when $T$ is small, the solution of the MAP deviates from the true solution, which effectively corresponds to $T=\infty$; when $T$ is large, the solution of the MAP is close to the true solution (Figure 1). Thus a sufficiently large $T$ is needed to capture the true MAP.

2.2 The Adaptive Minimum Action Method

The adaptive minimum action method was proposed to solve the clustering problem by numerically adjusting the image points $\{\varphi_{i}\}_{i=1}^{N}$ to be equally spaced with respect to arc-length. This method in [16] used the simplest discretized form of the integral (6) by the finite difference method and rectangle quadrature

[TABLE]

where $\Delta t_{i}=t_{i}-t_{i-1}$. Suppose at the $k_{th}$ step of some minimization method for (8), the image points are not equally spaced in arc-length, then the aMAM has a moving mesh and reparameterization step to redistribute the image points, based on the following relation between $t$ and a new parameter $\alpha$:

[TABLE]

where $\omega$ is the monitor function. Then the problem is actually solved based on the uniform partition $\{\alpha_{i}\}_{i=1}^{N}=i/N$.

The moving mesh and reparameterization step in [16] is implemented as follows: given the partition and image points $(t_{i}^{k},\varphi^{k}_{i}),\ i=1,\ldots,N$, at the $k_{th}$ iteration, the aMAM first finds the corresponding partition for the $\alpha$-parameter $\{\alpha_{i}^{k}\}_{i=1}^{N}$ by equation (9), and then finds a new partition $\{\hat{t}_{i}^{k}\}$ by interpolating $(\alpha_{i}^{k},t_{i}^{k})$ on the uniform partition $\{\hat{\alpha}_{i}^{k}=i/N\}$, and finally locates the new image points $\{\hat{\varphi}^{k}_{i}\}_{i=1}^{N}$ on numerical solution of the path from the new partition $\{\hat{t}^{k}_{i}\}_{i=1}^{N}$ by the cubic spline interpolation. Then the partition and image points $(\hat{t}_{i}^{k},\{\hat{\varphi}^{k}_{i}\}_{i=1}^{N})$ are used to solve the numerical scheme of the minimization problem (8) in the next step.

Equation (9) indicates that the partition is denser where the monitor function $\omega(t)$ is larger. Integrate (9), then we can see that

[TABLE]

where $\hat{K}=\int_{0}^{T}\omega(t)dt$. Here the grid size $\Delta\alpha$ for $\alpha$ is assumed constant. The aMAM in [16] selects the monitor function as

[TABLE]

with some constant $\tilde{K}>0$. Thus by solving equation (9) with a monitor function (11), the aMAM adjusts the partition of time parameter $t$ to be denser where the time derivative $|\varphi_{t}|$ is larger, $i.e.$, where we have fast dynamics. And the aMAM usually uses a large enough $\tilde{K}$ to allocate equally spaced image points. We can observe that when $\tilde{K}$ is large, equation (10) indicates that the arc length element, $\int_{t_{i}}^{t_{i+1}}\tilde{K}|\varphi_{t}|dt\approx\int_{t_{i}}^{t_{i+1}}\omega(t)dt=\hat{K}\Delta\alpha$, is approximately a constant. In this way, the aMAM capturers the geometric shape of the minimum action path efficiently and thus achieves better accuracy in most cases.

However, the singularity of the relationship between the increment of arc length ($d\varphi$) and the time variable ($dt$) can not be resolved, for the fact that $\frac{d\varphi}{dt}\rightarrow 0$ at the saddle point is the nature of the minimum action path, particularly for a large $T$ used in practice. But for a large $T$, a tangling phenomenon at the saddle may arise numerically in aMAM: the images on the path get locally trapped around the saddle by forming unpredictable zig-zags (the specific tangling profiles depend on initial condition and the frequency to apply the moving mesh step). And such tangling phenomena only exists in a local tiny region near the saddle point. Figure 2 illustrates this abnormal effect for the tunnel-diode model and the Maier-Stein model (Section 4). In these examples, we use the number of images as large as $N=2000$ for the large $T$ to show that increasing the number of images can not solve the problem. In fact, the more image points used, the more possible for the tangling to occur. It is noted from numerical experiments that the above aMAM in [16] is much more likely to generate such tangling than the original MAM without moving mesh in [2] (compare Figure 1(b) and 2(a)).

To illustrate the potential mechanism for the emergence and persistence of such tangling phenomena in the path optimization algorithms, we show in Figure 3(a) the error in calculating the time derivative $\varphi_{t}\approx\frac{\varphi_{i+1}-\varphi_{i}}{\Delta t_{i}}$ of a typical aMAM solution after a few numbers of iterations. Theoretically, $\left|\varphi_{t}\right|=\left|b(\varphi(t))\right|$ holds for the true optimal path; refer to [17, 9] as well as Section 3.1 for this conclusion. However, our plot clearly indicates that the approximation near the saddle point deviates more from the true value. This observation gives us a strong hint that the use of ${b(\varphi)}$ is more robust than $\varphi_{t}$ in our new method in the next section. Actually, in [29], the difference quantity $\left|\varphi_{t}(\cdot)\right|-\left|b(\varphi(\cdot))\right|$, is exactly the criteria for the posterior error estimate to implement their $h$-adaptivity. Next, we show the magnitude of the gradient of $S$ (i.e., the residual of E-L equation) in Figure 3(b). This figure suggests that near the saddle point, the sensitivity for the objective function $S$, $\delta S/\delta\varphi$, is quite small, while the error in path measured by the deviation to the true solution is relatively large near this saddle point. So, any type of local update near saddle point would not improve the objective function as effective as in other regions. In particular, during the moving mesh step, if the mesh redistribution or the interpolation of the path from the old partition to the new one accidentally triggers some slight tangling points around the saddle point region, then the gradient $\delta S/\delta\varphi$ does not feel such changes there to effectively untangle the bad points. This heuristically explains why in practice, once the tangling happens, it is almost impossible to get back in order again. In the original adaptive MAM [16], the reason why it is so easy to trigger the tangling is twofold: (1) the monitor function (11) used there is not robust because it is completely determined by $\varphi_{t}$, and the (relative) numerical error in the approximation of $\varphi_{t}$ is large at the saddle point. Therefore, the quality of the new partition is not locally very satisfying and deviates from the true arc-length parametrization near saddle points; (2) the cubic spline interpolation is dangerous: it either smooths out the corner (if no tangling) or magnify the effect once one or two image points start to oscillate. A simple procedure of adding more points locally apparently can not rescue this problem once the distribution of points has been locally ruined.

3 Methodology

The monitor function of moving mesh strategy in the aMAM scheme is problematic for the fact that it contains the numerical approximation of the derivative $\varphi_{t}$, $(\varphi_{i}-\varphi_{i-1})/\Delta t_{i}$. Based on the aMAM scheme, we propose an alternative choice of the monitor function $\omega$ in equation (9) as

[TABLE]

where $C=\int_{0}^{T}|b(\varphi(t^{\prime}))|^{r}dt^{\prime}$ is a normalizing constant, and $r$ is a positive constant. And the relation between the new variable $\alpha$ and time variable $t$ is

[TABLE]

With the boundary condition $\alpha(0)=0$, $\alpha(T)=1$, after integrating equation (12) on $t$, then we have that $\alpha$ is simply

[TABLE]

and $\alpha\in[0,1]$. For the system modeled by SDE (1), assume $\sigma(\varphi)=I$, $I$ is the identical matrix. The Euler-Lagrange equation associated with the minimization problem (6) is the following boundary value problem (see [2])

[TABLE]

We can obtain the minimum action path by numerically evolving the Euler-Lagrange equation (14) as the steepest descent dynamics of minimizing $S_{T}(\varphi)$. Denote $\varphi^{\prime}$ as the derivative of $\varphi$ with respect to $\alpha$, then we can observe that

[TABLE]

[TABLE]

Plug the two terms (15) and (16) into the Euler-Lagrange equation (14), we have the equivalent form of the Euler-Lagrange equation written in the new variable $\alpha\in[0,1]$:

[TABLE]

The action functional is calculated by

[TABLE]

Instead of solving the Euler-Lagrange equation (14) on time parameter $\{t_{i}\}_{i=1}^{N}$ with the moving mesh strategy in equation (12), we are actually solving its equivalent form (17) on $\{\alpha_{i}\}_{i=1}^{N}$.

3.1 The connection to gMAM

In (12), if $r=1$, the parametrization in $\alpha$ is just the arc-length parametrizaton, when $T$ is optimal ($T=\infty$). It is proven in the geometric minimum action method [17] and in [9] that for the optimal $T$, the solution to minimization problem (4), $\varphi^{*}$, satisfies the condition that the Hamiltonian $H(\varphi^{*},\varphi^{*}_{t}-b(\varphi^{*}))=0$. We then have

[TABLE]

i.e., the minimum action path $\varphi^{*}$ satisfies $|\varphi^{*}_{t}|=|b(\varphi^{*})|$. Then by $\eqref{q12}$, i.e., $\frac{d\alpha}{dt}=\frac{\left|b(\varphi^{*})\right|^{r}}{C}$, we have

[TABLE]

where $C=\int_{0}^{T}\left|b(\varphi^{*})\right|^{r}dt$. This can be understood as an $r$-dependent reparametrization of the path. $r=1$ corresponds to the arc-length parametrization. So $r=1$ gives $\left|\phi^{*}_{\alpha}\right|\equiv C$, i.e., $\alpha$ is the arc-length parameter. Thus, essentially, when $T$ is optimal and $r=1$, our method is equivalent to the gMAM (with the arc length parametrization), but we have more flexibility here. The Euler-Lagrangian equation in the gMAM has the similar form as (17) by replacing $\omega(t)$ by $\lambda(\alpha)=\frac{\left|b(\varphi)\right|}{\left|\varphi_{\alpha}\right|}$ which involves the derivative of $\varphi$. The monitor function $\omega={\left|b(\varphi)\right|^{r}}/{\int_{0}^{T}\left|b(\varphi)\right|^{r}dt}$ in our method needs the integral in time variable but without any derivatives.

3.2 The selection of $r$

We have shown that $r=1$ corresponds to the arc length parametrization. If $r$ is larger than $1$, then the equation (12) will favor the region where $\left|b\right|>1$ and less points (than the arc-length case $r=1$) will be placed around the saddle point where $\left|b\right|\ll 1$. So usually a large $r$ is not recommend to handle the non-smoothness near the saddle points. On the other hand, if $r=0$, then $\alpha=t/T$, this is just the (linearly rescaled) time parameter.

So, selecting $r$ within $(0,1]$ seems a good balance between the arc-length parametrization ($r=1$), which is good for capturing the geometric shape of the whole path, and the physical time parametrization ($r=0$), which favours the saddle points. Notice that at non-smooth saddle point, we need more image points to improve the accuracy of numerical solution. Using equally spaced image points is no longer a good choice. We can decrease the value of $r$ to make the image points denser at saddle points, but not too dense as for $r=0$. In the following Figure 4, we show the profiles of the different parametrizations with $r=1$ and $r=0.5$, respectively, of the (same) optimal path $\varphi(\cdot)=(x(\cdot),y(\cdot))$ for the Maier-Stein model. Note that $r=1$ actually gives the arc length parametrization in $s$. It is clearly seen that for instance, the $y$-component $y(s)$ (the dashed line) does not have a continuous derivative $dy/ds$ at the saddle point but if the $\alpha$-parametrization is associated with $r=0.5$, both $x(\alpha)$ and $y(\alpha)$ are smooth functions of $\alpha$. In fact, it is observed that the first order derivative $dx/d\alpha$ and $dy/d\alpha$ both vanish at the saddle point for $r=0.5$. The singularity is transferred to the relation between $\alpha$ (for $r=0.5$) and $s$ (for $r=1$): $d\alpha/ds\propto\left|b(\varphi)\right|^{0.5-1}=\left|b\right|^{-0.5}$. So in the right panel, at the saddle point where $b=0$, the derivative of $\alpha$ w.r.t. $s$ is infinitely large.

Next, we present some heuristic scaling analysis for the observations from the previous figure. Again we use $s$ to refer to the arc-length parameter ($r=1$) and $\alpha$ for the parametrization associated with a general $r>0$. The path is thus referred as $\phi(s)$ or $\varphi(\alpha)=\phi(s(\alpha))$. The non-smoothness in $s$ means at least for one component of the path $\phi$, say the first component $\phi^{1}(s)$, does not have a continuous derivative at $s_{*}$ which corresponds to the saddle point $\phi(s_{*})$, more precisely $\partial_{s}\phi^{1}(s_{*}-)\neq\partial_{s}\phi^{1}(s_{*}+).$ Therefore, for $s$ sufficiently close to $s_{*}$,

[TABLE]

or,

[TABLE]

where “sign” is the sign function. Note that we always have the norm of the vector $\left|\phi(s)-\phi(s_{*})\right|\propto\left|s-s_{*}\right|$ since $\left|\phi^{\prime}(s)\right|\equiv const$. Around the saddle point $\phi(s_{*})$, $b(\phi)=b(\phi)-b(\phi(s_{*}))\propto\left|\phi(s)-\phi(s_{*})\right|^{m}\propto\left|s-s_{*}\right|^{m}$ and $m=1$ if the Jacobi $\partial_{x}b(x)$ does not vanish. By the definition,

[TABLE]

Therefore,

[TABLE]

When $s\to s_{*}$, $\alpha\to\alpha_{*}$. So, the equivalent condition for the derivative of the component $\varphi^{1}$ w.r.t. $\alpha$ is continuous is that

[TABLE]

It easy to see that under this condition, the derivative $d\varphi^{1}/d\alpha$ actually vanishes at $\alpha=\alpha_{*}$. For the typical case of $m=1$, the condition becomes

[TABLE]

$\gamma\leq 1$ describes the type of the singularity of the cusp of the curve in one component. In general case of non-smoothness in multiple components with the corresponding $\left\{\gamma_{i}\right\}$ in each dimension, the condition for $r$ to eliminate the discontinuous $d\varphi/d\alpha(\alpha_{*})$ is

[TABLE]

In our two numerical examples here (and in most real applications, we guess), $\gamma=1$. So $r<1$ seems sufficient in theory. We numerically tested $r=1/2$ and $r=1$, and find that $r=1/2$ strictly less than one is indeed better for the error measured in $L_{\infty}$ norm. This observation is consistent with the empirical conclusion in the recent work [29] that the optimal distribution of the finite elements computed from the $h$-adaptivity is not equally spaced in arc-length, but is denser near the saddle point.

3.3 WENO Scheme for interpolation of the curve $\varphi$

To improve the accuracy of interpolations after the new partition is generated, the WENO method [25] is applied in the reparameterization steps. Next, we briefly present a basic review of this interpolation method. For a function $\varphi(\alpha)$, given the data points $(\alpha_{1},\ \varphi_{1})$, $(\alpha_{2},\ \varphi_{2})$, …, $(\alpha_{n},\ \varphi_{n})$, $i.e.$, $\varphi_{i}=\varphi(\alpha_{i})$, $i=0,1,2,...,n$, the Lagrangian interpolation function $\varphi^{*}(\alpha)$ is given as an approximation of $\varphi(\alpha)$:

[TABLE]

[TABLE]

To get the approximation of the function value at $\alpha_{*}$, where $(\alpha_{i-2}\leq\alpha_{*}\leq\alpha_{i+2})$, the following three interpolation can be used.

[TABLE]

[TABLE]

[TABLE]

These three interpolations are built in the local stencils $S_{1}=\{\alpha_{i-2},\alpha_{i-1},\alpha_{i}\}$, $S_{2}=\{\alpha_{i-1},\alpha_{i},\alpha_{i+1}\}$, and $S_{3}=\{\alpha_{i},\alpha_{i+1},\alpha_{i+2}\}$, respectively. If a larger stencil $S=\{\alpha_{i-2},\alpha_{i-1},\alpha_{i},\alpha_{i+1},\alpha_{i+2}\}$ is used, a fifth-order approximation is then obtained:

[TABLE]

Such a high order approximation $\tilde{\varphi}_{*}$ can be written as a linear convex combination of the three third-order approximations [27]:

[TABLE]

The WENO scheme chooses a convex combination of the three approximations $\varphi_{*}^{(1)},\varphi_{*}^{(2)},\varphi_{*}^{(3)}$ as:

[TABLE]

where the weights $\omega_{1},\omega_{2},\omega_{3}$ are given based on smooth indicator $\beta_{j}$ in each stencil $S_{j}$, $j=1,2,3$.

[TABLE]

The smoother the function $\varphi^{(j)}(\alpha)$ is in the target cell, the smaller this smoothness indicator $\beta_{j}$ is. The choice of $\omega_{j}$ is as follows.

[TABLE]

And in fact the weights might be negative in many cases. Then the above procedures can not provide us a stable scheme. One can apply the technique in [28] to deal with the negative weights.

3.4 Numerical Scheme of Improved aMAM for Non-smooth Path

The main concern of our problem is the case that $T$ is large. For a large time interval $T$, however, the use of (7) in time domain still has some issues near the saddle point for the reason that $\Delta t_{i}$ is still very large there, although it is relatively small compared with the total length $T$. So we prefer to working on the $\alpha$-domain $[0,1]$ to solve the Euler-Lagrangian equation (17) as a substitution of equation (7). Given the initial path in time discretization, $\varphi_{i}^{k}=\varphi(k\Delta\tau,\sum_{j=1}^{i}\Delta t_{j}^{k})$, $\Delta t_{j}^{k}=t^{k}_{j}-t^{k}_{j-1}$, where $k=1$ is the $1_{st}$ iteration, we can find the corresponding discrete $\alpha$ partition by equation (12). Then we reparametrize this path in $\alpha$ variable with equal grid size, still denoted as $\varphi_{i}^{k}=\varphi(k\Delta\tau,i\Delta\alpha)$ where $\Delta\alpha=1/N$. We then solve the following numerical scheme of the Euler-Lagrange equation (17) written in the variable $\alpha$.

[TABLE]

where $\omega_{\alpha}$ is calculated by $(\omega_{\alpha})_{i}^{k}=(\omega_{i+1}^{k}-\omega_{i-1}^{k})/2\Delta\alpha$. $\omega$ is defined in (12) where the constant $C=\int_{0}^{T}\left|b(\varphi)(t)\right|^{r}dt$ has to be calculated by using the numerical quadrature from the time partition grids. Both the Euler-Lagrange equation and the WENO scheme we worked on are associated with the parameter $\alpha$, but we keep updating the discretization with respect to time because it is needed for the computation of the constant $C_{k}$.

At the ${k}_{th}$ iteration, we need to numerically calculate the relation between $\alpha$ and $t$ to proceed the reparametrization and moving mesh. The disrectized form of equation (12) in the moving mesh method is

[TABLE]

where $C_{k}$ is found by calculating the integration $C_{k}=\int_{0}^{T}|b(\varphi^{k})|^{r}dt$. Notice there is singularity between $\alpha$ and $t$ where $b$ is zero, $i.e.$, when $\alpha_{i}$ is uniform, $\Delta t_{i}^{k}$ can be very large at the saddles where $|b(\varphi_{i}^{k})|\rightarrow 0$. Thus the integration (30) should be carefully calculated when $\Delta t_{i}^{k}$ is too large, for example, we can interpolate the discrete path with a very fine partition (for instance, break each interval into $M$ pieces) on interval $[t^{k}_{i-1},t^{k}_{i}]$ and then apply the quadrature scheme.

The reparameterization steps are as follows: given the partition and image points $(t_{i}^{k},\varphi^{k}_{i}),\ i=1,...,N$ at the $k_{th}$ iteration, we first find the corresponding parameter $\{\alpha_{i}^{k}\}_{i=1}^{N}$ by equation (30), and then find the new partition $\{\hat{t}_{i}^{k}\}$ and the new image points $\{\hat{\varphi}_{i}\}_{i=1}^{N}$ on the same path corresponding to the prescribed uniform partition $\{{\alpha}_{i}=i/N\}$. The interpolation between $\alpha$ and the path should be carefully done by the third-order WENO scheme since the path might be non-smooth. Then the new image points on the uniform $\alpha$-partition $(\alpha_{i},\{\hat{\varphi}_{i}\}_{i=1}^{N})$ are used to solve the numerical scheme of the Euler-Lagrange equation (29).

In summary, we present a new scheme for the calculation of non-smooth minimum action path, for a fixed and large $T$. Given the initial $\{(t_{i}^{0},\ \phi_{i}^{0})\}_{i=0,...,N}$, we update the points as the following algorithm (Algorithm 1). In all numerical examples we show here, the initial guess is always the straight line connecting two ends.

4 Numerical examples

4.1 Example 1. The tunnel-diode model

Consider random perturbations of the state model of a tunnel-diode circuit.

[TABLE]

where $h(x)=17.76x-103.79x^{2}+229.62x^{3}-226.31x^{4}+83.72x^{5}$, and $W_{1}$ and $W_{2}$ are independent Wiener processes. There are two stable fix points $Q_{1}=(0.0626,0.7582)$ and $Q_{2}=(0.8844,0.2103)$, and the saddle point is $Q_{3}=(0.2854,0.6098)$. We use our method to calculate the minimum action path, and compare it with the solution of the aMAM scheme. The initial path is just a straight line connecting $Q_{1}$ and $Q_{2}$.

4.1.1 Stopping criterion

We run the iteration algorithm until the difference

[TABLE]

is smaller than a threshold value $\delta(=10^{-10})$. And we also set the maximum for the total number of iterations, denoted by “max-it”. So, the algorithm stops when either the accuracy threshold is attained or the maximum number of steps is achieved.

4.1.2 Real solution

The numerical results show the solution of the MAP calculated from interval $[0,T]$ is extremely close to the true solution of the MAP at optimal $T=\infty$. Approximately, we take the minimum action path calculated by the gMAM scheme ($i.e.$ the solution of MAP when $T$ is optimal) with 10000 points as real solution. To get the high accuracy at the saddle point, we actually find the exact location of the saddle point first with the Newton method, then use gMAM to calculate the downhill and uphill pieces separately. Thus, the solution is the union of these two pieces.

4.1.3 Error

To quantify the error between the numerical path and the true path, we first define the distance from a point $x\in\mathbb{R}^{d}$ to a continuous path $\varphi=\varphi(\alpha)\subset\mathbb{R}^{d}$:

[TABLE]

where $|\cdot|$ is the standard Euclidian norm in $\mathbb{R}^{d}$. Then the Hausdorff metric between two continuous path $\varphi=\varphi(\alpha)$ and $\hat{\varphi}=\hat{\varphi}(\beta)$ (they may have different parametrization forms) is

[TABLE]

We call this $L_{\infty}$ error if $\varphi$ is the numerical solution and $\hat{\varphi}$ is the true solution. This $error_{\infty}$ describes the largest gap between the the numerical path and the true path. Likewise, to describe the “averaged” deviation, we define the $L_{2}$ error by the line integral in the following sense

[TABLE]

It is clear that these definitions of errors are free of parametrization. If either error is zero, then the two paths are two identical sets.

In practice, a numerical path $\varphi$ is represented by $m$ images, $\varphi_{1},\ \varphi_{2},...,\ \varphi_{m}$, and the true solution obtained from a fine partition is also represented by $K$ images $\hat{\varphi}_{j},~{}j=1,2,...K$, ($K>m$). Then the $L_{\infty}$ error in discrete form is

[TABLE]

And the $L_{2}$ error in discrete form is

[TABLE]

We first interpolate our numerical solution to $m$ equally spaced points $\varphi_{i},\ i=1,2,...m$ by linear interpolation, and compare with the real solution (obtained by the gMAM with an extraordinary large number of points) of $K$ equally spaced points $\hat{\varphi}_{j},~{}j=1,2,...K$. Here we select $K=m=10000$.

In comparison of the numerical results in Figure 5(a) and Figure 5(b), we find our method improves the efficiency of aMAM at saddles. The numerical illustration shows our method can solve the tangling problem at saddle and directly captures the whole minimum action path. Table 1 shows $L_{2}$ errors and $L_{\infty}$ errors of our numerical solution are both small after $1000$ iterations. Figure 6 shows the $L_{2}$ errors and $L_{\infty}$ errors both decay as the number of iteration increases. And our method achieves second order convergence in $L_{2}$ error and $L_{\infty}$ error (see Figure 7). The above results shows our method performs well for solving the minimum action path.

4.1.4 Effect of $r$

We can compare the result of $r=1$ and $r=0.5$. Figure 8 and Table 2 show, the effect of $r$ is to adjust the density of grid points around saddles. When $r=0.5$, the image points are denser at the saddle point than $r=1$, and the numerical solutions are more accurate.

4.2 Example 2. The Maier-Stein model.

Consider the following example of a diffusion process proposed by Maier and Stein [8]:

[TABLE]

where $W_{1}$ and $W_{2}$ are independent Wiener processes, and $\beta>0$ is a parameter. For all values of $\beta$ there are two stable equilibrium points $(\pm 1,0)$, and an unstable equilibrium point $(0,0)$. And here we select the parameter $\beta=10$. We calculate the minimum action path by our method, and compare it with the solution of the aMAM scheme. The initial path is simply chosen as $y=-0.5x^{2}+0.5$. The stopping criterion, real solution, and $L_{2}$, $L_{\infty}$ errors are defined as in Section (4.1.1-4.1.3), and $K=m=20000$.

Figure 9 shows that our solution captures the whole minimum action path better than the original aMAM, for our method avoids the tangling problem (Figure 9(a) vs 9(b)). Table 3 shows $L_{2}$ error and $L_{\infty}$ error of our solution are both small after $1000$ iterations. Figure 10 shows the $L_{2}$ errors and $L_{\infty}$ errors of our solution decay at each with the iterations. And Figure 11 shows second order convergence of our solution in $L_{2}$ error and $L_{\infty}$ error.

4.2.1 Effect of $r$

Also, to illustrate the effect of selecting different values of $r$, we list the result of $r=1$ and $r=0.5$ in Figure 12 and Table 4. From the results, we can see the image points are denser at the saddle point and the numerical errors are smaller when $r=0.5$. This result shows we can adjust the density of grid points by $r$ and thus improve the efficiency of the method.

5 Conclusion

From the numerical results, we can see that our algorithm is efficient in dealing with the path-tangling phenomenon at saddle points. Our scheme improves the original aMAM by excluding the effect of the large numerical errors caused by singularity of the relationship between arc-length and time variable in the calculation of minimum action path. This is mainly achieved with the help of a better choice of the monitor function in moving mesh strategy, which is free of derivative calculation and has a flexible choice of the value $r$. When the value of $r$ is smaller than the singularity index $\gamma$, the path parametrized in our new variable $\alpha$ becomes a smooth function in $C^{1}$ so that the corner problem is avoided. Our numerical results confirm that the most efficient parametrization in practical computation may not be the arc length parametrization ($r=1$), which is consist with the recent work of adaptive mesh refinement based on the posterior estimate in [29]. Our idea of using the drift term $b(\varphi)$ can be generalized to the construction of path parametrized in other situations like the gMAM or even the string method for some special needs. Moreover, with the help of WENO in doing reparameterization of non-smooth path, we can find more accurate numerical solutions. To achieve high order accuracy, one can consider high order numerical schemes of Euler-Lagrange equation (17) and do reparameterization with high order interpolation schemes.

Acknowledgments

The research of Xiang Zhou was supported by grants from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. CityU 11304314, 11337216, 11304715).