Asymptotic Behaviour of Time Stepping Methods for Phase Field Models

TL;DR

This paper analyzes the asymptotic behavior of various time stepping methods for phase field models like Allen Cahn and Cahn Hilliard, predicting optimal step sizes for metastable dynamics and confirming findings through computational experiments.

Contribution

It provides a formal asymptotic analysis of different time stepping methods, including energy stable and implicit schemes, for phase field models during metastable dynamics.

Findings

Energy stable methods often require more steps than others.

Backward Euler preserves energy decay and profile fidelity at larger steps.

Eyre type methods perform worse due to loss of profile fidelity.

Abstract

Adaptive time stepping methods for metastable dynamics of the Allen Cahn and Cahn Hilliard equations are investigated in the spatially continuous, semi-discrete setting. We analyse the performance of a number of first and second order methods, formally predicting step sizes required to satisfy specified local truncation error $\sigma$ in the limit of small order parameter $\epsilon \rightarrow 0$ during meta-stable dynamics. The formal predictions are made under stability assumptions that include the preservation of the asymptotic structure of the diffuse interface, a concept we call profile fidelity. In this setting, definite statements about the relative behaviour of time stepping methods can be made. Some methods, including all so-called energy stable methods but also some fully implicit methods, require asymptotically more time steps than others.The formal analysis is confirmed in computational studies. We observe that some provably energy stable methods popular in the literature perform worse than some more standard schemes. We show further that when Backward Euler is applied to meta-stable Allen Cahn dynamics, the energy decay and profile fidelity properties for these discretizations are preserved for much larger time steps than previous analysis would suggest. The results are established asymptotically for general interfaces, with a rigorous proof for radial interfaces. It is shown analytically and computationally that for most reaction terms, Eyre type time stepping performs asymptotically worse due to loss of profile fidelity.

Equations

$$u_t=\Delta u-\frac{1}{{ϵ}^2}f(u)$$

$$u_t=-ϵ\Delta \Delta u+\frac{1}{ϵ}\Delta f(u).$$

$$\mathcal{E}=\int \left(|\nabla u{|}^2/2+W(u)/{ϵ}^2\right)$$

$$\frac{u_{n+1}-u_n}{k_n}=\Delta u_{n+1}-\frac{1}{{ϵ}^2}f(u_{n+1}).$$

$$u-k\Delta u+\frac{k}{{ϵ}^2}f(u)=u_n$$

$$(1-k\Delta )w=-k\cdot \frac{f(u_1)-f(u_2)}{{ϵ}^2}=-\frac{k}{{ϵ}^2}\cdot f^{\prime }(s(x))w,$$

$$\left[1+\frac{k{f}^{\prime }(s)}{{ϵ}^2}-k\Delta \right]w=0,$$

$$\frac{1}{k}\int |u-u_n{|}^2+\frac{1}{2}\int \left(|\nabla u{|}^2-|\nabla u_n{|}^2+|\nabla u-\nabla u_n{|}^2\right)=-\frac{1}{{ϵ}^2}\left(f(u),(u-u_n)\right).$$

$$|F(u)-F(u_n)-f(u)(u-u_n)|=\left|{\int }_u^{u_n}{f}^{\prime }(s)(s-u_n)ds\right|\le \frac{f_{\infty }}{2}(u-u_n{)}^2.$$

$$(\frac{1}{k}-\frac{f_{\infty }}{2{ϵ}^2})\int |u-u_n{|}^2+E[u]-E[u_n]\le 0,$$

$$u-k\Delta u+\frac{k}{{ϵ}^2}{u}^3=u_n+\frac{k}{{ϵ}^2}{u}_n$$

$$u-k\Delta u+\frac{Mk}{{ϵ}^2}u=u_n-\frac{k}{{ϵ}^2}\left(u_n^3-(M+1)u_n\right)$$

$$u+k\Delta \Delta u-\frac{k}{{ϵ}^2}\Delta f(u)$$

$$u+k\Delta \Delta u-\frac{k}{{ϵ}^2}\Delta u^3$$

$$\mathcal{F}(u)=\Delta u-f(u)/{ϵ}^2$$

$$\mathcal{F}(u)=-ϵ\Delta \Delta u+\Delta f(u)/ϵ$$

$$u-\frac{k}{2}\mathcal{F}(u)u=u_n+\frac{k}{2}\mathcal{F}(u_n)$$

$$\frac{3u}{2}-k\mathcal{F}(u)=2u_n-\frac{1}{2}{u}_{n-1}.$$

$$(W(u)-W(u_n))/(u-u_n)$$

$$u_{\ast }-\alpha k\mathcal{F}(u_{\ast })=u^n$$

$$u-\alpha k\mathcal{F}(u)=u^n+(1-\alpha )k\mathcal{F}(u_{\ast })$$

$$u_p=u_n+\frac{k}{3}\left(\mathcal{F}(u_n)+4\mathcal{F}(u_{n+1})+\mathcal{F}(u_{n+2})\right)$$

$$\Vert u_{n+2}-u_p{\Vert }_{\infty }\le \sigma .$$

$$u\approx \pm 1.$$

$$u(x,t)\approx g(z)$$

$$u(x,t)\approx g(\mathrm{dist}(x,\Gamma )/ϵ)$$

$$u_t\approx g^{\prime }(\mathrm{dist}(x,\Gamma ))V/ϵ$$

$$\frac{{\partial }^{n}u}{\partial t^n}=O({ϵ}^{-n}).$$

$$k=O(\sqrt{\sigma }ϵ)\text{mbox}(BE)$$

$$u-u^n-k\Delta u+k\left[u^3-u\right]/{ϵ}^2$$

$$u-u^n-k\Delta u+k\left[u^3-u\right]/{ϵ}^2+k(u-u_n)/{ϵ}^2$$

$$k=O(\sqrt{\sigma }{ϵ}^{3/2})\text{mbox}(Eyre)$$

$$k^3{u}_{ttt}/12=O(k^3/{ϵ}^3)$$

Full text

Keywords: Allen-Cahn Equation, Allen-Cahn Equation, Phase Field Model, Time Stepping, Energy Stability.

1 Introduction

The mathematical literature for computational methods for Allen-Cahn (AC) dynamics [2], and its higher order relative Cahn Hilliard (CH) dynamics [4], is dominated by the proposal, use, and analysis of so-called energy stable schemes [20, 18, 23, 7]. AC and CH dynamics are gradient flows on an energy functional, and the solution should decrease that energy in time. Energy stable schemes guarantee that decrease no matter what time step is chosen. This is a desirable property not shared by standard fully implicit or semi-implicit (IMEX) time stepping methods. We will show in this work that some (but not all) fully implicit methods can outperform energy stable schemes when subject to fixed accuracy requirements. The recent article [24] gives especially clear evidence that when time steps are chosen appropriately, fully implicit methods are conditionally energy stable, and further that the large time steps allowed by energy stable schemes can come at the cost of significant loss of accuracy. We extend the commentary in [24] to show that in the metastable dynamic regime of AC and CH, some fully implicit methods can take optimally sized time steps. By optimal, we mean the asymptotically largest time steps as the order parameter $\epsilon\rightarrow 0$ that satisfy a given local error tolerance. Here, $\epsilon$ represents the width of interfacial layers in metastable dynamics and, like the authors of [24], we use the form of the equations scaled so that these dynamics transpire in an $O(1)$ time scale. When the dynamics are in this metastable regime, which dominates the time of typical phenomena of interest, definite statements about the behaviour of different time stepping methods can be made. This criteria does not take into account solver efficiency. However, we can make definite statements on how efficient solvers for nonlinear implicit time stepping need to be to outperform other methods.

A combination of asymptotic analysis and careful computational work backs up our claims. In addition, we present a rigorous result for implicit time stepping for meta-stable AC dynamics in radial geometry that shows that asymptotically larger time steps can be taken than previous analysis would suggest. These time steps preserve the diffuse interface structure (a property that we call proflie fidelity) and also the energy decay property of the equations. This result is shown for a class of reaction terms. An interesting result in Section 6.2 shows that Eyre-type time stepping can perform asymptotically worse with most reaction terms, while implicit time stepping has uniform asymptotic behaviour over a class of reaction terms. This was predicted by the analysis and confirmed computationally.

Our study focuses on pure materials science applications rather than the use of Cahn-Hilliard equations to track interfaces in so-called diffuse interface methods [25] in which the CH dynamics are coupled to other physics. We consider the simplest form of AC and CH dynamics, whose Gamma limit (as $\epsilon\rightarrow 0$) is well understood and use that well known structure to gain insight into the behaviour of the schemes. The authors believe that the insight gained from these studies will also apply to schemes used for other materials science models which are less well understood.

