Single Pass Computation of First Seismic Wave Travel Time in Three Dimensional Heterogeneous Media With General Anisotropy

Abstract

We present a numerical method for computing the first arrival travel-times of seismic waves in media defined by a general Hooke tensor, in contrast with previous methods which are limited to a specific subclass of anisotropic media, such as “tilted transversally isotropic” (TTI) media or “tilted orthorhombic” (TOR) media in Waheed et al. (G 80: 49-58, 2015), Bouteiller et al. (GJlI 212: 1498-1522, 2017). Our method proceeds in a single pass over the discretized domain, similar to the fast marching method, whereas existing methods for these types of anisotropy require multiple iterations, similar to the fast sweeping method. We introduce a new source factorization model, making it possible to achieve third-order accuracy in smooth media. We also validate our solver by comparing it with the solution of the elastic wave equation in a 3D medium with general anisotropy.

Appendices

Construction of the Synthetic Test

We describe the synthetic test used in § 4.1 to validate the convergence order of the proposed numerical scheme. For that purpose, we need to introduce the following notations.

Definition A.1

Let ${\varOmega } \subset {{\mathbb {R}}}^d$ be a domain, equipped with a metric $N_x(v)$ (resp. dual-metric $N_x^*(p)$), $x \in {\varOmega }$, $p,v \in {{\mathbb {R}}}^d$. Their pull-back by a diffeomorphism $\phi : {\widetilde{{\varOmega }}} \rightarrow {\varOmega }$, with Jacobian matrix denoted by ${\varPhi }$, is defined as

$$\begin{aligned} {\widetilde{N}}_x(v)&:= N_{\phi (x)}\big ({\varPhi }(x) v\big ),&\Big (\text {resp. } {\widetilde{N}}^*_x(p)&:= N^*_{\phi (x)} \big ({\varPhi }^{-T}(x) p\big ).\Big ) \end{aligned}$$

By construction, the geometrical quantities defined in § 2.1 and associated with the metrics $N_x$ and ${\widetilde{N}}_x$ are closely related: the path-length ${\widetilde{{{\mathcal {L}}}}} (\gamma ) = {{\mathcal {L}}}(\gamma \circ \phi )$, and distance ${{\tilde{d}}}(x,y) = d(\phi (x),\phi (y))$, where $x,y \in {{\tilde{{\varOmega }}}}$ and $\gamma : [0,1] \rightarrow {\varOmega }$. Likewise, if $u : {\varOmega } \rightarrow {{\mathbb {R}}}$ obeys the eikonal equation (2), then so does $u \circ \phi : {\widetilde{{\varOmega }}} \rightarrow {{\mathbb {R}}}$ with respect to the pulled-back dual-metric ${\widetilde{N}}_x^*$, with the appropriate seed point and boundary conditions. In our experiments, we use for simplicity a metric $N_x = N_c$ defined by a constant field of Hooke tensors c, and a star-shaped domain ${\varOmega }$ with respect to the origin, which is chosen as the seed point; the eikonal equation on ${\varOmega }$ (resp. ${\widetilde{{\varOmega }}}$) therefore admits the following explicit solution, as announced in (40):

$$\begin{aligned} u(x)&= N_c(x),&\Big (\text {resp. } {{\tilde{u}}}(x) = N_c\big (\phi (x)\big ).\Big ) \end{aligned}$$

In general, the pull-back of a metric defined by a Hooke tensor is not defined by a Hooke tensor alone, and one has in addition to keep track of the Jacobian matrix both symbolically and numerically. A special case of interest arises however for conformal transformations, for which the Jacobian is a scaled rotation, and which thus preserves the metric structure. More precisely, let $x\in {\widetilde{{\varOmega }}}$ be fixed, assume that $N_{\phi (x)}^*$ is defined as in (1) by a Hooke tensor c, and that ${\varPhi }(x) = \lambda R$ is the product of a scaling $\lambda >0$ and of a rotation R. Then ${\widetilde{N}}_x^*$ is defined by the Hooke tensor of components

$$\begin{aligned} {{\tilde{c}}}_{i'j'k'l'} = \lambda ^{-2} \sum _{i,j,k,l} c_{ijkl} R_{ii'} R_{jj'} R_{kk'} R_{ll'} . \end{aligned}$$

