Ordered Line Integral Methods for Computing the Quasi-potential

TL;DR

This paper introduces Ordered Line Integral Methods (OLIMs), a new family of algorithms for efficiently computing the quasi-potential in stochastic dynamical systems, offering significant speed and accuracy improvements over previous methods.

Contribution

The paper presents OLIMs, which improve upon the Ordered Upwind Method by optimizing updates and using quadrature rules, resulting in faster and more accurate quasi-potential computations.

Findings

OLIMs are 1.5 to 4 times faster than previous methods.

OLIMs achieve error reductions of two to three orders of magnitude.

Higher order quadrature rules significantly improve accuracy.

Abstract

The quasi-potential is a key function in the Large Deviation Theory. It characterizes the difficulty of the escape from the neighborhood of an attractor of a stochastic non-gradient dynamical system due to the influence of small white noise. It also gives an estimate of the invariant probability distribution in the neighborhood of the attractor up { to} the exponential order. We present a new family of methods for computing the quasi-potential on a regular mesh named the Ordered Line Integral Methods (OLIMs). In comparison with the first proposed quasi-potential finder based on the Ordered Upwind Method (OUM) (Cameron, 2012), the new methods are 1.5 to 4 times faster, can produce error two to three orders of magnitude smaller, and may exhibit faster convergence. Similar to the OUM, OLIMs employ the dynamical programming principle. Contrary to it, they (i) have an optimized strategy for the use of computationally expensive { triangle} updates leading to a notable speed-up, and (ii) directly solve local minimization problems using quadrature rules instead of solving the corresponding Hamilton-Jacobi-type equation by the first order finite difference upwind scheme. The OLIM with the right-hand quadrature rule is equivalent to OUM. The use of higher order quadrature rules in local minimization problems dramatically boosts up the accuracy of OLIMs. We offer a detailed discussion on the origin of numerical errors in OLIMs and propose rules-of-thumb for the choice of the important parameter, the update factor, in the OUM and OLIMs. Our results are supported by extensive numerical tests on two challenging 2D examples.

Figures

Tables

Table 1. Table 1: The least squares fits to the formulas E = C ​ N − q 𝐸 𝐶 superscript 𝑁 𝑞 E=CN^{-q} and T = C ​ N q 𝑇 𝐶 superscript 𝑁 𝑞 T=CN^{q} for the maximum absolute errors and the CPU times as functions of N 𝑁 N ( N = 2 p 𝑁 superscript 2 𝑝 N=2^{p} ) respectively for SDEs ( 24 ) (Columns 2 and 3) and ( 27 ) (Columns 4 and 5). The update factor K ​ ( p ) = p − 3 𝐾 𝑝 𝑝 3 K(p)=p-3 is chosen according to the Rule-of-Thumb for OUM and OLIM-R.

Method	Max Error	CPU time	Max Error	CPU time
OUM	$24.8\cdot N^{-0.860}$	$2.48\cdot {10}^{-6}\cdot N^{2.31}$	$19.1\cdot N^{-0.954}$	$3.52\cdot {10}^{-6}\cdot N^{2.30}$
OLIM-R	$24.8\cdot N^{-0.860}$	$0.817\cdot {10}^{-6}\cdot N^{2.26}$	$19.5\cdot N^{-0.957}$	$1.09\cdot {10}^{-6}\cdot N^{2.27}$
OLIM-MID	$47.5\cdot N^{-1.56}$	$1.55\cdot {10}^{-6}\cdot N^{2.25}$	$0.843\cdot N^{-0.944}$	$2.32\cdot {10}^{-6}\cdot N^{2.25}$
OLIM-TR	$50.3\cdot N^{-1.57}$	$1.60\cdot {10}^{-6}\cdot N^{2.26}$	$1.67\cdot N^{-1.02}$	$2.71\cdot {10}^{-6}\cdot N^{2.23}$
OLIM-SIM	$48.3\cdot N^{-1.56}$	$2.01\cdot {10}^{-6}\cdot N^{2.27}$	$0.923\cdot N^{-0.951}$	$3.87\cdot {10}^{-6}\cdot N^{2.23}$

Table 2. Table 2: Comparison of Maximal Errors, RMS errors and CPU times of the OUM and OLIM-R applied to SDE ( 24 ).

Method, $N$, $K$	Max Error	RMS error	CPU time, seconds
OUM, $N=512$
$K=3$	1.7669e-01	1.0440e-01	2.05
$K=5$	1.2133e-01	7.9878e-02	3.69
$K=7$	1.2058e-01	7.9659e-02	5.60
OLIM-R, $N=512$
$K=3$	1.8368e-01	1.0706e-01	0.55
$K=5$	1.2133e-01	7.9878e-02	0.88
$K=7$	1.2058e-01	7.9659e-02	1.26
OUM, $N=1024$
$K=4$	8.3905e-02	5.2161e-02	11.63
$K=6$	6.6225e-02	4.4102e-02	19.15
$K=8$	6.5912e-02	4.3992e-02	26.29
OLIM-R, $N=1024$
$K=4$	8.4836e-02	5.2346e-02	3.02
$K=6$	6.6225e-02	4.4102e-02	4.56
$K=8$	6.5912e-02	4.3992e-02	6.06
OUM, $N=2048$
$K=5$	4.1289e-02	2.6584e-02	62.80
$K=7$	3.5510e-02	2.3803e-02	94.16
$K=9$	3.5350e-02	2.3743e-02	127.43
OLIM-R, $N=2048$
$K=5$	4.1529e-02	2.6609e-02	15.04
$K=7$	3.5510e-02	2.3803e-02	21.55
$K=9$	3.5350e-02	2.3743e-02	28.60
OUM, $N=4096$
$K=5$	2.6959e-02	1.6629e-02	263.71
$K=7$	1.9089e-02	1.2776e-02	397.53
$K=9$	1.8656e-02	1.2599e-02	535.13
OLIM-R, $N=4096$
$K=5$	2.7204e-02	1.6702e-02	66.76
$K=7$	1.9051e-02	1.2769e-02	96.03
$K=9$	1.8656e-02	1.2599e-02	127.22

Table 3. Table 3: Comparison of Maximal Errors, RMS errors and CPU times of the OUM and OLIM-R applied to SDE ( 27 ).

Method, $N$, $K$	Max Error	RMS error	CPU time, seconds
OUM, $N=512$
$K=3$	5.1079e-02	2.1963e-02	2.71
$K=5$	5.0563e-02	2.1874e-02	5.05
$K=7$	5.0563e-02	2.2332e-02	7.71
OLIM-R, $N=512$
$K=3$	5.1925e-02	2.2371e-02	0.81
$K=5$	5.0850e-02	2.1940e-02	1.29
$K=7$	5.0845e-02	2.2371e-02	1.83
OUM, $N=1024$
$K=5$	2.5930e-02	1.1245e-02	20.50
$K=8$	2.5890e-02	1.1447e-02	36.65
$K=11$	2.5890e-02	1.1696e-02	53.57
OLIM-R, $N=1024$
$K=5$	2.6182e-02	1.1283e-02	5.36
$K=8$	2.6031e-02	1.1454e-02	8.71
$K=11$	2.6025e-02	1.1692e-02	12.21
OUM, $N=2048$
$K=5$	1.3478e-02	5.8271e-03	83.24
$K=10$	1.3204e-02	5.8695e-03	193.00
$K=15$	1.3189e-02	6.0301e-03	311.84
$K=20$	1.3190e-02	6.1811e-03	443.35
OLIM-R, $N=2048$
$K=5$	1.3656e-02	5.8558e-03	22.78
$K=10$	1.3215e-02	5.8607e-03	46.23
$K=15$	1.3200e-02	6.0178e-03	72.81
$K=20$	1.3201e-02	6.1618e-03	103.35
OUM, $N=4096$
$K=5$	7.2210e-03	3.0923e-03	341.70
$K=10$	6.6921e-03	2.9666e-03	791.00
$K=15$	6.6527e-03	3.0332e-03	1276.04
$K=20$	6.6506e-03	3.0975e-03	1811.45
OLIM-R, $N=4096$
$K=5$	7.3480e-03	3.1267e-03	101.85
$K=10$	6.7011e-03	2.9587e-03	210.50
$K=15$	6.6555e-03	3.0230e-03	326.48
$K=20$	6.6536e-03	3.0859e-03	463.64

Table 4. Table 4: The least squares fits to the formulas E = C ​ N − q 𝐸 𝐶 superscript 𝑁 𝑞 E=CN^{-q} for the maximum absolute errors and the RSM errors as functions of N 𝑁 N ( N = 2 p 𝑁 superscript 2 𝑝 N=2^{p} ) for SDEs ( 24 ) (Columns 2 and 3) and ( 27 ) (Columns 4 and 5). The update factor K ​ ( p ) = 10 + 4 ​ ( p − 7 ) 𝐾 𝑝 10 4 𝑝 7 K(p)=10+4(p-7) is chosen according to the Rule-of-Thumb for OLIM-MID, OLIM-TR, and OLIM-SIM.

Method	Max Error	RMS Error	Max Error	RMS Error
OLIM-MID	$0.817\cdot N^{-1.39}$	$0.705\cdot N^{-1.43}$	$2.47\cdot N^{-1.10}$	$5.85\cdot N^{-1.41}$
OLIM-TR	$1.31\cdot N^{-1.44}$	$1.16\cdot N^{-1.48}$	$1.61\cdot N^{-1.03}$	$0.646\cdot N^{-1.07}$
OLIM-SIM	$1.07\cdot N^{-1.42}$	$0.99\cdot N^{-1.46}$	$1.42\cdot N^{-1.02}$	$0.846\cdot N^{-1.15}$

Equations

$$d\mathbf{x}=\mathbf{b}(\mathbf{x})dt+\sqrt{ϵ}dW,\mathbf{x}\in {\mathbb{R}}^d$$