We consider a number of first and second order time stepping schemes: the energy stable Eyre’s method [9]; Backward (Implicit) Euler (BE) [10]; Trapezoidal Rule (TR) [10]; Second order Backward Differentiation Formula (BDF2) [10]; Secant [8]; standard semi-implicit (linear IMEX) methods of first and second order [3]; first and second order Scalar Auxiliary Methods (SAV) [19] for which a modified energy stability can be proved; and finally a second order Singular Diagonally Implicit Runge Kutta method with good stability properties (DIRK2) [10]. The resulting implicit systems are considered in the spatially continuous semi-discrete setting in a 2D periodic domain, with numerical validation done with a suitably refined Fourier spectral approximation. Time step schemes that result in nonlinear systems are solved with Newton’s iterations using the Preconditioned Conjugate Gradient Solver (PCG) developed in [6] at each iteration. Adaptive time stepping is done based on a user-specified local error tolerance $\sigma$. The variation of the number of time steps with $\epsilon$ for fixed $\sigma$ is predicted based on formal consideration of the local truncation error of the schemes in the metastable dynamics. The formal predictions are then validated in computational studies. With this criteria, first order BE performs better (asymptotically fewer time steps as $\epsilon\rightarrow 0$) than Eyre and first order IMEX and SAV. Second order TR and BDF2 perform better than Secant, DIRK2, and second order IMEX and SAV. The difference in both cases is asymptotically larger for CH than AC. These comparisons are also valid for computational time, using PCG counts as the measure, to similar accuracy. It is seen that optimal numbers of time steps are obtained when the dominant local truncation error is a higher order time derivative. This observation may have application in other systems with metastable dynamics. We observe that standard IMEX methods perform almost identically to SAV methods of the same order in the scenario we consider, at reduced computational cost.

It is observed that the global accuracy of BE is better than a naïve prediction based on the size of the local truncation error would suggest. A formal analysis of the scheme for the AC case shows that the dominant error made in one time step is asymptotically smaller than expected. This is due to a special structure of the local truncation error for BE, in which the asymptotically largest term lies in a strongly damped space.

We introduce the equations and numerical schemes in Section 2 with some introductory analysis. The scaling for AC and CH is chosen so that the metastable interface dynamics (approximate curvature motion for AC and Mullins-Sekerka flow [15] for CH) occurs in $O(1)$ time. In Section 3 we examine the metastable dynamics of the equations and make predictions for the behaviour of the time steps with $\epsilon$ and local error tolerance $\sigma$ under stability assumptions which are verified numerically in Section 4. We give an asymptotic analysis for the surprising accuracy and stability properties for BE with large time steps applied to AC in Section 5. In Section 6 we present the rigorous result for BE applied to AC with large time steps and also show the loss of profile fidelity for Eyre-type time stepping for most reaction terms. We end with a short discussion.

2 Equations and Schemes

We consider the simplest form of the AC dynamics for $u({\bf x},t)$ given by

[TABLE]

where $f(u)=u^{3}-u$ is the classical form of the reaction term. More general reaction terms are considered in Section 6. CH dynamics is described by a higher order partial differential equation

[TABLE]

For computational simplicity, we consider the two-dimensional (2D) cases of these equations in a doubly periodic cell $[0,2\pi]^{2}$. The time scaling in the equations above is chosen to give sharp interface (as $\epsilon\rightarrow 0$) motion in $O(1)$ time. The sharp interface limit yields curvature driven flow for AC and a nonlocal Mullins-Sekerka flow for CH [15]. Both types of dynamics have an associated energy functional

[TABLE]

where $W(u)=\frac{1}{4}(u^{2}-1)^{2}$ and the reaction term $f(u)=W^{\prime}(u)$. The energy $\mathcal{E}(t)$ is monotonic decreasing due to the gradient flow nature of the dynamics. For AC the gradient is in $L_{2}$ and for CH it is $H^{-1}$.

2.1 Time stepping

2.1.1 Backward Euler

We consider the simplest implicit scheme, first order Backward Euler (BE), also known as Implicit Euler. Applied to (1) keeping space continuous, we have

[TABLE]

where $u_{n}({\bf x})$ approximates the exact solution $u({\bf x},t_{n})$ and $k_{n}=t_{n+1}-t_{n}$ is the time step. We use the classical $f(u)=u^{3}-u$ as mentioned above. Dropping the subscript on the time step and the unknown solution at time level $n+1$ we have the nonlinear problem

[TABLE]

for $u$ given $u_{n}$.

Definition 1**.**

A time stepping scheme is said to have the energy decay property if ${\mathcal{E}}(u_{n+1})\leq{\mathcal{E}}(u_{n})$.

This property could be conditional on the choice of time step size. Additionally, it could depend on $u_{n}$. If a scheme has the energy decay property for any $u_{n}$ and $k$, the scheme is called unconditionally energy stable.

Theorem 1**.**

Consider (4), assume that $u_{n}\in H^{2}(\Omega)$ and $u_{n}$ takes values in $[-1,1]$, then there exists $u\in H^{2}(\Omega)$ that solves (4) with values in $[-1,1]$. Define $f_{\infty}:=\max\{|f^{\prime}(s)|,\,s\in[-1,1]\}$, then if $k\leq 2\epsilon^{2}/f_{\infty}$ the solution $u$ is unique and satisfies the energy decay property. Note that the energy stability result was established earlier in [24] with a different proof.

Proof.

The existence of $u$ follows from the standard method of sub-/super-solutions applied to comparison functions $-1$ and $+1$. To establish uniqueness, we assume $u_{1}$ and $u_{2}$ are solutions. Then their difference $w=u_{1}-u_{2}$ is a solution of

[TABLE]

where $s$ takes values between $u_{1}$ and $u_{2}$, and hence in $[-1,1]$. Isolating $w$ leads to the elliptic problem

[TABLE]

and if $k<\epsilon^{2}/f_{\infty}$ then the corresponding elliptic operator is strictly positive and $w$ is zero by the maximum principle. To establish energy decay, we take the inner product of (4) with the test function $u-u_{n}$:

[TABLE]

From the Fundamental Theorem of Calculus we develop the expansion,

[TABLE]

Using this relation to eliminate $f(u)$ yields the equality,

[TABLE]

which yields energy decay for $k<2\epsilon^{2}/f_{\infty}$. The Theorem is also true when homogeneous Neumann boundary conditions are specified. ∎

Thus we have existence of solutions to (4) for any time step size, and uniqueness and energy stability under the resitriction $k\leq 2\epsilon^{2}/f_{\infty}$. This is true for any $u_{n}$ under the restrictions of the Theorem. We shall see in Section 6 that asympoticaly larger time steps $k=o(\epsilon)$ can be taken when the dynamics are slow (interface motion) with locally unique, energy stable solutions. This is verified in computational tests.

2.1.2 Eyre’s Method

An alternative first order scheme to fully implicit BE was proposed by Eyre [9]:

[TABLE]

The scheme is derived conceptually by keeping a convex part of the reaction term $f(u)=u^{3}-u$ implicit and a concave part explicit. In this sense, it is an IMEX method but an unusual one since a nonlinear term is kept implicit and a linear term is handled explicitly. The method has appealing properties:

Theorem 2** (from [9]).**

The time step (5) has a unique solution $u$ for any $u_{n}$ and $k$ that is unconditionally energy stable.

Additional first order schemes considered are the SAV scheme [19] and a linear IMEX method [3]:

[TABLE]

with $M>0$, sometimes called a stabilization term. We take $M=2$ since that makes the left hand side a linearization about the far field values, but computational performance is relatively insensitive to $M$. The SAV scheme is energy stable with a modified energy. We use the same stabilization coefficient as above in the SAV scheme. There is a class of linearly implicit energy stable schemes [5, 12, 13] that require an asymptotically large stabilization term $O(\epsilon^{-p})$ with $p$ large and increasing from AC to CH and 2D to 3D for the analysis. These methods are theoretically interesting but are extremely inaccurate and not useful for practical applications. We have further discussion of these schemes in Remark 2.

All time stepping schemes can be applied to CH (2), with BE and Eyre shown below:

[TABLE]

In this case, BE is known to have unique solutions with the energy decay property when $k<\epsilon^{3}$ [24] and Eyre is unconditionally energy stable [9].

2.1.3 Second Order Schemes

We also consider the second order methods Trapezoidal Rule (TR), Secant (S) [8], Second Order Backward Differencing (BDF2), and Second Order Singular Diagonal Implicit Runge Kutta (DIRK2) [10] methods. These are described below for $u_{t}=\mathcal{F}(u)$ with

[TABLE]

With this notation:

[TABLE]

Secant is a variant of TR with the term $f(u)-f(u_{n})$ replaced by

[TABLE]

where $W$ is the energy term from (3). It is known to be conditionally energy stable [8]. For the simple form of $W$ we have taken, the expression above can be factored explicitly. DIRK2 is a two stage method