Another benefit of conformal transformations is that they leave invariant the length distortion and angular width of the metric, $\mu ({\widetilde{N}}_x) = \mu (N_{\phi (x)})$ and $\mu ({\widetilde{N}}_x) = \mu (N_{\phi (x)})$, see Definitions 2.2 and 2.4. Three dimensional conformal transformations include dilations, translations, rotations, the inversion $x\in {{\mathbb {R}}}^3 \setminus \{0\} \mapsto x/\Vert x\Vert ^2$, and compositions of these.

In our experiments we use a ”special conformal transformation”, see (40, right) and Fig. 3, which is smooth except for a singularity at $b/\Vert b\Vert ^2$, where $b \in {{\mathbb {R}}}^3$ is a parameter. It is obtained as the composition of an inversion, a translation by $-b$, and another inversion. More precisely, we choose $b = (1/6,1/9,1/18)$ and let ${\widetilde{{\varOmega }}} := ]-1,1[^3$ with seed at the origin, so that the singular point $b/\Vert b\Vert ^2 \notin {\widetilde{{\varOmega }}}$, and the image domain ${\varOmega } := \phi ({{\tilde{{\varOmega }}}})$ is star shaped with respect to the origin, see Fig. 3. Besides, we use the Hooke tensors of the olivine and mica as defined in Table 1, with a constant rotation of Euler axis (2, 1, 3) and angle $3\pi /5$.

Proof of Proposition 2.7

We estimate in this appendix the quantity ${\varTheta }(N)$, which measures the angular distortion associated with a norm N on ${{\mathbb {R}}}^d$, in terms of its length distortion, as announced in Proposition 2.7. Different proof techniques are used in the elliptic and anelliptic cases.

1.1 Anelliptic Norms

The announced estimate, established in Corollary B.2, follows from upper and lower bounds on the gradient of a norm, presented in the next lemma.

Lemma B.1

Let N be a norm on ${{\mathbb {R}}}^d$, differentiable at $v\in {{\mathbb {R}}}^d \setminus \{0\}$. Then

$$\begin{aligned} \mu _*(N) \le \Vert \nabla N(v)\Vert \le \mu ^*(N). \end{aligned}$$

Proof

Since N is 1-homogeneous, one has $\langle \nabla N(v),v\rangle = N(v)$ by Euler’s identity, and therefore

$$\begin{aligned} \mu _*(N) \le \frac{N(v)}{\Vert v\Vert } = \frac{\langle \nabla N(v),v\rangle }{\Vert v\Vert } \le \Vert \nabla N(v)\Vert . \end{aligned}$$

On the other hand, for any vector w, one obtains using successively the convexity of N and the triangular inequality

$$\begin{aligned} \langle \nabla N(v),w\rangle \le N(v+w) -N(v) \le N(w). \end{aligned}$$

Choosing $w:=\nabla N(v)$ yields the announced upper estimate and concludes the proof:

$$\begin{aligned} \mu ^*(N) \ge \frac{N(\nabla N(v))}{\Vert \nabla N(v)\Vert } \ge \frac{\langle \nabla N(v), \nabla N(v)\rangle }{\Vert \nabla N(v)\Vert } = \Vert \nabla N(v)\Vert . \end{aligned}$$

$\square$

Corollary B.2

For any norm N on ${{\mathbb {R}}}^d$, differentiable on ${{\mathbb {R}}}^d \setminus \{0\}$, one has $\mu (N) \cos {\varTheta }(N) \ge 1$.

Proof

Using Lemma B.1 we obtain as announced

$$\begin{aligned} \cos {\varTheta }(N) = \frac{\langle v,\nabla N(v)\rangle }{\Vert v\Vert \Vert \nabla N(v)\Vert } = \frac{N(v)}{\Vert v\Vert } \frac{1}{\Vert \nabla N(v)\Vert } \ge \mu _*(N)/\mu ^*(N) = 1/\mu (N). \end{aligned}$$

$\square$

1.2 Elliptic Norms

The announced estimate, established in Corollary B.4, follows from a classical inequality in analysis, which proof is recalled in the next lemma.

Lemma B.3

(Weighted Pólya-Szegö inequality) Let $(\lambda _i)_{i=1}^d$ be positive numbers, and $(\mu _i)_{i=1}^d$ be non-negative, let $\lambda _* = \min \{\lambda _1,\cdots , \lambda _d\}$ and let $\lambda ^* := \max \{\lambda _1,\cdots ,\lambda _d\}$. Then

$$\begin{aligned} \sqrt{\sum _{1 \le i \le d} \mu _i \lambda _i} \sqrt{\sum _{1 \le i \le d} \frac{\mu _i}{ \lambda _i}} \le \frac{1}{2} \left( \sqrt{\frac{\lambda ^*}{\lambda _*}}+\sqrt{\frac{\lambda _*}{\lambda ^*}}\right) \sum _{1 \le i \le d} \mu _i. \end{aligned}$$

Proof

Without loss of generality, we may assume that $\sum _{1 \le i \le d} \mu _i = 1$, and denote $E[\gamma ] := \sum _{1 \le i \le d} \mu _i \gamma _i$ for any sequence $(\gamma _i)_{i=1}^d$. Observing that $E[ (\lambda ^* - \lambda ) (1/\lambda _* - 1/\lambda )] \ge 0$, and developping, we obtain

$$\begin{aligned} \frac{\lambda ^*}{\lambda _*} + 1 \ge E\left[ \frac{\lambda ^*}{\lambda }\right] + E\left[ \frac{\lambda }{\lambda _*}\right] \ge 2 \sqrt{ E\left[ \frac{\lambda ^*}{\lambda }\right] E\left[ \frac{\lambda }{\lambda _*}\right] }. \end{aligned}$$

The second inequality follows from the arithmetic-geometric mean inequality $\frac{a+b}{2} \ge \sqrt{ab}$, $\forall a,b\ge 0$. The announced result follows. $\square$

Corollary B.4

For any elliptic norm N one has $\frac{1}{2} (\mu (N) + \mu (N)^{-1} ) \cos {\varTheta }(N) =1.$

Proof

Without loss of generality, up to a rotation, we may assume that for all $v\in {{\mathbb {R}}}^d$

$$\begin{aligned} N(v)^2&= \sum _{1\le i \le d} \lambda _i v_i^2,&\text {thus} \quad N(v) \nabla N(v)&= (\lambda _i v_i)_{i=1}^d, \end{aligned}$$

where $\lambda _1,\cdots ,\lambda _d>0$. Denote $\lambda _* = \min \{\lambda _1, \cdots , \lambda _d\}$ and $\lambda _* = \max \{\lambda _1, \cdots , \lambda _d\}$, so that $\mu (N) = \sqrt{\lambda ^*/\lambda _*}$. Then letting $\mu _i := \lambda _i v_i^2$ one obtains by Lemma B.3

$$\begin{aligned} \frac{\Vert \nabla N(v)\Vert \Vert v\Vert }{\langle \nabla N(v),v\rangle } = \frac{\sqrt{\sum _i \lambda _i^2 v_i^2} \sqrt{\sum _i v_i^2}}{\sum _i \lambda _i v_i^2} = \frac{\sqrt{\sum _i \mu _i \lambda _i} \sqrt{\sum _i \mu _i/\lambda _i}}{\sum _i \mu _i} \le \frac{1}{2} \left( \sqrt{\frac{\lambda ^*}{\lambda _*}}+\sqrt{\frac{\lambda _*}{\lambda ^*}}\right) . \end{aligned}$$

(41)

This proves $\frac{1}{2} (\mu (N) + \mu (N)^{-1} ) \cos {\varTheta }(N) \ge 1$. Adequately choosing v turns (41) into an equality, which concludes the proof. More precisely, we may assume without loss of generality that $\lambda _* = \lambda _1$ and $\lambda ^* = \lambda _2$, and then choose $v = (\sqrt{\lambda _2},\sqrt{\lambda _1},0,\cdots , 0)$. $\square$

Sequential Quadratically Constrained Programming

The numerical implementation of our eikonal equation solver involves the solution to optimization problems of the form

$$\begin{aligned} \max \{ \langle p,v\rangle ;\ f(p) \le 0\}, \end{aligned}$$

(42)

where f is smooth and strongly convex, and the vector v is fixed. They arise in the definition of the norm (25), which is used in the source factorization (20), as well as in evaluation of the update operator on vertices and edges (38), with an additional linear constraint in the latter case. In order to solve (42), we use an approach known as Sequential Quadratically Constrained Quadratic Programming (SQCQP) [21], which basic principle is to solve a sequence of simplified problems obtained by replacing the objective function and the constraints with their second-order Taylor expansion. We provide two basic results that are sufficient for our application, and refer to [21] for more details on this rich theory. Our first observation, which proof is left to the reader, is that the problem (42) has a closed form solution when f is a suitable quadratic function.

Lemma C.1

(Maximization of a linear function over an ellipsoid) Let $f : {{\mathbb {R}}}^d \rightarrow {{\mathbb {R}}}$ be a quadratic function such that the set $\{f < 0\}$ is a non-empty ellipsoid, and let $p,v \in {{\mathbb {R}}}^d$. Then

$$\begin{aligned} F(p) := p+M(p) (\lambda (p) v-V(p)) \end{aligned}$$

(43)

is the unique solution to (42), where

$$\begin{aligned} V(p)&:= \nabla f(p),&M(p)&:= (\nabla ^2 f(p))^{-1},&\lambda (p)&:= \sqrt{\frac{\langle V(p), M(p) V(p)\rangle - 2 f(p)}{\langle v,M(p) v\rangle }}. \end{aligned}$$

(44)

For convenience, the solution (43) is expressed in terms of the Taylor expansion of the quadratic function f at a given but arbitrary point p. Note however that, if f is a quadratic function as assumed in Lemma C.1, then the matrix M(p) in (44, center) is independent of p, and the value of F(p) is independent of p since it solves (42). The basic SQCQP framework consists in repeatedly evaluating (43) with a non-quadratic function f, thus generating a sequence of points $p_{n+1} = F(p_n)$, $n \ge 0$. This yields an iterative method for the optimization problem (42), enjoying a quadratic (Newton-like) local convergence rate, as shown in Proposition C.2. Variants of this method enjoy a global convergence guarantee [21] under suitable assumptions, but in our numerical experiments the basic method was adequate.

Proposition C.2

Let $f : {{\mathbb {R}}}^d \rightarrow {{\mathbb {R}}}$ be $C^3$ smooth, and let $v \in {{\mathbb {R}}}^d$. Assume that $p_* \in {{\mathbb {R}}}^d$ and $\lambda _*>0$ are such that

$$\begin{aligned} f(p_*)&=0&\nabla f(p_*)&=\lambda _* v&\nabla ^2 f(p_*) \succ 0. \end{aligned}$$

(45)

Then $p_*$ is an isolated local maximum for the optimization problem (42). In addition there exists a constant $C>0$ such that, for any $p_0 \in {{\mathbb {R}}}^d$ close enough to $p_*$, the sequence defined by $p_{n+1} = F(p_n)$, see (43), satisfies for all $n \ge 0$

$$\begin{aligned} \Vert p_n - p_* \Vert \le C^{-1} (C \Vert p_0 - p_*\Vert )^{2^n}. \end{aligned}$$

(46)

Sketch of proof

We recognize in (45) the second-order optimality conditions for the constrained optimization problem (42).

A first-order Taylor expansion shows that $\lambda (p_*+h) = \lambda _* + {{\mathcal {O}}}(h^2)$, and then $F(p_*+h) = p_* + {{\mathcal {O}}}(h^2)$. The estimate (46) follows by induction on $n\ge 0$. $\square$

Remark C.3

(Exponential transformation, and numerical stability) Assume that the constraint in (42) takes the form $g \le 0$, where $g = \exp (\alpha f)-1$ is a strongly convex function, defined in terms of a smooth (but non-convex) f and a positive constant $\alpha$. One can check that $(g,\nabla g, \nabla ^2 g)$ is positively proportional to

$$\begin{aligned} {{\tilde{f}}}&:= (1-\exp (-\alpha f)) /\alpha ,&\nabla f,&\nabla ^2 f + \alpha \nabla f \nabla f^T. \end{aligned}$$

Note also that f and ${{\tilde{f}}}$ vanish at the same points, and if $p_*$ is such a point then ${{\tilde{f}}}(p_*+h)=f(p_*+h) + {{\mathcal {O}}}(\Vert h\Vert ^2)$ for small h. In the sequential quadratic iterations, see Proposition C.2, one may thus replace $(g,\nabla g,\nabla ^2 g)$ with $(f,\nabla f,\nabla ^2 f + \alpha \nabla f \nabla f^T)$ and preserve the local quadratic convergence (46). This eliminates all exponentials, to the benefit of numerical stability.

Monotony and Causality in Fixed Point Problems

In this section, we review the properties of the numerical scheme considered in this paper, and derive the following guarantees: existence and uniqueness of a fixed point, convergence of an iterative method to find it, and validity of the fast marching method subject to an acuteness condition. We also discuss how these properties transfer to the source factored and high order scheme variants. Closely related arguments can be found in the literature devoted to semi-Lagrangian discretizations of the eikonal equation [2, 10, 38, 39, 53, 60]. We fix the grid scale $h>0$ in this appendix, and refer to [10] for a convergence analysis to the PDE solution as it is refined. Denote $X := {\varOmega }_h \setminus \{x_*\}$ the discretization set (4) minus the source point, and let ${{\mathbb {U}}}:= {{\mathbb {R}}}^X$ be the set of mappings from X to ${{\mathbb {R}}}$. Recall that the objective is to find $u \in {{\mathbb {U}}}$ such that for all $x \in X$

$$\begin{aligned} {\varLambda } u(x)&= u(x),&\text {where } {\varLambda } u(x)&:= \min _{y \in \partial {{\mathcal {V}}}_h^x} {{\,\mathrm{{{{\mathcal {I}}}}}\,}}_h^x u(y) + N_x(x-y), \end{aligned}$$

(47)

where ${{\mathcal {I}}}_h^x$ denotes the piecewise linear interpolation operator, on a polytope ${{\mathcal {V}}}_h^x$ enclosing x with vertices in the grid $h{{\mathbb {Z}}}^d$, see §2.2. By convention in (47, right), $u\in {{\mathbb {U}}}$ is extended to $h {{\mathbb {Z}}}^d \setminus X$ by $u(x_*)=0$ and $u=+\infty$ elsewhere. We make the following connectedness assumption: for any $x_0 \in X$, one can find $n\ge 1$ and $x_1,\cdots ,x_n\in X$, such that $x_{i+1}$ is a vertex of ${{\mathcal {V}}}_h^{x_i}$, for all $i<n$, and $x_*$ is a vertex of ${{\mathcal {V}}}_h^{x_n}$. Given $u,v\in {{\mathbb {U}}}$, the strict inequality “$u< v$” stands for “$\forall x \in X, u(x) < v(x)$”; and likewise for weak inequality $u\le v$. Given $u \in {{\mathbb {U}}}$ and $\tau \in {{\mathbb {R}}}$ we define

$$\begin{aligned} u^{\le \tau }(x) := {\left\{ \begin{array}{ll} u(x) &{} \text {if } u(x)\le \tau ,\\ +\infty &{} \text {else}. \end{array}\right. } \end{aligned}$$

Proposition D.1

The operator ${\varLambda } : {{\mathbb {U}}}\rightarrow {{\mathbb {U}}}$ defined by (47, right) is continuous and obeys the following properties, where $\delta _0,\delta _1$ are positive constants, and where $u,v \in {{\mathbb {U}}}$, and $s,t \ge 0$, $\tau \in {{\mathbb {R}}}$ are arbitrary

• Monotone: if $u\le v$ then ${\varLambda } u \le {\varLambda } v$.

• Subadditive: ${\varLambda }(u+t) \le {\varLambda } u + t$.

• $\delta _0$-submultiplicative: ${\varLambda }[(1+s) u] \le (1+s) {\varLambda } u - \delta _0 s$.

• Existence of a super-solution: there is ${\underline{u}} \in {{\mathbb {U}}}$ such that ${\varLambda } {\underline{u}} \le {\underline{u}}$.

If in addition ${\varTheta }(N_x,{{\mathcal {V}}}_h^x)<\pi /2$ for all $x \in X$, then the operator ${\varLambda }$ is also

• $\delta _1$-causal: if $u^{\le \tau } = v^{\le \tau }$ then $({\varLambda } u)^{\le \tau +\delta _1} = ({\varLambda } v)^{\le \tau +\delta _1}$.

Proof

The monotony of ${\varLambda }$ follows from the monotony of linear interpolation ${{\mathcal {I}}}_h^x$. Likewise, the subadditivitity of ${\varLambda }$ follows from the same property of ${{\mathcal {I}}}_h^x$ (actually ${\varLambda }(u+t) = {\varLambda } u + t$ at all points $x\in X$ for which the stencil ${{\mathcal {V}}}_h^x$ does not contain the source $x_*$). Submutiplicativity is established as follows, using the 1-homogeneity of the interpolation operator ${{\mathcal {I}}}_h^x$

$$\begin{aligned} \min _{y \in \partial {{\mathcal {V}}}_h^x} (1+s){{\,\mathrm{{{{\mathcal {I}}}}}\,}}_h^x u(y) + N_x(x-y) \le (1+s) \Big [\min _{y \in \partial {{\mathcal {V}}}_h^x} {{\,\mathrm{{{{\mathcal {I}}}}}\,}}_h^x u(y) + N_x(x-y)\Big ] - s \min _{y \in \partial {{\mathcal {V}}}_h^x} N_x(x-y), \end{aligned}$$

thus with $\delta _0 = \min \{N_x(x-y);\, x\in X,\, y \in \partial {{\mathcal {V}}}_h^x\}$. Consider the directed graph, with an edge (x, y) of length $N_x(x-y)$ whenever y is a vertex of ${{\mathcal {V}}}_h^x$. Then the distance from a given point $x_0\in X$ to the source $x_*$, denoted by ${\underline{u}}(x_0)$, is finite by assumption and obeys ${\varLambda } {\underline{u}} \le {\underline{u}}$. Finally, Proposition 2.3 establishes $\delta _1$-causality with $\delta _1$ the minimal value of $\Vert y-x\Vert \mu _*(N_x) \cos {\varTheta }(N_x,{{\mathcal {V}}}_h^x)$ among all $x\in X$ and all vertices y of the stencil ${{\mathcal {V}}}_h^x$. $\square$

In the remainder of this appendix, we do not use the specific form (47, right) of the operator ${\varLambda }$, but only the properties established in Proposition D.1. From monotony, subadditivity, and $\delta _0$-submultiplicativity, one derives the discrete comparison principle.

Proposition D.2

(Discrete comparison principle) Let $u, v \in {{\mathbb {U}}}$. If $u\le {\varLambda } u$ and ${\varLambda } v \le v$ then $u \le v$. In addition, if either inequality is strict then $u<v$.

Proof

Let $x \in X$ be such that $t:=u(x)-v(x)$ is maximal, so that $u \le v+t$ and $u(x) = v(x)+t$. Assuming that $t \ge 0$ we obtain $u(x)\le {\varLambda } u(x) \le {\varLambda } [v+t](x) \le {\varLambda } v(x) +t \le v(x)+t = u(x)$, by monotony and subadditivity. If either the first or last inequality is strict, we obtain a contradiction, thus $t<0$ and therefore $u<v$ as announced. Otherwise note that $v_\varepsilon := (1+\varepsilon ) v$ obeys $v_\varepsilon < {\varLambda } v_\varepsilon$ for any $\varepsilon >0$ by $\delta _0$-submultiplicativity, thus $u < v_\varepsilon$ by the previous argument, hence $u \le v$ by letting $\varepsilon \rightarrow 0$, which concludes the proof. $\square$

Using in addition the continuity of ${\varLambda }$ and the existence of a supersolution, one establishes that the fixed point problem (47, left) can be solved by iterating the operator. Finitely many iterations are sufficient if the operator is $\delta _1$-causal.

Proposition D.3

(Convergence of the global iterative method) The operator ${\varLambda }$ admits a unique fixed point ${\varvec{u}}$, and for any $u \in {{\mathbb {U}}}$ one has ${\varLambda }^n u \rightarrow {\varvec{u}}$ as $n \rightarrow \infty$. If in addition ${\varLambda }$ is $\delta _1$-causal and $u>0$, then ${\varLambda }^n u = {\varvec{u}}$ for all $n\ge \max ({\varvec{u}})/\delta _1$.

Proof

Proposition D.2 yields the uniqueness (but not the existence) of the fixed point ${\varvec{u}}$. The null function ${\overline{u}} = 0$ satisfies ${\varLambda } {\overline{u}} \ge \delta _0 \ge 0 = {\overline{u}}$, by $\delta _0$-submultiplicativity. Choose $t\ge 0$ sufficiently large so that ${\overline{v}}:={\overline{u}}-t \le u \le {\underline{u}}+t=:{\underline{v}}$, and note that ${\overline{v}} \le {\varLambda } {\overline{v}}$ and ${\varLambda } {\underline{v}} \le {\underline{v}}$ by subadditivity of ${\varLambda }$. Thus ${\overline{v}} \le \cdots \le {\varLambda }^n {\overline{v}} \le {\varLambda }^n u \le {\varLambda }^n {\underline{v}} \le \cdots \le {\underline{v}}$ by monotonicity of ${\varLambda }$, and induction on $n \ge 0$. By the monotone convergence theorem, ${\varLambda }^n {\overline{v}}$ and ${\varLambda }^n {\underline{v}}$ admit limits as $n \rightarrow \infty$. By continuity, these limits are fixed points of ${\varLambda }$, thus are equal to ${\varvec{u}}$ by uniqueness. By the squeeze theorem we obtain ${\varLambda }^n u \rightarrow {\varvec{u}}$ as announced.

Finally, assume that ${\varLambda }$ is $\delta _1$-causal, that $u>0$, and note that ${\varvec{u}}\ge {\varLambda } {\overline{u}}\ge \delta _0 > 0$. Then $u^{\le 0} = {\varvec{u}}^{\le 0}$, and thus by induction $({\varLambda }^n u)^{\le n \delta _1} = ({\varLambda }^n {\varvec{u}})^{\le n\delta _1} = {\varvec{u}}^{\le n \delta _1}$ for all $n \ge 0$. The result follows. $\square$

Global iteration is a poor way to allocate computational ressources in front propagation problems, and more efficient algorithms concentrate their efforts on a narrow band along the front. The convergence of iterative methods such as fast sweeping [47], the AGSI [10], or the FIM [24], follows from closely related arguments. The fast marching method Algorithm 1 [60] solves the fixed point problem (47, left) in finitely many steps with complexity ${{\mathcal {O}}}(N\ln N)$, see [41, Proposition A.2] for a proof based on the properties established in Proposition D.1, causality included. We next establish that the properties of Proposition D.1 are stable under perturbation.

Proposition D.4

(Operator perturbation) Let $\alpha _*,\alpha ^*\ge 0$, and for all $x\in X$ let $\alpha _x : X \rightarrow [-\alpha _*,\alpha ^*]$. Define ${{\tilde{{\varLambda }}}} : {{\mathbb {U}}}\rightarrow {{\mathbb {U}}}$ by ${{\tilde{{\varLambda }}}} u(x) := {\varLambda }[u+\alpha _x](x)$. Then ${{\tilde{{\varLambda }}}}$ is continuous, monotone, subadditive, is $(\delta _0-\alpha ^*)$-submultiplicative if $\delta _0> \alpha _*$, and admits the super-solution $(1+\alpha ^*/\delta _0){\underline{u}}$. If ${\varLambda }$ is $\delta _1$-causal with $\delta _1>\alpha _*$, then ${{\tilde{{\varLambda }}}}$ is $(\delta _1-\alpha _*)$-causal.

Proof

Fix $x\in X$, $u,v\in {{\mathbb {U}}}$, and $s,t\ge 0$. The continuity of ${{\tilde{{\varLambda }}}}$ immediately follows from the continuity of ${\varLambda }$. If $u\le v$, then $u+\alpha _x\le v+\alpha _x$, thus ${\varLambda }[u+\alpha _x] \le {\varLambda }[v+\alpha _x]$ since ${\varLambda }$ is monotone, therefore ${{\tilde{{\varLambda }}}}$ is monotone. One has ${\varLambda }[u+t+\alpha _x] \le {\varLambda }[u+\alpha _x] +t$ by subadditivity of ${\varLambda }$, thus ${{\tilde{{\varLambda }}}}$ is subadditive. One has ${\varLambda }[(1+s) u + \alpha _x] = {\varLambda }[(1+s)(u+\alpha _x)-s\alpha _x] \le {\varLambda }[(1+s)(u+\alpha _x)+ s\alpha _*] \le {\varLambda }[(1+s)(u+\alpha _x)]+s\alpha _* \le (1+s){\varLambda }[u+\alpha _x] -\delta _0 s + s\alpha _*$, using successively the monotony, subadditivity, and $\delta _0$-submultiplicativity of ${\varLambda }$, thus ${{\tilde{{\varLambda }}}}$ is $(\delta _0-\alpha _*)$-submultiplicative if $\delta _0>\alpha _*$. One has ${\varLambda }[(1+s){\underline{u}} + \alpha _x]\le {\varLambda }[(1+s) {\underline{u}}]+\alpha ^* \le (1+s){\varLambda } {\underline{u}}-\delta _0 s+\alpha ^*$, by subadditivity and submultiplicativity of ${\varLambda }$, thus choosing $s=\alpha ^*/\delta _0$ yields a supersolution of ${{\tilde{{\varLambda }}}}$. Finally, if ${\varLambda }$ is $\delta _1$-causal and $u^{\le \tau } = v^{\le \tau }$, then $(u+\alpha _x)^{\le \tau -\alpha _*} = (v+\alpha _x)^{\le \tau -\alpha _*}$ and therefore $({\varLambda }[u+\alpha _x])^{\le \tau -\alpha _*+\delta _1} \le ({\varLambda }[v+\alpha _x])^{\le \tau -\alpha _*+\delta _1}$, thus ${{\tilde{{\varLambda }}}}$ is $(\delta _1-\alpha _*)$-causal. $\square$

The first-order source factorization (21) falls in the framework of Proposition D.4, with $\alpha _x(y) := u_*(x)-u_*(y) + \langle \nabla u_*(x),y-x\rangle$, which satisfies $\alpha _x(y) = {{\mathcal {O}}}( h^2/ \Vert x-x_*\Vert )$, where $u_*$ is the source factor and $x_*$ is the source point. On the other hand, inspection of the proof of Proposition D.1 yields that $\delta _0 = {\hat{\delta }}_0 h$ and $\delta _1 = {\hat{\delta }}_1 h$ where ${\hat{\delta }}_0$ and ${\hat{\delta }}_1$ are independent of the grid scale h. Thus $\delta _0 \ge \Vert \alpha _x\Vert _\infty$ and $\delta _1\ge \Vert \alpha _x\Vert _\infty$ when h is sufficiently small (except for points x in a ball of radius ${{\mathcal {O}}}(h)$ around the source point $x_*$), and thus Proposition D.4 applies to the factored scheme.

The second and third-order schemes define perturbations (22) and (23) which depend on the unknown u, and thus do not fall in the framework of Proposition D.4. This is the reason why, following [51], we use them in a cautious way: only in the post-processing step of the fast marching method right before the accepted value is frozen⁴ see line 3 of Algorithm 1, and only if their magnitude does not exceed $C h ^2$ where C is an absolute constant. Together, these limitations ensure that the fast marching algorithm still terminates in a single pass over the domain, and produces an output obeying ${\varLambda } u = u + {{\mathcal {O}}}(h^2)$. Therefore $u = {\varvec{u}}+ {{\mathcal {O}}}(h)$ where ${\varvec{u}}$ is the solution of the original scheme, by Proposition D.5 below. In other words, we cannot prove that the high order variants of the scheme improve the solution accuracy, but at least they do not jeopardize first-order accuracy, and neither substantially increase computation time.

Proposition D.5

Let $u\in {{\mathbb {U}}}$ and let $k_*,k^*\ge 0$ be such that $k_* \le {\varLambda } u-u \le k^*$. Then $1-k^*/\delta _0 \le u/{\varvec{u}} \le 1+k^*/\delta _0$.

Proof

Let $s\ge 0$. Then ${\varLambda }[(1+s) u] \le (1+s) {\varLambda } u -\delta _0 s \le (1+s)(u+k^*)-\delta _0 s$, by $\delta _0$-submultiplicativity. Choosing $s = k^*/(\delta _0-k^*)$ yields ${\varLambda }[(1+s) u] \le (1+s) u$ and thus $(1+s)u \ge {\varvec{u}}$ by Proposition D.2. On the other hand, $(1+s){\varLambda }[u/(1+s)] \ge {\varLambda } u + s \delta _0 \ge u-k_*+s\delta _0$, again by $\delta _0$-submultiplicativity. Choosing $s=k_*/\delta _0$ yields ${\varLambda }[u/(1+s)] \ge u/(1+s)$ and thus $u/(1+s) \le {\varvec{u}}$ by Proposition D.2. The result follows. $\square$