Removing numerical dispersion from linear evolution equations

TL;DR

This paper introduces a novel method using Fourier integral operators to eliminate numerical dispersion errors in linear evolution equations, improving the accuracy of simulations over time.

Contribution

It presents a new approach employing time dispersion transforms to correct numerical errors caused by finite difference approximations in linear evolution equations.

Findings

The method effectively removes numerical dispersion in model equations.

It improves the accuracy of elastic and viscoelastic wave simulations.

The approach maintains correct evolution throughout the entire simulation lifespan.

Abstract

We describe a method for removing the numerical errors in the modeling of linear evolution equations that are caused by approximating the time derivative by a finite difference operator. The method is based on integral transforms realized as certain Fourier integral operators, called time dispersion transforms, and we prove that, under an assumption about the frequency content, it yields a solution with correct evolution throughout the entire lifespan. We demonstrate the method on a model equation as well as on the simulation of elastic and viscoelastic wave propagation.

Figures

Tables

Table 1. Table 1. The central finite difference weights used to compute the spatial first-order derivatives, truncated to 4 digits.

${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\alpha }_4$	${\alpha }_5$	${\alpha }_6$
1.2508	$-$0.1203	0.0321	$-$0.0101	0.0030	$-$0.0007

Equations

$$u^{\prime }(t)+Au(t)$$

$$u(t)$$

$$2\pi i\omega \hat{u}(\omega )+A\hat{u}(\omega )=\hat{f}(\omega ),$$

$$\frac{v(t+\Delta t)-v(t-\Delta t)}{2\Delta t}+Av(t)$$

$$v(t)$$

$$2\pi iq(\omega )\hat{v}(\omega )+A\hat{v}(\omega )=\hat{f}(\omega ),$$

$$q(\omega )=\frac{1}{2\pi i}\frac{e^{2\pi i\omega \Delta t}-e^{-2\pi i\omega \Delta t}}{2\Delta t}=\frac{\sin (2\pi \omega \Delta t)}{2\pi \Delta t}.$$

$$2\pi iq(\omega )\hat{u}(q(\omega ))+A\hat{u}(q(\omega ))=\hat{f}(q(\omega )).$$

$${\mathcal{P}}_{i}u(t,x)={\partial }_t^{n_i}{u}_i(t,x)+\sum _{j=0}^{n_i-1}\sum _{k=1}^K{\partial }_t^j{L}_{ijk}{u}_k(t,x),1\le i\le K,$$

$${\mathcal{P}}_{i}u(t,x)$$

$$u(t,x)$$

$${\partial }_t^{2}u(t,x)-c(x){\partial }_x^{2}u(t,x)=f(t,x),t>0,0\le x\le 1,$$

$${\partial }_t^{2}u(t,x_j)-c(x_j)\frac{u(t,x_{j-1})-2u(t,x_j)+u(t,x_{j+1})}{\Delta x^2}=f(t,x_j)$$

$$L_0=-\frac{1}{\Delta x^2}\left(\begin{array}{cccc}[r]c(x_1)& & & \\ & c(x_2)& & \\ & & \ddots & \\ & & & c(x_N)\end{array}\right)\left(\begin{array}{cccc}[r]-2& 1& & \\ 1& -2& 1& \\ & \ddots & \ddots & \ddots \\ & & 1& -2\end{array}\right).$$

$$t\mapsto \underset{x\in X}{\sup }|u_i(t,x)|\in H^{n_i}(\mathbb{R}).$$

$${\hat{u}}_i(\omega ,x)=\mathcal{F}(u_i(\cdot ,x))(\omega )={\int }_{-\infty }^{\infty }{e}^{-2\pi it\omega }{u}_i(t,x)dt.$$

$$\mathsf{D}{v}_i(t,x)=\sum _{n\in Z}{c}_{1,n}{v}_i(t+n\Delta t,x),$$

$${\mathsf{P}}_{i}v(t,x)={\mathsf{D}}^{n_i}{v}_i(t,x)+\sum _{j=0}^{n_i-1}\sum _{k=1}^K{\mathsf{D}}^j{L}_{ijk}{v}_k(t,x),1\le i\le K,$$

$${\mathsf{D}}^j=\underset{j}{\underbrace{\mathsf{D}\circ \dots \circ \mathsf{D}}},$$

$$\int_{-\infty}^{\infty}e^{-2\pi it\omega}\mathsf{D}^{j}v_{i}(t,x)\,dt=\bigg{(}\sum_{n}c_{1,n}e^{2\pi i\omega\Delta tn}\bigg{)}^{j}\widehat{v}_{i}(\omega,x).$$

$$\mathcal{F}\left({\partial }_t^j{u}_i(\cdot ,x)\right)(\omega )=(2\pi i\omega {)}^j{\hat{u}}_i(\omega ,x),$$

$$q(\omega )=\frac{1}{2\pi i}\sum _n{c}_{1,n}{e}^{2\pi in\omega \Delta t}$$

$$\mathcal{F}\left({\mathsf{D}}^j{v}_i(\cdot ,x)\right)(\omega )=(2\pi iq(\omega ){)}^j{\hat{v}}_i(\omega ,x).$$

$$q_0(\eta )=q(\eta /\Delta t)\Delta t=\frac{1}{2\pi i}\sum _n(c_{1,n}\Delta t)e^{2\pi in\eta }$$

$$\mathsf{D}v(t)=\frac{v(t+\Delta t)-v(t-\Delta t)}{2\Delta t},$$

$${\mathsf{D}}^{2}v(t)=\frac{v(t+\Delta t)+v(t-\Delta t)-2v(t)}{\Delta t^2},$$

$$\mathsf{D}v(t)=\frac{v(t+\Delta t/2)-v(t-\Delta t/2)}{\Delta t}.$$

$$\frac{u(t+\Delta t,x)-u(t,x)}{\Delta t}=\frac{u(t,x+\Delta x)-2u(t,x)+u(t,x+\Delta x)}{\Delta x^2}$$

$$\frac{{\mathcal{F}}_x(u(t+\Delta t,\cdot ))(\xi )-{\mathcal{F}}_x(u(t,\cdot ))(\xi )}{\Delta t}=\frac{e^{2\pi i\xi \Delta x}-2+e^{-2\pi i\xi \Delta x}}{\Delta x^2}{\mathcal{F}}_x(u(t,\cdot ))(\xi ).$$

$$\mathcal{F}_{x}(u(t+\Delta t,\cdot))(\xi)=\bigg{[}1+2\frac{\Delta t}{\Delta x^{2}}(\cos(2\pi\xi\Delta x)-1)\bigg{]}\mathcal{F}_{x}(u(t,\cdot))(\xi).$$

$$1+2\frac{\Delta t}{\Delta x^2}(\cos (2\pi \xi \Delta x)-1)=1+\Delta t(2\pi iq(\xi ){)}^2,$$

$$\mathcal{T}(f_i)(t,x)={\int }_{\Omega }{e}^{2\pi it\omega }{\hat{f}}_i(q(\omega ),x)d\omega .$$

Full text

1. Introduction

The difference between a continuous differential equation and its discretized counterpart is a source of numerical artifacts. Generally, the discretized system differs from the intended system in its dispersive and dissipative properties, so errors in the computation are referred to as numerical dispersion and numerical dissipation [23]. Here dispersion refers to a process in which energy separates into its component frequencies as the solution evolves, while dissipation refers to damping of energy during the evolution. Numerical dispersion thus refers to phase errors, while numerical dissipation refers to amplitude errors. The combined effect of the two numerical errors is sometimes described as numerical diffusion.

Numerical diffusion errors are typically studied through the local truncation error, i.e., the consistency between the discrete and continuous equation in terms of the discrete step size. If the method is stable, the Lax equivalence theorem [17] implies that the discretized equation converges to its continuous counterpart. As a consequence, the majority of the numerical methods for differential equations are designed with the intent of minimizing the local truncation error, with the expectation that the global error will then also be small. Examples are high-order accurate derivative schemes [11], [7] and high-order accurate integration schemes such as Runge-Kutta or ADER (Arbitrary high order schemes using DERivatives) [10], [26], [9], [18]. The high-order techniques typically lead to more accurate results compared to low-order methods, but come with a trade-off in increased computational cost.

In this paper we analyze numerical dispersion errors not through the local error but by comparing numerical solutions to true solutions. To illustrate, consider the evolution equation

[TABLE]

where $A$ is a linear operator independent of time $t$. Arguing formally we find by taking Fourier transforms with respect to $t$ of both sides of (1.1) that

[TABLE]

where $\widehat{f}(\omega)=\int_{-\infty}^{\infty}e^{-2\pi it\omega}f(t)\,dt$. Suppose now that we want to simulate the solution $u$ by obtaining a solution $v$ of the finite difference equation

[TABLE]

The finite difference approximation of the time derivative introduces an error which can be expressed by taking Fourier transforms and noting that $\widehat{v}$ will not satisfy (1.3) but instead

[TABLE]

where

[TABLE]

Such functions appear as (part of) the so-called amplification factor in von Neumann stability analysis of spatial finite difference operators [18], but here it is rotated from the spatial side into the time direction where it will be called a phase shift function. The connection is explored in Remark 2.4 below. As the name signifies, comparison between $q(\omega)$ and $\omega$ allows for a description of the numerical dispersion (i.e., phase) errors resulting from the finite difference approximation of the time derivative. We now make the simple observation that if (1.3) is evaluated at $q(\omega)$ instead of at $\omega$ then $\widehat{u}$ is seen to satisfy

[TABLE]

By comparing this equation with (1.6) we find that if $v$ is to correctly simulate the evolution of $u$ then the governing equation (1.4) for $v$ should be modified by replacing $f(t)$ with the function $g(t)$ that has Fourier transform $\widehat{g}(\omega)=\widehat{f}(q(\omega))$, while $u$ and $v$ should satisfy the relation $\widehat{v}(\omega)=\widehat{u}(q(\omega))$, if possible. Since $q$ is periodic, one obstruction is that unless (1.1)–(1.2) leads to a solution $u$ which is bandlimited in time, the relation $\widehat{v}(\omega)=\widehat{u}(q(\omega))$ will only capture part of the frequency content of the sought solution $u$. (Here $u$ is said to be bandlimited in time when $\widehat{u}$ has compact support.) As we shall see, this obstruction can be made negligible also for non-bandlimited $u$ as long as $\widehat{u}$ and $\widehat{f}$ has sufficiently rapid decay at infinity. The transform $f\mapsto g$ which modifies the finite difference equation is called the forward time dispersion transform (FTDT) and the transform $v\mapsto u$ which (ideally) turns a solution of the altered finite difference equation into a solution of the original evolution equation is called the inverse time dispersion transform (ITDT). In practice, there will typically be stability criteria placing an upper bound on the step size $\Delta t$. Throughout the paper we will assume that $\Delta t$ is chosen with such constraints in mind, but as we will demonstrate, one benefit of the proposed method is that it allows $\Delta t$ to be chosen at such upper bounds while still yielding very accurate results. This is true even for low order finite difference schemes, thus providing a new and computationally cheap way to obtain good simulation results.

The method of using dispersion transforms to correct for numerical errors caused by finite difference approximations of time derivatives was introduced in geophysics for the wave equation, see Stork [25] for the original approach. It was then further developed and improved by many authors [28], [8], [14], [15], [1], [16], [21] and seems to now have settled on the definitions of the transforms presented here. For a review of the development we refer to Koene et al. [15]. We mention that in the geophysical literature, usage of the FTDT and ITDT is usually described as applying a pre-computation and a post-computation filter, respectively. Section 2 presents a rigorous mathematical treatment of the method and shows using a unified approach that it has applications to a large class of linear evolution equations. After formally introducing the dispersion transforms (Definition 2.5) we use them to establish the correspondence between an evolution equation and its finite difference equivalent hinted at in the example above (Theorem 2.8), and show that they can be interpreted as Fourier integral operators. In §2.3 we provide discrete versions of the transforms and discuss implementation. Our main result shows that the proposed method of using the discrete dispersion transforms as pre- and post-computational filters yields a numerical solution that correctly models the desired evolution for any length of time (see Theorem 2.10 for the precise statement). Section 3 demonstrates the theoretic results by conducting numerical tests on a model equation where the solution obtained by the proposed method compares to the analytic solution with double precision accuracy (see Figure 2). The approximation error is small as long as the frequency content of the source term is negligible outside a window defined by the range of the phase function $q$ (see Figure 3). The size of this window becomes arbitrarily large as $\Delta t\to 0$.