$$U_A(\mathbf{x}):=\underset{\psi :\psi (0)\in A,\psi (L)=\mathbf{x}}{\inf }\left\{S(\psi ):={\int }_0^L\left[\Vert {\psi }^{\prime }\Vert \Vert \mathbf{b}(\psi )\Vert -{\psi }^{\prime }\cdot \mathbf{b}(\psi )\right]ds\right\}.$$

$$\Vert \nabla U_A(\mathbf{x}){\Vert }^2+2\mathbf{b}(\mathbf{x})\cdot \nabla U_A(\mathbf{x})=0,U_A(A)=0.$$

$$\frac{1}{-2\mathbf{b}(\mathbf{x})\cdot \hat{a}}\Vert \nabla U_A(\mathbf{x})\Vert =1,\hat{a}\equiv \frac{\nabla U_A(\mathbf{x})}{\Vert \nabla U_A(\mathbf{x})\Vert },U_A(A)=0,$$

$$\mathbf{b}(x)=\left[\begin{array}[]{cc}\partial_{x_{1}}b_{1}&\partial_{x_{2}}b_{1}\\ \partial_{x_{1}}b_{2}&\partial_{x_{2}}b_{2}\end{array}\right](\mathbf{x}-\mathbf{x}_{0})+O(\|\mathbf{x}-\mathbf{x}_{0}\|^{2})=:\mathbf{A}(\mathbf{x}-\mathbf{x}_{0})+O(\|\mathbf{x}-\mathbf{x}_{0}\|^{2})$$

$$U(\mathbf{x})$$

$$\mathcal{A}$$

$$\mathcal{B}$$

$$\mathcal{C}$$

$$\alpha$$

$$\mathbf{b}(\mathbf{x})=-\frac{\nabla U(\mathbf{x})}{2}+\mathbf{l}(\mathbf{x}),\nabla U(\mathbf{x})\cdot \mathbf{l}(\mathbf{x})=0.$$

$${\mathbf{b}(\mathbf{x})|}_C={\mathbf{l}(\mathbf{x})|}_C,$$

$$-\frac{\nabla U(\mathbf{x})}{2}=\mathbf{b}(\mathbf{x})-\frac{\mathbf{b}(\mathbf{x})\cdot \mathbf{b}({\mathrm{Proj}}_C\mathbf{x})}{\mathbf{b}({\mathrm{Proj}}_C\mathbf{x})\cdot \mathbf{b}({\mathrm{Proj}}_C\mathbf{x})}\mathbf{b}({\mathrm{Proj}}_C\mathbf{x}).$$

$$\mathsf{Neib}(C):=\bigcup _{0\le i<M}\mathsf{Rectangle}({\mathbf{x}}_i,{\mathbf{x}}_{\mathsf{next}(i)}).$$

$$U({\mathbf{x}}_f)=U({\mathbf{x}}_i)+{\int }_0^1\nabla U\left({\mathbf{x}}_i+t({\mathbf{x}}_f-{\mathbf{x}}_i)\right)\cdot ({\mathbf{x}}_f-{\mathbf{x}}_i)dt.$$

$${\mathbf{x}}^{\ast }={\mathbf{x}}_1+((\mathbf{x}-{\mathbf{x}}_1)\cdot \hat{\tau })\hat{\tau },\mathrm{where}\hat{\tau }:=\frac{{\mathbf{x}}_2-{\mathbf{x}}_1}{\Vert {\mathbf{x}}_2-{\mathbf{x}}_1\Vert }.$$

$$U(\mathbf{x})\approx \Vert \mathbf{x}-{\mathbf{x}}^{\ast }\Vert \left(4\Vert \nabla U({\mathbf{x}}_m)/2\Vert +\Vert \nabla U(\mathbf{x})/2\Vert \right)/3.$$

$${\mathsf{R}}_{1pt}({\mathbf{x}}_0,\mathbf{x})$$

$$U(\mathbf{x})$$

$${\mathsf{Q}}_{1pt}({\mathbf{x}}_0,\mathbf{x})=U({\mathbf{x}}_0)+\mathcal{Q}({\mathbf{x}}_0,\mathbf{x}),U(\mathbf{x})=\min \{{\mathsf{Q}}_{1pt}({\mathbf{x}}_0,\mathbf{x}),U(\mathbf{x})\}.$$

$${\mathsf{Q}}_{\Delta }({\mathbf{x}}_1,{\mathbf{x}}_0,\mathbf{x})$$

$$U(\mathbf{x})$$

$$g(s):=\frac{d}{ds}\left(sU({\mathbf{x}}_0)+(1-s)U({\mathbf{x}}_1)+\mathcal{Q}(s{\mathbf{x}}_0+(1-s){\mathbf{x}}_1,\mathbf{x})\right)=0,$$

$$\left[\begin{array}[]{c}u-u_{0}\\ u-u_{1}\end{array}\right]=\left[\begin{array}[]{c}(\mathbf{x}-\mathbf{x}_{0})^{T}\\ (\mathbf{x}-\mathbf{x}_{1})^{T}\end{array}\right]\nabla U,$$

$$\left[u-u_{0},u-u_{1}\right]P^{-T}P^{-1}\left[\begin{array}[]{c}u-u_{0}\\ u-u_{1}\end{array}\right]+2\mathbf{b}^{T}P^{-1}\left[\begin{array}[]{c}u-u_{0}\\ u-u_{1}\end{array}\right]=0.$$

$$\left[\begin{array}[]{c}c_{0}\\ c_{1}\end{array}\right]=P^{-T}\left(\mathbf{b}(\mathbf{x})+P^{-1}\left[\begin{array}[]{c}u-u_{0}\\ u-u_{1}\end{array}\right]\right)$$

$$u=\underset{s\in [0,1]}{\min }\left[s{u}_0+(1-s)u_1+\Vert \mathbf{b}\Vert \Vert \mathbf{x}-{\mathbf{x}}_s\Vert -\mathbf{b}\cdot (\mathbf{x}-{\mathbf{x}}_s)\right],$$

$$\text{where}{\mathbf{x}}_s=s{\mathbf{x}}_0+(1-s){\mathbf{x}}_1.$$

$$U(\mathbf{x})=\underset{t\in \mathbb{R}}{\min }[U_0+f(t)],$$

$$f(t):=|\mathbf{b}|\frac{h}{\cos (\theta (t))}-|\mathbf{b}|\frac{h}{\cos (\theta (t))}\cos (\pi -(\alpha +\theta )).$$

$${\mathbf{x}}_0=\arg \underset{\mathbf{y}\in \mathsf{Accepted}\mathsf{Front},\Vert \mathbf{x}-\mathbf{y}\Vert \le Kh}{\min }{\mathsf{Q}}_{1pt}(\mathbf{y},\mathbf{x}).$$

$$U(\mathbf{x})=\min \{\underset{{\mathbf{x}}_1\in \mathsf{Accepted}\mathsf{Front},{\mathbf{x}}_1\in \mathcal{N}({\mathbf{x}}_0)}{\min }{\mathsf{Q}}_{\Delta }({\mathbf{x}}_1,{\mathbf{x}}_0,\mathbf{x}),U(\mathbf{x})\}.$$

$$d{x}_1$$

$$d{x}_2$$

Full text

Keywords: Quasi-potential; Ordered Line Integral Methods; Ordered Upwind Method; Accuracy; CPU time; Update radius; Hierarchical Update Strategy

MSC 2010 Numbers: 65N99; 49L20; 60J60

1 Introduction

The quasi-potential is a key concept of Large Deviation Theory [6]. Let us consider a nongradient stochastic differential equation (SDE) of the form

[TABLE]

where $\mathbf{b}(\mathbf{x})$ is a twice continuously differentiable vector field, $W$ is the standard $d$-dimensional Brownian motion, and $\epsilon$ is a small parameter. We assume that the corresponding deterministic system $\dot{\mathbf{x}}=\mathbf{b}(\mathbf{x})$ has a finite number of attractors. Let $A\subset\mathbb{R}^{d}$ be an attractor and $\mathbf{x}$ be a point. The quasi-potential at $\mathbf{x}$ with respect to the attractor $A$ can be defined as (details are provided in Appendix A)

[TABLE]

In Eq. (2), $\|\cdot\|$ is the 2-norm, the functional $S(\psi)$ is the geometric action [8, 9], $\psi(s)$ is an absolutely continuous path from $A$ to $\mathbf{x}$ parametrized by its arc length (i.e., $\|\psi^{\prime}\|=1$), $L$ is its length which can be infinite. According to the Large Deviation Theory, the expected first passage time from a small neighborhood of the attractor $A$ to a small neighborhood of the point $\mathbf{x}$ lying in its basin of attraction is logarithmically equivalent to $\exp(U_{A}(\mathbf{x})/\epsilon)$ [6]. The curve corresponding to the maximum likelihood path to $\mathbf{x}$ from $A$ is the global minimizer of the geometric action $S(\psi)$ [6, 8, 9]. The invariant probability distribution in the basin of attraction of $A$ is logarithmically equivalent to $\exp(-U_{A}(\mathbf{x})/\epsilon)$ [6]. A nice and visual account on the significance of the quasi-potential as a tool for quantifying the stability of attractors is found in [14].

Nongradient SDEs of the form (1) often arise in modeling of ecological and biological systems. A number of models of population dynamics were analyzed in [14] using the quasi-potential. The quasi-potential for a genetic switching model was constructed in [12].

Therefore, the quasi-potential is of crucial importance for the quantification of the dynamics of systems evolving according to SDE (1) with a small noise. Unfortunately, it is not readily available for nongradient SDEs and can be found analytically only in special cases. For example, analytic formulas for the quasi-potential are available for linear SDEs [1, 3]. For general nonlinear SDEs, finding the quasi-potential is a difficult task. In any bounded neighborhood of an attractor $A$ of $\dot{\mathbf{x}}=\mathbf{b}(\mathbf{x})$, the quasi-potential is a Lipschitz-continuous but not necessarily differentiable function as it is for the Maier-Stein [13, 1] model. One can show that the quasi-potential $U_{A}(\mathbf{x})$, i.e., the solution of the functional minimization problem (2), is also a viscosity solution of a Hamilton-Jacobi-type PDE [1]