[TABLE]

with $\alpha=1-1/\sqrt{2}$. Both DIRK2 and BDF2 are A-stable, and so preferable to TR and Secant from the perspective of stiff ODE solver theory [10]. A second order linear IMEX method (SBDF2 [3]) and two variants of second order SAV methods based on BDF2 are also considered.

2.2 Spatial discretization and solution procedure

The current work concentrates on the time stepping errors, and it is convenient to consider the semi-discrete, spatially continuous approximation. This idealization is approximated well by the Fourier spectral spatial discretization. The computational results shown have sufficient spatial resolution that spatial errors do not affect the results in the digits shown. The computations are done in a full 2D setting, rather than in a reduced dimensional radial setting as could be done, in order to give PCG iteration counts for the nonlinear time stepping methods that have meaning for more general computations. Note that the PCG counts are independent of spatial resolution when the problem is resolved.

2.3 Error estimation and adaptive time stepping

We perform two time steps of the same size $k$ in order to use a specialized predictor $u_{p}$ for $u_{n+2}$.

[TABLE]

where $\mathcal{F}(u)=\Delta u-f(u)/\epsilon^{2}$ for AC and $-\epsilon\Delta\Delta u+\Delta f(u)/\epsilon$ for CH as above. Time step sizes are adjusted so that

[TABLE]

The predictor $u_{p}$ is formally one order more accurate than the numerical approximation $u_{n+2}$ from time stepping, up to fifth order. The predictor has an inherent dominant local error $k^{5}u_{ttttt}/90$ that is a pure time derivative of $u$, which is shown below in Section 3.1.3 to be a desirable property.

For the one step methods, the time step is adjusted adaptively to maintain a local error below $\sigma$ as described in [6]. For BDF2 and its linear variants, time steps are only adjusted by a factor of two. When time steps are reduced (using Hermite cubic interpolation for the restart value) or increased, four time steps are taken before checking the local error to allow relaxation of the initial error layer.

3 Local Truncation Errors in Metastable Dynamics

3.1 Metastable dynamics

In our formulation, it is known that after a short time $O(\epsilon^{2})$ solutions to AC tend to interfaces between regions of solution near the equilibrium values,

[TABLE]

These interfaces have width $\epsilon$ and move approximately with curvature motion. We refer to this dynamics as metastable or slow, even though with the particular time scaling we have chosen it occurs in in $O(1)$ time. For the majority of the time, the solution will be in this regime, so we concentrate now on the expected and observed behaviour of time stepping in this setting. With the choice of $f(u)=u^{3}-u$, we have

[TABLE]

with $g(z):=\tanh(z/\sqrt{2})$ and $z={\rm dist}(x,\Gamma)/\epsilon$, where $\Gamma$ is the approximate interface with arc length parameter $s$ moving with curvature motion (normal velocity equal to curvature). We fix its location at the $u=0$ level set. The local coordinates $(s,z)$ are shown in Figure 1. This structural result on the metastable solution can be obtained with formal asymptotics. In the outer asymptotic region for AC the solution takes the form $u=\pm 1$ to all orders. Curvature motion as the limit $\epsilon\rightarrow 0$ has been proven rigorously [16, 1].

CH has the same metastable solution structure (8) with normal interface velocity given by Mullins-Sekerka flow, in $O(1)$ time in our scaling (2). We refer the reader to the review article [17] for details.

From (8), we see that time and space derivatives are large near the interface. Starting with

[TABLE]

we can take a time derivative to obtain:

[TABLE]

where $V$ is the normal velocity at the point on $\Gamma$ closest to $x$. Formally taking higher derivatives in this pattern yields:

[TABLE]

This is used to analyze the truncation error of the time stepping schemes.

3.1.1 Predicted time step sizes for AC

A standard strategy for adaptive time stepping is to have a user specified local error tolerance of $\sigma$. The error for each time step is estimated and the time step adjusted so that there is an estimated error in that single time step less than $\sigma$. It is known that the dominant local truncation error for BE is $k^{2}u_{tt}/2$ which in metastable dynamics is $O(k^{2}/\epsilon^{2})$ from (9). The local truncation error restriction then requires time steps of size

[TABLE]

We now proceed to determine the expected behaviour of time steps with $\epsilon$ and $\sigma$ from the other schemes. We can write the BE scheme (4) and Eyre’s scheme (5) for AC in an instructive way

[TABLE]

Knowing that the truncation error for BE is $O(k^{2}/\epsilon^{2})$ we see that the truncation error for the Eyre scheme is dominated by the last term in its expression above, which has leading order $k^{2}u_{t}/\epsilon^{2}=O(k^{2}/\epsilon^{3})$. Our time step prediction in this case is

[TABLE]

Thus, the advantage of the Eyre scheme to be able to take large time steps and remain energy stable is never realized if accurate computational results are required. Reference [24] has an alternate way to view the loss of accuracy that does not highlight this asymptotic difference. The first order IMEX and SAV schemes have the same asymptotic behaviour as Eyre.

Remark 1**.**

The formal local error analysis above relies on the stability of the schemes in metastable dynamics under the resulting time step restrictions. More than simple stability, the analysis requires that the time stepping preserves the asymptotic structure of the diffuse interface. This is the concept we have named profile fidelity. All predicitions described in this section lead to time stepping that preserves profile fidelity for the classical choice of $f(u)=u^{3}-u$. We observe the predicted time step behaviour in $\epsilon$ and $\sigma$ computationally. In Section 6.2 we show that for (most) other reaction terms, Eyre time stepping loses profile fidelity for time steps $k=O(\epsilon^{3/2})$ and in these cases, $k=O(\epsilon^{2})$ is needed for accuracy.

Remark 2**.**

The first order, linearly implicit energy stable scheme for 2D AC is analyzed in [5]. The analysis requires a stabilization term of order $\epsilon^{-2}|\ln\epsilon|$. If such a scheme were implemented, the time steps required for a local error tolerance of $\sigma$ would be $k=O(\sqrt{\sigma}\epsilon^{5/2}/|\ln\epsilon|)$, prohibitively small for practical computation.

We can determine the dominant term in the local truncation errors of the second order schemes applied to AC:

[TABLE]

We consider two second order SAV variants based on how an extrapolated approximation is computed. If the extrapolated value of $u^{n+1}$ is taken as $2u^{n}-u^{n-1}$ the scheme (referred to as SAV2-A) behaves similarly to SBDF2. If the extrapolated value is computed with a first order linear IMEX scheme as suggested in [19] (referred to as SAV2-B), the scheme has a local truncation error of order $k^{3}/\epsilon^{5}$. The results are summarized in Table 1. It is clear that BE takes asymptotically (as $\epsilon\rightarrow 0$) fewer time steps than Eyre, although they are both first order in time step size. TR and BDF2 take asymptotically fewer time steps than Secant, DIRK2, SBDF2, SAV2-A and SAV2-B although they are all second order methods. The computations in Section 4 below show that these time step estimates correspond to real computational behaviour.

Remark 3**.**

We predict the number $M$ of time steps in Tables 1 and 2 and how it varies with $\epsilon$ and $\sigma$. As shown in Figure 3 we are also predicting how a profile of time steps $k(t)$ behaves with $\epsilon$ and $\sigma$.

3.1.2 Predicted time step sizes for CH

The same local truncation analysis can be done for the CH in the metastable regime where the solution has the same interface structure (8) with the interface $\Gamma$ moving approximately with Mullins-Sekerka flow in $O(1)$ time. BE, TR, BDF2, and SBDF2 have the same error expressions as above, but Eyre, Secant and DIRK2 have local truncation errors when applied to CH listed below:

[TABLE]

where we have used the fact that the Laplacian $\Delta$ increases the size of terms by $1/\epsilon^{2}$ near the interface. The first order IMEX and SAV schemes have the same asymptotic behaviour as Eyre. SAV2-A behaves similarly to SBDF2 as before, with SAV2-B worse by a power of $\epsilon$ as for the AC case above. The results are summarized in Table 2. The predictions in this table are validated in the numerical experiments in the next section. Although the methods all have the formal order of accuracy in terms of time step size, the behaviour as $\epsilon\rightarrow 0$ varies significantly. Note that the gap between BE and the other first order schemes, and between TR/BDF2 and Secant/DIRK2/SBDF2/SAV2-A is wider for CH dynamics than it was for AC.

3.1.3 Discussion: the source of increased local error