One shortcoming of the method is the need to store the solution for all moments in time, which for two- or three-dimensional problems may require a very large memory capacity. However, even though simulations might be carried out over a global domain the solution is often only desired at a small subset to which the corrections by means of the dispersion transforms can be restricted, at a cost that is far smaller than running the simulation with a globally large accuracy. We demonstrate such an application in Section 4, where we perform elastic and viscoelastic wave simulation for the earth’s subsurface, but only use the dispersion transforms to correct the resulting ground motion at certain points on the boundary (the source and receiver positions). The viscoelastic simulations show that the method can deal even with dissipative wave physics while still yielding highly accurate solutions.

The paper is concluded with three appendices. In Appendix A we have gathered results of tangential or supplementary nature referenced in the main text. In Appendix B one can find the implementation of the finite difference scheme used in the viscoelastic wave simulations. Finally, in Appendix C we provide codes for implementing the dispersion transforms in MATLAB.

2. Numerical dispersion in evolution equations

Let $X\subset\mathbb{R}^{d}$ and $u=(u_{1},\ldots,u_{K})$ be a vector valued function of $(t,x)\in\mathbb{R}\times X$. Introduce the $K\times K$ system of differential operators

[TABLE]

where $\partial_{t}^{j}u_{k}(t,x)=\partial^{j}u_{k}(t,x)/\partial t^{j}$ and the $L_{ijk}$ are linear spatial operators depending on $x\in X$ but independent of time $t$ so that $\partial_{t}$ and $L_{ijk}$ commute. Consider the evolution equation in $\mathbb{R}\times X$ given by

[TABLE]

Since the system is translation invariant in $t$ we may without loss of generality assume that $t_{0}=0$ below. We will assume that the problem is well posed and that, depending on the spatial operators $L_{ijk}$, appropriate spatial conditions are imposed to ensure a unique solution. For an extensive background on partial differential equations we refer to Hörmander [12] and Evans [5].

When solving (2.1)–(2.2) by means of finite difference methods, numerical dispersion errors inevitably occur as a result of approximating the time derivatives with finite differences. The purpose of this paper is to establish a method by which to alter the chosen finite difference system and capture the correct time evolution of the solution $u$ to (2.1)–(2.2).

In this work, the exact structure of the spatial operators will not be essential. In the applications we have in mind, each $L_{ijk}$ will typically be a differential operator in $x$ with coefficients that depend on $x$ but not on $t$, or (going one step further in the direction of obtaining a numerical solution) each $L_{ijk}$ will be the result of discretizing such a differential operator by means of some numerical scheme. To isolate the effects of time dispersion we will in the latter case assume that any resulting space dispersion errors are essentially fully decoupled from the time dispersion errors, and that appropriate stability criteria that may govern the possible size of the time step $\Delta t$ are satisfied. We will only be concerned with finite difference schemes which are numerically stable and depend continuously on the initial data. A comprehensive treatment of finite difference methods can be found in LeVeque [18].

Example 2.1**.**

As a prototype, consider the scalar, one-dimensional acoustic wave equation in heterogeneous media:

[TABLE]

where $u$ is the acoustic pressure and the propagation velocity $c(x)$ depends on position $x$. Assume Dirichlet boundary values $u(t,0)=u(t,1)=0$. This is an example of (2.1) with $K=1$ and $\mathscr{P}u$ of the form $\mathscr{P}u=\partial_{t}^{2}u+L_{0}u+\partial_{t}L_{1}u$, where $L_{0}=-c(x)\partial_{x}^{2}$ and $L_{1}\equiv 0$. On the other hand, if we want to solve this equation numerically, we may choose to discretize the equation in $x$ by sampling at $x=x_{j}$ where $x_{j}=j\Delta x$ for $j=0,1,\ldots,N+1$, and replace $\partial_{x}^{2}$ with a finite difference, say

[TABLE]

for $1\leq j\leq N$, with the boundary values dictating $u(t,x_{0})=u(t,x_{N+1})=0$. This is also an example of (2.1) with $K=1$ and $\mathscr{P}u$ of the form $\mathscr{P}u=\partial_{t}^{2}u+L_{0}u+\partial_{t}L_{1}u$, where the unknown $u$ is now the $t$-dependent vector $u(t)=(u(t,x_{1}),\ldots,u(t,x_{N}))$ and while we still have $L_{1}\equiv 0$, $L_{0}$ is now a tridiagonal matrix which factorizes as

[TABLE]

Whenever we discuss solving this type of equation by means of a finite difference scheme in time of step size $\Delta t$, we will always assume that the corresponding Courant-Friedrich-Lewy (CFL) condition is satisfied by $\Delta t$ in terms of the interval length $\Delta x$, so that the resulting equation is numerically stable.

Note that because of the initial condition (2.2), the source terms $f_{i}$ in (2.1) will have to vanish identically for $t\leq 0$. We will in addition assume they have sufficient decay as $t\to\infty$ and that each $L_{ijk}$ and $f_{i}$ are regular enough for (2.1)–(2.2) to admit a strong solution $u_{i}$, integrable with respect to $t$, such that for each $1\leq i\leq K$

[TABLE]

Here, $H^{s}=W^{s,2}$ is the usual $L^{2}$ Sobolev space of order $s$. In particular, the partial Fourier transforms $\widehat{u}_{i}(\omega,x)$ and $\widehat{f}_{i}(\omega,x)$ are assumed to be well defined and square integrable, where

[TABLE]

We remark that (2.3) does not automatically hold in general – in the example (1.1)–(1.2) discussed in the introduction it would depend on properties of the operator $A$. In the applications we focus on in this paper, however, the energy is expected to dissipate in the region of interest, due e.g. to damping or the geometry of the spatial domain, which makes (2.3) natural.

Remark*.*

Realizations of the Cauchy problem (2.1)–(2.2) can for example be found in initial value problems for:

• (1)

Ordinary differential equations with constant coefficients.

• (2)

Heat equations, linear parabolic equations.

• (3)

Wave equations, linearly damped wave equations, Maxwell’s equations, linear elasticity.

• (4)

Visco-acoustic and viscoelastic equations solved via memory variables (see §4.1).

• (5)

Strictly hyperbolic pseudodifferential equations, Tricomi equations.

2.1. Finite difference system

Let $\mathsf{D}$ denote the finite difference operator

[TABLE]

where $Z$ is a finite set usually consisting of a subset of the integers or half-integers, and the coefficients $c_{1,n}$ are chosen so that $\mathsf{D}$ becomes an approximation of the first order time derivative $\partial_{t}$. Introduce the finite difference operators

[TABLE]

corresponding to the differential operators $\mathscr{P}_{i}$ discussed above, obtained by approximating the time derivatives by means of $\mathsf{D}$. Note that the spatial operators $L_{ijk}$ thus are the same in $\mathscr{P}_{i}$ and $\mathsf{P}_{i}$. Here and for the majority of the paper $\mathsf{D}^{j}$ denotes the composition

[TABLE]

which means in particular that the same scheme $\mathsf{D}$ is assumed to be used as a basis for $\mathsf{D}^{j}$ in each of the $K$ operators $\mathsf{P}_{i}$. The case of non-matching finite difference schemes is discussed briefly on page Nonmatching finite difference schemes below and again in §A.3 in the appendix.

Taking a partial Fourier transform of (2.4) we observe that

[TABLE]

In view of this identity and the fact that

[TABLE]

we define a phase shift function $q$ as

[TABLE]

so that

[TABLE]

We will assume that $c_{1,n}$ is chosen in such a way that $q(\omega)$ is real-valued and invertible as a mapping $q:\Omega\to q(\Omega)$ for some subset $\Omega=\Omega(\Delta t)\subset\mathbb{R}$. For a comment on the case when $q$ is not real-valued (which happens e.g., in the case of a forward Euler scheme), see the remark on page Backward and forward type schemes. Note also that with respect to the normalized variable $\omega\Delta t$, the right-hand side of (2.7) is invertible for all $\omega\Delta t$ belonging to some fixed, $\Delta t$-independent set. In fact, under the natural assumption that $c_{1,n}\Delta t$ is independent of $\Delta t$, it follows that

[TABLE]

is a trigonometric polynomial independent of $\Delta t$.

Example 2.2**.**

Let $\mathsf{D}$ be given by (2.4), where the index $n$ ranges over the integers, and choose coefficients $c_{1,\pm 1}=\pm 1/(2\Delta t)$ and $c_{1,n}=0$ for all other values of $n$. Then $\mathsf{D}$ is the central difference operator

[TABLE]

and $q(\omega)=\sin(2\pi\omega\Delta t)/2\pi\Delta t$. It follows that $q$ is invertible for $\omega\in\Omega$ where $\Omega=[-\frac{1}{4\Delta t},\frac{1}{4\Delta t}]$. In other words, $q$ is invertible when the normalized variable $\omega\Delta t$ satisfies $\lvert\omega\Delta t\rvert\leq 1/4$. Moreover, $q_{0}(\eta)=\sin(2\pi\eta)/2\pi$.

Example 2.3**.**

Let $\mathsf{D}^{2}$ be the approximation of a second order derivative given by

[TABLE]

compare with the approximation of $\partial_{x}^{2}$ in Example 2.1. Then $\mathsf{D}^{2}$ is equivalent to two applications of the central difference operator from Example 2.2 at half the step size, namely

[TABLE]

It follows that $q(\omega)=\sin(\pi\omega\Delta t)/\pi\Delta t$ for $\omega\in\Omega=[-\frac{1}{2\Delta t},\frac{1}{2\Delta t}]$ and $q_{0}(\eta)=\sin(\pi\eta)/\pi$. This finite difference operator appears in connection with certain leapfrog schemes. We mention that $q$ can also be found by applying a Fourier transform to both sides of (2.9) followed by easy calculations, compare with Remark 2.4 below.

Remark 2.4*.*

As mentioned in the introduction, there is a connection between the phase shift function and the so-called amplification factor appearing in von Neumann stability analysis; in fact, their definitions use the same idea although it is applied to the time domain for the phase shift function whereas it is applied to the spatial domain for the amplification factor. To see this, consider for example the one-dimensional heat equation $\partial_{t}u(t,x)=\partial_{x}^{2}u(t,x)$. Suppose we discretize the equation as

[TABLE]

where the time derivative is approximated using a forward Euler scheme and the second order spatial derivative is approximated using the finite difference scheme in Example 2.3 but now in $x$ instead of $t$. By taking a partial Fourier transform in $x$ of both sides we may write this as

[TABLE]

Rearranging terms and using Euler’s formula gives

[TABLE]

The function in brackets is called the amplification factor at wave number $\xi$ and it depends on the choice of spatial discretization scheme. By the double angle formula we see that for this particular choice it satisfies

[TABLE]

where $q(\xi)=\sin(\pi\xi\Delta x)/\pi\Delta x$ is the phase shift function from Example 2.3, now defined with respect to $\Delta x$ and evaluated at wave number $\xi$ instead of $\Delta t$ and frequency $\omega$. (In von Neumann analysis one usually considers the case of a plane wave $u(t,x)=e^{2\pi ix\xi}$ at fixed time $t$ and assumes that $u(t+\Delta t,x)=g(\xi)u(t,x)$ for some amplification factor $g$ which is found by inserting these expressions in (2.10), see e.g., LeVeque [18, Ch. 9]; the result is the same.)

2.2. Time dispersion transforms