[TABLE]

Eq. (3) has at least two viscosity solutions due to the fact that $A$ is an attractor and $\nabla U(\mathbf{x})=0$ for all $\mathbf{x}\in A$: $U_{A}(\mathbf{x})\equiv 0$ and $U_{A}(\mathbf{x})$ is the quasi-potential (i.e, the solution of Eq. (2)). We refer an interested reader to Refs. [5] and [11] where the concept of viscosity solution was introduced and the questions of existence and uniqueness of solution of Eq. (3) were investigated respectively.

1.1 Background

Minimum Action Paths (MAPs), i.e., minimizers of the geometric action $S(\psi)$ ((Eq. (2)) can be found numerically by path-based methods. The Geometric Minimum Action Method (GMAM) [8, 9] and the Adaptive Minimum Action Method (AMAM) [21, 22] iteratively update paths connecting given initial and final points starting from user-provided initial guesses so that the paths approach a local minimizer of the geometric action $S(\psi)$ (Eq. (2)) and a local minimizer of the original Freidlin-Wentzell action (Eq. (A-1)) respectively. The advantage of the path-based methods is that they are suitable for high- and even infinite-dimensional systems, i.e., stochastic partial differential equations (SPDEs) (e.g., [8, 7]). However, the found action minimizer is biased by the initial path and hence might not be the global minimizer, and it can be inexact due to slow convergence of the iterative process and numerical effects in the case if the actual MAP exhibits complex behavior. If the global action minimizer is found, the quasi-potential can be calculated along it.

Finding the quasi-potential in the whole region surrounding the attractor has important advantages over the search for MAPs. First, the quasi-potential allows us to estimate the invariant probability density near the attractor up to exponential order [6, 1]. Second, suppose a system has more than two attractors. Once the quasi-potential is computed with respect to an attractor $A$ in a large enough region, the most likely escape set $E$ from the basin of attraction of $A$ is automatically detected as the quasi-potential is constant along any trajectory running from $E$ to another attractor. The escape set $E$ can be a saddle point, a hyperbolic periodic orbit, or a more complex set of points in the phase space. Once the quasi-potential is computed, the global minimizer of the geometric action connecting the attractor $A$ with the escape set $E$ is guaranteed to be found by straightforward numerical integration [1] (see Section 3, Fig. 4 and its caption). On the contrary, path-based methods might have wrong endpoints, e.g., the most likely escape from $A$ is to the attractor $B$, while a path-based method is set up to seek action minimizers connecting attractors $A$ and $C$. Furthermore, even if the endpoints are identified correctly, path-based methods might converge to local minimizers.

On the other hand, it is clear that the computation of the quasi-potential in whole regions of the phase space is limited to low-dimensional systems.

To the best of our knowledge, the first attempt to compute the quasi-potential on a regular mesh was undertaken in [1]. It was done for 2D systems for both possible types of attractors, asymptotically stable equilibrium and stable limit cycle. The proposed numerical technique was an adjustment of the Ordered Upwind Method (OUM) [17, 18] for the case where the anisotropy coefficient $\Upsilon$ was unbounded. (What is the anisotropy coefficient will be explained right below.) The original OUM was designed for solving Hamilton-Jacobi equations of the form $\frac{1}{F\left(\mathbf{x},\hat{a}\right)}\|\nabla u(\mathbf{x})\|=1$, where $\hat{a}\equiv\frac{\nabla u}{\|\nabla u\|}$, with the boundary condition $u(A)=0$ where $A$ can be a point or a curve. The function $F(\mathbf{x},\hat{a})$, the speed of front propagation in the normal direction, was assumed to be bounded: $0<F_{\min}\leq F(\mathbf{x},\hat{a})\leq F_{\max}<\infty$. The Hamiltonian $\|\nabla u\|/F(\mathbf{x},\hat{a})$ was assumed to be convex, Lipschitz-continuous, and homogeneous of degree one in $\nabla u$. The anisotropy coefficient $\Upsilon$ is the ratio of the maximal and the minimal values of the speed function: $\Upsilon:=F_{\max}/F_{\min}$. One can think of $u(\mathbf{x})$ being the minimal possible travel time from $A$ to the point $\mathbf{x}$ in the following associated control problem. Suppose an astronaut on Mars (a 2D manifold with no roads) needs to minimize his traveltime from $A$ to $\mathbf{x}$. He can pick the direction of motion $\hat{q}$ of his rover at every moment of time. The speed of the rover is the function $f(\mathbf{y},\hat{q})$ of the position $\mathbf{y}$ and the direction of motion $\hat{q}$. The speed functions $f(\mathbf{y},\hat{q})$ and $F(\mathbf{x},\hat{a})$ relate via a Legendre transform. We refer an interested reader to Ref. [18].

The OUM [17, 18] employs the dynamical programming principle, i.e., at each step it solves the driver’s optimal control program by means of a finite difference upwind scheme for $\frac{1}{F\left(\mathbf{x},\hat{a}\right)}\|\nabla u(\mathbf{x})\|=1$. The role of the anisotropy coefficient $\Upsilon$ is very important. The value of $u$ at a mesh point $\mathbf{x}$ can be updated only from the mesh points lying within the distance of $\Upsilon h$ from it ($h=\max\{h_{1},h_{2}\}$, $h_{1}$, $h_{2}$ are the mesh steps in $x_{1}$ and $x_{2}$ respectively).

Eq. (3), the Hamilton-Jacobi equation for the quasi-potential, can be put in the form $\frac{1}{F\left(\mathbf{x},\hat{a}\right)}\|\nabla u(\mathbf{x})\|=1$:

[TABLE]

Hence, the front speed function is $F\left(\mathbf{x},\hat{a}\right)=\frac{1}{-2\mathbf{b}(\mathbf{x})\cdot\hat{a}}$. The corresponding speed function $f$ is the integrand in Eq. (2). It is clear that the speed function $F\left(\mathbf{x},\hat{a}\right)=\frac{1}{-2\mathbf{b}(\mathbf{x})\cdot\hat{a}}$ is unbounded at any direction $\hat{a}$ at the points $\mathbf{x}$ where $\mathbf{b}(\mathbf{x})=\mathbf{0}$, and at the directions $\hat{a}$ normal to $\mathbf{b}(\mathbf{x})$ at any point $\mathbf{x}$ where $\mathbf{b}(\mathbf{x})\neq\mathbf{0}$. This difficulty was overcome in [1] by the introduction of the update factor $K$ ($K$ is a positive integer) such that the value of $U_{A}$ at a mesh point $\mathbf{x}$ could be updated only from the mesh points lying within the distance of $Kh$ (i.e., within the update radius $Kh$) from it. The additional numerical error due to this adjustment was investigated [1] and shown to decay quadratically for a fixed update radius with mesh refinement. The OUM for finding the quasi-potential was incorporated into an R-package QPot available at cran.r-project.org [15, 16].

1.2 A brief summary of main results

In this paper, we present a new family of methods for computing the quasi-potential on a regular rectangular mesh. Like the OUM, these methods are based on the dynamical programming principle. Contrary to OUM, they do not use the upwind scheme for solving local minimization problems. Moreover, they completely abandon Hamilton-Jacobi Eq. (3) and refer only to the minimization problem (2). At each step, the local minimization problem is solved on the set of local straight line paths of length at most $Kh$111 Daisy Dahiya: ddahiya@math.umd.edu; Maria Cameron cameron@math.umd.edu;

Department of Mathematics, University of Maryland, College Park, MD 20742.

The Lagrangian $L(\psi,\psi^{\prime}):=\|\psi^{\prime}\|\|\mathbf{b}(\psi)\|-\psi^{\prime}\cdot\mathbf{b}(\psi)$ is integrated along them using quadrature rules. Therefore, we name this family of methods the Ordered Line Integral Methods. The names of the methods in this family reflect which quadrature rule is employed:

[TABLE]

We prove a theorem showing that the solution of the minimization problem in OLIM-R method is equivalent to the solution of the upwind scheme in the OUM. Our extensive numerical experiments with the other OLIMs show that they are significantly more accurate than the OUM. Our least squares fits to the error formula $E=Ch^{q}$ indicate that their error constants are up to 100 times smaller and the exponents $q$ are larger than those for OUM and OLIM-R. To make sure, OLIMs like OUM are at most first order due to the use of linear interpolation, but the major portion of the numerical error, the quadrature rule error, decays at least quadratically with the mesh refinement, which affects the exponents $q$ (makes them greater than 1) obtained by the least squares fits. Furthermore, we propose a CPU-time-saving implementation for OLIMs that reduces the number of calls of computationally expensive triangle update. As a result, the OLIM-R is about four times faster than the OUM. The other OLIMs are also faster than the OUM by some more modest factors. The graphs of the CPU times versus errors eloquently display that the OLIMs with second and higher order quadrature rules produce at least as accurate solutions as the OLIM-R (and hence the OUM) at smaller by several orders of magnitude CPU times.

So far, we have implemented the OLIMs in 2D. Our C codes OLIM_righthand.c, OLIM_midpoint.c, OLIM_trapezoid.c and OLIM_simpson.c are posted on M. Cameron’s web site [2]. A promotion of OLIMs to 3D in underway and will be reported elsewhere in the future.

The rest of the paper is organized as follows. The OLIMs are described in Section 2. The results of our numerical tests are presented in Section 3. The origin of numerical errors in the OLIMs are investigated in Section 4. The results of this work are summarized in Section 5.

2 Ordered Line Integral Methods

Throughout the rest of the paper, we will assume that the attractor $A$ of the system under consideration is fixed and omit the subscript $A$ in the notation for the quasi-potential with respect to $A$: $U_{A}(\mathbf{x})\equiv U(\mathbf{x})$. We consider a regular rectangular mesh with mesh steps $h_{1}$ and $h_{2}$ in $x_{1}$ and $x_{2}$ directions respectively, and set $h:=\max\{h_{1},h_{2}\}$. The mesh defines the sets of nearest neighbors for each mesh points. Every inner mesh point $P$ has eight nearest neighbors surrounding it as shown in Fig. 1

Similar to as it is done in the OUM [17, 18, 1], the mesh points are divided into the following four categories.

Unknown points: the points where the solution $U$ has not been computed yet, and none of its nearest neighbors is Accepted or Accepted Front.

• 1.

Considered: the points that have Accepted Front nearest neighbors. Tentative values of $U$ that might change as the algorithm proceeds, are available at them.

• 2.

Accepted Front: the points at which $U$ has been computed and no longer can be updated, and they have at least one Considered nearest neighbor.

• 3.

Accepted: the points at which $U$ has been computed and no longer can be updated, and they have only Accepted and/or Accepted Front nearest neighbors.

The general outline of the OLIMs coincides with the one of the OUM.

**Initialization

**Start with all mesh points being Unknown. Compute tentative values of $U$ at mesh points near the attractor $A$ and make them Considered.

**The main body

while** $\{$ the boundary of the mesh is not reached $\}$ and $\{$ the set of Considered points is not empty $\}$

       1: Make the Considered point $\mathbf{x}$ with the smallest tentative value of $U$ Accepted Front.

       2: Make all Accepted Front nearest neighbors of $\mathbf{x}$ that no longer have Considered nearest neighbors Accepted.

       3: Update all Considered points within the distance $Kh$ from $\mathbf{x}$ using $\mathbf{x}$ and maybe its nearest Accepted Front neighbors.

       4: Make all Unknown nearest neighbors of $\mathbf{x}$ Considered and compute tentative values of $U$ at them using the Accepted Front points lying within the distance $Kh$ from them.

end while

In the rest of this section, we will elaborate the initialization and steps 3 and 4 of the while-cycle.

2.1 Initialization

The way the initialization is performed is important for the accuracy of the computation of the quasi-potential. The two types of attractors in 2D, the asymptotically stable equilibrium point and the asymptotically stable limit cycle, correspond to the boundary condition $U(A)=0$ given at an initial point or curve respectively. The initialization procedure depends on whether the computation starts from the initial point or the initial curve.

2.1.1 Initialization from the initial point

Let $\mathbf{x}_{0}$ be an asymptotically stable equilibrium of $\dot{\mathbf{x}}=\mathbf{b}(\mathbf{x})$. Since $\mathbf{b}(\mathbf{x})=[b_{1}(\mathbf{x}),b_{2}(\mathbf{x})]^{T}$ is twice continuously differentiable and $\mathbf{b}(\mathbf{x}_{0})=\mathbf{0}$, we get

[TABLE]

in the neighborhood of $\mathbf{x}_{0}$. We calculate the quasi-potential for the linear approximation $\mathbf{b}(\mathbf{x})\approx\mathbf{A}(\mathbf{x}-\mathbf{x}_{0})$ at the four mesh points surrounding the point $\mathbf{x}_{0}$ (Fig. 2, Left) using the analytic formula [1] for the quasi-potential for a linear SDE with an asymptotically stable equilibrium. Let $A_{11}$, $A_{12}$, $A_{21}$, and $A_{22}$ be the entries of the matrix $\mathbf{A}$ in Eq. (5). Then the quasi-potential in the neighborhood of $\mathbf{x}_{0}$ is approximated by

[TABLE]

These four mesh points surrounding the point $\mathbf{x}_{0}$ become Considered.

2.1.2 Initialization from the initial curve

The initialization of the mesh points in the neighborhood of the initial curve $C$ is based on the observation that the gradient of the quasi-potential vanishes at $C$ (see Fig. 2, Center and Right). Assuming that $C$ is a smooth curve, there is a neighborhood $\mathcal{N}(C)$ of $C$ where the gradient of the quasi-potential $\nabla U(\mathbf{x})$ is a smooth vector field. Hence, in $\mathcal{N}(C)$, the following decomposition into two smooth vector fields takes place:

[TABLE]

One can easily deduce this decomposition using Eq. (3). Since the rotational component $\mathbf{l}(\mathbf{x})$ is smooth, its direction in $\mathcal{N}(C)$ changes continuously and can be approximated by the direction of $\mathbf{l}$ at the orthogonal projection ${\rm Proj}_{C}\mathbf{x}$ of $\mathbf{x}$ onto the curve $C$. Taking into account that $\nabla U(\mathbf{x})$ vanishes at $C$ and therefore

[TABLE]

we approximate the gradient of the quasi-potential in $\mathcal{N}(C)$ by

[TABLE]

The second term in the right-hand side of Eq. (12) is the projection of $\mathbf{b}(\mathbf{x})$ onto $\mathbf{b}({\rm Proj}_{C}\mathbf{x})$. The curve $C$ is provided by the user by a set of points $\{\mathbf{x}_{i}\}_{i=0}^{M-1}$. It is assumed that $\{\mathbf{x}_{i},\mathbf{x}_{i+1}\}$, $i=0,\ldots,M-2$, and $\{\mathbf{x}_{M-1},\mathbf{x}_{0}\}$ are consecutive points along the curve. We define ${\sf next}(i)=i+1$ for $0\leq i\leq M-2$, and ${\sf next}(M-1)=0$. For each pair of the consecutive points along the curve $(\mathbf{x}_{i},\mathbf{x}_{{\sf next}(i)})$, we define the smallest rectangle containing $\mathbf{x}_{i}$ and $\mathbf{x}_{{\sf next}(i)}$ with sides lying on the mesh lines and denote it by ${\sf Rectangle}(\mathbf{x}_{i},\mathbf{x}_{{\sf next}(i)})$. A collection of such rectangles is shown in Fig. 2, Center, by red contours. We will initialize all mesh points lying in the union of the rectangles

[TABLE]

Each mesh point $\mathbf{x}$ located within ${\sf Neib}(C)$ is initialized using Simpson’s quadrature rule as follows. First, observe that for any pair of points $\mathbf{x}_{i}$ and $\mathbf{x}_{f}$ we have

[TABLE]

Let us take $\mathbf{x}_{f}\equiv\mathbf{x}$, the point to be initialized, and $\mathbf{x}_{i}={\rm Proj}_{C}\mathbf{x}$, the projection of $\mathbf{x}$ onto the curve $C$. Note that $U({\rm Proj}_{C}\mathbf{x})=0$ and $\nabla U({\rm Proj}_{C}\mathbf{x})=0$. An approximation $\mathbf{x}^{\ast}$ to ${\rm Proj}_{C}\mathbf{x}$ is constructed as follows. Let $\mathbf{x}_{1}$ be the closest point of the curve $C$ to the point $\mathbf{x}$, and the point $\mathbf{x}_{2}$ be the neighboring point of $\mathbf{x}_{1}$ along the curve $C$ such that the angle between the vectors $\mathbf{x}-\mathbf{x}_{1}$ and $\mathbf{x}_{2}-\mathbf{x}_{1}$ does not exceed $\pi/2$ (Fig. 2, Right). The projection ${\rm Proj}_{C}(\mathbf{x})$ is approximated by the orthogonal projection of $\mathbf{x}$ onto the line passing through $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$:

[TABLE]

To use Simpson’s rule for computing the integral in Eq. (14), we need to evaluate $U(\mathbf{x}^{\ast})$, $\nabla U(\mathbf{x}^{\ast})$, $\nabla U(\mathbf{x})$, and $\nabla U(\mathbf{x}_{m})$, where $\mathbf{x}_{m}$ is the midpoint between $\mathbf{x}$ and $\mathbf{x}^{\ast}$. Since $U=0$ and $\nabla U=0$ on $C$, $U(\mathbf{x}^{\ast})\approx 0$ and $\nabla U(\mathbf{x}^{\ast})\approx 0$. Next, we note that on the line segment $[\mathbf{x}^{\ast},\mathbf{x}]$, $\nabla U$ is nearly parallel to $(\mathbf{x}-\mathbf{x}^{\ast})$, hence $\nabla U\cdot(\mathbf{x}-\mathbf{x}^{\ast})\approx\|\nabla U\|\|\mathbf{x}-\mathbf{x}^{\ast}\|$. Finally, we approximate $\|\nabla U(\mathbf{x}_{m})\|$ and $\|\nabla U(\mathbf{x})\|$ using Eq. (12) and get

[TABLE]

All initialized mesh points become Considered.

2.2 The one-point updates and the triangle updates

The OUM [17, 18, 1] involves two types of updates: the one-point update and the triangle update. The OUM one-point update applied to (3) reads

[TABLE]

where $\mathbf{x}$ is the mesh point being updated and $\mathbf{x}_{0}$ is an Accepted Front point within the distance $Kh$ from $\mathbf{x}_{0}$. Eq. (16) is merely the integrand in Eq. (2) integrated along the line segment from $\mathbf{x}_{0}$ to $\mathbf{x}$ using the right-hand rectangle quadrature rule. The OUM triangle update is done by the upwind finite difference scheme applied to PDE (3). Its details are worked out in Section 2.3 below.

In the proposed OLIMs, the finite difference scheme as well as PDE (3) are completely abandoned. The one-point update and the triangle update are done as follows. Let $\mathcal{Q}$ be a basic quadrature rule, i.e., right-hand rectangle, midpoint, trapezoid, or Simpson’s. The OLIM one-point update results from the application of $\mathcal{Q}$ for the integration of the integrand of Eq. (2) along the line segment $[\mathbf{x}_{0},\mathbf{x}]$:

[TABLE]

I.e., we compute the value ${\sf Q}_{1pt}(\mathbf{x}_{0},\mathbf{x})$ and replace the current value of $U$ at $\mathbf{x}$ with ${\sf Q}_{1pt}(\mathbf{x}_{0},\mathbf{x})$ if and only if it is smaller than the current value. The OLIM triangle update at the mesh point $\mathbf{x}$ from the triangle $(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})$, where the points $\mathbf{x}_{1}$ and $\mathbf{x}_{0}$ are assumed to be nearest neighbors, results from solving the minimization problem

