Parareal Convergence for Oscillatory PDEs with Finite Time-scale Separation

TL;DR

This paper extends the convergence analysis of the Parareal method for oscillatory PDEs to cases with finite time-scale separation, providing practical guidelines for parameter selection.

Contribution

It proves convergence of the Parareal method for finite time-scale separation and introduces an optimization-based approach for choosing averaging parameters.

Findings

Convergence is achievable with finite time-scale separation.

An optimization problem determines the averaging window and time step.

Guidelines for parameter selection improve practical implementation.

Abstract

A variant of the Parareal method for highly oscillatory systems of PDEs was proposed by Haut and Wingate (2014). In that work they proved superlinear conver- gence of the method in the limit of infinite time scale separation. Their coarse solver features a coordinate transformation and a fast-wave averag- ing method inspired by analysis of multiple scales PDEs and is integrated using an HMM-type method. However, for many physical applications the timescale separation is finite, not infinite. In this paper we prove con- vergence for finite timescale separaration by extending the error bound on the coarse propagator to this case. We show that convergence requires the solution of an optimization problem that involves the averaging win- dow interval, the time step, and the parameters in the problem. We also propose a method for choosing the averaging window relative to the time step based as a function of the finite frequencies inherent in the problem.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\frac{1}{\epsilon }\mathcal{L}\mathbf{u}+\mathcal{N}(\mathbf{u},\mathbf{u})=0,$$