Let $\mathbf{1}_{\Omega}$ denote the characteristic function of a set $\Omega$, so that $\mathbf{1}_{\Omega}(\omega)=1$ if $\omega\in\Omega$ and $\mathbf{1}_{\Omega}(\omega)=0$ if $\omega\notin\Omega$. Based on the previous discussion we will henceforth assume that the function $q$ introduced above is restricted to the largest subset $\Omega=\Omega(\Delta t)$ of its domain of definition containing the origin where the mapping $q:\Omega\to q(\Omega)$ is invertible. The inverse function $q^{-1}$ is assumed to be defined on $q(\Omega)$. We also assume that $\Omega$ and $q(\Omega)$ both exhaust $\mathbb{R}$ in the limit as $\Delta t\to 0$.

Definition 2.5**.**

Let $f_{i}(t,x)$ be a function integrable in $t$. Given a finite difference operator $\mathsf{D}$, let $q$ be the corresponding phase shift function in (2.7). Define the forward time dispersion transform (FTDT) of $f_{i}(t,x)$ as

[TABLE]

Define the inverse time dispersion transform (ITDT) of $f_{i}(t,x)$ by

[TABLE]

The definition extends in the natural way to distributions with well-defined Fourier transforms which are integrable on $q(\Omega)$ and $\Omega$. For example, since the Dirac measure $\delta(t)$ has Fourier transform $\widehat{\delta}(\omega)\equiv 1$ the FTDT of $\delta(t)$ is

[TABLE]

Example 2.6**.**

Let $q(\omega)=\sin(2\pi\omega\Delta t)/2\pi\Delta t$ for $\omega\in[-\frac{1}{4\Delta t},\frac{1}{4\Delta t}]$, so that $q$ is the phase shift function corresponding to the finite difference operator in Example 2.2. Then

[TABLE]

where $\operatorname{sinc}(t)=\sin(t)/t$ is the sinc function.

For future purposes we record the fact that

[TABLE]

which follows by a straightforward change of variable. Similarly, we also have

[TABLE]

Finally, note that

[TABLE]

which together with a straightforward calculation shows that

[TABLE]

In other words, $\mathcal{I}(\mathcal{T}(f_{i}))$ does not equal $f_{i}$, but the bandlimited version of $f_{i}$ with frequency support contained in the range of $q$. Thus, if $f_{i}$ is already bandlimited then $\mathcal{I}(\mathcal{T}(f_{i}))=f_{i}$ for sufficiently small $\Delta t$. As we will demonstrate, the effects of the approximation $\mathcal{I}(\mathcal{T}(f_{i}))\approx f_{i}$ are negligible also for non-bandlimited functions with sufficient decay at infinity as long as $\Delta t$ is chosen sufficiently small (so that $q(\Omega)$ is large enough). This latter situation is analyzed in depth in what follows, in particular in Section 3.

Example 2.7**.**

Consider again the case when $q(\omega)=\sin(2\pi\omega\Delta t)/2\pi\Delta t$ for $\omega\in[-\frac{1}{4\Delta t},\frac{1}{4\Delta t}]$. Let $f$ be a bandlimited function so that $f(t)=\mathcal{F}^{-1}(\mathbf{1}_{[-B,B]}\widehat{f})(t)$ for some minimal number $B$ (the bandwidth). By the Nyquist-Shannon sampling theorem, the sampling rate necessary to accurately represent $f$ is $f_{s}>2B$. However, in order to utilize the entire frequency content of $f$ when computing the forward dispersion transform $\mathcal{T}(f)$, the sampling rate has to be doubled since $[-B,B]\subset\Omega$ if and only if $B<\tfrac{1}{4\Delta t}$, i.e., $f_{s}=1/\Delta t>4B$. Furthermore, the sampling rate has to be effectively tripled in order for $\mathcal{I}(\mathcal{T}(f))$ to equal $f$, since $[-B,B]\subset q(\Omega)$ if and only if

[TABLE]

i.e., $f_{s}=1/\Delta t>2B\pi$. These drawbacks can sometimes be removed, respectively improved, by using a staggered grid provided that the original equation (2.1) permits such leapfrog discretization schemes (compare with Example 2.3). See Section 4 and Appendix B for an example of such an implementation.

We shall now examine the applications of Definition 2.5 for evolution equation modeling. The following theorem establishes a correspondence between the system of differential equations (2.1) and its counterpart obtained by approximating time derivatives with finite differences. In particular, it shows how to use the FTDT to modify the source function in (2.1) when changing to a system of finite difference equations, and how the ITDT turns a function satisfying said finite difference system into an exact solution of the original evolution equation.

Theorem 2.8**.**

Let $u=(u_{1},\ldots,u_{K})$ be a solution to the evolution equation (2.1). Set $g_{i}=\mathcal{T}(f_{i})$ and $v_{i}=\mathcal{T}(u_{i})$. Then $v=(v_{1},\ldots,v_{K})$ satisfies the finite difference system

[TABLE]

for each value of $t$, where $\mathsf{P}_{i}$ is the finite difference operator given by (2.5), obtained from $\mathscr{P}_{i}$ by approximating time derivatives with finite differences. Conversely, suppose that $v=(v_{1},\ldots,v_{K})$ is a function satisfying (2.16) for all $t$, where each $g_{i}$ and $v_{i}$ is integrable in $t\in\mathbb{R}$. Set $f_{i}=\mathcal{I}(g_{i})$ and $u_{i}=\mathcal{I}(v_{i})$. Then $u=(u_{1},\ldots,u_{K})$ satisfies (2.1) for all $t$.

Note that in neither of the two statements of the theorem is the function $v$ obtained by solving the finite difference system (2.16), which comes with considerations of stability when choosing step size $\Delta t$. Obtaining $v$ by numerically solving (2.16) is of course the ultimate goal, but this first requires a discussion about discrete versions of the dispersion transforms and is postponed until subsection 2.3. There we also address the issue of interpolating a discrete solution to make it possible to verify that it satisfies the continuous equation (2.1), see Theorem 2.9.

Proof.

First note that applying the Fourier transform to (2.1) and evaluating at $q(\omega)$ gives

[TABLE]

Next, using the definition of $\mathsf{P}_{i}$ together with (2.6)–(2.7) we get

[TABLE]

by the Fourier inversion formula. Now substitute $\widehat{v}_{i}(\omega,x)=\mathbf{1}_{\Omega}(\omega)\widehat{u}_{i}(q(\omega),x)$ and use (2.17) to obtain $\mathsf{P}_{i}v_{i}(t,x)=\int_{\Omega}\widehat{f}_{i}(q(\omega),x)e^{2\pi it\omega}\,d\omega$. By construction, the right-hand side equals $g_{i}(t,x)$, which proves the first part of the theorem.

To prove the converse statement, suppose that $v$ is some function satisfying (2.16) for all $t\in\mathbb{R}$. Since the $v_{i}$ are integrable we may apply the Fourier transform to both sides of (2.16). Doing so we find by inspecting (2.5) and using (2.6)–(2.7) that

[TABLE]

Setting $u_{i}=\mathcal{I}(v_{i})$ we see by (2.13) that $u_{i}(t,x)=\int_{\Omega}e^{2\pi itq(\omega)}\widehat{v}_{i}(\omega,x)q^{\prime}(\omega)\,d\omega$. Applying $\mathscr{P}_{i}$ to $u=(u_{1},\ldots,u_{K})$ and differentiating under the integral sign shows that $\mathscr{P}_{i}u(t,x)$ is equal to

[TABLE]

By (2.18) we conclude that $\mathscr{P}_{i}u(t,x)=\int_{\Omega}e^{2\pi itq(\omega)}\widehat{g}_{i}(\omega,x)q^{\prime}(\omega)\,d\omega=\mathcal{I}(g_{i})(t,x)$, where the last identity follows by (2.13). This completes the proof. ∎

We conclude this subsection with a few general remarks.

Fourier integral operators*.*

Close inspection of (2.13) and (2.14) using the normalized phase shift function $q_{0}(\eta)$ defined in (2.8) shows that the dispersion transforms $\mathcal{I}$ and $\mathcal{T}$ can be formally interpreted as Fourier integral operators depending on a small semiclassical parameter $h=\Delta t$ (see Appendix A.1). As such, they are associated with a canonical map $\chi$ and its inverse $\chi^{-1}$ acting on phase space via

[TABLE]

The physical meaning of this is well understood in terms of dynamics of wave packets [6]. We provide a detailed presentation in §A.1.1, briefly summarized as follows: let $(t_{0},\eta_{0})$ be a point in phase space and consider a Gaussian wave packet defined by

[TABLE]

When $t\neq t_{0}$ we have $\varphi_{(t_{0},\eta_{0})}(t)=O(\Delta t^{\infty})$ as $\Delta t\to 0$, where $O(\Delta t^{\infty})$ means $O(\Delta t^{N})$ for all $N>0$. Similarly, the semiclassical (i.e., scaled) Fourier transform

[TABLE]

is $O(\Delta t^{\infty})$ as $\Delta t\to 0$ if $\eta\neq\eta_{0}$. Such a function is said to be microlocally small outside $\{(t_{0},\eta_{0})\}$. By Proposition A.2, the ITDT of $\varphi_{(t_{0},\eta_{0})}$ is microlocally small outside

[TABLE]

and the FTDT of $\varphi_{(t_{0},\eta_{0})}$ is microlocally small outside

[TABLE]

Thus in this sense, as $\Delta t\to 0$, the ITDT of $\varphi_{(t_{0},\eta_{0})}$ behaves like the wave packet $\varphi_{\chi(t_{0},\eta_{0})}$ and the FTDT of $\varphi_{(t_{0},\eta_{0})}$ behaves like the wave packet $\varphi_{\chi^{-1}(t_{0},\eta_{0})}$. This phenomenon is illustrated in Figure 1.

We also mention that by using arguments similar to those in the proof of Theorem 2.8, it is straightforward to check that $\mathcal{T}\mathscr{P}_{i}=\mathsf{P}_{i}\mathcal{T}$. Viewing the dispersion transforms as Fourier integral operators, the proof of the first statement in Theorem 2.8 would then proceed by simply noting that, by assumption, $v_{i}=\mathcal{T}u_{i}$ and $g_{i}=\mathcal{T}f_{i}$, so

[TABLE]

The converse statement of Theorem 2.8 can be proved in a similar way. In the sequel we shall continue to prefer elementary proofs using explicit formulas instead of relying on the framework of microlocal analysis. However, this interpretation does succinctly highlight the obstruction caused by allowing time-dependent coefficients in (2.1), see the discussion in Appendix A.1.

Initial conditions*.*

These are not mentioned in Theorem 2.8. Since the finite difference schemes we have considered so far are multistep methods, initial data for $t\leq 0$ is required to get started. The natural choice is to impose the same initial conditions for (2.16) as for (2.1); this is also motivated by the fact that for $t\leq 0$, dispersion should not yet have started to affect the numerical solution. However, suppose that $u_{i}$ is the solution to (2.1) with initial condition (2.2), and let $v_{i}$ be a function satisfying (2.16) with the same initial condition. Then $v(t,x)\equiv 0$ for $t\leq 0$, but due to the non-local nature of the dispersion transforms, this is not true for $\mathcal{T}(u_{i})$ which introduces an approximation error between $\mathcal{T}(u_{i})$ and $v_{i}$. On the other hand, according to Lemma A.3 the error is small and controlled by the time-step size $\Delta t$ (see §A.2 for precise statements). Assuming that solutions of (2.16) depend continuously on the initial data we thus conclude by Theorem 2.8 that $\mathcal{T}(u_{i})$ will continue to stay close to $v_{i}$ for all time. The introduction of this error is also mitigated by the fact that when both dispersion transforms are used together in a modeling scenario as pre- and post-filters, then a reverse error is introduced during the post-filtering process. This is given credence by the numerical results in Sections 3 and 4. In view hereof we will from now on always assume that unless stated to the contrary, initial conditions are given by (2.2).

Backward and forward type schemes*.*

In this case, the phase shift function $q$ will not be real-valued in general. Still, under certain conditions one can define a version of the FTDT and ITDT, although this requires sufficiently fast (exponential) decay of the solutions $t\mapsto u_{i}(t)$ for the definitions above to make sense. This happens, e.g., if $u_{i}(t,x)\equiv 0$ for $t<0$, $\lvert u_{i}(t,x)\rvert\leq Ce^{-2\pi\alpha t}$ for some constants $C$ and $\alpha$, while $\operatorname{Im}q(\omega)<\alpha$ for $\omega\in\Omega$, for then the integral