[TABLE]

where the vector field $\mathbf{b}$ is approximated with a linear field within the triangle $(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})$ for the purpose of the application of the quadrature rule. The minimization problem (18) is solved by taking the derivative of the function to be minimized, setting it to zero, i.e.,

[TABLE]

and then trying to find a root of $g$ in the interval $[0,1]$. For this purpose, we use the hybrid secant/bisection method [20, 19].

Depending on the quadrature rule used, we denote the methods OLIM-R, OLIM-MID, OLIM-TR, and OLIM-SIM (see Section 1). The details of the triangle update for each OLIM are worked out in Appendix B.

2.3 Equivalence of update rules in the OUM and OLIM-R

The one-point update in OLIM-R and the OUM are done according to Eq. (16). In this section, we show that the triangle updates in OLIM-R and the OUM, if successful, give identical results.

The triangle update in the OUM [1] is done as it is proposed in [18]. The quasi-potential $U(\mathbf{x})$ is assumed to be linear within the triangle $(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})$ and have Accepted values $u_{1}$ and $u_{0}$ at $\mathbf{x}_{1}$ and $\mathbf{x}_{0}$ respectively (see Fig. 13). Its value at $\mathbf{x}$ is to be found. The finite-difference scheme proposed in [17, 18] applied to Eq. (3) is derived from the observation that