In the metastable regime, the two terms in AC and CH (diffusion and nonlinear reaction) are both large but approximately cancel to give the slow dynamics. The methods with asymptotically (as $\epsilon\rightarrow 0$) small local errors (BE, TR, BDF2) have dominant truncation errors that are pure time derivatives of the solution, which inherit this high order cancellation. The other methods which have large local errors have truncation errors that involve the reaction term individually. This imbalance amplifies the size of the error. As an example, DIRK2 applied to $u_{t}=\mathcal{F}(u)$ has an error proportional to $\mathcal{F}^{\prime\prime}u_{t}^{2}$. From this discussion, we believe the ranking of the schemes in this work will also apply to other nonlinear problems with metastable dynamics.

4 Computational Results

4.1 Allen Cahn

We take initial conditions in the form of a radial front

[TABLE]

and compute with $\epsilon$ = 0.2, 0.1, 0,05 and 0.025. The benchmark for accuracy is the time $T$ at which the value at the domain centre $(\pi,\pi)$ changes from negative to positive. Except for the exponentially small (in $\epsilon$) derivative discontinuities at the periodic boundaries, the dynamics approximate the sharp interface limit of curvature motion of a circle. The expectation from asymptotic analysis of the sharp interface limit is that

[TABLE]

This is confirmed by the numerical solutions below. Some snapshots of the dynamics are shown in Figure 2. A video of the dynamics is also available [21].

4.1.1 First order methods

The PCG approach is known to have bounded condition number under the scaling $k=C\epsilon^{2}$ for BE with $C<1$ [24] and we observe good behaviour in the example below even with $C>1$ in the metastable regime. It is observed computationally in this work that the PCG for Eyre’s method is independent of $k$ and $\epsilon$ although the authors are not aware of a proof in the literature.

Results of the numerical experiments in which $\sigma$ and $\epsilon$ were varied for BE and Eyre are shown in Tables 3 and 4. Spatial errors do not affect the digits shown in any of the computational results in this paper.

Table 3 validates the second order $O(k^{2})$ local truncation error since the number of time steps was predicted to be $M=O(1/\sqrt{\sigma})$ for both methods with $\epsilon$ constant, noting that $\sqrt{10}\approx 3.16$. Such results for other schemes and for the CH benchmark problem below are not shown, but verify the formal accuracy of the schemes. Table 4 validates the prediction of $M=O(1/\epsilon)$ for BE and $M=O(1/\epsilon^{3/2})$ for Eyre with $\sigma$ constant, noting that $2^{3/2}\approx 2.83$. Both tables validate the prediction that for the same local tolerance $\sigma$, Eyre involves more computational work than BE and gives less accurate answers. CG counts for both methods are small as expected. You see (unexpectedly) that the final accuracy of BE does not seem to degrade as $\epsilon\rightarrow 0$ for fixed $\sigma$. This is discussed in Section 5 below. Although BE does not guarantee energy stability, no step accepted by the local error tolerance exhibited an energy increase.

For completeness, we show the time step sizes as a function of time for BE in Figure 3 with $\epsilon$ and $\sigma$ varied. As mentioned in Remark 3 our predictions for the behaviour of the time steps sizes $k$ as $\epsilon$ and $\sigma$ are varied describe a profile $k(t)$.

We repeat the $\epsilon\rightarrow 0$ study for IMEX1 and SAV1 in Table 5. These methods require a fixed number of FFT calculations per time step to invert the constant coefficient linear implicit aspect of the schemes, with SAV1 requiring four times as many solves as IMEX1. It is seen that IMEX1 behaves almost identically to SAV1 and both are superior to Eyre’s method when computational cost is considered. In the context of this study, there is no benefit from the theoretical guarantees of energy stable schemes and BE is the optimal (with our asymptotic definition) first order scheme with IMEX1 the runner up. This will remain true for other nonlinear solver strategies for BE as long as they require fewer than $O(1/\sqrt{\epsilon})$ iterations when adaptive time steps are taken.

Remark 4**.**

Note that for the BE computation for $\epsilon=0.025$ we can still get reasonable accuracy taking $\sigma=10^{-2}$. In this case, the maximum value of $k/\epsilon^{2}$ is 14.6. Clearly, the theory which guarantees existence of solutions and energy decay for $k<\epsilon^{2}$ [24] can be improved for metastable dynamics. This is explored in the analysis in Sections 5 and 6 below.

4.1.2 Second Order Methods

The CG counts of all the nonlinear second order methods are relatively insensitive to $\epsilon$, similar to the first order methods shown above. We show the number of time steps used for the seven methods in Table 6, for $\sigma=10^{-4}$ fixed and $\epsilon$ varied. The superiority of TR and BDF2 is clearly seen with $M=O(1/\epsilon)$, compared to $M=O(1/\epsilon^{4/3})$ (noting that $2^{4/3}\approx 2.52$) for Secant, DIRK2, SBDF2, SAV2-A and $M=O(1/\epsilon^{5/3})$ (noting that $2^{5/3}\approx 3.18$) for SAV-B as predicted above. The pattern in the number of time steps for the multi-step methods is a bit rougher due to the strict criteria we have used for adaptive time step change. As above, we see no benefit from the theoretical guarantees of energy stable schemes. Fully implicit methods TR and BDF2 are asymptotically optimal in terms of the number of time steps and are computationally optimal if the solvers require fewer than $O(1/\sqrt[3]{\epsilon})$ iterations when adaptive time steps are taken (which appears to be the case with the Newton PCG solver we used). SBDF2 is the runner up and notably it is comparable to the fully implicit DIRK2 method but does not have the overhead of a nonlinear solve.

It is interesting to note that the slight change in the extrapolation procedure in the SAV2 schemes makes such a difference to their asymptotic performance. It is confirmation that merely considering the order of time stepping scheme and its theoretical energy stability properties is not the whole story.

4.2 Cahn Hilliard

For the initial conditions we take

[TABLE]

with $r=\sqrt{(x-\pi)^{2}+(y-\pi)^{2}}$ and compute with $\epsilon$ = 0.2, 0.1, 0,05 and 0.025. As before, the benchmark is the time $T$ at which the value at the domain centre $(\pi,\pi)$ changes from negative to positive. The dynamics are shown in Figure 4 and a video of the dynamics is also available [22] .

4.2.1 First order methods

Results of the numerical experiments in which $\epsilon$ is varied for the first order methods are shown in Table 7. These validate the prediction of $M=O(1/\epsilon)$ for BE and $M=O(1/\epsilon^{2})$ for Eyre and IMEX1 with $\sigma$ constant. As for the AC case, SAV1 behaves similarly to IMEX1 at increased computational cost. For CH, the implicit problem for BE is more difficult to solve as $\epsilon\rightarrow 0$ with fixed $\sigma$, but it is still more accurate than Eyre stepping for equivalent computational cost. It will be asymptotically more efficient as long as the solution strategy for the nonlinear problem requires fewer than $O(1/\epsilon)$ iterations with adaptive time stepping. As with AC, we see that BE does not suffer from global accuracy decrease as $\epsilon\rightarrow 0$.

4.2.2 Second order methods

The CG counts for the second order methods behave like those of BE with $\epsilon$ as shown above. We show the number of time steps used for the four methods in Table 8, for $\sigma=10^{-4}$ fixed and $\epsilon$ varied. The superiority of TR and BDF2 is clearly seen, consistent with $M=O(1/\epsilon)$ , compared to $M=O(1/\epsilon^{5/3})$ (noting that $2^{5/3}\approx 3.17$) for Secant, DIRK2, and SBDF2 as predicted above. Results for SAV2-A are comparable to those for SBDF2. Again, the implications for the asymptotic computational superiority of fully implicit TR and BDF2 under the assumption of sufficient solver efficiency are clear.

5 Asymptotic Analysis of Properties of BE AC Solutions

The results in Table 4 present the accuracy for BE applied to AC with fixed local error tolerance $\sigma=10^{-4}$ under various values of $\epsilon$. It is remarkable the accuracy in the benchmark time does not degrade as $\epsilon\rightarrow 0$. This is unexpected, as a naïve prediction would be that the final accuracy scaled like $M\sigma=O(\sqrt{\sigma}/\epsilon)$ where $M$ is the number of time steps. It is clear that the resulting solution accuracy for the schemes under specified local error tolerance is a nontrivial question.

We present below the asymptotic analysis of a fully implicit BE time step (4) in two dimensions assuming the solution is in the meta-stable regime. That is, $u_{n}$ is approximately described as a curve ${\bf x}_{n}(s)$ parametrized by arc length with normal $\hat{n}$, dressed with the heteroclinic profile (8). We take the scaling $k=c\epsilon$ with $c$ independent of $\epsilon$, both sufficiently small depending only on the curve ${\bf x}_{n}$. We consider the formal asymptotics for the implicit time step $u$ of (4) in this setting, anticipating that $u$ will have the same local dependence $u(s,z)$. Using

[TABLE]

where $\kappa$ is the curvature of the interface, we find at leading order $O(\epsilon^{-1})$ that $u$ has the same homoclinic structure around a new curve ${\bf x}(s)$. That is,