[TABLE]

defining $\widehat{u}_{i}(q(\omega),x)$ is absolutely convergent for all $\omega\in\Omega$. In particular, $v_{i}(t,x)=\mathcal{T}(u_{i}(\cdot,x))(t)$ is well defined. If the $f_{i}$ satisfy similar decay conditions then the first statement in Theorem 2.8 immediately generalizes to cover this situation. The converse statement can also be generalized using similar adjustments. We leave it to the interested reader to fill in the details.

Nonmatching finite difference schemes*.*

Due to the coupled nature of (2.1) it was essential in the proof of Theorem 2.8 that the same finite difference approximation of the time derivative was used for all involved operators $\mathscr{P}_{i}$. As soon as this is not the case, the result ceases to hold without appropriate modifications. For comparison, one such example of nonmatching finite difference schemes is provided in §A.3.

2.3. Discrete transforms

Theorem 2.8 shows how to use the FTDT to compensate for numerical dispersion when passing from a continuous equation to a finite difference equation, and that the ITDT should be used afterwards to turn a solution of the finite difference equation into a solution of the desired continuous equation. However, Theorem 2.8 failed to address how to apply the ITDT to a solution obtained by actually solving a finite difference equation (modified using the FTDT). For this, we must introduce suitable discrete versions of the transforms. In the process, we will obtain a methodology for correctly simulating the solution to an evolution equation of type (2.1).

We will demonstrate how to simulate the solution for $0\leq t\leq T$, where $T>0$ is the desired lifespan. Discretizing the equations in time and correcting for time dispersion leads us to solve the difference system (2.16). Suppose therefore that $v_{i}$ is a computed solution to (2.16), with known values $v_{i}(t_{n},x)$ at times $t_{n}=n\Delta t$, $0\leq n<T/\Delta t$. (We describe below how to compute the right-hand side of (2.16) using discrete sums.) We assume that $T/\Delta t=N$ for some integer $N$, so that

[TABLE]

and denote by $S$ the set of sampling points

[TABLE]

We will assume that $\Delta t$ is chosen small enough depending on the spatial operators $L_{ijk}$ so that the resulting difference system (2.16) is numerically stable.

We begin with a general discussion and let $f(t)$ be a function of $t\in[0,T]$ with known values at the points in $S$. Let $\omega_{m}\in\Omega$ and introduce

[TABLE]

This is a Riemann sum of the integral $\int_{0}^{T}f(t)e^{-2\pi it\omega_{m}}dt$ and as such an approximation of $\widehat{f}(\omega_{m})$ provided $f$ vanishes outside $[0,T]$. Inspecting the definition (2.13) of $\mathcal{I}(f)$ we then choose a partition of $\Omega$ and define a function of the continuous variable $t\in[0,T]$ via

[TABLE]

Here, $\Delta\omega$ is the distance between two consecutive points $\omega_{m+1}$ and $\omega_{m}$ in the partition. The formula is thus a Riemann sum of the integral $\int_{\Omega}\widehat{f}(\omega)e^{2\pi iq(\omega)t}q^{\prime}(\omega)\,d\omega$. In view of (2.13), this is clearly a discrete representation of the ITDT defined in §2.2. Its usage allows for modeling the desired solution of (2.1) with correct evolution in time.

Theorem 2.9**.**

Let $\Delta t=T/N$ and $v=(v_{1},\ldots,v_{K})$ be a solution of (2.16) computed at times $t_{n}=n\Delta t$ for $0\leq n<N$. Define $u_{i}(t,x)=\operatorname{\mathcal{I}_{\mathrm{disc}}}(v_{i}(\cdot,x))(t)$ and $f_{i}(t,x)=\operatorname{\mathcal{I}_{\mathrm{disc}}}(g_{i}(\cdot,x))(t)$. Then $u=(u_{1},\ldots,u_{K})$ solves (2.1) for $0<t<T$.

Proof.

In the proof we let $x$ be fixed and suppress it from the notation. If $f$ is a function sampled on $S$ and $\mathsf{D}$ is given by (2.4) then a simple calculation shows that

[TABLE]

The second factor on the right is identified as $2\pi iq(\omega_{m})$, with $q$ given by (2.7). We record the fact that if $v_{i}$ solves (2.16) then $a_{m}(\mathsf{P}_{i}v)=a_{m}(g_{i})$, which in view of (2.22) means that

[TABLE]

Next, inserting the definition of

[TABLE]

into (2.1) and differentiating we get

[TABLE]

In view of (2.23) we conclude that

[TABLE]

By definition, the right-hand side is equal to $f_{i}(t)$, which completes the proof. ∎

Having verified that $\operatorname{\mathcal{I}_{\mathrm{disc}}}(v_{i})(t)$ evolves correctly in time, we now discuss how the transform $\operatorname{\mathcal{I}_{\mathrm{disc}}}$ acts on arbitrary vectors in a discrete setting. Given a solution $v_{i}$ to (2.16) with known values $v_{i}(t_{n},x)$ at times $t_{n}=n\Delta t$, $0\leq n<T/\Delta t$, we first construct the function $\operatorname{\mathcal{I}_{\mathrm{disc}}}(v_{i}(\cdot,x))(t)$ as above. To obtain a function sampled on $S$ we simply evaluate $\operatorname{\mathcal{I}_{\mathrm{disc}}}(v_{i}(\cdot,x))(t)$ at the points $t=k\Delta t$, $k=0,\ldots,N-1$. This immediately generalizes to an arbitrary vector of length $N$: given any vector $(f_{0},\ldots,f_{N-1})$ we define its inverse time dispersion transform by

[TABLE]

for $k=0,\ldots,N-1$.

We now describe how to compute the FTDT (of, e.g., the right-hand side of (2.1)) using discrete sums. For any function $f$ sampled on $S$ we define a modified version of the samples $a_{m}(f)$ by

[TABLE]

where the frequencies $\omega_{m}$ are as above. Thus $b_{m}(f)$ is defined by replacing $\omega_{m}$ by $q(\omega_{m})$ in the definition of $a_{m}(f)$. Next, define a function of the continuous variable $t\in[0,T]$ via

[TABLE]

which in view of the previous discussion is a Riemann sum of the integral defining $\mathcal{T}(f)(t)$. To obtain a function sampled on $S$ we evaluate $\operatorname{\mathcal{T}_{\mathrm{disc}}}(f)(t)$ at the points $t=k\Delta t$, $k=0,\ldots,N-1$. Finally, to define the FTDT of a vector we identify $f(n\Delta t)$, $n=0,\ldots,N-1$, with a vector $f=(f_{0},\ldots,f_{N-1})$ and define the forward time dispersion transform of $(f_{n})$ as

[TABLE]

As with the inverse time dispersion transform, this immediately generalizes to an arbitrary vector of length $N$. Given any vector $(f_{0},\ldots,f_{N-1})$ we thus define its forward time dispersion transform by

[TABLE]

Combined with Theorem 2.9, the previous discussion yields the main result of this section.

Theorem 2.10**.**

Let $\Delta t=T/N$. Given a source function $f_{i}^{\mathrm{orig}}$ of (2.1), set $g_{i}=\operatorname{\mathcal{T}_{\mathrm{disc}}}(f_{i}^{\mathrm{orig}})$ and let $v=(v_{1},\ldots,v_{K})$ be a solution of (2.16) computed at times $t_{n}=n\Delta t$ for $0\leq n<N$. Define $u_{i}(t,x)=\operatorname{\mathcal{I}_{\mathrm{disc}}}(v_{i}(\cdot,x))(t)$. Then, for $0<t<T$, $u=(u_{1},\ldots,u_{K})$ solves (2.1) with $f_{i}^{\mathrm{orig}}$ replaced by $\operatorname{\mathcal{I}_{\mathrm{disc}}}(\operatorname{\mathcal{T}_{\mathrm{disc}}}(f_{i}^{\mathrm{orig}}(\cdot,x)))(t)$.

We stress that, as mentioned in §2.2, the composition of the FTDT and ITDT is not the identity mapping since $\mathcal{I}(\mathcal{T}(f))$ is a bandlimited version of $f$ with frequency support contained in $q(\Omega)$. In particular, suppose we want to simulate a solution to (2.1) with source term $f_{i}^{\mathrm{orig}}$. To do so, the method prescribed by Theorem 2.10 is to compute $g_{i}=\operatorname{\mathcal{T}_{\mathrm{disc}}}(f_{i}^{\mathrm{orig}})$, i.e., the (discrete) FTDT of the source term, and solve the discretized equation (2.16) with source term $g_{i}$. If $v_{i}$ is the obtained solution, the theorem implies that $u_{i}=\operatorname{\mathcal{I}_{\mathrm{disc}}}(v_{i})$ simulates the evolution of the solution to the original equation (2.1) but with source term

[TABLE]

which is an approximation of the bandlimited version of $f_{i}^{\mathrm{orig}}$ with frequency support contained in $q(\Omega)$. If the frequency content of $f_{i}^{\mathrm{orig}}$ is negligible outside a compact set then one can make sure to capture its most relevant features by choosing $\Delta t$ sufficiently small. This follows since $q(\Omega)\to\mathbb{R}$ as $\Delta t\to 0$ and is investigated in detail in Section 3 below (see Figure 3).

Remark*.*

Note that a priori, the vector $(f_{n})$ in (2.24) and (2.26) should be a vector representing a function sampled on $S$. If not, interpreting the continuous FTDT and ITDT as Riemann sums lead to different discrete formulas since the range of the index $n$ changes. Note also that although the evolution equation (2.1) is translation invariant, the FTDT and ITDT transforms are not. In particular, we have

[TABLE]

in general. However, this is not a problem since we do in fact have

[TABLE]

as a consequence of (2.15) (in analogy with the Fourier inversion formula). Thus, when solving a Cauchy problem on, say, $[t_{0},T+t_{0}]$, one can still apply (2.24) and (2.26) to a vector representing a function sampled on $[t_{0},T+t_{0}]$. Heuristically, this amounts to the same as translating the original equation to $[0,T]$, applying the transforms there, and translating back. In view of the discussion preceding this remark one is then simulating a solution to an evolution equation with source term

[TABLE]

2.4. Fast implementation

In practice, the formulas (2.24) and (2.26) can often be simplified once a specific choice of phase shift function $q$ is made. Specifically,

• •

using the normalized variable $\omega\Delta t$ can allow the formulas to be interpreted as discrete Fourier transforms which can be implemented using the fast Fourier transform (FFT), and

• •

using symmetry properties of $q$ and $\Omega$ can allow for more efficient algorithms.

Both situations are showcased in the following example.

Example 2.11**.**

Let $\mathsf{D}$ be the central difference operator from Example 2.2,

[TABLE]

Then $q(\omega)=\sin(2\pi\omega\Delta t)/2\pi\Delta t$. Assume as above that $\Delta t=T/N$ where $N$ is the number of sampling points in the time domain, and $T$ is the desired lifespan of the solution. Then $\Omega=\{\omega:\lvert\omega\rvert\leq 1/4\Delta t\}$. To avoid cumbersome notation we will assume that $N$ is even so that $N/2$ is an integer. Inspecting (2.24) we see that we can compute the inner sum by means of the discrete Fourier transform by choosing $\omega_{m}$ appropriately. We pick

[TABLE]

so that $\omega_{m}\in\Omega$ when $m=-N/2,\ldots,N/2$. Substitution into (2.24) gives after cancellations that

[TABLE]