[TABLE]

where $u=U(\mathbf{x})$, $u_{0}=U(\mathbf{x}_{0})$, and $u_{1}=U(\mathbf{x}_{1})$ are the values of $U$ at the points $\mathbf{x}$, $\mathbf{x}_{0}$, and $\mathbf{x}_{1}$ respectively. The vector field $\mathbf{b}$ within the triangle $(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})$ is approximated by the constant field $\mathbf{b}\equiv\mathbf{b}(\mathbf{x})$ ($\mathbf{b}$ evaluated at $\mathbf{x}$). Throughout the rest of this section, we will omit the argument of $\mathbf{b}$. Denoting the matrix on the right-hand side of Eq. (19) by $P$ and plugging $\nabla U$ into Eq. (3), we obtain the following quadratic equation for $u$:

[TABLE]

If Eq. (20) has a solution $u$ satisfying $u\geq\min\{u_{0},u_{1}\}$, it is subjected to the consistency test that the characteristic passing through the point $\mathbf{x}$ crosses the segment $[\mathbf{x}_{1},\mathbf{x}_{0}]$. The direction of the characteristic is given by $\mathbf{b}(\mathbf{x})+\nabla U$ [1]. The condition that it crosses the segment $[\mathbf{x}_{1},\mathbf{x}_{0}]$ is equivalent to the fact that $\mathbf{b}(\mathbf{x})+\nabla U=c_{0}(\mathbf{x}-\mathbf{x}_{0})+c_{1}(\mathbf{x}-\mathbf{x}_{1})$ is a linear combination with nonnegative coefficients $c_{0}$, $c_{1}$, i.e., the vector

[TABLE]

has nonnegative entries. If $u$ has passed the consistency test then $U(\mathbf{x})=\min\{u,U(\mathbf{x})\}$.

OLIM-R performs the triangle update by solving the following minimization problem

[TABLE]

Like Eq. (20), Eq. (22) is set up under the assumption that the function $U(\mathbf{x})$ is linear within the triangle $(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})$.

The solution by the finite difference scheme (20) and the solution of the minimization (22) are equivalent in the sense specified by the following theorem.

Theorem 1**.**

Let $u$ be a solution of Eq. (20) such that the coefficients $c_{0}$ and $c_{1}$ in Eq. (21) are positive, i.e., the characteristic passing through the point $\mathbf{x}$ crosses the open interval $(\mathbf{x}_{1},\mathbf{x}_{0})$. Then $u$ is also the solution of the minimization problem (22) and the corresponding minimizer $s^{\ast}$ satisfies $0<s^{\ast}<1$.

Conversely, let $u$ be the solution of the minimization problem (22) with the corresponding minimizer $s^{\ast}$ satisfying $0<s^{\ast}<1$. Then $u$ is also a solution of Eq. (20) such that the characteristic passing through the point $\mathbf{x}$ crosses the open interval $(\mathbf{x}_{1},\mathbf{x}_{0})$.

The proof of Theorem 1 is given in Appendix C.

2.4 Reducing CPU time: a hierarchical update strategy

Clearly, the one-point updates require significantly fewer floating point operations than the triangle updates in both the OUM and the OLIMs. Therefore, we want to develop of a strategy for reducing the number of calls for the triangle update. Suppose that the mesh is fine enough so that $\mathbf{b}$ is approximately constant within the update radius $Kh$, and level sets of $U$ can be approximated by line segments within the distance $Kh$ from the mesh point $\mathbf{x}$ to be updated.

Motivated by this, we consider the following idealized situation. Let $\mathbf{b}$ be constant and the level set $U=U_{0}$ be a straight line. Let $\mathbf{x}$ be a point at the distance $h$ from this level set as shown in Fig. 3 (a,Top). Then

[TABLE]

where $t$ is the axis along the level set, and $f(t)$ is the update function

[TABLE]

The angles $\alpha\in[0,\pi/2)$ and $\theta\in(-\pi/2,\pi/2)$ are defined as shown in Fig. 3 (a, Top). Differentiating $f$ with respect to $\theta$, setting its derivative to zero, and then returning to the variable $t$, we find that $f$ has a single minimum on $-\infty<t<\infty$ achieved at $t=h\tan(\alpha)$ (see Fig. 3 (a, Bottom)).

Now we address the case where the quasi-potential $U$ is not differentiable along some curve, and hence its level sets have kinks. Typically, the vector field at the kink curve is directed approximately toward the kink. For example, see the quasi-potential computed in Ref. [1] for the Maier-Stein model [13]. To account for such a case, we consider the following idealized situation. Let $\mathbf{b}$ be constant, and the level set $U=U_{0}$ consists of two line segments. Let $t$ be the arclength parameter along the level set, $t=0$ at the kink, and the angles $\alpha$ and $\theta$ be defined as shown in Fig. 3 (b, Top). Let the point $\mathbf{x}$ lie on the kink line. The update function $f(t)$ for this case is obtained from the one in Eq. (23) using an appropriate reflection. Its graphs for three values of $\alpha$ are shown in Fig. 3 (b, Bottom). It has two local minima at $t=\pm 2h\tan(\alpha)$.

These considerations suggest that the update values of $U$ at a Considered point $\mathbf{x}$ lying near the Accepted Front as functions of an arclength parameter along the level set of $U$ approximating the Accepted Front will behave similar to the graphs in Fig. 3(a, Bottom, or b, Bottom), i.e., have a single local minimum, or two equal local minima. Therefore, we have rationalized the update procedures as follows. When step 4 of the main body of the algorithm is performed, i.e., an Unknown nearest neighbor $\mathbf{x}$ of the new Accepted Front point $\mathbf{x}^{\ast}$ becomes Considered, we compute tentative values at $\mathbf{x}$ from every Accepted Front point at distance at most $Kh$ from $\mathbf{x}$ using only the cheap one-point update. Suppose the minimal tentative value at $\mathbf{x}$ by the one-point update has been computed from the Accepted Front point $\mathbf{x}_{0}$, i.e.,

[TABLE]

Only then the triangle updates ${\sf Q}_{\Delta}(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})$ are called for all Accepted Front nearest neighbors $\mathbf{x}_{1}$ of $\mathbf{x}_{0}$ (if any). Then

[TABLE]

The symbol $\mathcal{N}(\mathbf{x}_{0})$ denotes the set of the nearest neighbors of $\mathbf{x}_{0}$ (eight nearest neighbors for every inner point as shown in Fig. 1). We will refer to this update strategy as the hierarchical update strategy.

Our numerical experiments in Section 3 below indicate that this strategy reduces CPU time in OLIM-R in comparison with the OUM by the factor of approximately 4, while the numerical errors in them either coincide, or the ones in OLIM-R exceed the ones in the OUM by less than 1%.

3 Numerical tests

We compare performances of the OLIMs (OLIM-R, OLIM-MID, OLIM-TR, and OLIM-SIM) and the original OUM-based [1] quasi-potential solver on two examples for which analytic formulas for the quasi-potential are available. These two examples are quite challenging from the computational point of view due to the large rotational components of $\mathbf{b}$ in comparison with $\nabla U$ and large curvatures of their MAPs (Minimum Action Paths). Furthermore, in the second example, the quasi-potential grows as a fourth degree polynomial.

A linear SDE. For the linear SDE [1]

[TABLE]

the quasi-potential with respect to the origin, the asymptotically stable equilibrium of the corresponding deterministic system, is the quadratic function

[TABLE]

The parameter $a$ is the quotient of the magnitudes of the rotational and the potential components ( $\mathbf{l}$ and $-\tfrac{1}{2}\nabla U$ respectively) of the vector field $\mathbf{b}$:

[TABLE]

We set $a=10$, i.e., the rotational component exceeds the potential component in magnitude by the factor of 10 that makes the computation of the quasi-potential challenging. The computational domain for this example is the square $[-1\leq x_{1}\leq 1]\times[-1\leq x_{2}\leq 1]$. The Minimum Action Path for SDE (24) arriving at the point $(0,0.9)$ is shown in Fig. 4(a).

An SDE with a limit cycle. The unit circle $C:=\{(x_{1},x_{2})\in\mathbb{R}^{2}~{}|~{}x_{1}^{2}+x_{2}^{2}=1\}$ is the asymptotically stable limit cycle of the deterministic system corresponding to the SDE [10, 1]