[TABLE]

with $g(z)=\tanh(z/\sqrt{2})$ and where we have changed coordinates to $(s,z)$ with

[TABLE]

based on the curve ${\bf x}(s)$ after the implicit time step. In the language of Remark 1 we predict that the scheme preserves profile fidelity and show below that this is asymptotically consistent. In (10), $v(z,s)$ is the correction to the leading order solution. We will identify the size and structure of this term below.

We take

[TABLE]

where $\rho$ is the average normal speed through the time step. Recalling that $k=c\epsilon$ and the spatial scaling of $z$, we have

[TABLE]

A diagram is shown in Figure 5. Note that the variation in normal direction appears in higher order asymptotic terms, so it is consistent in what follows to use the same $\hat{n}$ as normal direction for both curves, i.e. the same “$z$”.

Considering now the next order term $O(1)$ in (4) with the forms (10) and (12):

[TABLE]

where $\mathcal{L}:=\partial^{2}/\partial z^{2}+f^{\prime}(g)\cdot$ and we have used the smallness of $c$ for the cubic Taylor approximation of $g(z)-g_{n}(z)$ on the right hand side. We consider (13) at each $s$ in the $L^{2}({\mathbb{R}})$ orthogonal decomposition of $G:={\rm span}\{g^{\prime}(z)\}$ and $G^{\perp}$. Note that $g^{\prime\prime}\in G^{\perp}$ (this does not depend on the specific reaction term $f=u^{3}-u$ chosen here) and $\mathcal{L}$ has $G$ as its kernel and has bounded inverse on $G^{\perp}$ from standard Fredhold theory [11]. Thus we have $\rho=\kappa+O(c^{2})$ and $v=O(c)$ in $G^{\perp}$. Careful examination of these results shows that the errors in $G^{\perp}$ do not accumulate and are globally of size $O(c\epsilon)=O(k)$ and so decrease as $\epsilon\rightarrow 0$. Global errors in interface position after $O(1/k)$ time steps are of size $c^{2}$, independent of $\epsilon$. Global solution value errors due to the position error have size $O(c^{2}/\epsilon)=O(k^{2}/\epsilon^{3})$ and so it is seen that BE behaves like a second order method in this scaling. This explains the unexpected accuracy in AC BE computations as $\epsilon\rightarrow 0$.

Remark 5**.**

Note that the error estimator (7) uses $\mathcal{F}(u)$ which sees the undamped dominant truncation error term, which is why the number of time steps behaves with $\epsilon$ in the manner predicted in Section 3.1.1. Thus for BE applied to AC in the metastable regime, the estimator asymptotically over-estimates the local errors actually made.

The formal asymptotic results can also be used to show that the implicit time steps in this scaling lead to energy decrease. Neglecting the $O(c^{2})$ terms in the interface motion, we have from (11)

[TABLE]

Using the identities for arc length parametrized curves $|{\bf x}_{s}|=1$, $\kappa\hat{n}={\bf x}_{ss}$ and ${\bf x}_{s}\cdot{\bf x}_{sss}=-\kappa^{2}$ it follows by taking the $s$ derivative of the equation above and the dot product with ${\bf x}_{s}$ at each $s$ that

[TABLE]

This shows that the metastable curve at time $n$ is longer than at the next step $n+1$. Since the energy $\mathcal{E}$ is proportional to curve length to highest order in the metastable regime [14], we have shown formally that implicit time stepping for AC has the energy decay property under this time step scaling. Large, accurate, fully implicit time steps can be taken in computations validated in Section 4.

In the next section we show a closely related rigorous result in a radial geometry. The main result is in Proposition 1. A key ingredient is an identification of a dominant term in the space $G$ that represents the interface motion, separate from heavily damped terms in the perpendicular space, as shown here. Care must be taken to control the size of terms which are formally neglected in this asymptotic analysis.

6 Rigorous Radial Analysis of AC With BE and Eyre Time Stepping

We derive rigorous asymptotic evolution of a radially symmetric profile for BE and first order Eyre-type methods for the Allen-Cahn equation in $\mathbb{R}^{2}$. Extensions to radial profiles in $\mathbb{R}^{d}$ is immediate. More precisely we consider a splitting $f=f_{+}-f_{-}$ and study the iterative scheme

[TABLE]

For simplicity, we assume that $f$ is smooth, odd about $u=0$, has precisely three simple zeros at $u=\pm 1$ and at $u=0$, and tends to $\pm\infty$ as $u\to\pm\infty$. This includes the classical choice of $f(u)=u^{3}-u$ but we consider other reaction terms in this class since Eyre’s method can have quite different behaviour as shown in Section 6.2. The BE scheme corresponds to the choice $f_{-}\equiv 0$ while Eyre-type schemes take $f_{+}^{\prime},f_{-}^{\prime}\geq 0$. We pose the problem on the affine space

[TABLE]

with $u_{n}\in Y$ as a given. The assumptions on $f$ imply that the continuous 1D Allen-Cahn equation has a steady state solution

[TABLE]

which is heteroclinic to $\pm 1$; that is $g\to\pm 1$ as $z\to\pm\infty.$ Considering $R>1$, we modify this $g$ at the order $O(e^{-1/\epsilon})$ so that $g^{\prime}=0$ on $(-\infty,-1/\epsilon)$ for $\epsilon\ll 1$. This introduces exponentially small residuals in the sequel that have no impact upon the salient results of our analysis.

We introduce $z=\frac{r-R}{\varepsilon}$, the weighted inner product

[TABLE]

and the associated spaces $L^{2}_{R}$ and $H^{1}_{R}$. We rewrite the iterative equation as

[TABLE]

on the domain $z\in(-R/\epsilon,\infty).$ We decompose $u_{n}$ and $u$ as

[TABLE]

where $R_{n}$ and $v_{n}$ are taken as given and $R$ and $v$ are to be determined. The profile associated to $u_{n}$ is denoted $g_{n}$ and observe that it admits the expansion

[TABLE]

In the sequel we will enforce the orthogonality conditions

[TABLE]

and denote the corresponding subspaces of $L^{2}_{R}$ by $X^{\perp}$ and $X_{n}^{\perp}$ respectively with the associated orthogonal projections $\Pi$ and $\Pi_{n}$.

At this point the analysis of the implicit and Eyre-type schemes diverges sufficiently that we approach them distinctly.

6.1 Backward Euler estimates

For BE we take $f_{-}^{\prime}\equiv 0$, $f=f_{+}$, and write the iterative map as

[TABLE]

where we have introduced the linear operator

[TABLE]

and the nonlinearity

[TABLE]

The operator $L$ is self-adjoint in the weighted inner product for which the eigenvalue problem takes the form

[TABLE]

subject to $\partial_{z}\psi(-R/\epsilon)=0$ and $\psi\to 0$ as $z\to\infty.$ Since the profile $g$ solves (14), it will be useful to compare $L$ to the simpler operator

[TABLE]

arising as the linearization of (14) about $g$ in $L^{2}(\mathbb{R})$. The operator $L_{0}$ is self-adjoint on $L^{2}(\mathbb{R})$, and since $g$ is heteroclinic with $g^{\prime}>0$, the Sturm-Liouville theory on $L^{2}(\mathbb{R})$ implies that $L_{0}$ has a simple, ground-state eigenvalue at $\lambda=0$ with eigenfunction $g^{\prime}$ and the remainder of the spectrum of $L_{0}$ is strictly positive, in particular $L_{0}$ is uniformly coercive on the space $\{g^{\prime}\}^{\perp}_{L^{2}(\mathbb{R})}$. While $L$ does not generically have a kernel, it does have an eigenspace with a small associated eigenvalue. However, for $\epsilon$ sufficiently small, it inherits the coercivity of $L_{0}$.

Lemma 1**.**

Fix $\epsilon_{0}>0$ sufficiently small, then there exists $\alpha>0$, independent of $R\geq 1$ and of $\epsilon\in(0,\epsilon_{0})$, such that

[TABLE]

*for all $v\in H^{1}_{R}$ satisfying $\langle v,g^{\prime}\rangle_{R}=0.$ *

Proof.

We defer the proof of $L^{2}_{R}$ coercivity to the appendix. To extend coercivity to $H^{1}_{R}$ we observe that

[TABLE]

so that for any $t\in(0,1)$ we may write

[TABLE]

where we have introduced $\tilde{\alpha}:=\alpha/(1+\alpha+\|f^{\prime}(g)\|_{\infty})>0.$ Dropping the tilde, we have (20) with $\alpha$ independent of $\epsilon>0$ and $R>1$. ∎