As stated, the inner sum is the value at $m$ of the discrete Fourier transform of $\tilde{f}$, where $\tilde{f}$ is $f=(f_{n})$ zeropadded to twice the length (i.e., the $m$:th Fourier mode of the vector $(f_{0},\ldots,f_{N-1},0,\ldots,0)$ of length $2N$), and can be computed, e.g., using the FFT. The outer sum is the value at $k$ of a modified discrete inverse Fourier transform (truncated to use only the Fourier modes for $-N/2\leq m\leq N/2$ instead of the full range $-N\leq m\leq N-1$). If discrete transforms of numerous samples are to be computed, it is advantageous to interpret (2.27) as a linear map acting on the vector $(f_{n})$ and compute the corresponding matrix. The cost of this operation scales as $O(N^{2}+N\log N)$. Details for implementation in MATLAB can be found in Appendix C.1.

In a similar manner we find by substituting the expression for $\omega_{m}$ into (2.26) that

[TABLE]

Here, the inner sum is a modified discrete Fourier transform while the outer sum is a truncated discrete inverse Fourier transform at $k$. The outer sum can be computed using the inverse FFT. We also observe that if $f=(f_{n})_{n=0}^{N-1}$ is a vector with real entries, then the inner sum in $(\operatorname{\mathcal{T}_{\mathrm{disc}}}(f_{n}))_{k}$ equals $b_{m}(f)/\Delta t$ in the notation above, where $b_{-m}(f)=\overline{b_{m}(f)}$ and bar denotes complex conjugation. This is a consequence of the fact that sine is an odd function. Similarly, $a_{-m}(f)=\overline{a_{m}(f)}$, and these symmetries can be used for a more efficient implementation. See Appendix C for details. Note that both (2.27) and (2.28) only contain frequency content up to a quarter of the sampling rate, i.e., up to half the Nyquist (folding) frequency. This situation is avoided when a leapfrog scheme can be employed, see Appendix C.2.

Remark*.*

An alternative definition of $\operatorname{\mathcal{I}_{\mathrm{disc}}}$ found in the literature [15] is obtained by using Riemann sums to approximate (2.12) instead of (2.13). One such example is

[TABLE]

where the $\xi_{m}$ are points evenly distributed in $q(\Omega)$ and the first factor is the distance between consecutive points $\xi_{m+1}$ and $\xi_{m}$. For implementation using the discrete Fourier transform, a natural option is to choose $\xi_{m}$ so that $e^{2\pi i\xi_{m}k\Delta t}=e^{2\pi imk/2N}$ for those $m$ for which $\xi_{m}\in q(\Omega)$. Then $\xi_{m}=\omega_{m}$ for $m$ in a subset of $[-N,N-1]$, and the formula above reduces to

[TABLE]

(The absence of the factor $q^{\prime}$ found in (2.24) is explained by the relation

[TABLE]

where $\widetilde{\omega}_{m}$ is the preimage of $\xi_{m}\in q(\Omega)$.) It is easy to see that with $q$ as in Example 2.11 this results in

[TABLE]

where $M$ is the largest integer such that $M\leq N/\pi$. Here, the inner sum is a modified discrete Fourier transform while the outer sum can be computed using the inverse FFT.

3. Numerical simulations of a model equation

Here we propose to examine the accuracy of the method by solving a family of ordinary differential equations with known analytic solutions and comparing the resulting numerical solutions, corrected to account for dispersion, with simulations of the analytic expressions. To describe the method’s inherent limitation that comes from restricting the frequency support of the adjusted source term (see the discussion after Theorem 2.10), we shall perform tests with source terms of varying frequency support. We consider the simple model

[TABLE]

where the source $f$ is a modulated Gaussian window function given by

[TABLE]

This is the probability density function of a normal distribution with mean $\mu$ and variance $\sigma^{2}$, modulated by the factor $\exp\left(2\pi ia(t-\mu)\right)$ with modulation parameter $a$ controlling the location of the frequency support of $f$.111In contrast to the wave packets discussed in the remark on page Fourier integral operators and in Appendix A.1, the parameter $a$ is a priori independent of $\Delta t$. In addition, in line with the conventions of probability theory, the factor of normalization has been taken here with respect to the usual $L^{1}$ norm instead of the $L^{2}$ norm. Because of the initial condition in (3.1), the source term should vanish for $t\leq 0$. Moreover, in most applications that we have in mind, the energy is assumed to have dissipated by the end of the experiment. For these reasons we will center $f$ at (say) $t=5$ by taking $\mu=5$, take $\sigma$ so small that $f(t)$ is (practically) zero for $t\leq 0$, and run the experiment well past $t=10$. In particular, if $H$ is the Heaviside function then we will not distinguish between the functions $f(t)$ and $H(t)f(t)$ in what follows. Taking Fourier transforms we see that if $u$ is a solution of (3.1) then

[TABLE]

where we identify the first factor on the right as the Fourier transform of $t\mapsto h(t)=e^{-t}H(t)$. By the Fourier inversion formula it follows that $u$ is given by the convolution

[TABLE]

Applying the methodology presented in Section 2 we shall compare a sample of $u(t)$ for $0\leq t\leq 20$ with a numerically computed solution using the time dispersion transforms. To this end, we consider

[TABLE]

where $\mathsf{D}$ is the central difference operator (appearing in Example 2.11) given by

[TABLE]

and $g=(g_{k})$ is the FTDT of $f$, i.e.,

[TABLE]

compare with (2.28). The sample $g$ is computed using the implementation of the FTDT described in Appendix C.1. After solving (3.4) we finally compute the ITDT of $v=(v_{n})$ using formula (2.27), again implemented as described in Appendix C.1. To minimize potential wraparound effects resulting from using the dispersion transforms (inherited from the FFT and the inverse FFT) on a modulated source function, we will solve the difference equation for $0\leq t\leq 20$ and apply a tapered cosine window, affecting the final sample points when $18\leq t\leq 20$, before computing the ITDT of the result. For transparency, we include plots obtained both with and without this taper.

3.1. Varying the modulation

In Figure 2a we display the analytic solution $u(t)$ computed using (3.3) and sampled at $t=n\Delta t$ with $\Delta t=0.02\,\mathrm{s}$. The source function $f$ was chosen to have mean $\mu=5\,\mathrm{s}$, variance $\sigma^{2}=0.1$ and modulation $a=0$. We furthermore show the numerical approximation of $u(t)$ and its error due to the standard central finite difference scheme and the forward Euler scheme. Finally, we use the time dispersion transform method to compute the solution, and show the difference between $u(t)$ and $\operatorname{\mathcal{I}_{\mathrm{disc}}}(v)(t)$ with and without using a taper. The numerical results are computed on a desktop with Intel Xeon CPU E5-1650 v3 $@$ $3.50\,\mathrm{GHz}$, running MATLAB 2017. It takes 0.083 seconds to compute and apply the FTDT to the source function; 0.00049 seconds for the 1000 time integration steps; and 0.077 seconds to compute and apply the ITDT to the solution vector. The time dispersion transform method outperforms the standard schemes by at least 9 orders of magnitude – when the taper is used we even obtain accuracy up to an order of $10^{-15}$ on the range $t\in[0,18]$.

Figure 2b shows the result of adding a modulation by changing $a=0$ to $a=4$. We see that the Fourier support of the source function $f$ still sits comfortably within the critical frequency set $q(\Omega)$, which for $\Delta t=0.02\,\mathrm{s}$ is given by $q(\Omega)=\{\omega:\lvert\omega\rvert\leq 25/\pi\}$ with $\omega$ measured in $\mathrm{Hz}$. The method continues to perform remarkably well, particularly in comparison with the forward Euler and central finite difference schemes. The computation time is identical to the previous case.

It should be mentioned that the discrete Fourier transform (and its implementation, the FFT) can be viewed as a trapezoidal rule quadrature applied to the Fourier transform integrals. For bandlimited functions, the approximation order is superalgebraic (meaning infinite approximation order) once the sampling is finer than what the bandlimitation prescribes. Implementing the transforms using the FFT and its inverse therefore makes the proposed method especially well suited for the type of essentially bandlimited source functions considered here and helps explain the high quality of these results.

3.2. Varying the frequency support

In Figure 3a we have tried to break the method by setting $a=7.5$. We see that a part of the Fourier support of the source function $f$ is now outside the critical frequency set $q(\Omega)=\{\omega:\lvert\omega\rvert\leq 25/\pi\}$ and the reconstruction of the analytic solution is quite poor. This is in part due to the strong oscillations of $f$; we see in Figure 3a that the Euler scheme also completely breaks down. However, we see that adding a taper results in partial recovery. The computation time is identical to the previous two cases.

As explained in the paragraph following (2.26), we can improve the recovery by decreasing $\Delta t$, thus making sure that $q(\Omega)$ is large enough to capture the most relevant frequency content of $f$. The result of taking $\Delta t=0.01\,\mathrm{s}$ and keeping all other parameters the same can be seen in Figure 3b. Again, our proposed method performs at least 8 orders of magnitude better than the standard schemes. It takes 0.215 seconds to compute and apply the FTDT to the source function; 0.0175 seconds to compute the 2000 time integration steps; and 0.240 seconds to compute and apply the ITDT to the solution.

4. Viscoelastic wave simulation

Here we test the accuracy of the method on seismic wave simulation, with the purpose of supplementing the previous section with a more realistic situation in which an analytic expression for the sought solution is not available. For comparison we provide simulations both of non-dissipative (elastic) as well as dissipative (viscoelastic) waves.

4.1. The viscoelastic equations

A common approach to model wave propagation in anelastic media exhibiting both elastic and viscous behavior is to use viscoelastic theory [24]. For elastic media, it is common to use the analogy of a spring to model the medium under strain (by assuming that Hooke’s law of proportionality between force and displacement holds). This is a linearized assumption that holds well for small displacements such as those found for seismic waves. In tensor notation, each component of the stress tensor will then be a linear combination of all components of the strain tensor. For viscoelastic media, the spring analogy is replaced by the analogy of a spring and dashpot in series, in parallel with another spring. The resulting model is known as a standard linear solid, and several standard linear solids can be connected in parallel to emulate a desired viscoelastic behavior. In a viscoelastic medium, then, each component of the stress tensor will be a linear combination of the entire history of the strain tensor, rather than only its current value. As it is memory-intensive to store the entire strain history, the system of equations is typically reformulated with the use of two constitutive relations, one that expresses the stress as a function of strain and a so-called memory variable, another that expresses the memory variable in terms of the strain and the memory variable itself. Below we recall the resulting equations in two and three dimensions; however our findings also apply to the one-dimensional case. For details of the derivation we refer to [24], [2] and [3].

Wave propagation in a viscoelastic medium with $N$ standard linear solids can be described by Newton’s second law

[TABLE]

together with the stress-strain relation

[TABLE]

with the so-called memory equations:

[TABLE]

In the equations above and throughout this section we employ Einstein notation and sum over repeated indices. The meaning of the symbols is as follows:

• $\rho$

is the density,

• $\sigma_{ij}$

denotes the $ij$:th component of the stress tensor ($i,j=1,2$ in two dimensions and $i,j=1,2,3$ in three dimensions),

• $v_{i}$

denote the components of the particle velocity,

• $f_{i}$

are the components of external body force,

• $r_{ijn}$

are the $N$ memory variables $(n=1,\ldots,N$) corresponding to the stress tensor $\sigma_{ij}$,

• $\tau_{\sigma n}$

is the stress relaxation time of the $n$:th standard linear solid for both pressure and shear waves,

• $\tau^{p}$, $\tau^{s}$

define the level of dissipation for pressure and shear waves, respectively,

• $\pi$

denotes the relaxation modulus corresponding to pressure waves analogous to $\lambda+2\mu$ in the elastic case, where $\lambda$ and $\mu$ are the Lamé parameters,

• $\mu$

is the relaxation modulus corresponding to shear waves and is the analog of the Lamé parameter $\mu$ in the elastic case.

The equations describe the propagation of mechanical waves through a solid, in terms of strains (particle displacements away from their resting position) and stress disturbances (away from the reference stress states), as a function of spatially varying material properties. In geophysics, these equations can be used to model the propagation of seismic waves through realistic dissipative Earth models. The viscoelastic behavior is governed by the $2+N$ parameters $\tau^{p}$, $\tau^{s}$ and $\tau_{\sigma n}$, $1\leq n\leq N$, which all depend on the spatial variable $x$. Here $\tau^{p}$ and $\tau^{s}$ are computed using a variable so-called $Q$ model consisting of one component $Q_{p}$ for the pressure wave and one component $Q_{s}$ for the transverse shear wave via