[TABLE]

The quasi-potential with respect to the unit circle is given by the quartic function111Actually, in our codes, the maximal update length for the one-point update is $Kh$, while it is $Kh+\sqrt{h_{1}^{2}+h_{2}^{2}}$ for the triangle update.

[TABLE]

11footnotetext: There is an error in Eq. (89) in [1]. It should be $U=\tfrac{1}{2}(r^{2}-1)^{2}$.

The decomposition of the vector field into the potential and the rotational components is

[TABLE]

We set the computational domain for this example to be the square $[-2\leq x_{1}\leq 2]\times[-2\leq x_{2}\leq 2]$. The Minimum Action Paths for SDE (24) arriving at the points $(0,1.9)$ and $(0,0.1)$ are shown in Fig. 4(b).

We test the OLIMs and the OUM on SDEs (24) and (27). The mesh sizes are $N\times N$ where $N=2^{p}$, $p=7,8,\ldots,12$. The update factor $K$ varies from $K=1$ to $K=50$. We have measured the maximum absolute error, the RMS error, and the CPU time for each mesh size $N$ and for each value of $K$.

3.1 Dependence of the accuracy on the update factor $K$

The choice of the optimal value of the update factor $K$ is a subtle issue. On one hand, if the Minimum Action Paths (MAPs) for a considered SDE would be straight lines, large $K$ would enable all mesh points to be updated from the correct triangles and hence enhance the accuracy. However, typically, MAPs are not straight lines. Furthermore, large $K$ leads to numerical integration by simple quadrature rules along long linear segments and hence increases the integration error. Finally, large $K$ increases the CPU time. As a result, the optimal $K$ should be not too large and not too small.

The results of our measurements are shown in Fig. 5. For all methods, and all mesh sizes $2^{p}\times 2^{p}$, we plot the maximal absolute error versus $K$ for SDEs (24) and (27). Note that the optimal values of $K$ are larger for SDE (24) than those for SDE (27).

Based on our plots, we propose the following Rules-of-Thumb for choosing $K$. The points of the graphs corresponding to the proposed Rules-of-Thumb are marked on the graphs by large dots.

The Rule-of-Thumb for the OUM and OLIM-R. For an $N\times N$ mesh where $2^{7}\leq N\leq 2^{12}$, pick

[TABLE]

The Rule-of-Thumb for OLIM-MID, OLIM-TR, and OLIM-SIM. For an $N\times N$ mesh where $2^{7}\leq N\leq 2^{12}$, pick

[TABLE]

3.2 Comparison of OLIM-MID, OLIM-TR, and OLIM-SIM to OLIM-R and the OUM

Our comparison of OLIM-MID, OLIM-TR, and OLIM-SIM to OLIM-R and the OUM is unambiguously in favor of the former ones. The comparison is conducted using the update factor $K(p)=p-3$ according to the Rule-of-Thumb for the OUM and OLIM-R, i.e., in the way benefiting the OUM and OLIM-R rather than the OLIMs with higher order quadrature rules.

The graphs of the maximum absolute error and the CPU time as functions of $N$ respectively are shown in Fig. 6 (a,b,c,d). The least squares fits to the formulas $E=CN^{-q}$ and $T=CN^{q}$ for the maximum absolute errors and the CPU times respectively are given in Table 1. Fig. 6 (a,b) and Table 1 (Columns 2 and 4) show that OLIM-MID, OLIM-TR, and OLIM-SIM are 100 to 1000 times more accurate then OLIM-R and the OUM for values of $K$ optimized for the latter methods. Fig. 6 (c,d) and Table 1 (Columns 3 and 5) show that OLIM-MID, OLIM-TR, and OLIM-SIM are faster than the OUM. OLIM-MID is faster than the OUM at least by the factor of 1.5. OLIM-R is faster than the OUM at least by the factor of 3. The CPU time (in seconds) versus the update factor $K$ for $N=1024$ is plotted for the OLIMs and the OUM in Fig. 7. These plots illustrate the advantage of our time-saving update strategy.

The CPU time as the function of the maximum absolute error is plotted in Fig. 6. It is clear that OLIM-MID, OLIM-TR, and OLIM-SIM are significantly better methods than OLIM-R and the OUM in terms of the balance between the accuracy and the CPU time.

3.3 Effects of the hierarchical update strategy on numerical errors in OLIM-R

The one-point update and the triangle update in OLIM-R and the OUM are equivalent. Therefore, if the CPU-saving hierarchical update strategy would not be implemented in OLIM-R, the numerical errors produced by these methods would coincide (in the exact arithmetics). However, the hierarchical update combined with numerical errors might lead to non-identical numerical solutions by the OUM and OLIM-R. Figs. 6(a) and (b) show that the maximal errors in these two methods are very close (the curves visually coincide). More detailed data for the maximal and RMS errors and CPU times in the OUM and OLIM-R are displayed in Tables 2 and 3 for SDEs (24) and (27) respectively. They indicate that the hierarchical update, in some cases, might increase the numerical error, but this increase is negligible. On the other hand, the CPU times in OLIM-R are approximately 4 times smaller.

3.4 Comparison of OLIM-MID, OLIM-TR, and OLIM-SIM to each other

The comparison of OLIM-MID, OLIM-TR, and OLIM-SIM is conducted using $K(p)=10+4(p-7)$, i.e., according to the Rule-of-Thumb for OLIM-MID, OLIM-TR, and OLIM-SIM. Fig. 8 shows that all these three methods are quite close in accuracy. Fig. 7 indicates that OLIM-SIM is somewhat slower than OLIM-MID and OLIM-TR. The least squares fits to the formulas $E=CN^{-q}$ for the maximum absolute errors and the RSM errors as functions of $N$ ($N=2^{p}$) are given in Table 4. OLIM-MID gives the best results on SDE (24) while OLIM-TR and OLIM-SIM challenge it on SDE (27) for rough meshes.

We favor OLIM-MID method. For fine meshes, it has the best balance between the accuracy and the CPU time among all methods considered in this work.

The error plots for each of the methods for SDEs (24) and (27) for the $256\times 256$ meshes and $K=11$ (corresponds to the Rule-of-Thumb) are shown in Fig. 9. We will return to the discussion on the error distributions in Section 4 below. Now we just note that these error plots are consistent with the results of Section 4. As we will show, the midpoint quadrature rule tends to underestimate the line integral, the trapezoid rule tends to exaggerate it, and the Simpson rule can make errors of either sign.

4 Origin of errors in the OLIMs

In this section, we discuss various factors which contribute to the error of the numerical solution by the OLIMs.

4.1 Errors of quadrature rules

In this Section, we discuss quadrature rule errors applied to the integral along a straight line segment connecting the points $\mathbf{x}_{0}$ and $\mathbf{x}_{1}$

[TABLE]

where

[TABLE]

Let $\mathbf{b}_{0}$ be the vector field at point $\mathbf{x}_{0}$. Assuming that $l$ is sufficiently small, we approximate $\mathbf{b}(\mathbf{x}(t))$ by its Taylor expansion around $\mathbf{x}_{0}$:

[TABLE]

where

[TABLE]

Then the integral in Eq. (32) becomes

[TABLE]

Approximating it using the quadrature rule $\mathcal{Q}$ with the error $E_{\mathcal{Q}}$, we obtain

[TABLE]

The errors of the right-hand, midpoint and trapezoid quadrature rules involve $f^{\prime}$, $f^{\prime\prime}$, and $f^{\prime\prime}$ (see e.g. [4]). For $f(t)$ given by Eq. (33), one can easily evaluate these derivatives at $t=0$. The resulting error estimates for the line integrals are

[TABLE]

where $\mathcal{R}$, $\mathcal{M}$, and $\mathcal{T}$ denote the right-hand, the midpoint, and the trapezoid basic quadrature rules respectively, and

[TABLE]

Eqs. (35) and (38) show that the integration error in OLIM-R is $O(l^{2})$ and can be of either sign. Eqs. (36) and (37) indicate that the integration errors in OLIM-MID and OLIM-TR are $O(l^{3})$. The first term of $f^{\prime\prime}(0)$ (Eq. (39)) is due to the linear part of $\mathbf{b}$. It is nonnegative due to the Schwarz inequality. The second term of $f^{\prime\prime}(0)$ is due to the nonlinearity of $\mathbf{b}$. It can have an arbitrary sign. Therefore, the contribution to the error due to the linear part of $\mathbf{b}$ is nonpositive for the midpoint rule, and nonnegative for the trapezoid rule. This is consistent with the error plots in Fig. 9(a) and (c) for the linear SDE.

Simpson’s quadrature rule has an error $O(l^{5})$:

[TABLE]

It is easy to calculate the contribution to $f^{(4)}(0)$ due to the linear part of $\mathbf{b}$:

[TABLE]

It can have an arbitrary sign. It is clear that the integration error in OLIM-SIM is small, and the total error of the numerical solution by OLIM-SIM, which is comparable to those by OLIM-MID and OLIM-TR (see Sec. 3.4), is due to the other factors discussed in Sections 4.3 and 4.2.

The lengths $l$ in Eqs. (35)-(37) and (40) can be recorded while the numerical solution is computed. For each mesh point $\mathbf{x}$ we define the update length as the distance

[TABLE]

for the one-point update, and

[TABLE]

for the triangle update, where $s^{\ast}$ is the solution of the corresponding minimization problem in Eq. (18). We update $l(\mathbf{x})$ every time when the value of the quasi-potential at $\mathbf{x}$ is updated. The update lengths $l$ in the numerical solutions by OLIM-MID on $256\times 256$ mesh with $K=11$ for SDEs (24) and (27) are shown in Fig. 10 (a) and (b) respectively. The OLIMs are designed so that the update length $l$ can be at most $(K+\sqrt{2})h$ which is $0.097$ and $0.194$ for the cases in Figs. 10 (a) and (b) respectively. For SDE (24), the update length is close to its maximal value at a significant fraction of mesh points. For SDE (27), it is close to its maximum in the outer neighborhood of the unit circle from which the computation starts. For the other OLIMs, the update lengths are similar. Therefore, the reduction of the integration error over rather long line segments by the use of second and higher order quadrature rules significantly improves the accuracy. This is exactly what we observe in Figs. 5 and 6.