We assume throughout our analysis that $\|v\|_{H^{1}_{R}}$ and $\|v_{n}\|_{H^{1}_{R}}$ are uniformly bounded by $\delta\ll 1$. Returning to (17), we denote its right-hand side as ${\cal F}_{\rm BE}$. To have the inversion of the operator on the left-hand side be contractive the term ${\cal F}_{\rm BE}$ must be approximately orthogonal to the small eigenspace of $L$. As Lemma 1 shows it is sufficient to be $L_{R}$-orthogonal to $g^{\prime}$, the kernel of $L_{0}.$ To this end we determine $R=\hat{R}_{\rm BE}(v,v_{n},R_{n})$ such that ${\cal F}_{\rm BE}\in X^{\perp}$, or equivalently

[TABLE]

Assuming this condition has been enforced we introduce

[TABLE]

and may rewrite the BE iteration in the equivalent formulation

[TABLE]

The key step is the introduction of the operator $\Pi$, the orthogonal projection onto $X^{\perp}$. This is redundant when ${\cal F}_{\rm BE}\in X^{\perp}$, but preserves contractivity for choices of $(v,v_{n},R_{-}R_{n})$ when it is not. Our goal is to show the function ${\cal G}_{\rm BE}$ is a contraction mapping and to develop asymptotic formula for $R$ and $v$.

Lemma 2**.**

The function $R=\hat{R}_{\rm BE}$ satisfies the implicit relation

[TABLE]

where

[TABLE]

Moreover we have the Lipshitz estimate

[TABLE]

so long as $k\delta^{2}\ll\epsilon^{2}.$

Proof.

Due to parity considerations, we remark that $\|g^{\prime}\|_{R}^{2}=R\|g^{\prime}\|^{2}_{L^{2}(\mathbb{R})}$, up to exponentially small terms. For brevity, and as an element of foreshadowing, we approximate $(R-R_{n})$ by $k$ in the $O$-error terms. We address the terms in ${\cal F}_{\rm BE}$ and derive the following elementary estimates,

[TABLE]

For this scheme, ${\cal F}_{\rm BE}$ depends upon $v$ only through ${\cal N}$. Collecting terms in the orthogonality condition that are linear in $R-R_{n}$ and identifying relevant higher order terms yields the relation

[TABLE]

where $b_{1}$ is given in (24). Under the assumptions on $k$ and $\delta$ we have the leading order result $R-R_{n}=-k/R$. Substituting this relation into (30) yields the result (23).

To obtain the Lipschitz estimate we observe from the estimates above that

[TABLE]

The nonlinearity satisfies the Lipschitz properties

[TABLE]

while

[TABLE]

Adding and subtracting $\langle{\cal N}(\tilde{v}),g^{\prime}\rangle_{R}$ and using (29), we arrive at the estimates

[TABLE]

Imposing the condition $k\delta^{2}\ll\epsilon^{2}$ yields (25). ∎

To establish bounds on the map ${\cal G}_{\rm BE}$ defined in (22) we apply $M$ to both sides of the relation and take the $L^{2}_{R}$ inner product with respect to ${\cal G}_{\rm BE}$. Using the coercivity estimate (20) we find

[TABLE]

Taking $v,v_{n}\in B_{H^{1}_{R}}(\delta)$ for $\delta\ll 1$ and recalling that the projection $\Pi$ crucially cancels the leading order term in $g-g_{n}$, we estimate

[TABLE]

For the BE system we examine distinguished limits $k=\epsilon^{s}$, for $s\in(1,2)$, which we call the large time-step regime, for which the $H^{1}_{R}$ term is dominant on the left-hand side of (31). We drop the $L^{2}_{R}$ term to find,

[TABLE]

Taking $\delta=\epsilon^{s^{\prime}}$ for any $s^{\prime}>\max\{s/2,2(s-1)\}$ then we determine that

[TABLE]

for $\epsilon$ sufficiently small. In particular, since $s>\max\{s/2,2(s-1)\}$ in the large time-stepping regime, we may take $\delta=k=\epsilon^{s}$, so that, viewing ${\cal G}_{\rm BE}$ as a map on $(v,v_{n})$, we have ${\cal G}_{\rm BE}:B_{H^{1}_{R}}(k)\times B_{H^{1}_{R}}(k)\mapsto B_{H^{1}_{R}}(k),$ for all $s$ in the large time-step regime.

Proposition 1**.**

Fix $1<s<2$, then in the distinguished limit $k=\epsilon^{s}$ the function ${\cal G}_{\rm BE}$ defined in (22) with $R:=R_{n+1}=\hat{R}_{\rm BE}(v;v_{n},R_{n})$ maps $B_{H^{1}_{R}}(k)\times B_{H^{1}_{R}}(k)$ into $B_{H^{1}_{R}}(k)$ and is a strict contraction. In particular it has a unique solution $v\in B_{H^{1}_{R}}(k)$, denoted by $v_{n+1}$ which satisfies

[TABLE]

In particular there exists $c>0$ such that for all $v_{0}\in B_{H^{1}_{R}}(ck)$ and $R_{0}>1$ the sequence $\{(v_{n},R_{n})\}_{n=1}^{N}$ satisfies $v_{n}\in B_{H^{1}_{R}}(ck)$ while $\{R_{n}\}_{n=0}^{N}$ satisfies the backwards Euler iteration

[TABLE]

where $b_{1}>0$ is given by (24). Here $N$ is the iteration number such that $R_{N}>1$ and $R_{N+1}<1.$

Proof.

We have established the mapping property. To establish the contractivity we must control the impact of $f$ upon the projection $\Pi$ through the motion of the front $R$. We assume that $v,\tilde{v},v_{n}\in B_{H^{1}_{R}}(k)$ and denote $R=R(v)$ and $\tilde{R}=R(\tilde{v})$, with the associated front profiles denoted by $g$ and $\tilde{g}$. The estimate (25) establishes that $\hat{R}_{\rm BE}$ is Lipschitz with constant $ck\delta/\epsilon$, which in the the large time-step regime reduces to $ck^{2}/\epsilon$. Following the proof of (25) we find that

[TABLE]

In the large time-step regime, using (20) we deduce the bound

[TABLE]

We wish to obtain a bound on the difference of ${\cal G}_{\rm BE}$ at two values of $v$:

[TABLE]

We first bound the difference

[TABLE]

The analysis is complicated by the fact that $M$ is only uniformly invertible on the range of $\Pi$. To factor these projected inverses we act with $M$, observing

[TABLE]

where we used that fact that $\tilde{M}\tilde{M}^{-1}\tilde{\Pi}=\tilde{\Pi}$ and hence $\tilde{M}\tilde{M}^{-1}\tilde{\Pi}+(I-\tilde{\Pi})=I.$ Since the right-hand side of (36) lies in the range of $\Pi$ we may invert boundedly,

[TABLE]

To recover the whole ${\rm g}_{\rm BE}$ we act with $(I-\Pi)$ on (35) obtaining

[TABLE]

Adding (38) to (37) yields a regularized expression that accounts for the shifts in the projections

[TABLE]

The operators $M^{-1}\Pi$ and $\tilde{M}^{-1}\tilde{\Pi}$ are bounded using (34), while

[TABLE]

where $\|\cdot\|_{R*}$ denotes the operator norm from $L^{2}_{R}$ into itself. The projections satisfy

[TABLE]

Applying these estimates to (39) and using (31) to estimate $\Pi{\cal F}_{\rm BE}$ we obtain

[TABLE]

Finally we write

[TABLE]

and using (33), (40) estimate

[TABLE]

which is contractive so long as $k\ll\epsilon^{\frac{2}{3}}$ which holds with the large time-step regime.

Within the large time-step regime the leading order iteration (30) simplifies as $k\ll k^{2}/\epsilon^{2}$ and the dominant correction is given by the $b_{1}$ term. To compare to standard notation we rewrite the regime as $\epsilon^{2}\ll k=\delta\ll 1$ and replace the internal parameter $\delta$ with $k$, the result is the large time-step interation (32). ∎

6.2 Eyre-type iterations

For an Eyre iteration the map (15) takes the form

[TABLE]

where we have introduced the Eyre linear operator

[TABLE]

the explicit-term residual

[TABLE]

and the nonlinearity

[TABLE]

which we further decompose into implicit and explicit parts

[TABLE]

The operator $L_{+}$ is self-adjoint in the weighted inner product for which the eigenvalue problem takes the form

[TABLE]

subject to $\partial_{z}\psi(-R/\epsilon)=0$ and $\psi\to 0$ as $z\to\infty.$ The coercivity estimate is substantially simpler than for BE as the operator $L_{+}$ is strictly positive without constraint.

Lemma 3**.**

There exists $\alpha_{+}>0$, independent of $R\geq 1$, such that

[TABLE]

*for all $v\in H^{1}_{R}$. *

Proof.

Since $f_{+}^{\prime}\geq 0$ the normalized ground-state eigenfunction $\psi_{0}$ of $L_{+}$, satisfies