[TABLE]

$Q$ is a quality factor that approximately measures the amount of energy dissipation [3], with $Q_{p},Q_{s}\to\infty$ corresponding to the elastic, undamped case.

We end by noting that (4.1)–(4.3) constitutes a system of equations of the form (2.1). Indeed, considering the two-dimensional (2D) case, denote $v_{1},v_{2}$ by $u_{1},u_{2}$, denote the three distinct $\sigma_{ij}$, $i,j=1,2$, by $u_{3},u_{4},u_{5}$, and the $3N$ distinct memory variables $r_{ijn}$ by $u_{6},\ldots,u_{3N+5}$. Finally, we note by inspecting the equations above that the spatial operators involved are linear and independent of $t$.

4.2. Model introduction

We apply the theory of the previous section on a viscoelastic wave modeling example, using the leapfrog scheme described in detail in Appendix B to solve (4.1)–(4.3). Since leapfrog schemes are the simplest energy-conserving integrators [9] this is a natural choice in order to avoid numerical dissipation errors and thus isolate the effects caused by numerical dispersion errors.

We use the open-source 2D modeling engine SOFI2D developed by Bohlen et al. [4] to perform the 2D viscoelastic simulation. Time is discretized into steps of constant length $\Delta t$. Similarly, continuous space is discretized into a 2D grid with spacing of $\Delta x$ and $\Delta z$ in the $x$ and $z$ directions. The wave equation is then solved using staggering of quantites in space, and using the leapfrog method to integrate the wave equation in time [27]. The spatial derivatives are efficiently approximated with a central 1D finite difference stencil of half-order 6:

[TABLE]

for which the weights are given in Table 1.

These weights correspond to an equiripple (minimax) filter that keeps the group-velocity error of the first-order derivative approximation confined to within 0.1%. Such ‘optimal’ finite difference coefficients are customary in geophysical finite difference modeling [22], see e.g. [13] for the design procedure. We state the Courant-Friedrichs-Lewy (CFL) condition that ensures stability of the 2D simulation given the second-order accurate integration of the equations in time, as a function of the chosen discretizations and maximum velocity encountered in the simulation:

[TABLE]

We will see that the maximum velocity present in the model is $4700\,\mathrm{m/s}$, and we choose a spatial discretization of $\Delta x=\Delta z=12.5\,\mathrm{m}$. The maximum stable time-step then follows as $\Delta t=1.3\,\mathrm{ms}$. We choose this as the ‘coarse’ time-step. We can compare this coarse solution against an additional ‘fine’ simulation, which uses a time-step of $\Delta t=0.013\,\mathrm{ms}$, which we consider to be the reference solution for our purposes.

The model used for the simulation is the Marmousi 2 model [19], which provides a density model ($\rho$), a model for compressional wave velocities ($v_{p}$) and transverse shear velocities ($v_{s}$) which reflect the instantaneous elastic deformation modes for (4.1)–(4.3):

[TABLE]

The models for $\rho$, $v_{p}$ and $v_{s}$ are shown in Figure 4a. For the viscoelastic modeling we furthermore create the $Q$ model used in (4.4) by smoothing the $v_{s}$ and $v_{p}$ models and normalizing them to a maximum $Q$ of 350, as shown in Figure 4b. Additionally shown in the figures are the source location at $(x,z)=(0,25)$ and a series of recorders along the entire upper model boundary at $(x,z)=(n\cdot 12.5,62.5)$ for $n=-679,\dots,679$, with the coordinates in meters. One specific recorder at $(x,z)=(4625,62.5)$ is highlighted as an arbitrarily shown recorder that will be zoomed in upon in the results.

The source-time function of the model is a typical seismic source wavelet, described as a $15\,\mathrm{Hz}$ peak frequency Ricker wavelet with a time-delay of $0.15\,\mathrm{s}$:

[TABLE]

The source is injected as an explosive source that radiates equally in all directions. The recorders along the upper model boundary record the pressure variations (the diagonal stress components $\sigma_{xx}+\sigma_{zz}$) as a function of time at every $\Delta t$ simulated.

4.3. The elastic model results

We first test the theory in an elastic case, in which we take $Q\to\infty$ so there is no damping, and with $N=0$ so there is no relaxation mechanism at all. The evolution of the wave equation is then computed from time 0 to $5\,\mathrm{s}$ with three model runs:

1. (1)

using a coarse time-step of $\Delta t=1.3\,\mathrm{ms}$ without correcting the source or receiver time-series, 2. (2)

using a coarse time-step of $\Delta t=1.3\,\mathrm{ms}$, but using the FTDT to correct the source injection time-series (4.5), and the ITDT to correct the recorded time-series, 3. (3)

using a fine time-step of $\Delta t=0.013\,\mathrm{ms}$ as a reference solution.

As the implemented wave equation solver scales linearly in time (keeping the spatial discretization $\Delta x=\Delta z=12.5\,\mathrm{m}$ in place), the third simulation thus costs 100 times more computational time. In this elastic instance, it takes 50 seconds to compute simulations (1) and (2), but 75 minutes for the simulation (3). Application in simulation (2) of the FTDT on the source-time function takes half a second to compute, and applying the ITDT to all 1359 recorded signals takes seven seconds in total. This is, essentially, of negligible cost compared to the fine simulation (3).

After finishing all the computations, we subsample the fine simulation to be able to compare the results sample-by-sample. The computed result is then shown in Figure 5. We show a zoom on a single recording (its location is denoted in cyan in Figure 4 but was chosen arbitrarily). It is clearly visible that the coarse simulation (1) created a recording that differs starkly from that made within the fine simulation (3). Conversely, after applying the time dispersion transforms on the source-time function and the receiver time-function, we obtain a solution that follows the correct phase and amplitude of the fine simulation. The two images below the graph in Figure 5 subtract the results of the coarse simulations from the fine simulation – confirming that the correction procedure in this paper removes the dispersion error effectively for all recordings. The sum of all 1359 root mean square (RMS) errors along all traces is 1634 for the coarse case, and 1.6 for the simulation with the proposed time dispersion transforms – the error energy is thus reduced by a factor 1018. The remaining errors seem to be of localized impact only, affecting strong peaks and throughs in the time-series, but do not seem to accumulate over time.

4.4. The viscoelastic model results

The viscoelastic model uses three relaxation mechanisms ($N=3$) to model the spatially heterogeneous $Q$ model. Apart from these changes to the physical model, we proceed in exactly the same way as in the previous example. The computed result is displayed in Figure 6. The amplitudes in this model decrease with time as exemplified by the now 10 times smaller amplitudes in the graph compared to the elastic case, due to the damping. Like the elastic case before, the coarse simulation with $\Delta t=1.3\,\mathrm{ms}$ differs significantly from the fine simulation with $\Delta t=0.013\,\mathrm{ms}$. Conversely, applying the proposed corrections to the coarse simulation creates an adequate fit to the fine simulation. The computational time of all simulations is roughly doubled compared to the elastic simulation at 100 seconds for the coarse simulations and over 2 hours for the viscoelastic simulation. The FTDT and ITDT are still applied to the source-time function in half a second and to all 1359 traces in 7 seconds, respectively. The sum of all 1359 RMS errors along all traces is 212 for the coarse case, and 1.47 for the coarse case using the proposed transforms. Again, the error energy is reduced (now by a factor of 144) at very little additional cost. Again, these errors do not accumulate for longer simulation times.

Acknowledgements

Jens Wittsten was supported by the Swedish Research Council grants 2015-03780 and 2019-04878. Erik Koene was supported by SNF grant 2-77220-15. We gratefully acknowledge the work of Thomas Bohlen, Denise De Nil, Daniel Köhn and Stefan Jetschny on the open-source code SOFI2D used in the viscoelastic simulations. We wish to thank Christof Stork for interesting discussions. We also thank the referee for valuable comments and suggestions.

Appendix A Auxiliary results

A.1. Fourier integral operators

Here we show how the dispersion transforms can be naturally understood as Fourier integral operators (FIO). In this context it will be convenient to view $\Delta t$ as a small (semiclassical) parameter $h>0$. We define the semiclassical Fourier transform of a function $f(t)$ by

[TABLE]

so that $\widehat{f}(\omega)=h^{\frac{1}{2}}\mathcal{F}_{h}(f)(\omega h)$. The Fourier inversion formula then takes the form

[TABLE]

Standard references for semiclassical analysis are Martinez [20] and Zworski [29].

Recall the normalized phase shift function $q_{0}$ introduced in (2.8) which satisfies $q_{0}(\omega h)/h=q(\omega)$, and define $\Omega_{0}$ by $\eta=\omega h\in\Omega_{0}$ if and only if $\omega\in\Omega$. By changing variables in (2.13) it is easy to see that

[TABLE]

Note that $q^{\prime}(\omega)=q_{0}^{\prime}(\omega h)$ and $q^{-1}(\omega)=q_{0}^{-1}(\omega h)/h$, which implies that $q^{\prime}(q^{-1}(\omega))=q_{0}^{\prime}(q_{0}^{-1}(\omega h))$, so making the change of variables $\eta=\omega h$ in (2.11) similarly gives

[TABLE]

If $f$ is a function whose semiclassical Fourier transform has support contained in a set $U$ we will say that $f$ is $h$-bandlimited in $U$. Applying the ITDT to a function already $h$-bandlimited in $\Omega_{0}$ can be naturally understood, in view of (A.1), as the action of a semiclassical FIO (call it $A$) which in $\mathbb{R}\times\Omega_{0}\subset T^{\ast}(\mathbb{R})$ has phase function $\varphi(t,\eta)=tq_{0}(\eta)$ and symbol $a(t,\eta)=q_{0}^{\prime}(\eta)$. $A$ is associated to the canonical transformation locally given by

[TABLE]

i.e., by (2.19). Similarly, $\mathcal{T}$ is a semiclassical FIO (call it $B$) which in $\mathbb{R}\times q_{0}(\Omega_{0})$ has phase function $\psi(t,\eta)=tq_{0}^{-1}(\eta)$ and symbol $b(t,\eta)=1/q_{0}^{\prime}(q_{0}^{-1}(\eta))$. $B$ is associated to the inverse map $\chi^{-1}$. The composition $BA$ acts as the identity operator on functions $h$-bandlimited in $\Omega_{0}$. The composition $AB$ acts as the identity operator on functions $h$-bandlimited in $q_{0}(\Omega_{0})$. Furthermore, using the Fourier inversion formula and arguments similar to those in the proof of Theorem 2.8, it is straightforward to check that $B\mathscr{P}_{i}=\mathsf{P}_{i}B$, so that, as operators acting on functions $h$-bandlimited in $\Omega_{0}$, $\mathsf{P}_{i}=B\mathscr{P}_{i}A$.

Note that the previous discussion can also be had in the framework of microlocal analysis for fixed $\Delta t$, i.e., without viewing $\Delta t$ as a semiclassical parameter. Our choice was made in preparation for §A.1.1 below. If one instead takes the other viewpoint and repeats the arguments above one finds that the dispersion transforms are realized as FIOs associated to the canonical map

[TABLE]

and its inverse. This gives a formula for the appropriate discrete operator to be used for given choice of discrete approximation $\mathsf{D}$ of the time derivative, even in the case when time-dependent coefficients are allowed in (2.1). In fact, if $q$ is the corresponding phase shift function, then the previous paragraph shows that $\mathscr{P}_{i}$ should be replaced by (a discretized version of)

[TABLE]

By Egorov’s theorem, this operator is a pseudodifferential operator with an integral representation

[TABLE]

where, with abuse of notation, the symbol $Q_{i}(t,\omega)$ is a function defined on phase space (omitting all dependence on the spatial variable $x$). Assuming that $\mathscr{P}_{i}=p_{in_{i}}(t)\partial_{t}^{n_{i}}$ plus lower order terms, the principal symbol $\sigma(Q_{i})$ of $Q_{i}$ is given by

[TABLE]