4.2 Errors due to the curvatures of the level sets of the quasi-potential and MAPs

The OLIMs involve two contributions to their numerical errors due to the approximations of curved segments with straight line segments (Fig. 11). First, in both, the one-point and the triangle updates, curved segments of MAPs of lengths at most $Kh(1+o(1))$ and $(Kh+\sqrt{h_{1}^{2}+h_{2}^{2}})(1+o(1))$ respectively are approximated by straight line segments. Second, in the triangle update, curved segments of length at most $\sqrt{h_{1}^{2}+h_{2}^{2}}(1+o(1))$ along which the quasi-potential changes linearly, are approximated by straight line segments.

Let $l$ be the length of the line segment in the one-point update or in the triangle update approximating the MAP segment. Suppose that there is a MAP connecting the endpoints of this line segment, and its curvature is $\kappa>0$. Assume that $l$ is small and $\kappa$ is constant. Then the length of the MAP segment is

[TABLE]

OLIM-MID and OLIM-SIM involve the evaluation of $\mathbf{b}$ at the midpoint of the line segment. The midpoint of the MAP segment is at distance $d$ from it:

[TABLE]

The norm of the difference of $\mathbf{b}$ evaluated at these two midpoints is proportional to $d$. The curvature of MAPs can blow up near the equilibria of $\mathbf{b}$. As an example, the graphs of the curvature of the MAPs versus the quasi-potential along them for SDEs (24) and (27) are plotted in Fig. 12.

The bases of the triangles $[\mathbf{x}_{1},\mathbf{x}_{0}]$ in the triangle update are at most of length $h\sqrt{2}$. In the triangle update, the quasi-potential along the line segment $[\mathbf{x}_{1},\mathbf{x}_{0}]$ is assumed to change linearly, while the actual curve along which the quasi-potential changes linearly, is not a straight line in general. For example, if $U(\mathbf{x}_{1})=U(\mathbf{x}_{0})$, such a curve is the corresponding level set. For SDEs (24) and (27), such level sets are ellipses and circles respectively. Let $\mu>0$ be the curvature of the curve connecting $\mathbf{x}_{1}$ and $\mathbf{x}_{0}$ along which the quasi-potential is a linear function. Then the length of the segment of MAP connecting the point $\mathbf{x}$ to be updated with this curve is shorter (as in Fig. 11) or longer than $l_{MAP}$ in Eq. (42) by at most

[TABLE]

In summary, the local error of the triangle update due to approximating curved segments with line segments is bounded by

[TABLE]

where $J$ is the Jacobian matrix of $\mathbf{b}$. Eq. (45) explains relatively large numerical errors near the origin by OLIM-MID, OLIM-TR, and OLIM-SIM applied to SDE (27) (see Fig. 9 (b), (d), and (f)).

The comparison of $E_{curve}$ and the integration errors of the OLIMs (Eqs. (35), (36), (37), and (40)), combined with the results of our numerical tests in Section 3, show that the OLIMs are at most first order accurate regardless of the quadrature rule used. However, for a wide range of reasonable mesh sizes, OLIM-MID, OLIM-TR, and OLIM-SIM might appear to have superlinear convergence (see the least squares fits in Table 4) due to the change of the relative magnitudes of error terms of different orders. On the contrary, the OLIM-R exhibits sublinear convergence (see the least squares fits in Table 1).

4.3 Finite update radius

The direct application of the OUM developed in [17, 18] to the Hamilton-Jacobi PDE (3) for the quasi-potential would require an infinite update radius. The OUM adjusted for computing the quasi-potential [1] as well as the OLIMs set the update radius to $Kh$ where the update factor $K$ is a finite positive integer, and $h=\max\{h_{1},h_{2}\}$ is the mesh step. The extra error due to the finite update radius may appear in the case where the angle between the normal to the level set of $U$ approximating the Accepted Front and the vector field $\mathbf{b}$ at the point $\mathbf{x}$ to be updated is close to $\pi/2$ (the angle $\alpha$ in Fig. 3 (a, Top)). This extra error was quantified in [1] and shown not to exceed $0.5\|\mathbf{b}(\mathbf{x})\|h(K^{-1}+O(K^{-3}))$ for sufficiently small $h$.

To minimize this extra error, one might be tempted to use large values of the update factor $K$. However, there is a trade-off. The quadrature errors are approximately proportional to $h^{p}K^{p}$ where $p=2,3,$ or 5 depending on the quadrature rule used (see Section 4.1), and the errors due to the curvature of MAPs and level sets of the quasi-potential are approximately $C_{1}h^{2}K^{2}+C_{2}h^{3}K^{3}$. Hence, the optimal $K$ is the minimizer of the local error function of the form

[TABLE]

where the coefficients $A$, $B$, $C$ and $D$ depend on the vector field $\mathbf{b}$ and its derivatives, on the direction and the curvature of the MAPs, and on the curvature of the level sets of the quasi-potential. They change from one mesh point to another, and can have different orders of magnitude. Differentiating $E_{h}$ with respect to $K$ we get the following equation for the optimal $K$:

[TABLE]

In two simple cases, where $B\ll Ch$ or $Ch\ll B$, the optimal $K$ is $O(h^{-1/3})$ or $O(h^{-1/2})$ respectively, i.e., $O(N^{1/3})$ or $O(N^{1/2})$. In the general case, the dependence of the optimal $K$ on $N$ is more complicated, however, it is clear that $K$ must grow with $N$ but not faster that $N^{1/2}$. Our numerical experiments suggest simple Rules-of-Thumb for choosing $K$ (see Section 3.1). It turns out that $K$ does not need to be very large.

5 Conclusion

We have introduced the family of Ordered Line Integral Methods (OLIMs) for computing the quasi-potential on a rectangular mesh. The update rules in OLIM-R, employing the right-hand rectangle quadrature rule, are equivalent to those in the OUM. OLIM-MID, OLIM-TR, and OLIM-SIM (employing the midpoint, the trapezoid, and Simpson’s quadrature rules) converge faster and admit errors two to three orders of magnitude smaller than OLIM-R and the OUM (See Fig. 5). Nevertheless, asymptotically, the OLIMs, like the OUM, are at most first order accurate due to the use of linear interpolation. While the use of second or higher order quadrature rules requires the use of a nonlinear solver and hence makes the triangle update more expensive than the one in the OUM, the proposed hierarchical update strategy for reducing the number of calls for the triangle update makes the OLIMs faster than the OUM. Due to this strategy, OLIM-R is about four times faster than the OUM. OLIM-MID and OLIM-TR are about 1.5 times faster than the OUM for the optimal $K$’s and about three times faster for $K=50$ (see Fig. 7). Our C codes OLIM_righthand.c, OLIM_midpoint.c, OLIM_trapezoid.c and OLIM_simpson.c implementing the corresponding OLIMs are available on M. Cameron’s website [2].

We have investigated the dependence of the error upon the update factor $K$ and concluded that, while large $K$ nearly eliminates the error in the direction of the MAP due to insufficient update radius in the case of a slowly changing vector field $\mathbf{b}$, it increases the integration error and the error due to the curvature of the MAP. Based on our study of the relationship between the numerical error and the update factor $K$, we have proposed Rules-of-Thumb for choosing $K$.

Our comparison of the OLIMs with second or higher order quadrature rules shows that the best one of them in terms of the balance between the accuracy and the CPU time is achieved by OLIM-MID.

The present work is focused on the 2D case with isotropic diffusion. The OLIMs can be extended to the case with anisotropic diffusion and to 3D. We will report these developments in the future.

Acknowledgements

We thank Professor A. Vladimirsky for a valuable discussion. This work was supported in part by the NSF grant DMS1554907.

Appendix A. The Freidlin-Wentzell action vs the Geometric action

The Freidlin-Wentzell action functional for SDE (1) is defined on the set of absolutely continuous paths $\phi(t)$ by [6]

[TABLE]

The original definition of the quasi-potential [6] with respect to a compact set $A$ (an attractor of $\dot{\mathbf{x}}=\mathbf{b}(\mathbf{x})$) at a point $\mathbf{x}$ is

[TABLE]

The minimization with respect to the travel-time $T$ can be performed analytically [6, 8, 9] resulting at the geometric action $S(\psi)$. Let $\phi(t)$ be a fixed absolutely continuous path $\phi(t)$. Expanding $\|\cdot\|^{2}$ in Eq. (A-1) and using the inequality $y^{2}+z^{2}\geq 2yz$ for all nonnegative real numbers $y$ and $z$, we get:

[TABLE]

The inequality in Eq. (A-3) becomes an equality if and only if $\|\dot{\phi}\|=\|\mathbf{b}(\phi)\|$. Let $\chi$ be the path obtained from $\phi$ by a reparametrization such that $\|\dot{\chi}\|=\|\mathbf{b}(\chi)\|$. Then

[TABLE]

Note that $T_{\chi}$ can be infinite. The integral in right-hand side of Eq. (A-5) is invariant with respect to the parametrization of the path $\chi$. Hence, we can pick the most convenient one, for example, the arclength parametrization, and denote the reparametrized path by $\psi$. Hence,

[TABLE]

where $L$ is the length of the paths $\chi$ and $\psi$ (corresponding to the same curve). For computation of the quasi-potential, it is more convenient to deal with the geometric action $S(\psi)$ than with the Freidlin-Wentzell action $S_{T}(\phi)$.

Appendix B. The triangle updates for the OLIMs

**OLIM-R

**OLIM-R performs the triangle update by solving the following minimization problem

[TABLE]