[TABLE]

Since the ground-state eigenvalue is strictly positive, this establishes the $L^{2}_{R}$ coercivity of $L_{+}$ with $\alpha_{+}=\lambda_{0}^{+}$. The $H^{1}_{R}$ coercivity follows as in Lemma 1. ∎

We assume throughout our analysis that $\|v\|_{H^{1}_{R}}$ and $\|v_{n}\|_{H^{1}_{R}}$ are uniformly bounded by $\delta\ll 1$. We denote the right-hand side of (41) by ${\cal F}_{\rm E}$ and introduce

[TABLE]

which is strictly contractive on the full space $L^{2}_{R}$, and re-write the Eyre iteration as

[TABLE]

For the Eyre iteration the role of the projection $\Pi$ is diminished as $M_{+}$ is contractive without it. Our goal is to show the existence of a map $R=\hat{R}_{\rm E}(v,v_{n},R_{n})$, for which

[TABLE]

to establish the contractive mapping properties of ${\cal G}_{\rm E}$, and to develop asymptotic formula for $R$ and $v$. We do this in the long time-stepping regime, $k\gg\epsilon^{2}$, which has no lower bound for the Eyre scheme.

Lemma 4**.**

Assume $k\gg\epsilon^{2}$. There exists a smooth function $\hat{R}_{\rm E}:B_{H^{1}_{R}}(\delta)\times B_{H^{1}_{R}}(\delta)\times\mathbb{R}\mapsto\mathbb{R}$ such that the profile $g=g(z;R)$ satisfies (45). The function $R=\hat{R}_{\rm E}$ satisfies the implicit relation

[TABLE]

where we have introduced the leading order Eyre time constant

[TABLE]

when $f_{-}^{\prime}\not\equiv 0$. Moreover we have the Lipshitz estimate

[TABLE]

so long as $\delta\ll 1$.

Proof.

Due to parity considerations, we remark that $\|g^{\prime}\|_{R}^{2}=R\|g^{\prime}\|^{2}_{L^{2}(\mathbb{R})}$, up to exponentially small terms. For brevity, and as an element of foreshadowing, we approximate $(R-R_{n})$ by $\epsilon^{2}$ in the $O$-error terms. Addressing the terms in ${\cal F}_{\rm BE}$ one by one, we record

[TABLE]

With these reductions we can simplify the orthgonality condition, identifying terms that are linear in $R-R_{n}$ and most relevant higher order terms. The result is the balance

[TABLE]

where $c_{-}$, introduced in (47) is positive since $f_{-}^{\prime}\geq 0$ by assumption. The largest terms and error terms come from the residual, and we kept the lower order constant on the left-hand side to emphasize that in the long time-stepping regime, the residual dominates the natural time-step term. Indeed, the iteration is independent of step size, $k$, given at leading order by (54).

To obtain the Lipshitz estimate we observe from the bounds above that dependence of $\hat{R}_{\rm E}$ on $v$ arises from the balance of the linear $R-R_{n}$ term in the residual against the nonlinearity. Since both these terms are multiplied by $k/\epsilon^{2}$ this factor cancels and we have the balance

[TABLE]

The nonlinearity satisfies the Lipshitz properties

[TABLE]

while

[TABLE]

Adding and subtracting $\langle{\cal N}(\tilde{v},v_{n}),g^{\prime}\rangle_{R}$ and using (53), we arrive at the estimates

[TABLE]

Imposing the condition $\delta\ll 1$ yields (48). ∎

We may now establish the main result on the Eyre sequence.

Proposition 2**.**

There exists $c>0$ such that for any $k\gg\epsilon^{2}$ the function ${\cal G}_{\rm E}$ defined in (44) with $R:=\hat{R}_{\rm E}(v;v_{n},R_{n})$ maps $B_{H^{1}_{R}}(c\epsilon)\times B_{H^{1}_{R}}(c\epsilon)$ into $B_{H^{1}_{R}}({\cal G}_{\rm E}(0,v_{n}),\epsilon^{2})$ and is a strict contraction, satisfying

[TABLE]

In particular ${\cal G}_{\rm E}$ has a unique fixed point in that set, which we denote $v_{n+1}$. Moreover, if the Eyre balance parameter

[TABLE]

then for any $\rho\in(0,1)$ there exists $c>0$ such that for all $v_{0}\in B_{H^{1}_{R}}(c\epsilon)$ and $R_{0}>1$, the sequence $\{(v_{n},R_{n})\}_{n=1}^{N}$ satisfies $v_{n}\in B_{H^{1}_{R}}(c\epsilon)$ and

[TABLE]

where the Eyre number, $c_{E}$, is defined by

[TABLE]

*where $K_{+}>0$ is defined in (65) and $K_{+}L_{+}^{-1}\Pi>0$ is self-adjoint. *

Proof.

To establish the contractivity of ${\cal G}_{\rm E}$ we follow the arguments for backward Euler, sketching only the differences. We introduce

[TABLE]

and derive the expression

[TABLE]

The operators $M_{+}^{-1}\Pi$ and $\tilde{M}_{+}^{-1}\tilde{\Pi}$ are bounded as $L_{+}$ has no small eigenvalues. Using (48) we estimate

[TABLE]

Similarly the projections satisfy

[TABLE]

Applying these estimates to (58,59) and following the proof of (48) to estimate $\Pi{\cal F}_{\rm E}$ we obtain

[TABLE]

Finally we write

[TABLE]

and estimate the ${\cal F}_{\rm E}$ term from which the dominant contribution comes from the residual

[TABLE]

where we used (48) in the last inequality. In particular we deduce that

[TABLE]

Combining (60), (62) and (61), imposing $\delta=\epsilon$, and using $k\gg\epsilon^{2}$ we arrive at strict contractivity on $B_{H^{1}_{R}}(c\delta)$ for any fixed $c>0$.

To establish bounds on the the fixed point $v_{n+1}$ of ${\cal G}_{\rm E}(\cdot;v_{n})$ we observe from (43) that in the large time-stepping regime

[TABLE]

Using this result we expand

[TABLE]

Inverting $M_{+}$ we find, at leading order

[TABLE]

In particular we deduce that

[TABLE]

Arguing inductively, since the Eyre balance parameter $\gamma<1$ and the functions $\|L_{+}^{-1}\Pi f_{-}^{\prime}(g)g^{\prime}\|_{H^{1}_{R}}$ are uniformly bounded for all $R\geq 1$, we deduce that if $\delta:=\|v_{0}\|_{H^{1}_{R}}=O(\epsilon)$ then the sequences $\{(R-R_{n})\epsilon^{-2}\}_{0}^{N}$ and $\{\epsilon^{-1}\|v_{n}\|_{H^{1}_{R}}\}_{0}^{N}$ are uniformly bounded, independent of $\epsilon\ll 1$ and $k\gg\epsilon^{2}$ for all $n\leq N$ so long as $R_{n}>1$ for all $n=0,\ldots,N.$

To improve this bound we require Lipschitz estimates on the $v_{n}$ component of ${\cal G}_{\rm E}$. To this end we find

[TABLE]

Here we introduce the quasi-steady parameter $\rho\in(0,1).$ Since $|R_{n}-R_{m}|=O(\epsilon^{2-\rho})$ for $|n-m|\leq N_{\rho}:=\leq\epsilon^{-\rho}$ we infer that

[TABLE]

for all such $n$ and $m$. For $n>N_{\rho}$ we define the quasi-equilibrium

[TABLE]

where $E_{n}$ is the $R=R_{n}$ translate of

[TABLE]

Here the self-adjoint operator

[TABLE]

is well defined since $\|L_{+}\Pi\circ f_{-}^{\prime}(g)\|_{H^{1}_{R}*}=\gamma<1$ by assumption. Using the Lipschitz property (6.2) of ${\cal G}_{\rm E}$ and the quasi equilibrium relation

[TABLE]

we deduce that

[TABLE]

for $k=n-m_{\epsilon},\ldots,n$. Since $\gamma^{N_{\rho}}\ll\epsilon$ we deduce from an inductive argument that

[TABLE]

for all $n>N_{s}.$ Inserting this result in (46) we arrive at the leading order Eyre iteration (57). ∎

Remark 6**.**

There are two examples of particular relevance

[TABLE]

with the decomposition $f_{+}=(1+\beta)u^{3}$ and $f_{-}=u+\beta u^{3}$ for $\beta>0$. The choice $\beta=0$ is classical and very degenerate, as in this case $f_{-}^{\prime}(u)=1$ and the corresponding Eyre balance parameter $\gamma$, defined in (56) is zero, and the Eyre number, (57) is 1. In this case it is possible to rewrite Eyre’s method as backward Euler with a rescaled time. In particular the slow convergence to equilibrium will not be in evidence. For larger values of $\beta$ the balance parameter increases from zero and the Eyre number decreases from 1. As the balance parameter increases through $1$ we anticipate enhanced slowing of the front profile as the Eyre number tends to zero. The choice of non-zero $\beta$ can be viewed as spurious, a deliberate attempt to foul the method. A more robust example of non-zero balance arises naturally through the model