The lower order terms of $Q_{i}$ can be expressed in terms of $\chi_{q}$ and derivatives of the symbol of $\mathscr{P}_{i}$. However, due to the dependence of $\omega$ for example in $p_{in_{i}}(t/q^{\prime}(\omega))$, the two factors on the right cannot be separated in such a way that the operator $\mathcal{T}\mathscr{P}_{i}\mathcal{I}$ is directly realized as a finite difference operator. Investigating the case of time-dependent coefficients is therefore beyond the scope of the current paper and will not be pursued further here.

A.1.1. Dynamics of wave packets

For $(x_{0},\xi_{0})\in T^{\ast}(\mathbb{R})$, a function of the form

[TABLE]

will be called a Gaussian wave packet. Here $\varphi_{(x_{0},\xi_{0})}$ has been normalized with respect to the usual inner product in $L^{2}(\mathbb{R})$. We see that when $h\ll 1$, $\varphi_{(x_{0},\xi_{0})}(t)=O(h^{\infty})$ is negligible if $t\neq x_{0}$, where $O(h^{\infty})$ means $O(h^{N})$ for all $N>0$. Similarly,

[TABLE]

is negligible if $\eta\neq\xi_{0}$. These notions are combined in the following definition.

Definition A.1**.**

Let $u=u(h)$, $0<h\leq h_{0}$, be a family of functions in $L^{2}(\mathbb{R})$. We say that $u$ is microlocally small near $(x_{0},\xi_{0})\in T^{\ast}(\mathbb{R})$ if the inner product

[TABLE]

uniformly for $(x,\xi)$ in a neighborhood of $(x_{0},\xi_{0})$. The complement of such points $(x_{0},\xi_{0})$ is called the semiclassical wavefront set of $u$, denoted $\mathrm{WF}_{h}(u)$.

Another common name for the semiclassical wavefront set is frequency set, usually denoted $\mathrm{FS}(u)$. For other equivalent definitions, including those employing the Fourier-Bros-Iagolnitzer (FBI) transform we refer to Martinez [20] and Zworski [29]. Our presentation is inspired by Faure [6].

As alluded to above, $\mathrm{WF}_{h}(\varphi_{(x_{0},\xi_{0})})=\{(x_{0},\xi_{0})\}$ which is made evident by computing the inner product

[TABLE]

The following result describes how the wavefront set of a Gaussian wave packet is affected by the dispersion transforms.

Proposition A.2**.**

Let $\varphi_{(x_{0},\xi_{0})}$ be a Gaussian wave packet, and let $\chi$ be the canonical map given by (2.19). Then

[TABLE]

Proof.

We will prove the first identity; the proof of the second is similar and is left to the reader. Changing variables in (A.1) we see that

[TABLE]

so an application of the Plancharel formula gives

[TABLE]

In view of (A.3), $(\mathcal{I}(\varphi_{(x_{0},\xi_{0})}),\varphi_{(x,\xi)})_{L^{2}(\mathbb{R})}$ is therefore equal to the integral

[TABLE]

Due to the quadratic terms in the exponential, this is clearly $O(h^{\infty})$ in the semiclassical limit $h\to 0$ if there is no $\eta_{0}$ such that $\eta_{0}-\xi=q_{0}^{-1}(\eta_{0})-\xi_{0}=0$, i.e., such that $\xi=\eta_{0}=q_{0}(\xi_{0})$. Writing the remaining oscillatory factors as $e^{2\pi i\phi(\eta)/h}$ with $\phi(\eta)=x\eta-x_{0}q_{0}^{-1}(\eta)$, it follows from the principle of non-stationary phase that the integral is also $O(h^{\infty})$ unless $\phi^{\prime}(\eta_{0})=0$, see e.g., [29, Lemma 3.10]. But $\phi^{\prime}(\eta_{0})=0$ implies that

[TABLE]

so $\mathrm{WF}_{h}(\mathcal{I}(\varphi_{(x_{0},\xi_{0})}))\subset\{(x_{0}/q_{0}^{\prime}(\xi_{0}),q_{0}(\xi_{0}))\}$. That we in fact have equality follows by applying the method of stationary phase to (A.4) with $(x,\xi)$ determined above. The phase function

[TABLE]

satisfies $\operatorname{Im}\Phi\geq 0$ and has a unique critical point at $\eta_{0}=q_{0}(\xi_{0})$, which is non-degenerate since

[TABLE]

Multiplying $e^{2\pi i\Phi/h}$ by a cutoff function which is identically $1$ near $\eta_{0}$ thus gives

[TABLE]

see Hörmander [12, Theorem 7.7.5]. Since $\lvert e^{2\pi i\Phi(\eta_{0})/h}\rvert=1$, the right-hand side is not smaller than $O(h^{\frac{1}{2}})$ as $h\to 0$, which yields the result. ∎

A.2. Initial conditions

Let $u_{i}$ be the solution to (2.1) with initial conditions (2.2), so $u_{i}(t,x)\equiv 0$ for $t\leq 0$. By virtue of Theorem 2.8, $\mathcal{T}(u_{i})$ satisfies the corresponding finite difference equation modified to account for time dispersion, namely (2.16). Suppose we want to use (2.16) to obtain a sampling of $\mathcal{T}(u_{i})$. This would require adding initial conditions to (2.16), and it is natural to use the same initial conditions (2.2) as before. Let therefore $v_{i}$ be a function satisfying (2.16) with initial conditions (2.2), so that $v_{i}(t,x)\equiv 0$ for $t\leq 0$. Then because the equation but not the initial condition for $v_{i}$ has been modified using the FTDT, this introduces an approximation error between $\mathcal{T}(u_{i})$ and $v_{i}$. Indeed, $v_{i}(t,x)\equiv 0$ for $t\leq 0$ while we have, by the Fourier inversion formula and the definition of $\mathcal{T}(u_{i})$,

[TABLE]

which does not equal

[TABLE]

and thus $\mathcal{T}(u_{i})(t,x)\neq u_{i}(t,x)\equiv 0$ for $t\leq 0$. We see that the “correct” finite difference initial value problem to solve in order to obtain a sampling of $\mathcal{T}(u_{i})$ would be (2.16) with initial condition given by (A.5) for $t\leq 0$. Since this would involve using knowledge of $u_{i}$ this is obviously not possible in practice. On the other hand, since $\mathcal{T}(u_{i})$ and $v_{i}$ have the same evolution in time according to Theorem 2.8, $\mathcal{T}(u_{i})$ will continue to stay close to $v_{i}$ for all $t$ as long as $\mathcal{T}(u_{i})(t,x)$ is approximately identically zero for $t\leq 0$. The accuracy of approximation is ensured by the following lemma, expressed in terms of derivatives $\partial_{t}^{j}\mathcal{T}(u_{i})$ of orders $j$ for which the regularity assumption (2.3) is valid.

To keep the presentation general, we make the assumptions that

• •

$\mathbf{1}_{\Omega}(\omega)$ converges pointwise to 1 as $\Delta t\to 0$, i.e., $\Omega\subset\mathbb{R}$ exhausts $\mathbb{R}$ in the limit as $\Delta t\to 0$,

• •

$q(\omega)$ converges pointwise to $\omega$ as $\Delta t\to 0$,

• •

$\lvert q(\omega)\rvert\geq c\lvert\omega\rvert$ for $\omega\in\Omega$ where $c$ is a real-valued constant independent of $\Delta t$.

We also assume that $\Omega$ has a decomposition $\Omega=\Omega_{\mathrm{inn}}\cup\Omega_{\mathrm{out}}$ consisting of an inner and outer region where $\Omega_{\mathrm{inn}}\to\mathbb{R}$ as $\Delta t\to 0$, such that for some real-valued constants $C_{1},C_{2},C_{3}$ independent of $\Delta t$,

• •

$q^{\prime}(\omega)\geq C_{1}$ if $\omega\in\Omega_{\mathrm{inn}}$,

• •

$\lvert\omega\rvert\geq C_{2}/\Delta t$ if $\omega\in\Omega_{\mathrm{out}}$,

• •

$\Omega_{\mathrm{out}}$ has Lebesgue measure $\lvert\Omega_{\mathrm{out}}\rvert\leq C_{3}/\Delta t$.

To illustrate, if $q(\omega)$ is the function described in Example 2.2 then these assumptions are satisfied with $c=2/\pi$, $\Omega_{\mathrm{inn}}=\{\omega:\lvert\omega\rvert\leq(8\Delta t)^{-1}\}$ and $C_{1}=1/\sqrt{2}$, $C_{2}=1/8$, $C_{3}=1/4$.

Lemma A.3**.**

Let $u_{i}$ solve (2.1)–(2.2). Then

[TABLE]

and the convergence is uniform with respect to $x\in X$ in the sense of (2.3). The rate of convergence depends on the definition of the phase shift function $q$ in (2.7).

Proof.

Inspecting the definitions and recalling that $\Omega=\Omega_{\mathrm{inn}}\cup\Omega_{\mathrm{out}}$ we see in view of the Fourier inversion formula that the result is proved by showing

[TABLE]

with uniform convergence in $t$. Indeed, the left-hand side is the limit of $\partial_{t}^{j}\mathcal{T}(u_{i})(t,x)$ as $\Delta t\to 0$ and the right-hand side is equal to $\partial_{t}^{j}u_{i}(t,x)$ which is zero for $t\leq 0$ by assumption. Note that (A.6) is essentially a consequence of the Lebesgue dominated convergence theorem. For the benefit of the reader we include the details.

Before treating each integral on the left separately we make two observations. First, (2.3) implies

[TABLE]

which means that $\lvert\widehat{u}_{i}(\omega,x)\rvert\leq g_{i}(\omega,x)^{1/2}(1+\lvert\omega\rvert^{2})^{-n_{i}/2}$ for some integrable function $\omega\mapsto g_{i}(\omega,x)$, where $g_{i}(\omega,x)\to 0$ as $\lvert\omega\rvert\to\infty$ by the Riemann-Lebesgue lemma. Second, by assumption we have $\lvert q(\omega)\rvert\geq c\lvert\omega\rvert$ for $\omega\in\Omega$ with $c$ independent of $\Delta t$, so

[TABLE]

We begin by treating the first integral on the left of (A.6). Recall that by our standing assumptions, $u_{i}(t,x)$ is integrable in $t$ which means that $\widehat{u}_{i}(\omega,x)$ is continuous in $\omega$, while $\mathbf{1}_{\Omega_{\mathrm{inn}}}(\omega)\to 1$ and $q(\omega)\to\omega$ pointwise as $\Delta t\to 0$. Next, note that

[TABLE]

since $\xi\in q(\Omega_{\mathrm{inn}})$ implies that $\omega=q^{-1}(\xi)\in\Omega_{\mathrm{inn}}$ for which we have $q^{\prime}(\omega)\geq C_{1}$ by assumption. Since $C_{1}$ is independent of $\Delta t$ and the right-most integral is convergent by (A.7), Lebesgue’s dominated convergence theorem together with (A.8) implies that

[TABLE]

uniformly in $t$.

To treat the second integral on the left of (A.6), note that

[TABLE]

for $0\leq j\leq n_{i}-1$. Since $\lvert\omega\rvert\geq C_{1}/\Delta t$ when $\omega\in\Omega_{\mathrm{out}}$ and $\lvert\Omega_{\mathrm{out}}\rvert\leq C_{2}/\Delta t$, it is then easy to see that

[TABLE]

with $g_{i}$ and $c$ independent of $\Delta t$. Since $g_{i}(\omega,x)\to 0$ as $\lvert\omega\rvert\to\infty$ it follows that the supremum above tends to 0 as $\Delta t\to 0$. In view of (A.8) we conclude that

[TABLE]

uniformly in $t$. This proves (A.6). From the proof it is clear that the convergence is uniform with respect to $x\in X$ in the sense of (2.3), and that the rate of convergence depends on $q(\omega)$. ∎

A.3. Non-matching finite difference schemes

Here we briefly discuss what can happen if non-matching finite difference approximations are used to define $\mathsf{D}^{j}$ in (2.5). To highlight the effect we choose a simple prototype of the evolution equation (2.1): Let $u=(u_{1},u_{2})$ and consider the system