Taking the derivative of the function to be minimized

[TABLE]

with respect to $s$ and setting it to zero, we obtain the following equation for $s$:

[TABLE]

Regrouping terms and taking squares we obtain the following quadratic equation for $s$:

[TABLE]

We solve Eq. (B-3), select its root $s^{\ast}$, if any, on the interval $[0,1]$, and verify that it is also the root of Eq. (B-2). In the case of success, the triangle update returns

[TABLE]

Otherwise, it returns ${\sf Q}_{\Delta}(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})=+\infty$.

**OLIM-MID

**OLIM-MID performs the triangle update by solving the following minimization problem

[TABLE]

Taking the derivative of

[TABLE]

with respect to $s$ and setting it to zero, we obtain the following equation for $s$:

[TABLE]

The hybrid nonlinear solver [20, 19] is used for finding a root $s^{\ast}$ of Eq. (B-8) in the interval $[0,1]$. In the case of success, the triangle update returns

[TABLE]

Otherwise, it returns ${\sf Q}_{\Delta}(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})=+\infty$.

**OLIM-TR

**OLIM-TR performs the triangle update by solving the following minimization problem

[TABLE]

Taking the derivative of

[TABLE]

with respect to $s$ and setting it to zero we obtain the following equation for $s$:

[TABLE]

The hybrid nonlinear solver [20, 19] is used for finding a root $s^{\ast}$ of Eq. (B-10) in the interval $[0,1]$. In the case of success, the triangle update returns

[TABLE]

Otherwise, it returns ${\sf Q}_{\Delta}(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})=+\infty$.

**OLIM-SIM

**OLIM-SIM performs the triangle update by solving the following minimization problem

[TABLE]

Taking the derivative of

[TABLE]

with respect to $s$ and setting it to zero, we obtain the following equation for $s$:

[TABLE]

The hybrid nonlinear solver [20, 19] is used for finding a root $s^{\ast}$ of Eq. (B-12) in the interval $[0,1]$. In the case of success, the triangle update returns

[TABLE]

Otherwise, it returns ${\sf Q}_{\Delta}(\mathbf{x}_{1},\mathbf{x}_{0},\mathbf{x})=+\infty$.

Appendix C. Proof of Theorem 1

Proof.

Without the loss of generality we assume that $\mathbf{x}_{1}$ is the origin.

Step 1. Show that $u$ is a solution of Eq. (20) if and only if $u-u_{1}=\|\mathbf{x}\|(U_{\xi}\cos(\alpha)+U_{\eta}\sin(\alpha))$ where (see Fig. 13) $\alpha$ ($0<\alpha<\pi$) is the angle between the vectors $\mathbf{x}_{0}$ and $\mathbf{x}$, $U_{\xi}=\|\mathbf{x}_{0}\|^{-1}(u_{0}-u_{1})$, and $U_{\eta}$ is a solution of

[TABLE]

which is Eq. (3) written in the $(\xi,\eta)$-coordinates at the point $\mathbf{x}$.

First observe that both Eqs. (20) and (22) are invariant with respect to translations. Therefore, we shift $\mathbf{x}_{1}$ to the origin as shown in Fig. 13 without changing their solutions.

Second, Eq. (20) is invariant with respect to orthogonal transformations. Indeed, the multiplication of $\mathbf{x}$ and $\mathbf{x}_{0}$ by an orthogonal matrix $O$ converts Eq. (19) to

[TABLE]

Hence the matrix $P$ in Eq. (20) changes to $PO^{T}$ and $\mathbf{b}$ becomes $O\mathbf{b}$ leading to the equation

[TABLE]

which is equivalent to Eq. (20). Hence, we apply an orthogonal transformation to map the original coordinate system onto the $(\xi,\eta)$ system in which $\mathbf{x}_{0}$ lies on the positive $\xi$-semiaxis and the $\eta$-coordinate of $\mathbf{x}$ is positive:

[TABLE]

where $\alpha$ ($0<\alpha<\pi$) is the angle between vectors $\mathbf{x}_{0}$ and $\mathbf{x}$ as shown in Fig. 13.

Finally, if $u$ is a solution of Eq. (20) then

[TABLE]

Hence, if $u$ is the solution of Eq. (20), then $U_{\xi}$ is exactly $(u_{0}-u_{1})/\|\mathbf{x}_{0}\|$ which shows that it is independent of $u$. Hence Eq. (20) can be rewritten as an equation Eq. (C-1) for $U_{\eta}$.

Step 2. Find geometric conditions guaranteeing the existence of solution(s) of Eq. (C-1) satisfying the consistency check and determine the selection rule if it has two solutions.

Eq. (3) implies that $\nabla U$ is orthogonal to $2\mathbf{b}+\nabla U$. Therefore, the locus of the vectors $\nabla U$ satisfying Eq. (3) is the circle [1] shown in Fig. 13. This circle passes through the origin and has center at the end of the vector $-\mathbf{b}$ originating from the origin. Since $\|\nabla U\|^{2}=U_{\xi}^{2}+U_{\eta}^{2}$, Eq. (3) has a solution if and only if the line normal to the $\xi$-axis and passing through the point $(U_{\xi},0)$ (the red dashed line in Fig. 13) intersects the circle. The MAP is collinear to the vector $\mathbf{b}+\nabla U$ [1]. The consistency condition requires that the MAP passing through the point $\mathbf{x}$ crosses the interval $[\mathbf{x}_{1},\mathbf{x}_{0}]$. This means that the angle between the vector $\mathbf{b}+\nabla U$ and the positive $\xi$-semiaxis should be not less than the angle $\alpha$ between the vector $\mathbf{x}-\mathbf{x}_{1}\equiv\mathbf{x}$ and the positive $\xi$-semiaxis, and not greater than the angle between the vector $\mathbf{x}-\mathbf{x}_{0}$ and the positive $\xi$-semiaxis. Drawing rays parallel to $\mathbf{x}$ and $\mathbf{x}-\mathbf{x}_{0}$ from the center of the circle and then dropping normals from their intersections with the circle to the $\xi$-axis as shown in Fig. 13, we obtain the interval on the $\xi$-axis where $U_{\xi}$ should belong in order to make the solution $U_{\eta}$ of Eq. (C-1) satisfy the consistency condition. This interval bounded by the endpoints of the thin brown and green-blue dashed lines in Fig. 13. Note that the consistency condition can be satisfied only by the larger root of Eq. (C-1), i.e., we should select the root

[TABLE]

Step 3. Find the solution of the minimization problem (22) and show that, if the minimizer $s^{\ast}\in(0,1)$ then it coincides with $u=\|\mathbf{x}\|(U_{\xi}\cos(\alpha)+U_{\eta}\sin(\alpha))$, where $U_{\xi}=(u_{0}-u_{1})/\|\mathbf{x}\|$ and $U_{\eta}$ is given by Eq. (C-4).

Consider the function to be minimized in Eq. (22) rewritten for $\mathbf{x}_{1}$ shifted to the origin:

[TABLE]

where $\gamma$ is the angle between the vectors $\mathbf{b}$ and $\mathbf{x}-s\mathbf{x}_{0}$. The point $s\mathbf{x}_{0}$, and hence the value of $s$, is uniquely determined by the angle $\gamma$ (Fig. 13):

[TABLE]

Moreover, since $\cot(\gamma-\beta)$ is a monotone function on the interval $0<\gamma-\beta<\pi$, there is a one-to-one correspondence between $-\infty<s<\infty$ and $\beta<\gamma<\beta+\pi$. Therefore, the function $f(s)=:F(\gamma(s))$, $\beta<\gamma<\beta+\pi$, where

[TABLE]

If $(s^{\ast},f(s^{\ast}))$ is a minimum of $f(s)$, then there is a unique minimum $(\gamma^{\ast},F(\gamma^{\ast})=f(s^{\ast}))$ of $F(\gamma)$.

Let us minimize $F(\gamma)$. Its derivative is given by:

[TABLE]

Setting it to zero, cancelling the positive constant $\|\mathbf{x}\|\sin(\alpha)/\sin^{2}(\gamma-\beta)$, regrouping the terms, and applying trigonometric formulas, we obtain the following equation for $\gamma$:

[TABLE]

Hence, the optimal angle $\gamma$ satisfies:

[TABLE]

Let us denote by $\gamma^{\ast}$ the solution of Eq. (C-9) lying in the interval $(\beta,\beta+\pi)$. To check whether $\gamma^{\ast}$ is a maximizer or a minimizer, we evaluate the second derivative of $F(\gamma)$ at $\gamma^{\ast}$ and find:

[TABLE]

as the angle $\gamma-\beta\in(0,\pi)$ by construction. Hence the optimal $\gamma$ is the minimizer of $F$. Next, we recall Eq. (3): $\|\nabla U\|+2\mathbf{b}\cdot\nabla U=0$. Adding $\|\mathbf{b}\|^{2}$ to both sides, we obtain $\|\nabla U+\mathbf{b}\|^{2}=\|\mathbf{b}\|^{2}$. Then Eq. (C-9) and the equality $\|\nabla U+\mathbf{b}\|=\|\mathbf{b}\|$ imply

[TABLE]

Therefore,

[TABLE]

On the other hand, from Eq. (C-9) we obtain:

[TABLE]

Plugging Eq. (C-13) into Eq. (C-12) we get

[TABLE]

which coincides with Eq. (C-4).

Finally, the solution of the minimization problem (22)

[TABLE]

is achieved either at $s^{\ast}$ if $0\leq s^{\ast}\leq 1$, or at the endpoints $s=0$ or $s=1$. Hence, if $0<s^{\ast}<1$, then the solution of the minimization problem (22) coincides with the one of the finite difference scheme (20), and the latter meets the consistency conditions. Conversely, the solution of the finite difference scheme (20) satisfying the consistency conditions coincides with the one of the minimization problem (22), and the corresponding minimizer $s^{\ast}\in[0,1]$.

∎