[TABLE]

with $\beta\geq 1.$ This suggests the optimal decomposition $f_{+}=u^{5}$ and $f_{-}=\beta u^{3}$. Here, unambiguously, increasing $\beta$ increases the balance parameter and will lead to non-trivial enhanced slow-down with potential instability as $\gamma$ increases through $1.$ These analytic predictions are validated in a computational study below.

Remark 7**.**

To leading order, in the large time-stepping regime $k\gg\epsilon^{2}$, the Eyre iteration recovers backward Euler with the substitution $k\mapsto c_{E}\epsilon^{2}.$ This reduces to the exact result for the case $f(u)=u^{3}-u$ and $f_{-}(u)=u$, for which $f_{-}^{\prime}=1$, as the Eyre constant reduces to $1$ since $\Pi f_{-}(g)g^{\prime}=\Pi g^{\prime}=0.$

The strong contractivity of ${\cal G}_{\rm E}$ with respect to $v$, given in (55), arises from the strong convexity with respect to $v$, but the slow evolution and marginal convergence to the quasi-equilibrium, given in (6.2) arises from the balance between the implicit and explicit terms. The parameter $\gamma$ measures this balance, with the quasi-equilibrium structure lost as $\gamma$ increases towards 1. Indeed, since $\|K_{+}\|_{H^{1}_{R}*}\sim(1-\gamma)^{-1}$, the Eyre constant will generically tend to zero as $\gamma\to 1$.

6.2.1 Computational Validation of Remark 6

We perform computations for AC with the non-classical $f(u)=u^{5}-u^{3}$ (which also leads to meta-stable dynamics of curvature motion) using the same initial conditions and accuracy criteria as described in Section 4.1. BE performs almost identically to the results shown in Tables 3 and 4 for the classical $f(u)=u^{3}-u$ in terms of accuracy and variation of time steps with $\epsilon$ and $\sigma$. This matches the theory in Section 6.1 which can be summarized as BE has profile fidelity when $k=o(\epsilon)$.

When Eyre’s method is applied to the dynamics with $f(u)=u^{5}-u^{3}$, with the natural splitting suggested in Remark 6, profile fidelity is lost as predicted. The formal prediction of $k=O(\epsilon^{3/2})$ which was seen computationally for $f(u)=u^{3}-u$ in Table 4 is not observed for $f(u)=u^{5}-u^{3}$. Rather, we see $k=O(\epsilon^{2})$ as predicted by the theory in the previous section. The numerical results are shown in Table 9.

7 Summary and Future Work

We have identified the time step scaling for several first and second order schemes for AC and CH under the restriction of fixed local truncation error, $\sigma$. In particular, we derive the asymptotic behaviour of time-step number with $\sigma$ and asymptotic parameter $\epsilon$ during meta-stable dynamics. These predictions are made under the assumption that the time steps preserve the asymptotic structure of the diffuse interface, a concept we refer to as profile fidelity. The predictions are verified in numerical experiments. We see that methods whose dominant local truncation error can be expressed as a pure time derivative have optimal asymptotic performance in this particular limit. BE, TR, and BDF2 all have this desirable property. We believe these methods will also have superior performance for other problems with metastable dynamics. Our numerical results show that BE performs better than expected and we have shown an explanation of this behaviour with formal asymptotics.

The optimal fully implicit methods asymptotically computationally outperform all linearly implicit methods in the limit we consider. We present precise criteria on the computational cost of nonlinear solvers for this comparison. The provably energy stable first and second order SAV schemes had higher computational cost than standard IMEX methods for similar results. As a final result, we present a rigorous proof that large time steps with fully implicit BE can be taken with locally unique solutions that are energy stable. This is done for the 2D radial AC equation in meta-stable dynamics. Eyre-type iteration is also considered in this analytic framework, and it is shown that in general this approach loses profile stability unless very small time steps are taken.

Our work gives strong evidence that some fully implicit schemes for phase field models should be given more consideration and that provably energy stable schemes are inaccurate and give no benefits during meta-stable dynamics. Perhaps a hybrid scheme that switches between the two approaches depending on the dynamic regime should be considered.

Extending the analysis to the non-radial case and to CH is an interesting question. We observed that the question of global accuracy is not trivial in Section 5 and should be considered for other schemes. Accurate local error estimation for these problems is another interesting question to pursue.

Appendix: $L_{R}^{2}$ coercivity of $L$

Here we show the technical argument for Lemma 1. For $u,v\in H^{1}(\mathbb{R})$ we define the inner product

[TABLE]

with the standard norms $L^{2}_{\ell}$ and $H^{1}_{\ell}$ while $L^{2}_{\ell^{c}}$ is defined in $\mathbb{R}\backslash[-\ell,\ell].$ Let $L_{0}$ be as defined in (19).

Lemma 5**.**

Fix $\ell_{0}>0$ sufficiently large there exists $\alpha>0$ such that for all $\ell>\ell_{0}$

[TABLE]

for all $u\in H^{1}_{\ell}\cap L^{2}(\mathbb{R})$ satisfying $\langle u,g^{\prime}\rangle_{L^{2}(\mathbb{R})}=0$ and $1=\|u\|_{L^{2}_{\ell}}\geq\|u\|_{L^{2}_{\ell^{c}}}.$

Proof.

Let $\phi$ be the minimizer of $\langle L_{0}u,u\rangle_{\ell}$ over $H^{1}(\mathbb{R})$ subject to $\|u\|_{L^{2}_{\ell}}=1$ and the full-line orthogonality $\langle u,\psi_{0}\rangle_{L^{2}(\mathbb{R})}=0$. By scaling, the minima is attained with $\|\phi\|_{L^{2}_{\ell}}=\frac{1}{2}$ and satisfies

[TABLE]

subject to Neumann boundary conditions $\phi_{x}(\pm\ell)=0$, in addition to the full line orthogonality condition. The operator $L_{0}$ on the truncated domain has eigenvalues $\lambda_{0}^{\ell}<\lambda_{1}^{\ell}<...$ which are $O(e^{-d\ell})$ far away from the eigenvalues of $L_{0}$ on the full line. In particular $\lambda_{0}^{\ell}$ may be negative, but the rest are uniformly positive. In $L^{2}_{\ell}$ we partition $\phi=\beta\psi_{0}^{\ell}+\phi^{\perp}$, where $\psi_{0}^{\ell}$ is the $L^{2}_{\ell}$ ground state of $L_{0}$ and $\phi^{\perp}\in L^{2}_{\ell}$ is $L^{2}_{\ell}$ orthogonal to $\psi_{0}^{\ell}.$ Then we have

[TABLE]

On the other hand the orthogonality condition implies that

[TABLE]

where the subscript $\ell^{c}$ denotes integration over $\mathbb{R}\backslash[-\ell,\ell]$ with the corresponding norms. In particular we deduce that

[TABLE]

Since $g^{\prime}$ decays exponentially at $\pm\infty$, is complementary norm is exponentially small in $\ell.$ From orthogonality of $\psi_{0}^{\ell}$ and $\phi^{\perp}$ we have

[TABLE]

or equivalently

[TABLE]

and taking $\ell$ large enough we use these bound in (67) to show that $\alpha$ is exponentially close to $\lambda_{1}^{\ell}>0.$ ∎

To complete the proof of Lemma 1 we take $\ell$ sufficiently large to apply Lemma 5 and then take $\epsilon$ sufficiently small that $\epsilon|z|\leq\epsilon\ell\ll 1$. Under these conditions $L^{2}_{\ell}$ and $L^{2}_{R}(-\ell,\ell)$ are equivalent norms, uniformly in $\epsilon$, and we have uniform $L^{2}_{R}(-\ell,\ell)$ coercivity of $L$. Conversely, $L$ is clearly $L^{2}_{R}$ coercive on $[-R/\epsilon,\infty)\backslash[-\ell,\ell]$ since $f^{\prime}(g)$ is strictly positive there. Clearly $L$ is uniformly coercive on function with more than half their $L^{2}_{R}$ mass in $[-R/\epsilon,\infty)\backslash[-\ell,\ell]$. The $g^{\prime}$ orthogonality condition implies approximate orthogonality to $\psi_{0}$ for $\ell$ large. From these we deduce the full $L^{2}_{R}$-coercivity of $L$ over $X$.

Acknowledgements

KP recognizes support from the NSF DMS under award 1813203. BW acknowledges support from an NSERC Canada grant. XC is supported by a UBC International Doctoral Fellowship.