$${\mathbf{U}}_n^k={\overline{\phi }}_{\Delta T}({\mathbf{U}}_{n-1}^k)+({\phi }_{\Delta T}({\mathbf{U}}_{n-1}^{k-1}-{\overline{\phi }}_{\Delta T}({\mathbf{U}}_{n-1}^{k-1})),k=1,2,\dots$$

$$\begin{array}{c}\frac{\partial \overline{\mathbf{u}}(\mathbf{x},t)}{\partial t}+\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{\int }_0^{\tau }{e}^{s\mathcal{L}}\mathcal{N}(e^{-s\mathcal{L}}\overline{\mathbf{u}}(\mathbf{x},t),e^{-s\mathcal{L}}\overline{\mathbf{u}}(\mathbf{x},t))\mathrm{d}s=0,\\ {\overline{\mathbf{u}}(\mathbf{x},t)|}_{t=0}={\mathbf{u}}^0(\mathbf{x}),\end{array}$$

$${\mathbf{u}}^0(x,t,\tau )=e^{-\tau \mathcal{L}}\bar{\mathbf{u}}(\mathbf{x},t).$$

$$\frac{\partial \overline{\mathbf{u}}}{\partial t}+e^{t\mathcal{L}/\epsilon }\overline{\mathcal{N}}(e^{-t\mathcal{L}/\epsilon }\overline{\mathbf{u}},e^{-t\mathcal{L}/\epsilon }\overline{\mathbf{u}})=0,$$

$$\overline{\mathcal{N}}(\overline{\mathbf{u}}(t))$$

$$\approx \frac{1}{\overline{M}}\sum _{m=0}^{\overline{M}-1}\rho \left(\frac{s_m}{T_0}\right)e^{s_m\mathcal{L}}\mathcal{N}(e^{-s_m\mathcal{L}}\overline{\mathbf{u}}(t))$$

$$\frac{\partial \overline{\mathbf{u}}}{\partial t}=-\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{\int }_0^{\tau }\sum _{\mathbf{k}\in {\mathbb{Z}}^2}\sum _{\alpha =-1}^1\left[\sum _{\mathbf{k}={\mathbf{k}}_1+{\mathbf{k}}_2}\sum _{{\alpha }_1,{\alpha }_2}{\sigma }_{{\mathbf{k}}_1}^{{\alpha }_1}{\sigma }_{{\mathbf{k}}_2}^{{\alpha }_2}{C}_{\mathbf{k},{\mathbf{k}}_1,{\mathbf{k}}_2}^{\alpha ,{\alpha }_1{\alpha }_2}{e}^{i(\mathbf{k}\cdot \mathbf{x})-i{\Omega }_{\mathbf{k},{\mathbf{k}}_1,{\mathbf{k}}_2}^{\alpha ,{\alpha }_1,{\alpha }_2}s/\epsilon }\right]{\mathbf{r}}_{\mathbf{k}}^{\alpha }\mathrm{d}s.$$

$${\mathcal{S}}_{\mathbf{k},\alpha }=\{({\mathbf{k}}_1,{\mathbf{k}}_2,{\alpha }_1,{\alpha }_2):\mathbf{k}={\mathbf{k}}_1+{\mathbf{k}}_2,{\omega }_{\mathbf{k}}^{\alpha }={\omega }_{{\mathbf{k}}_1}^{{\alpha }_1}+{\omega }_{{\mathbf{k}}_2}^{{\alpha }_2}\}.$$

$$\frac{\partial {\sigma }_{\mathbf{k}}^{\alpha }}{\partial t}+\sum _{{\mathcal{S}}_{\mathbf{k},\alpha }}{\sigma }_{{\mathbf{k}}_1}^{{\alpha }_1}{\sigma }_{{\mathbf{k}}_2}^{{\alpha }_2}{C}_{{\mathbf{k}}_1,{\mathbf{k}}_2,\mathbf{k}}^{{\alpha }_1,{\alpha }_2,\alpha }=0,$$

$$e^{s\mathcal{L}/\epsilon }\mathcal{N}(e^{-s\mathcal{L}/\epsilon }\overline{\mathbf{u}}(t),e^{-s\mathcal{L}/\epsilon }\overline{\mathbf{u}}(t))$$

$$=\sum _{{\mathcal{S}}_{\mathbf{k},\alpha }}{\mathcal{N}}_n(\overline{\mathbf{u}}(t))+\sum _{\beta =1}^{\infty }\left(\sum _{{\mathcal{S}}_{\mathbf{k},\alpha }^{{ϵ}_{\beta }}}{e}^{i{\lambda }_{n}s}{\mathcal{N}}_n(\overline{\mathbf{u}}(t))\right)$$

$${\mathcal{S}}_{\mathbf{k},\alpha }^{{ϵ}_{\beta }}=\left\{({\mathbf{k}}_1,{\mathbf{k}}_2,{\alpha }_1,{\alpha }_2):\mathbf{k}={\mathbf{k}}_1+{\mathbf{k}}_2,{ϵ}_{\beta -1}<\frac{1}{\epsilon }|{\omega }_{\mathbf{k}}^{\alpha }-{\omega }_{{\mathbf{k}}_1}^{{\alpha }_1}+{\omega }_{{\mathbf{k}}_2}^{{\alpha }_2}|\le {ϵ}_{\beta }\right\},$$

$${\mathbf{v}}_t(t)=e^{t\mathcal{L}/\epsilon }N\left(e^{-t\mathcal{L}/\epsilon }\mathbf{v},e^{-t\mathcal{L}/\epsilon }\mathbf{v}\right),t\in [0,\Delta T],$$

$$\tilde{\mathbf{v}}\left(\tau \right)=\mathbf{v}\left(t\right).$$

$${\partial }_t\mathbf{v}\left(t\right)={\partial }_t\tilde{\mathbf{v}}\left(t/(\epsilon \Delta T)\right)=\frac{1}{\epsilon \Delta T}{\partial }_{\tau }\tilde{\mathbf{v}}\left(t/(\epsilon \Delta T)\right)=\frac{1}{\epsilon \Delta T}{\partial }_{\tau }\tilde{\mathbf{v}}\left(\tau \right).$$

$${\partial }_{\tau }\tilde{\mathbf{v}}\left(\tau \right)=\epsilon \Delta T{e}^{\tau \Delta T\mathcal{L}}N\left(e^{-\tau \Delta T\mathcal{L}}\tilde{\mathbf{v}}(\tau ),e^{-\tau \Delta T\mathcal{L}}\tilde{\mathbf{v}}\left(\tau \right)\right).$$

$${\mathbf{x}}_t=\epsilon \mathbf{f}\left(\mathbf{x},t\right).$$

$${\mathbf{y}}_t=\epsilon \overline{\mathbf{f}}\left(\mathbf{y},t\right),$$

$$\overline{\mathbf{f}}\left(\mathbf{x},t\right)=\frac{1}{\eta }{\int }_0^{\eta }\rho \left(\frac{s}{\eta }\right)\mathbf{f}\left(\mathbf{y},t+s\right)ds,$$

$$\Vert \mathbf{f}\left({\mathbf{x}}_1,t\right)-\overline{\mathbf{f}}\left({\mathbf{x}}_2,t\right)\Vert \le \beta \Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert .$$

$$M=\underset{\mathbf{x}\in D}{\sup }\underset{0\le t\le L}{\sup }\Vert \mathbf{f}\left(\mathbf{x},t\right)\Vert .$$

$$\Vert \mathbf{x}-\mathbf{y}\Vert \le M\left(1+\frac{1}{2}\beta \epsilon \right)\epsilon \Delta T\eta ,$$

$$\Vert {\partial }_{\mathbf{y}}\mathbf{f}(\mathbf{y},t)\Vert \le M_1,\mathbf{y}(t)\in D\subset {\mathbb{R}}^n.$$

$$\underset{j}{\max }\Vert \frac{{\partial }^j\mathbf{f}}{\partial {\mathbf{y}}_k^j}\Vert \le M,1\le k\le n,0\le j\le p.$$

$$\Vert \mathbf{y}(t)-{\mathbf{y}}_{\Delta T}(t)\Vert \le CM\epsilon \Delta T^3\underset{x\in \mathbb{R}}{\max }\left({\lambda }_n^2\frac{1}{\eta }{\int }_0^{\eta }\rho \left(\frac{s}{\eta }\right)e^{i{\lambda }_{n}s}ds\right),$$

$$\Vert \mathbf{y}(t)-{\mathbf{y}}_{\Delta T}(t)\Vert \le C_t{\left(\Delta T\right)}^{p+1}\underset{t}{\max }{\Vert \frac{{\mathrm{d}}^{p+1}\mathbf{y}}{\mathrm{d}{t}^{p+1}}\left(t\right)\Vert }_2,$$

$$\frac{\mathrm{d}{y}_j}{\mathrm{d}t}=\epsilon \frac{1}{\eta }{\int }_0^{\eta }\rho \left(\frac{s}{\eta }\right)\sum _n\Delta T{e}^{i\Delta T{\lambda }_n(t+s)}{\mathbf{N}}_{n,j}(\mathbf{y})\mathrm{d}s,$$

$$\begin{gathered} \frac{{\mathrm{d}}^3{y}_j(t)}{\mathrm{d}{t}^3}=\epsilon \frac{1}{\eta }{\int }_0^{\eta }\rho \left(\frac{s}{\eta }\right)(i^2\Delta T^3\sum _n{\lambda }_n^2{e}^{i\Delta T{\lambda }_n(t+s)}{\mathbf{N}}_{n,j}+ \\ 2i\Delta T^2\sum _n\sum _k{\lambda }_n{e}^{i\Delta T{\lambda }_n(t+s)}\frac{\partial {\mathbf{N}}_{n,j}(\mathbf{y})}{\partial y_k}\frac{\mathrm{d}{y}_k(t)}{\mathrm{d}t}+ \\ \Delta T\sum _n\sum _{k,l}{e}^{i\Delta T{\lambda }_n(t+s)}\frac{{\partial }^2{\mathbf{N}}_{n,j}(\mathbf{y})}{\partial y_k\partial y_l}\frac{\mathrm{d}{y}_k(t)}{\mathrm{d}t}\frac{\mathrm{d}{y}_l(t)}{\mathrm{d}t}+ \\ \Delta T\sum _n\sum _k{e}^{i\Delta T{\lambda }_n(t+s)}\frac{\partial {\mathbf{N}}_{n,j}(\mathbf{y})}{\partial y_k}\frac{{\mathrm{d}}^2{y}_k(t)}{\mathrm{d}{t}^2})\mathrm{d}s. \end{gathered}$$

$$P(\eta )=\frac{1}{\eta }{\int }_0^{\eta }\rho \left(\frac{s}{\eta }\right)e^{i\Delta T{\lambda }_{n}s}\mathrm{d}s;e_n(t)=e^{i\Delta T{\lambda }_{n}t},$$

$$\frac{\mathrm{d}^{3}y_{j}(t)}{\mathrm{d}t^{3}}=\varepsilon\Delta T^{3}P(\eta)\sum_{n}\left[-\lambda_{n}^{2}e_{n}\mathbf{N}_{n,j}+\left(2i\varepsilon\sum_{k}\lambda_{n}e_{n}\frac{\partial\mathbf{N}_{n,j}}{\partial y_{k}}\right)\left(\sum_{n^{\prime}}e_{n^{\prime}}\mathbf{N}_{n^{\prime},j}\right)+\right.\\ \left.\varepsilon^{2}\sum_{k,l}e_{n}\frac{\partial^{2}\mathbf{N}_{n,j}}{\partial y_{k}\partial y_{l}}\left(\sum_{n^{\prime}}e_{n^{\prime}}\mathbf{N}_{n^{\prime},j}\right)\left(\sum_{n^{\prime\prime}}e_{n^{\prime\prime}}\mathbf{N}_{n^{\prime\prime},j}\right)+\right.\\ \left.\sum_{k^{\prime}}e_{n}\frac{\partial\mathbf{N}_{n,j}}{\partial y_{k^{\prime}}}\Biggl{(}\varepsilon\sum_{n^{\prime}}\lambda_{n^{\prime}}e_{n^{\prime}}\mathbf{N}_{n^{\prime},j}+\varepsilon^{2}\sum_{n^{\prime\prime}}e_{n^{\prime\prime}}\mathbf{N}_{n^{\prime\prime},j}\sum_{n^{\prime\prime\prime}}\sum_{l^{\prime}}e_{n^{\prime\prime\prime}}\frac{\partial\mathbf{N}_{n^{\prime\prime\prime},j}}{\partial y_{l^{\prime}}}\Biggr{)}\right],$$

Full text

1 Introduction

A variation of the Parareal method[28][29] for highly oscillatory systems of equations was proposed in [21]. The method constructs the coarse solver based on a coordinate transformation and fast-wave averaging motivated by multiscale analysis of PDEs [7, 26, 33] and performs the integration using techniques from the Heterogeneous Multiscale Method [14]. They proved that the method provides parallel speedups[21] in the limit of infinite time scale separation. Finite time scale separation is an important case to understand for physical applications because many physical phenomena, such as those occurring in numerical weather prediction, have finite frequencies inherent in the problem, e.g. Earth’s rotation rate is finite. In this paper we extend the work of [21] by showing that rapid convergence of the method is also possible for finite timescale separation by proving error bounds on the coarse solver in the case of finite timescale separation.

Examples of applications of the Parareal algorithm being applied to parabolic PDEs include simulations of financial markets (i.e. the Black-Scholes equation for an American put [5] and a nonlinear parabolic evolutionary equation via the finite element method [22]. Hyperbolic systems solved with Parareal include simulation of molecular dynamics [3], fluid/structure interaction [15], solution of the Navier-Stokes equations [16], and reservoir modelling [20]. In all of these applications, the degree of oscillatory stiffness was not sufficient to impede convergence, but it is known (cf. [1]) to be an issue for the Parareal method.

There have been several modifications to the Parareal method which apply to highly oscillatory systems which assume that the system may be separated into fast and slow variables. In terms of ODEs, [27] have proposed a multiscale method for singularly perturbed ODEs where the fast dynamics are dissipative. Ariel et al (2016)[1] propose a method for highly oscillatory ODEs which is multiscale in nature but does not require explicit knowledge of the fast and slow variables. Gander and Hairer (2014) [19] suggest Parareal methods for Hamiltonian dynamics. Approaches using symplectic integrators with applications to molecular dynamics are presented in, for example, [2] and [6]. Finally, [21] proposed a method which is motivated by a asymptotic solutions for fast singular limits of nonlinear evolutionary PDEs. It is an extension of this method which we study here and which we refer to as Asymptotic Parallel-in-Time (APinT). It takes its name from the modified coarse solver which is inspired by methods used in the asymptotic analysis of PDEs.

In this work we are primarily interested in oscillatory stiffness. Oscillatory stiffness places a restriction on the convergence of the Parareal method due to accuracy and stability limitations it places on the timestep size. We consider oscillatory stiffness to be a phenomenon arising from the presence of rapid oscillations which restricts the coarse timestep, with the degree of oscillatory stiffness being the degree to which the timestep is restricted. As discussed by Higham and Trefethen[23], stiffness is a transient phenomenon involving finite time intervals. This is important here as in this paper we are concerned with mitigating oscillatory stiffness over the interval of the coarse timestep.

We consider as a model equation a PDE of the form:

[TABLE]

where $\mathbf{u}$ is the vector of unknowns, $\mathcal{L}$ is a skew-Hermitian linear operator with purely imaginary eigenvalues, and $\mathcal{N}(\cdot,\cdot)$ is a nonlinear operator which is assumed to be of quadratic type. We further assume that the solution, $\mathbf{u}\in L^{2}$ and that we may approximate eq. 1 as a finite system of ODEs. The linear term then induces temporal oscillations on an $\mathcal{O}(\varepsilon)$ timescale, which can require the use of prohibitively small timesteps for standard numerical integrators if $\varepsilon$ is small. If accuracy is required, as is necessary for the Parareal method, implicit methods must also use a small time step.

The contribution of this work lies in an improved error estimate of the asymptotically motivated coarse solver which permits a mathematical description of the relationship between the time stepping error and the time-averaging of the the nonlinear operator. This understanding leads to an improved mathematical understanding of the convergence of the APinT method. With this improved notion of accuracy, we are then able to prove that the APinT algorithm converges for finite time scale separation. This is an advance on the previously shown case in the limit of $\varepsilon\to 0$.

The slow solution relies on an averaged version of eq. 1, with the average taken over an infinite window in the limit as $\varepsilon\to 0$. This is estimated numerically by a finite sum over a sufficiently large window. It has been shown[12][14][21] that the length required for this window may be reduced through the use of a smooth kernel of integration. In the small-$\varepsilon$ limit, we find that this method provides a convergent algorithm. We have also found that for finite $\varepsilon$ the averaging window may be chosen such that the slow solution is sufficiently accurate that the Parareal method remains convergent, as shown in fig. 1 and discussed in section 4.3.

In the next section we give an overview of the Parareal method to set the context of the work. In section 3 we discuss the slow solution which is found by fast-wave averaging and which forms the coarse propagator of the APinT method. We will discuss the implications of the existence of near resonances on the slow solution. With that in mind, we proceed in section 4 to prove the error bounds and therefore the convergence of the APinT method. Finally, in section 5 we will show numerical experiments on the one-dimensional rotating shallow water equations and discuss some particularities of solving these equations with the APinT method.

2 The Parareal Algorithm

In this section we briefly review the Parareal method proposed by Lions et al.[28], and further expanded upon by Maday and Turinici[29]. They proposed a generalisation of the concept of domain decomposition to the temporal domain, in which a ‘coarse’ approximation to the solution is computed which is then refined, parallel in time, by the ‘fine’ timestep. The solution has been shown to iteratively converge to the fine solution[18]. In practice, this method requires that the coarse timestepping method permits large timesteps, that it be inexpensive to compute and sufficiently accurate that the method converges quickly. In general, the maximum timestep is $\mathcal{O}(\varepsilon)$, so in the case of $\varepsilon=\mathcal{O}(1)$, i.e. the less-stiff case, Parareal may be applied without any modifications. The insight of [21] was that a slow solution based on a coordinate transformation and a time average over the fast waves in the nonlinear operator provides a convergent and efficiently-computable coarse approximation. In fact, they showed that under suitable assumptions of smoothness, superlinear convergence is obtained as $\varepsilon\to 0$.

Also related to this work is that of [30] who paid particular attention to the parallel implementation. They note that the coarse solver may employ a coarser timestep [28], a coarser discretisation [16], and/or a simpler physical model [29]. The APinT coarse solver as presented here us a combination of the first and third of these.

A sketch of the algorithm is in fig. 2. We assume for the sake of simplicity that we are interested in solving eq. 1 on the interval $t\in[0,1]$. Let $\varphi_{t}(\mathbf{u}_{0})$ denote the evolution operator associated with eq. 1 such that $\mathbf{u}(t)=\varphi_{t}(\mathbf{u}_{0})$ solves the full equation. Similarly $\overline{\varphi}_{t}(\mathbf{u}_{0})$ solves the averaged equations.

We then divide the time domain into $N$ finite subintervals, $[n\Delta T,(n+1)\Delta T]$, where $n=0,\ldots N-1$. The Parareal algorithm begins with a coarse solve and then proceeds by computing approximations to the solution, $\mathbf{U}_{n}^{k}$, iteratively, following:

[TABLE]

Here, since the quantities $\mathbf{U}_{n-1}^{k-1}$ in the difference $(\varphi_{\Delta T}(\mathbf{U}^{k-1}_{n-1}-\overline{\varphi}_{\Delta T}(\mathbf{U}^{k-1}_{n-1}))$ are already computed at iteration $k$, the difference can be computed in parallel for all $n$. Since the computation of $\overline{\varphi}_{\Delta T}(\mathbf{U}_{n-1}^{k})$ is cheap, the overall computation is quick in a parallel sense if the iterates converge quickly.

3 The Slow Solution

Our interest in solving general PDEs which arise in physical modelling, in particular those of weather and climate, requires that we confront the problem of oscillatory stiffness over the interval of a coarse time step. In this section we shall describe the mathematical roots of the slow solution, provide a short description of its historical context, and review its discretisation as a coarse solver for the Parareal algorithm. As the convergence of the coarse solver for the APinT algorithm depends on the quality of the approximation, we will then investigate the numerical behaviour of this solver before providing an improved proof of the performance of the slow solver.

As an example of an application where the actual physics was thought to be asymptotic, but later shown not to be we look to the field of numerical weather prediction. Historically, the physical notion of ‘slow’ dynamics, called Quasi-Geostrophic (QG) equations was a major advance in understanding weather. This insight was due to Charney[8] who derived the ‘slow’ equations. These reduced equations allowed the fast waves, which cause the oscillatory stiffness, to be filtered while still resolving the large-scale motions of the fluid. The work was later expanded upon, leading to what is generally recognised as the first successful numerical weather prediction [9] [10]. Since those early days, the reduced equations have been rigorously shown to hold asymptotically in the limit of $\varepsilon\rightarrow 0$[13]. In contrast to these results, modern weather prediction has found that the reduced equations are not accurate enough to be predictive and therefore they rely on numerical approximations of the full equations of motion[11]. As such, we are confronted with the problem that at least some of the oscillations matter even for the large scale flow and so we must find some way to resolve the fast waves in order to capture the full dynamics.

To address this problem for the Parareal method, [21] constructed a numerical approximation to the governing equations eq. 1 based on the consideration of fast singular limits[33]. One of the conclusions of this theory for when the linear operator, $\mathcal{L}$ is skew-Hermitian is that though the leading order dynamics is not slow itself, the slow dynamics evolve independently of the fast, while the fast dynamics are ‘swept’ by the slow [33]. For realistic weather, we expect $\varepsilon\sim 0.01$ to $\varepsilon\sim 0.1$ [34]. Since $\varepsilon$ is finite for realistic cases, the timescale separation is also finite. Working in the limit of small $\varepsilon$ and applying the method of multiple scales, an averaged equation for equations of the form eq. 1 was found by [13].

For a slow timescale, $t$, and a fast timescale, $\tau$, they showed that averaged equations in the asymptotic limit of $\varepsilon\to 0$ must satisfy:

[TABLE]

where $\bar{\mathbf{u}}$ denotes the averaged $\mathbf{u}$ and where the integral is taken over the nonlinear operator, not the solution itself, and where there is a mapping by the exponential of the linear operator between the averaged and ‘full’ solutions:

[TABLE]

The above condition, studied in detail by [7], [26], and [33], motivates our averaging, described below. We follow [21] and write the averaged equation in the following form:

[TABLE]

from which the full solution may be obtained through the application of the matrix exponential as in eq. 4. The averaged equation, eq. 5, has lost the factor of $1/\varepsilon$ from the full equation, eq. 1, the main source of the oscillations, although another derivative will regain this term and so the oscillations have not been entirely eliminated.

In order to use this equation as a coarse solver for Parareal with finite timescale separation, [21] retreated from the asymptotic limit by taking the integral over the nonlinear operator in eq. 3 over a finite time averaging window, rather than the infinite limit associated with $\epsilon\rightarrow 0$. This integral is approximated numerically by using a smooth kernel, $\rho(s)$, $0\leq s\leq 1$ which is chosen such that the length $T_{0}$ of the time window for the averaging is as small as possible, and approximate the averaged nonlinear operator eq. 3 as:

[TABLE]

The Heterogeneous Multiscale Method (HMM) [14] is then applied to the slow equation by computing the averages numerically, as in eq. 6. In practice, the width of the averaging window may be freely chosen, as illustrated in fig. 3. The choice of this window has a significant effect on the convergence of the method, as illustrated in fig. 1. As the computational cost of Parareal is proportional to the number of iterations required, an optimally-chosen window is necessary. In addition to a new error estimate our insight into the coarse error allows us to choose the optimal averaging window length as a function of the timescale separation.

3.1 Triad Resonances

In order to permit a long coarse timestep, the coarse solver proposed by [21] and described in section 3 filters the nonlinear operator. The effect of this is a change in the content of the nonlinear triad interactions which is a function of both the degree of near resonance of the interaction and the length of the averaging window. As we discuss in section 4, the extent to which near-resonant sets are retained or rejected has an important impact on the convergence of the Parareal method. Therefore, in this section we review nonlinear triad resonances for the model problem that we use in our tests in section 5.2, noting that a similar approach applies to all systems of the form eq. 1.

Systems governed by a quadratic nonlinearity with dispersive waves exhibit triad resonances[13][17][24]. Since we are motivated by geophysical modelling, we consider the general averaged equation for the rotating shallow water equations, which are commonly used as a test case for geophysical solvers. Following the notation of [13], we decompose the right-hand side of eq. 3 in terms of its basis of eigenvectors and write:

[TABLE]

where $\Omega_{\mathbf{k},\mathbf{k}{1},\mathbf{k}{2}}^{\alpha,\alpha_{1},\alpha_{2}}=\omega_{\mathbf{k}_{1}}^{\alpha_{1}}+\omega_{\mathbf{k}_{2}}^{\alpha_{2}}-\omega_{\mathbf{k}}^{\alpha}$, $\alpha=-1,0,1$ refers to the different branches of the eigenvalues, $\mathbf{k}$, $\mathbf{k}_{1}$, and $\mathbf{k}_{2}$ are the wavenumbers, $\omega_{\mathbf{k}}^{\alpha}$ is the dispersion relation at a given $\alpha$ and wavenumber, $\sigma$ denotes the Fourier coefficient in this basis, $\mathbf{r}_{\mathbf{k}}^{\alpha}$ is the right eigenvector of the linear operator, and $C_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}}^{\alpha_{1},\alpha_{2},\alpha}$ is an interaction coefficient [31].

In the asymptotic case, the limit $\varepsilon\to 0$, by the orthogonality of the Fourier series, the only waves which remain after the wave averaging procedure (i.e. where $\tau\to\infty$) are the direct three-wave resonances (cf. [13], [31], [35]), i.e. the elements of the resonant set, $\mathcal{S}_{\mathbf{k},\alpha}$ i.e.:

[TABLE]

The wave-averaged solution then follows:

[TABLE]

It is this three-wave resonance condition from which the behaviour of the averaging kernel can be understood. In the limit as $\tau\to\infty$, only the direct resonances should remain. Because we have finite timescale separation, this integral is approximated over a finite averaging window which must be large enough to filter the non-resonant triad.

As shown in [21], using a finite time-averaging window and numerically integrating with respect to a smooth, finitely supported kernel permits this algorithm to result in a convergent Parareal algorithm in the limit of small $\varepsilon$.

For finite $\varepsilon$ we take a finite average and so the solution set is larger than the direct resonant set and this has an important effect on the convergence of Parareal. To better explain the finite-$\varepsilon$ case, we define concentric shells of near-resonances, i.e. we rewrite the triad-based form eq. 7 as:

[TABLE]

where $\mathcal{S}_{\mathbf{k},\alpha}^{\epsilon_{\beta}}$, $\beta=1,2,\ldots$ refers to a near-resonant set, i.e.:

[TABLE]

where $\epsilon_{0}=0$ by definition. The direct-resonant set results in a solution consisting of only the slow dynamics of the system – which was shown in [13] to be equivalent to the reduced equations for this system. Again, and as we will see in section 4, the extent to which the near-resonant sets are retained and rejected by the averaging procedure is fundamental to the convergence of the APinT variation of the Parareal method.

4 Error Bounds

Now that we have discussed the key elements of the algorithm, we are in a position to discuss and prove convergence for the case when $\varepsilon$ is finite. We shall first construct an improved error estimate for the coarse solution, and then use that result to prove the convergence of the Parareal method.

4.1 A bound on the errors due to time-stepping and time-averaging in the coarse solver

In this section we employ the idea of near-resonant sets to extend the existing proof of APinT convergence[21] to the case of finite $\varepsilon$. As has been demonstrated above (cf. fig. 1), the choice of the averaging window width, $\eta$, has a profound effect on the convergence of the method. While the choice of $\eta$ is well-understood for the limit of small $\varepsilon$[21], we show here that $\eta$ may be similarly chosen to provide convergence for $\varepsilon$ up to $\mathcal{O}(1)$ for an appropriate coarse timestep. We first reduce eq. 1 to a standard form for ODEs. Following section 3, we write:

[TABLE]

i.e. we are interested in the solution over a $\Delta T$ timescale. Let $\tau=t/(\varepsilon\Delta T)$, and so $\tilde{\mathbf{v}}\left(\tau\right)$, defined on the interval $\left[0,1/\varepsilon\right]$,

[TABLE]

Then differentiation gives,

[TABLE]

Upon this substitution into the coarse solver eq. 5 over the discrete time interval eq. 13, we arrive at the desired form which permits us to use the framework given in [32] where they have derived bounds for averaging methods. The aim of this is to modify and reapply their result for the error bound due to averaging, which holds on a general dynamical system of finite ODEs. This averaging error is one of the two major sources of error in the timestepping of the coarse solver. We then write the coarse solver in the form:

[TABLE]

While our interest is in solving PDEs describing physical systems, in practice we employ a Fourier spectral method, which has the effect of treating the PDE as a finite-dimensional system of ODEs. In this paper we show that the APinT method is convergent for finite systems of ODEs. This gives us access to the machinery of the numerical analysis of ODEs and averaging methods, following [32]. Let $\mathbf{x}$ solve the governing equations when they are written as a system of ODEs, i.e. in the form shown in eq. 16. For example, in the numerical experiments given in section 5.2, $\mathbf{x}$ is the Fourier solution. Then we may write:

[TABLE]

Similarly, we consider the coarse solver eq. 16 written as a system of ODEs. Let $\mathbf{y}$ solve this averaged form of eq. 16, i.e.:

[TABLE]

where the averaging follows directly from the averaged equation, eq. 6 and is written:

[TABLE]

where $\eta$ denotes the finite length of the averaging window.

Lemma 1**.**

Considering the initial value problems in $\mathbf{x}$ and $\mathbf{y}$ as stated above where $\mathbf{f}$ is $\mathbb{R}^{n}\times\mathbb{R}$ Lipschitz continuous with constant $\beta$ in $\mathbf{x}$ on $D\subset\mathbb{R}^{n}$ and $t$ on an $\mathcal{O}(1)$ timescale, i.e. for all $\mathbf{x}_{1},\mathbf{x}_{2}\in D$, $\beta$ is such that:

[TABLE]

Let:

[TABLE]

Then we can bound the difference between the exact solution $\mathbf{x}$ and the averaged solution $\mathbf{y}$ as:

[TABLE]

The above lemma follows from a modification of Lemma 4.2.8 in [32] in order to include the kernel of integration (cf. appendix A). We have here bounded the error over an $\mathcal{O}(1)$ time interval instead of $\mathcal{O}(1/\varepsilon)$ so that the rate of convergence at different degrees of scale separation may be more easily compared, as in practice we are interested in simulations over fixed timescales. Taking the unmodified lemma provides a slightly different result as it gives the averaging error over a simulation time which scales with $\varepsilon$. Due to the numerical nature of the proof here, the appropriate timescale is over a coarse timestep.

eq. 22 places a bound on the error committed by averaging over the fast waves, independent of the numerical methods used for spatial or temporal discretisation. Next, we consider the error arising from the numerical approximation of eq. 18. In doing so, we will need to assume bounds on $\mathbf{f}$ in the region of phase space where $\mathbf{y}$ exists. Then assume that

[TABLE]

We assume that such a bound exists for higher spatial derivatives of $\mathbf{f}$, such that

[TABLE]

Lemma 2**.**

Denote the numerical approximation to the averaged solution $\mathbf{y}(t)$ with timestep $\Delta T$ by a second-order timestepping method as $\mathbf{y}_{\Delta T}(t)$. Assume that $\mathbf{y}(t)=\varepsilon\mathbf{\overline{f}}(\mathbf{x},t)$ and that $\mathbf{\overline{f}}\in D$ as in eq. 23. Assume that integration is performed with respect to a smooth kernel, $\rho(\cdot)$, and let $\lambda_{n}$ denote the $n$-th near resonant triad (cf. section 3.1). Then the local time-stepping error of a second order time-stepping scheme applied to eq. 18 satisfies:

[TABLE]

for some constant, $C\in\mathbb{R}<\infty$ and where $M$ is the bound over the nonlinear operator as given in eq. 22.

Proof.

The timestepping error of a $p$-th order scheme is bounded by [25]:

[TABLE]

First, decompose $\mathbf{f}$ in terms of its basis of eigenvectors as discussed in section 3.1. As with eq. 7, we may write the solution as a sum of ODEs, each for a specific resonant nearness, $\lambda_{n}$. Then for the $j$-th component of $\mathbf{y}$, we write

[TABLE]

where the nearness of the resonances in any particular ODE is exposed through the eigenvalue sum, $\lambda_{n}$, in the exponent and where the subscript $,j$ denotes the $j$-th component and not a derivative, as it would with Einstein’s notation. We then seek the third time derivative, which is found to be

[TABLE]

This is then the right-hand side which is integrated with respect to the smooth kernel. The magnitude of the near-resonant triad, $\lambda_{n}$, now presents itself as a multiplier on the complex exponential. It is then clear that it is this value, which is zero for direct resonances but becomes large in general, which is the source of numerical stiffness. For convenience, we introduce

[TABLE]

then

[TABLE]

In bounding the timestepping error, we are interested in the norm of this quantity. Recalling that we are working with a finite-dimensional system of ODEs and applying the triangle and Cauchy-Schwarz inequalities we find that

[TABLE]

Now, as $\mathbf{N}$ and all of its spatial derivatives up to and including $p=2$ are bounded by $M$ by eq. 23, we write

[TABLE]

where $C_{f}$ is a positive constant. We will now assume that $|\lambda_{n}|\neq 0$ as we are interested in the sup-norm of these values, which is nonzero when near-resonances are included. The directly resonant case has been treated by [21]. Then we must consider two possibilities. Firstly, if $|\lambda_{n}|\leq 1$, then we define some constant, $K_{1}$,

[TABLE]

If $|\lambda_{n}|>1$, the binomial theorem yields:

[TABLE]

And then we may write

[TABLE]

As for the Rotating Shallow Water Equations there must always be a value of $\lambda_{n}$ which is strictly greater than one, we shall assume that it is the second value which is the maximum. We now let $C=C_{t}K$. Finally, we bound the nonlinear term in the same fashion as eq. 22, where the fact that:

[TABLE]

completes the proof by providing an upper bound for the nonlinear operator as in eq. 22. This provides a bound for the error due to timestepping which does not depend directly on the solution, but rather on the general properties of the nonlinearity, in particular the triadic interactions. ∎

With these two lemmas describing our primary sources of error, we will seek a bound on the error in the APinT algorithm and use this to show convergence. From eq. 25, it follows that the timestepping error depends on:

[TABLE]

We define the following term which describes the filtering, independent of the gain due to the scale separation and the coarse timestep, and which is the key insight into understanding how to regularise an oscillatory problem over a finite time interval.

[TABLE]

$\Lambda(\eta)$ provides a measure of the extent to which the averaging integral mitigates the numerical stiffness. Recall that when the maximum $\lambda_{n}$ is large, as it is for highly oscillatory problems, it contributes to large gradients on the right-hand side requiring a small numerical timestep. In contrast, the integral component tends to zero as $\lambda_{n}$ gets large, and does so superlinearly because of the smooth kernel, $\rho(s)$ [21]. This term is then where we see precisely how the averaging procedure filters the fast oscillations, causing $\Lambda(\eta)$ to achieve a lower magnitude than $\lambda_{n}^{2}$ does on its own and therefore reducing the numerical stiffness.

In seeking a bound on the error in the timestepping, it is necessary to bound this term for some particular averaging kernel, $\rho(s)$. The choice of averaging kernel affects the error bounds through this function. $\Lambda(\eta)$ is bounded and tends rapidly to zero as $\eta\to\infty$ (q.v. section 5.3).

With this in mind, we now prove theorem 1 which bounds the error committed by the coarse timestepping as compared to the fine. This will later allow us to prove error bounds on APinT subject to finite timescale separation.

Theorem 1**.**

Let $\Delta T$ denote the coarse timestep for a second order numerical method. We assume a finite scale separation on the order of $\varepsilon$. For an averaging window of length $\eta$, the total error in the coarse timestepping for the APinT algorithm is bounded by:

[TABLE]

where $M$ is the sup-norm over the nonlinear operator as in lemmas 22 and 25 and $C_{0}$, $C_{1}$, and $D_{1}$ are finite constants.

Proof.

By the triangle inequality, we may write:

[TABLE]

eq. 22 is used to bound the first term, i.e.:

[TABLE]

Applying eq. 25 and eq. 37 to the second term yields:

[TABLE]

where $\Lambda(\eta)$ is bounded independently of $\lambda_{n}$ for any averaging window length, $\eta$. Combining the bounds in equations eq. 39 and eq. 41 gives the theorem as desired.

∎

4.2 Proof of Parareal Convergence

We may now derive error bounds for the Parareal iteration on finite systems of ODEs, given in eq. 2. Using our improved error bound for the coarse solver which holds for finite $\varepsilon$, we modify the proof given in [21], which held only as $\varepsilon\to 0$. For consistency we define several operators following [21]. Let $\tilde{\varphi}_{\Delta T}(\cdot)$ be the evolution operator associated with numerically solving the slow equation using an $\mathcal{O}(p)$ method, such that $\tilde{\varphi}_{\Delta T}(\cdot)$ is a numerical approximation of $\overline{\varphi}_{\Delta T}(\cdot)$. Furthermore, let $\varphi_{\Delta T}(\cdot)$ denote the evolution operator for the fine solution. We then define:

[TABLE]

Then, as in [4], [21], and section 4.1 we make the following assumptions:

The operators $\varphi(\cdot)$ and $\overline{\varphi}(\cdot)$ are uniformly bounded for $0\leq t\leq 1$:

[TABLE] 2. 2.

The averaging method is accurate in the sense that:

[TABLE] 3. 3.

The averaged evolution operator satisfies:

[TABLE]

and the numerical approximation to the evolution equation satisfies:

[TABLE] 4. 4.

Following eq. 22 and eq. 25 and eq. 43, the error operators satisfy:

[TABLE]

and:

[TABLE]

We have now quantified the major sources of error in the coarse timestepping which will affect the convergence of Parareal. The following proof of the convergence follows directly from these bounds.

Theorem 2**.**

Subject to the above assumptions, the error, $\mathbf{u}(T_{n})-\mathbf{U}_{n}^{k}$, after the $k$-th Parareal iteration is bounded by:

[TABLE]

Proof.

The proof is by induction on $k$. When $k=0$:

[TABLE]

where we have used eq. 44, which bounds the error induced by the averaging procedure, to bound the first term and eq. 25, which governs the timestepping error, for the second. Now assume that:

[TABLE]

We may then write the Parareal iteration, eq. 2 in the following form, using eq. 47 and eq. 48:

[TABLE]

By directly substituting equations eq. 46, eq. 47, and eq. 48, we have:

[TABLE]

Finally, application of the discrete Gronwall inequality gives:

[TABLE]

∎

eq. 49 is one of the key contributions of this work. Using the understanding of near-resonance and the result of theorem 1, it generalises the proof given by [21] of convergence for the asymptotic limit as $\varepsilon\to 0$ to finite $\varepsilon$. This is a significant improvement as for many physical applications such as weather and climate modelling $\varepsilon$ remains finite. As the averaging window length, $\eta$, may be freely chosen we may select an optimal $\eta$ for a wide range of $\varepsilon$ subject to the other constants and choice of $\Delta T$ such that the method is convergent. We discuss this in the next section.

4.3 Convergence for any $\varepsilon$

Given eq. 49 we are in finally in a position to discuss convergence for any timescale separation. For the APinT algorithm to converge, we require that:

[TABLE]

We are then left with the problem of choosing an appropriate averaging window length, $\eta$, depending on the degree of scale separation, $\epsilon$, and the filtered contribution of the triads, $\Lambda(\eta)$. In the interest of demonstrating that one exists, we assume the scaling (for example):

[TABLE]

We then have:

[TABLE]

as $\varepsilon\to 0$, our error also decreases for any value of the power $s$. $\Lambda\left(\frac{\Delta T}{\varepsilon^{s}}\right)$ is bounded, so as $\varepsilon\to 1$, all terms remain bounded and we may choose our coarse timestep accordingly to ensure convergence. This means that the method proposed here may be applied across the full range of $\varepsilon\in(0,1]$ with only a change of averaging window length, which allows convergence for physical problems where the time scale separation may change throughout the computation. This is in contrast to the proof in the limit [21] which proved convergence only for $\varepsilon\to 0$.

5 The One-Dimensional Rotating Shallow Water Equations

We now consider an example, using the one-dimensional rotating shallow water equation as a test-case, as did [21]. Let the unknown vector be:

[TABLE]

We then write the linear and nonlinear operators in the full model eq. 1 as:

[TABLE]

for some constant, $F\in\mathbb{R}$. The corresponding eigenvalues are:

[TABLE]

In general, as $\varepsilon\to 0$ we expect that $\Lambda(\eta)\to 0$ as well due to cancellation of oscillations in the integral[13]. As discussed in section 4, oscillatory stiffness arises due to the magnitude of the gain term outside of the integral, which is large for highly oscillatory systems. The integral itself, however, is bounded from above by one, and achieves this value only for directly resonant triads (cf. section 3.1), where the gain is zero. As the distance of resonance (i.e. the magnitude of $|\omega_{k}^{\alpha}-\omega_{k_{1}}^{\alpha_{1}}-\omega_{k_{2}}^{\alpha_{2}}|$) increases, the integral tends to zero as well (and does so faster with larger $\eta$).

The choice is then for a given degree of scale separation, $\varepsilon$, to choose an $\eta$ which mitigates the stiffness sufficiently to allow the necessary coarse timestep, while retaining as much fidelity to the full equations as possible (cf. eq. 22, where the averaging error is proportional to $\eta$). We shall deal with the practical implications of the form of $\Lambda(\eta)$ in section 5.3.

5.1 The Optimally-Averaged Slow and Fast Solutions

In order to illustrate the slow averaged solution over which the timestepping is performed and its relation to the full solution, fig. 4 compares the slow, full and true solutions for the stiff case where $\varepsilon=0.01$. The spatio-temporal oscillations are very rapid in the stiff case, which is the source of the timestep limitation by the CFL condition. However, the slow solution over which the timestepping is performed lacks these rapid oscillations and so permits the large timestep.

Spatio-temporal oscillations from an initially stationary Gaussian height field are shown on a domain which is spatially periodic, i.e. the top and bottom boundaries of the plots wrap around. The decay of the height field into waves travelling in opposite directions is visible in fig. 4. The optimal averaging window (q.v. section 5.3) was applied, and convergence to single precision was obtained in six iterations.

5.2 Numerical Results on the 1-D RSWE

In this section we present numerical results for the one-dimensional rotating shallow water equations which build on those presented in [21]. Figure 5 shows the norm of the coarse error, i.e. $\|\mathbf{x}(t)-\mathbf{y}(t)\|_{2}$ computed relative to the fine timestep versus the width of the averaging window, $T_{0}$, which denotes the numerical choice of $\eta$, where $0\leq t\leq 1$, $\delta t=2\mathrm{e}{-4}$, and the spatial resolution is $N_{x}=64$. This spatial resolution and fine timestepping regime were found to be within the asymptotic range of the timestepping. A second-order Strang splitting method was used for both the coarse and fine solves. The initial flow was stationary with a Gaussian height field.

For the smallest $\varepsilon$ the asymptotic behaviour is well approximated, as the fidelity of the coarse timestepping increases as the averaging window increases. This is consistent with the behaviour described in section 3.1, where the theory predicts that $\eta\to\infty$ as $\varepsilon\to 0$. However, for larger $\varepsilon$ such as the two cases shown, there is a clear optimal size for the averaging window to take, i.e. the minimum in the red and green curves in fig. 5. The location but not the magnitude of this point is predicted by eq. 62.

Theorem 1 states that we should expect that outside of the small-$\varepsilon$ limit the iterative error should decrease with $k$ as in eq. 49 and exhibit a minimum where the sum of the timestepping and averaging errors is smallest. In the case of $\varepsilon=0.01$, which is near the limit as $\varepsilon\to 0$, we expect that for a large enough averaging window we will have optimal convergence, with no improvement in solution quality for a larger averaging window. The numerical results are then consistent with the theory developed in this paper.

In practice, we seek a choice for $T_{0}$ for which the solution is non-stiff on an $\mathcal{O}(\Delta T)$ interval, and therefore as $\Delta T$ increases, so must the averaging window. Similarly, the oscillatory stiffness is proportional to $1/\varepsilon$, and so as $\varepsilon\to 0$ it is necessary to choose a longer averaging window, and therefore to apply stronger smoothing to the solution.

Comparing fig. 5 to fig. 6, which shows the iterative error in the APinT method after three iterations for the same parameters, the direct computation of the coarse timestepping error provides good qualitative agreement with the optimal choice of $\eta$ for the different values of $\varepsilon$. This is in direct agreement with the prediction of eq. 49.

As $T_{0}$ is taken smaller, instability is observed for all $\epsilon$. This corresponds to the explicit CFL limit being violated, as reducing the length of the averaging window increases the maximum wave speed in the solution. For large $T_{0}$, the iterative error roughly stabilises for finite $\varepsilon$. In the limit as $\eta$ is taken very large, the coarse timestepping corresponds to an incorrect equation (e.g. as in the QG equations for our example) being solved in a numerically stable fashion. The difference in the coarse and fine equations is sufficient to inhibit convergence, but does not violate the timestepping limit.

5.3 Optimal Averaging for the 1-D RSWE

It was shown in section 4.3 that it is possible to choose the averaging window in such a way as to ensure convergence. Beyond doing this, we may choose the window optimally to obtain the fastest possible convergence (cf. fig. 6).

The reason we are able to describe the qualitative behaviour this way is a direct result of the Parareal algorithm and the sources of error present in it. The Parareal method consists of an initial approximation to the solution performed by the coarse solver which is accurate to within the coarse error predicted by theorem 1. This is followed by a series of parallel-in-time corrections which converge to the full solution, derived from the difference between the coarse and fine solutions. The closer the initial approximation is to the solution, the less correction is required to converge.

The optimal choice of $\eta$ may be written as an optimisation problem:

[TABLE]

for some as-yet unknown constants $C_{1}$, $C_{2}$, and $C_{3}$. It is here that the fact that both the timestepping and averaging errors are bounded proportionally to the norm of the nonlinear term, $M$, becomes serendipitous, as this constant may be somewhat non-optimal in practice. In seeking the location, but not the magnitude, of the minimum coarse error, the bound on the norm of the nonlinear operator plays no role. Seeking stationary points with respect to $\eta$, this then requires us to find $\eta$ such that:

[TABLE]

The result of eq. 59 is used to choose the optimal averaging window. This result captures the relationship of parameters such as the timestep and the scale separation on the optimal averaging, but relies on several unknown constants. If these constants $C_{n}$ were known, the optimal averaging window could be determined computationally. Equation 59 would then provide an approximation to the optimal window. Given some initial data of the type shown in fig. 6, these constants may be fit by least-squares. Doing so fits the known trend to the known data, and permits the optimal averaging window to be recomputed ‘on the fly’ in a computation.

Certain practical issues arise in the computation of $\eta$. Firstly, the computation of $\frac{\mathrm{d}\Lambda}{\mathrm{d}\eta}$ requires all triads to be investigated, i.e. the maximum is taken over the set of all near-resonant sets. Doing so is computationally expensive, although if this computation were to be performed infrequently the cost could be negligible compared to the simulation cost. Additionally, finding $\eta$ requires solving a transcendental equation in at least two variables ($\eta$, $\varepsilon$), both for the initial fitting of constants, and for the optimisation on the fly. We therefore propose a simpler model based on the behaviour of $\Lambda(s)$.

Restricting ourselves for this example to a Gaussian kernel, we may consider the asymptotic behaviour of the kernel as $\lambda_{n}$ is large. This gives:

[TABLE]

We then multiply our approximation by $1=\eta^{2}/\eta^{2}$, to obtain:

[TABLE]

for some constant, $D_{1}$, since $x^{2}e^{-x^{2}}$ is bounded independently of $x$ (cf. 7). We then replace our first term in eq. 58 and seek fixed points corresponding to the minimum error. This yields:

[TABLE]

This equation provides an estimate for the optimal averaging window length, $\eta_{\mathrm{opt}}$, in terms of the computational parameters and the empirically-fit constants. This result is consistent with eq. 49, as it exhibits a clear minimum for $\mathcal{O}(1)$ values of $\varepsilon$, with the optimal averaging window increasing as $\varepsilon\to 0$, as the asymptotic theory predicts. Both approximations are shown in fig. 8 for a set of minima extracted from a series of runs of the algorithm.

The full model given in eq. 59 provides a much closer approximation both to the behaviour for $\varepsilon=1$ and as $\varepsilon\to 0$, and as an actual fit to the points. It does this, however, at the cost of several orders of magnitude more computational difficulty. The simple model of eq. 62, on the other hand, provides a reasonable approximation to the error as a function of $\varepsilon$, but has the disadvantage of poorly resolving the trend in the limit as $\varepsilon\to 0$. While the simple prediction underestimates the optimal as $\varepsilon\to 0$, the behaviour in this range is well-understood (cf. [21], [1]) and so a hybrid model may easily be applied in practice.

6 Conclusion

We have investigated the convergence of a Parareal method using the APinT coarse solver, which provides a technique by which oscillatory-stiff equations may be solved with the Parareal method. The convergence of this method is due to the averaging applied to the coarse solution, which filters the fast waves and mitigates the oscillatory stiffness present in many of the equations of mathematical physics. This averaging must be performed over the entirety of the nonlinear operator due to the role the direct and near-resonances play in the oscillatory stiffness of the system.

By describing the error of the coarse solver in terms of the interplay between the average over the rapid oscillations and the timestepping, we show the method converges for finite scale separation, significantly extending the domain of applicability for this method.

We have shown here that this method is convergent across a wide range of scale separation, which is an improvement on the prior result[21] which held only in the small-$\varepsilon$ limit. Further, in section 5.3 we considered both a full and a reduced model to predict the optimal averaging window in practical codes.

Appendix A Proof of eq. 22

Consider

[TABLE]

and its averaged version

[TABLE]

where with $\eta^{\prime}=\varepsilon h\eta$,

[TABLE]

Then

[TABLE]

To prove this, change variables: $\tau=t/\left(h\varepsilon\right)$ and

[TABLE]

[TABLE]

Then

[TABLE]

Thus,

[TABLE]

and, with $h\tau^{\prime}=s/\varepsilon$, $ds=\varepsilon hd\tau^{\prime}$,

[TABLE]

Thus,

[TABLE]

Define

[TABLE]

Lemma 3**.**

If $\phi\left(t\right)$ is Lipschitz-continuous with Lipschitz-constant $\lambda$. Then

[TABLE]

where

[TABLE]

Proof.

Using that

[TABLE]

we have that

[TABLE]

∎

Lemma 4**.**

Consider

[TABLE]

with $\mathbf{f}$ continuous in each argument. Also assume that

[TABLE]

and

[TABLE]

Then defining

[TABLE]

we have that

[TABLE]

Proof.

We calculate that

[TABLE]

where

[TABLE]

and

[TABLE]

In the last inequality, we used that $0\leq t\leq\varepsilon^{-1}$.∎

Consider

[TABLE]

with the same assumptions as in the previous lemmas. Let

[TABLE]

Then

[TABLE]

Note that

[TABLE]

By lemma 4,

[TABLE]

where

[TABLE]

Therefore,

[TABLE]

where

[TABLE]

Also, since

[TABLE]

we have that

[TABLE]

Finally, by Gronwall’s inequality,

[TABLE]

Acknowledgments

We would like the acknowledge the support of the University of Exeter and Los Alamos National Laboratory.