[TABLE]

where the $L_{ij}$ are linear spatial operators independent of $t$ and $C$ is a constant. We will as before let $\mathsf{D}$ denote a scheme satisfying all prior assumptions and use it to model the time derivative in (A.9), but we assume that (A.10) is an auxiliary equation and allow a different scheme to model its time derivative. (An example is provided by the usage of memory variables in (4.2)–(4.3), where in large-scale supercomputer global seismological simulations it can be desirable that the stress and displacement are updated at every step in time, but the memory variables only every 4 steps in time, which saves computational costs without being dramatically worse in performance.) Denote it by

[TABLE]

and define

[TABLE]

The next result explains how the discretized equations should be modified in order to accurately model the evolution of a solution to (A.9)–(A.10).

Proposition A.4**.**

Let $u=(u_{1},u_{2})$ be a solution of (A.9)–(A.10). Set

[TABLE]

and define $v_{i}=\mathcal{T}(u_{i})$ for $i=1,2$. Then $v=(v_{1},v_{2})$ solves

[TABLE]

for each value of $t$, where $\ast$ denotes time convolution.

Proof.

Note that $G$ is well defined since $q$ is real-valued and the integration domain is a compact set. Also,

[TABLE]

Fix $x$ and suppress it from the notation, and let $L_{j}$ be the $1\times 2$ system $L_{j}=(L_{j1},L_{j2})$, $j=1,2$. Using the Fourier inversion formula and the definition of $v_{1}$ we have

[TABLE]

Taking a Fourier transform of (A.9) and evaluating at $q(\omega)$ shows that

[TABLE]

so $v_{1}$ solves (A.11).

To prove that $v_{2}$ solves (A.12), we observe that

[TABLE]

This formula is easily obtained by taking a Fourier transform of (A.10), solving for $\widehat{u}_{2}$ and evaluating the result at $q(\omega)$. It follows that

[TABLE]

Since $\widehat{v}_{2}(\omega)=\mathbf{1}_{\Omega}(\omega)\widehat{u}_{2}(q(\omega))$ we find in view of (A.13) and (A.14) that

[TABLE]

Since $\widehat{u}_{1}(q(\omega))\mathbf{1}_{\Omega}(\omega)=\widehat{v}_{1}(\omega)$, this is equivalent to (A.12) by virtue of the Fourier inversion formula. ∎

Proposition A.4 shows that the price one has to pay for using different finite difference schemes to approximate the time derivatives in (A.9) and in (A.10), is the appearance of a convolution in (A.12). Ignoring the convolution results in an approximation of the desired evolution that can be estimated in terms of the amount by which (A.13) differs from 1.

Note that if the constant $C$ in (A.12) is replaced by a spatial operator $L_{22}$ which, while independent of $t$, is not simply constant in $x$ then the previous result has to be modified accordingly. By minor changes, the proof of Proposition A.4 shows that the result remains valid if $G$ is replaced by

[TABLE]

and (A.12) is replaced by

[TABLE]

We remark that $G$ is well defined due to the assumption that $\lvert q(\omega)\rvert\geq c\lvert\omega\rvert$ for $\omega\in\Omega$.

Appendix B Viscoelastic finite difference equations

Here we discuss the removal of time dispersion from 2D and 3D viscoelastic finite difference modeling for a specific leapfrog scheme developed in [24] (see Bohlen [3] for an explicit implementation). Recall (4.1)–(4.3). The time derivative of a function is approximated by

[TABLE]

In this case, the phase shift function $q(\omega)$ is found to be

[TABLE]

which is invertible for $\omega\in\Omega$ where $\Omega=[-\frac{1}{2\Delta t},\frac{1}{2\Delta t}]$, see Example 2.3. Here the upper limit $1/2\Delta t$ coincides with the Nyquist frequency which is an improvement compared to finite difference scheme employed in Section 3. The drawback is the need to use a time average of the memory variables as described below.

Let $\mathsf{M}f(t,x)$ denote the time average

[TABLE]

of a function $f(t,x)$. Equations (4.1)–(4.3) are discretized in time by

[TABLE]

together with

[TABLE]

and

[TABLE]

If in addition the spatial derivatives are discretized using a fourth-order staggered forward operator and backward operator as is done by Bohlen [3], one arrives at the discrete equations [3, (A.3)–(A.17)] after some straightforward calculations. (In this paper we have extended it to a twelfth-order scheme.) We then have the following.

Theorem B.1**.**

Let $v_{i}$ and $\sigma_{ij}$ solve (4.1)–(4.2), with memory variables solving (4.3). Define $V_{i}=\mathcal{T}(v_{i})$, $\Sigma_{ij}=\mathcal{T}(\sigma_{ij})$ and $g_{i}=\mathcal{T}(f_{i})$. Then for each value of $t$, $V_{i}$ and $\Sigma_{ij}$ solve (B.4) exactly and (B.5) approximately, where the $R_{ijn}$ are exact solutions to (B.6). In the numerical simulations of Section 4, the approximation error is $O(\Delta t^{2})$.

Before the proof we recall from §4.1 that, according to Theorem 2.8, the functions $V_{i}$, $\Sigma_{ij}$ and $R_{ijn}$ are exact solutions to the equations obtained by removing all occurrences of the averaging operator $\mathsf{M}$ from (B.4)–(B.6).

Proof.

We will keep $x$ fixed and omit it from the notation. We prove the statements when $i\neq j$, the other case being similar. We first observe that for a function $f(t)$ we have

[TABLE]

We also record the fact that if $v_{i}$ and $\sigma_{ij}$ solve (4.1)–(4.3) then

[TABLE]

where

[TABLE]

which follows from (4.3) and a straightforward computation.

By definition, $\widehat{V}_{i}(\omega)=\mathbf{1}_{\Omega}(\omega)\widehat{v}_{i}(q(\omega))$. Similar formulas hold for $\Sigma_{ij}$ and $g_{i}$. Hence,

[TABLE]

by the Fourier inversion formula. Using (B.8) evaluated at $q(\omega)$ instead of $\omega$ we see that the right-hand side is equal to $\int_{\Omega}\widehat{f}_{i}(q(\omega))e^{2\pi it\omega}\,d\omega=g_{i}(t)$, which proves that (B.4) holds.

Next, write

[TABLE]

Applying (B.9) evaluated at $q(\omega)$ instead of $\omega$ we get

[TABLE]

We now take a Fourier transform in $t$ of (B.6). Using (B.7), elementary computations show that

[TABLE]

for $\omega\in\Omega$, where the last identity follows by inserting the definition of $V_{i}$ and inspecting (B.10). Using (B.7) again it is straightforward to check that

[TABLE]

where the last factor is uniformly bounded for $\omega\in\Omega$, and the second factor is $O(\Delta t^{2})$ when $\omega$ is restricted to a bounded, $\Delta t$-independent set. In the simulations in Section 4 it turns out that $\widehat{R}_{ijn}(\omega)$ is indeed supported in a $\Delta t$-independent set, see Figure 7. In view of (B.11) we thus conclude that

[TABLE]

The result now follows by applying the Fourier inversion formula to the integral on the right. ∎

Naturally, there are also versions of Theorems 2.9 and 2.10 corresponding to Theorem B.1, as well as a version of the converse statement in Theorem 2.8. We leave for the reader to fill in the details.

Appendix C Implementation

Here we show how to implement the discrete dispersion transforms in MATLAB in two specific cases, namely the finite difference scheme from Example 2.11 that is used in the numerical simulations of Section 3, and the leapfrog scheme from Appendix B that is used in the viscoelastic simulations of Section 4. The interested reader should then be able to adapt the procedure to other cases without difficulty.

C.1. Central difference scheme

Consider the finite difference operator from Example 2.11 and recall (2.27). We see by inspection that we can view $(\operatorname{\mathcal{I}_{\mathrm{disc}}}(f_{n}))_{k}$ in matrix terms as row $k+1$ of a matrix $A$ applied to $\tilde{f}=(f_{0},\ldots,f_{N-1},0,\ldots,0)$, where $A=(a_{(k+1)(n+1)})_{k,n=0}^{2N-1}$ is the matrix with element

[TABLE]

at position $(k+1)(n+1)$. We may view the sum as ranging over $-N\leq m\leq N-1$ with terms for $-N/2+1\leq\lvert m\rvert\leq N$ being zero; in particular the term for $m=-N$ is zero, and so would the term with $m=N$ be. Changing variables $m\mapsto-m$ we thus see that this is the inverse discrete Fourier transform of $m\mapsto g_{k}(m)$ evaluated at $n$, where

[TABLE]

$A$ can therefore be computed by applying the inverse FFT to the matrix with columns $g_{k+1}$ and taking real transpose, e.g., via

where we also take advantage of conjugate symmetry. The last line truncates the matrix so it can be applied directly to the original vector $f$ without having to zeropad the sample manually as this is already built in.

Similarly, by inspecting (2.28) we see that we can view $(\operatorname{\mathcal{T}_{\mathrm{disc}}}(f_{n}))_{k}$ as row $k+1$ of a matrix $B$ applied to $\tilde{f}=(f_{0},\ldots,f_{N-1},0,\ldots,0)$, where $B=(b_{(k+1)(n+1)})_{k,n=0}^{2N-1}$ is the matrix with element

[TABLE]

at position $(k+1)(n+1)$. As before we may view the sum as ranging over $-N\leq m\leq N-1$ with terms for $-N/2+1\leq\lvert m\rvert\leq N$ being zero. We thus see that this is the inverse discrete Fourier transform of $m\mapsto h_{n}(m)$ evaluated at $k$, where

[TABLE]

$B$ can therefore be computed (without taking transpose) by applying the inverse FFT to the matrix with columns $h_{n+1}$, e.g., via

As before, the last line truncates the matrix thus avoiding the need to zeropad the sample $f$ manually. The matrices $A$ and $B$ are depicted in Figure 8.

C.2. Leapfrog scheme

Consider now the leapfrog scheme from Appendix B and recall that in this case, the phase shift function is $q(\omega)=\sin(\pi\omega\Delta t)/\pi\Delta t$ by (B.2), so $q$ is invertible for $\omega\in\Omega$ where $\Omega=[-\frac{1}{2\Delta t},\frac{1}{2\Delta t}]$. We remark that this is not the same as replacing $\Delta t$ by $\Delta t/2$ in the previous subsection since the time-step size for each individual function is in fact $\Delta t$. Repeating the arguments in Example 2.11 for this choice of $q$ we find by inserting

[TABLE]

into (2.24) that

[TABLE]

Here the inner sum is the value at $m$ of the discrete Fourier transform of the vector $(f_{0},\ldots,f_{N-1})$ zeropadded to twice the length. The outer sum is the value at $k$ of a modified discrete inverse Fourier transform (without truncation). This is the image of $f$ under the action of row $k+1$ of the matrix $A=(a_{(k+1)(n+1)})_{k,n=0}^{2N-1}$ with element

[TABLE]

at position $(k+1)(n+1)$. Changing variables $m\mapsto-m$ we identify this as the inverse discrete Fourier transform of $m\mapsto g_{k}(m)$ evaluated at $n$, where

[TABLE]

$A$ can be computed by applying the inverse FFT to the matrix with columns $g_{k+1}$ and taking real transpose via

Next, by inserting the expression for $\omega_{m}$ (shifted one index) into (2.26) we find that

[TABLE]

Thus we can view $(\operatorname{\mathcal{T}_{\mathrm{disc}}}(f_{n}))_{k}$ as row $k+1$ of a matrix $B$ applied to $(f_{0},\ldots,f_{N-1})$ zeropadded to twice the length, where $B=(b_{(k+1)(n+1)})_{k,n=0}^{2N-1}$ has element

[TABLE]

at position $(k+1)(n+1)$. This is the inverse discrete Fourier transform of $m\mapsto h_{n}(m)$ evaluated at $k$, where

[TABLE]

$B$ can be computed by applying the inverse FFT to the matrix with columns $h_{n+1}$ via