Provably Stable Full-Spectrum Dispersion Relation Preserving Schemes

by   Christopher Williams, et al.

The dispersion error is often the dominant error for computed solutions of wave propagation problems with high-frequency components. In this paper, we define and give explicit examples of α-dispersion-relation-preserving schemes. These are dual-pair finite-difference schemes for systems of hyperbolic partial differential equations which preserve the dispersion-relation of the continuous problem uniformly to an α%-error tolerance. We give a general framework to design provably stable finite difference operators that preserve the dispersion relation for hyperbolic systems such as the elastic wave equation. The operators we derive here can resolve the highest frequency (π-mode) present on any equidistant grid at a tolerance of 5% error. This significantly improves on the current standard that have a tolerance of 100 % error.



page 1

page 2

page 3

page 4


Upwind Summation By Parts Finite Difference Methods for Large Scale Elastic Wave Simulations In Complex Geometries

High-order accurate summation-by-parts (SBP) finite difference (FD) meth...

Sharp stability for finite difference approximations of hyperbolic equations with boundary conditions

In this article, we consider a class of finite rank perturbations of Toe...

Spatially-optimized finite-difference schemes for numerical dispersion suppression in seismic applications

Propagation characteristics of a wave are defined by the dispersion rela...

Free-stream preserving finite difference schemes for ideal magnetohydrodynamics on curvilinear meshes

In this paper, a high order free-stream preserving finite difference wei...

Approximating moving point sources in hyperbolic partial differential equations

We consider point sources in hyperbolic equations discretized by finite ...

Computing semigroups with error control

We develop an algorithm that computes strongly continuous semigroups on ...

Annihilation Operators for Exponential Spaces in Subdivision

We investigate properties of differential and difference operators annih...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

It is well known that the dispersion error of a computational scheme has a great impact on the accuracy of a numerical simulation for hyperbolic partial differential equations (PDEs). Various efforts have been made to minimise the dispersion error for computational schemes [tam1993dispersion, linders2015uniformly], but all of these efforts have been unable to resolve high-frequency wave modes, such as the highest-frequency mode on a computational grid called the -mode. These attempts have also exclusively used central finite difference (FD) stencils on co-located grids, as this has generally been accepted as necessary for a stable numerical scheme. Recently, the dual-pairing summation-by-parts (SBP) FD framework [Mattsson2017, DovgilovichSofronov2015] has shown that this is not necessarily true. In some applications, staggered FD methods [Yee1968, Virieux1986] on Cartesian meshes can minimise numerical dispersion errors. However, the design of high order accurate and stable staggered FD methods in complex geometries is a challenge. See some progress in this direction [OssianAnders2020, Ossian2021].

In this paper, we exploit the dual-pairing SBP FD framework on co-located grids to design FD schemes which can resolve high frequencies to an acceptable error tolerance for use in application. The dual-pairing SBP FD framework introduces additional degrees of freedom that can be tuned to diminish numerical dispersion errors. In doing so we minimise two major numerical issues currently present in computational wave simulation, these are: the presence of spurious high-frequency wave modes; and, the reduced grid-refinement needed for the same error bounds by a factor of more than two in one-spatial dimension. The result of this improvement is the absence of computationally fatal spurious wave modes in numerically computed solutions, and an efficiency increase that is exponential with the dimension of the problem. For instance, in three-dimensions for a time-dependent problem our schemes achieve the same error bounds as traditional methods with

times less computational effort.

The structure of this paper is as follows. First we recall the traditional SBP [kreiss1974finite, strand1994summation] and dual-pairing SBP [Mattsson2017, DovgilovichSofronov2015] framework. Here we identify what is necessary for our scheme to be provably stable. Next we define what an -dispersion-relation-preserving scheme is within this framework. We then turn to the construction of a -DRP scheme, giving a methodology for constructing upwind dual-pairing schemes that are able to resolve high-frequencies to error. We find that this problem can be reduced to two optimisation problems. The first optimisation pertains to the design of an interior stencil that is a non-convex smooth optimisation problem with a convex relaxation. The second relates to the boundary stencils needed, and can be posed as a non-smooth convex minimisation problem. We then compute explicit examples and present the numerical error properties of these schemes. Small test examples are presented and compared to standard methods.

2 The Dual-Pairing SBP FD Framework

Previous work on FD methods has primarily focused on central stencil schemes, with carefully chosen boundary stencils so that the SBP property holds. It was thought that central FD stencils was necessary for provably stable schemes for general systems. Recently, progression on skewed stencil schemes in the SBP framework has shown that non-central stencils may also be stable. In this section we review what is necessary for a skewed stencil scheme to maintain stability.

2.1 Notations

The following definitions are required for our following study. We will always use the variable to be the spatial variable and work with the closed interval of length . Let with the standard inner product


Let be the number of grid points on a uniformly spaced grid,


Here are the grid-points and is the spatial step. For a given

we define the grid vector



For a function we use the bold character to denote the grid function, that is the restriction of on the grid,


We introduce the diagonal and positive definite matrix


and define the discrete inner product


We note that with the non-dimensional constants being the weights of a composite quadrature rule and is the spatial step. The following vectors are of repeated interest:


We will also use the sequence spaces , those are the space of square summable sequences and doubly infinite bounded sequences. That is


2.2 Measure Approximation

Traditionally it has been the case that derivative operators have been approximated by FD quotients as an approximation of the limit given by a classical derivative. Finite element methods were then formulated to solve PDE’s in their weak form, and provided a strong theoretical foundation for numerical methods. The SBP FD framework [kreiss1974finite, strand1994summation] was then formulated as to mimic the theoretical guarantees of spectral methods for finite difference schemes. It was found that SBP methods also produced quadrature rules for numerical integration. Here, we generalise this idea by starting with an approximating discrete measure, and extrapolating differential operators in a flexible framework which allows for other critical components of the PDE to be approximated.

We will always use an atomic measure to approximate Lebesgue measure. For simplicity we will restrict our attention to the interval and denote the Borel measurable subsets to be .

Definition 2.1.

An atomic measure is a measure of the form


where are quadrature weights and are the quadrature points.

A family of atomic measures converges weakly to Lebesgue measure if for all continuous functions on we have


See that


In particular, we may form the inner-product


So we may view the functions as discrete evaluations, or the measures as discrete approximations to the continuous setting. We make this clear as to inherit the theoretical framework from the finite element setting for finite difference operators. To gain this in our dual-pairing SBP framework we demand two conditions on our approximating measure , these are:

  1. in the weak sense; and,

  2. there exists linear operators so that


On the right hand side, is a boundary integral operator given through . This operator is identically zero on the interior of our domain granting the negative duality


when restricted to the interior of the domain. There are three primary advantages to using the dual-pairing SBP framework in the finite difference setting, these are:

  1. numerical stability is preserved from the continuous setting;

  2. weak derivatives are approximated instead of strong derivatives in traditional finite difference schemes; and,

  3. the stencils can be designed to mimic the dispersion relation in the continuous setting.

The first and second properties come from the general dual-pairing framework, the last is the primary topic of this paper for the design of such schemes. We briefly show the second property now.

Proposition 2.2.

Let and , then as we have that .


From the assumptions above we have


The right hand side converges to zero as tends to infinity. This is because


strongly as is smooth and is of at least first order. Write where minimally , then the convergence is given from the weak convergence of as is necessarily absolutely continuous. ∎

Thus even when the test function space comprises of smooth functions, the finite-difference operators may yield legitimate approximations to weak derivatives. In summary, the dual-pairing framework allows for weak-derivative approximation and the preservation of the summation-by-parts.

2.3 Discrete Dual-Pairing SBP Operators

The matrix encodes the 1D boundary integral of the interval . Again, throughout our study, the matrix is always a diagonal positive definite matrix. We now may state the integration-by-parts and the summation-by-parts properties.

For we have the integration-by-parts property


Let so that approximates the first derivative grid function . The summation-by-parts property can be stated as


for suitable grid functions . We may now separate the criteria of differentiation on a test function and the SBP criteria through considering the matrix equation


This is the standard SBP framework [kreiss1974finite, strand1994summation]. We will analyse the extension where we have two differential operators , that obey on test function vectors


and the modified SBP property


that translates to the matrix equation


As is standard notation, define the matrices and . To simplify finding these operators we will put the additional constraint


Note that


We will assume that the stencil width for is for . When the matrix is diagonal, it is necessarily true that for that is identical to the identity operator for with away from one or . For the points such that we call the position an interior point of the stencil. The compliment of this set is the boundary points of the operator. Throughout this study, we will use polynomials to be the test functions considered. We may now formally define a dual-pairing SBP framework.

Definition 2.3.

A dual-pairing SBP operator of accuracy is the 3-tuple where is a diagonal and positive matrix, and


for polynomials of degree or less on the interior points and polynomials of degree on the boundary and the involved matrices obey


The high frequency modes (in the neighborhood of the -modes) for both traditional stencils and the stencils we consider carries the majority of the numerical error. It is useful to consider some form of artificial dissipation for high-frequency components. For such scenarios we consider upwind dual-pairing SBP operators.

Definition 2.4.

A upwind dual-pairing SBP operator is a dual-pairing SBP operator that additionally satisfies


is negative semi-definite.

Before we can define a dispersion-relation-preserving-scheme we must recall some basic properties of the dispersion relation for linear hyperbolic problems.

2.4 -Dispersion Relation Preserving Schemes

In this section we recall what the dispersion relation is for a simple hyperbolic system. We then define -dispersion-relation-preserving schemes (-DRP schemes).

Consider the simplest hyperbolic system


Scaling constants and tensor products may be used to generalise our method to higher dimensional linear hyperbolic systems, but here we will only consider this simple case. Using the notation

, take the Fourier solution basis of the form


Here is the angular frequency, is the spatial wave number and is the constant polarisation vector. For the continuous model problem (28), the wave-mode belongs to the solution if and only if


This is the dispersion relation for this model problem (28).

Remark 2.5.

We can take square roots in (28) yielding the linear dispersion relation


For the current study, we are only interested in analysing the dispersion relation for the interior of a given stencil.

We consider the -periodic functions in interval and the uniform discretisation (2) of the interval.

Note that due to periodicity of the grid functions we must have

Definition 2.6.

Let be a dual-pairing SBP operator, and let be the interior stencil for . Let be the space of doubly infinite bounded sequences of real numbers, then . These operators are the interior upwind operators.

Recall for interior points , where is the -th unit vector, and so for the interior operators we only need


for the summation-by-parts criteria.

Now let us consider the semi-discrete counterpart


where we have replaced the spatial derivatives with the finite difference operators .

We consider the interior stencils only, thus

with the consistency requirements

Here, is the uniform grid spacing, are the non-dimensional constant coefficients defining the upwind finite difference stencils. Note that for .

Let be the numerical dispersion relation for a given discrete operator for the semi-discrete approximation (34). For a given dual-paring SBP scheme we may find an explicit expression for the discrete dispersion relation.

Let be an internal upwind pairing, then introduce the numerical wave number and define the discrete Fourier symbols


such that



Lemma 2.7.

Consider the semi-discrete approximation (34) with . The semi-discrete discrete dispersion relation for (34) is given by


where is numerical frequency and is the numerical wave number.


Inserting the plane waves

in the semi-discretised PDE (34), we have


and the solvability condition


Re-arrangement of the terms to find and using , as defined in (35) gives the result. ∎

We define the numerical dispersion relation


Similarly, for the traditional SBP operator we have the numerical dispersion relation


where are the non-dimensional constant coefficients of the finite difference operator, with , , and satisfying the consistency requirements

By comparing the numerical dispersion relation , defined in (40)–(41) and the continuous dispersion relation , defined in (30), we can determine the numerical dispersion error. Define the point-wise dispersion error


We are interested in schemes that are accurate across the whole spectrum , so we define the maximal relative dispersion error to be


Note that for the traditional SBP operators based on central finite difference scheme, from (41) we have . Therefore, all central stencil schemes have a maximal dispersion error of , and are completely unable to resolve the highest frequency modes on a given computational grid [tam1993dispersion, linders2015uniformly].

We are interested in constructing schemes whose dispersion error for high frequencies is acceptable for most application problems. This motivates our definition for a -DRP scheme.

Definition 2.8.

Let be an interior dual-pairing SBP operator. We call an -DRP operator if the maximal relative dispersion error is bounded by . When a dual-pairing SBP operator has a -DRP operator on its interior, we call it a -DRP dual pairing SBP operator.

For many applications, the tolerance is an acceptable error. Colloquially, we will call a scheme DRP-preserving if it is a -DRP scheme with . In the next section we show how we can construct such a scheme.

3 Construction of -Dispersion Relation Preserving Schemes

In this section we show how to construct DRP-schemes which are provably stable through the dual-pairing SBP framework. Previously work [tam1993dispersion, linders2015uniformly] has looked at optimising the dispersion relation for central stencil schemes. All such attempts are -DRP schemes that support spurious wave modes in their solution, being unable to resolve wave numbers higher than . Due to the dual-pairing framework, we may resolve the full spectrum to error.

3.1 Interior Stencil

First we focus on designing the interior stencil of our scheme. Recall we are interested in the minimisation


for a discrete dispersion relation from a family of candidate models and all . Ideally, we would consider the infinity norm of this quantity but for computational feasibility we will consider the least-squares minimisation on the whole spectrum. To the authors knowledge, this is the first successful attempt to optimise the dispersion error for the whole spectrum. Note that continuous angular frequency and the numerical counterpart take the sign of its argument . We consider the following least squares optimisation


for an appropriate family of candidate models . It is important that our approximant is exact around , as this determines the error convergence with grid refinement.

As stated above, for convenience we will work with quadratic positive functions , , however, using the following Lemma the analysis holds also for , .

Lemma 3.1.

Let and , that take the sign of its argument , then we have


First by the reverse triangle inequality and non-negativity,


Now we have


So we gain


but as we are compactly supported:


If we now consider to be a convex set, then both




are strictly convex and share the same minimiser if it is in the candidate model set due to Lemma 3.1.

Now the function to be approximated, , is smooth around zero. To make our optimisation computationally feasible, we will place both and our elements of into orthonormal cosine bases. That is


where are the Fourier coefficients. This is a natural choice for the function we are aiming to approximate. The coefficients for the continuous dispersion relation have a closed form expression, namely

The coefficients for the numerical dispersion relation depend on the finite difference stencils and can be computed using a symbolic mathematical software, such as Maple. The coefficients form a finite square-summable sequence, and vanishes after a finite , that is

Lemma 3.2.

For positive functions ,


where the square summable space of sequences whose entries are zero after entries, are such that


with the projection of into and


Applying Lemma 3.2 and Parseval’s formula yields


Now by Bessel’s equality for the level truncation of the sequence and its compliment we have


For the entry in decays like as the extended periodic function on is once weakly differentiable, and so the square sum of these terms behaves like completing the claim. ∎

Last, it would be beneficial if we could find a parameterised family that forms a linearly independent basis set. For instance if is of the form


for , then the functions form a linearly independent basis for the family


Further if and as , . This gives a approximation space for , but it is still of great importance that the approximant is uniformly a good approximant to around zero. For our stencils, we fix so for each member of the generated uniformly approximates around , then use the remaining free-parameters to approximate the function in the higher-frequency region.

All that is left is to find a family which obeys our given assumptions. Let be the central difference stencil of order , then following [Mattsson2017] we may construct


where is the relevant scaling constant with being the standard first order backward and forward FD operators, see [gustafsson1995time], defined by

This produces a finite-difference approximant of order , taking the negative transpose of gives an order approximant that obeys


The co-efficients for these operators are given in Table 5.



This family is naturally parameterised by the parameters . The following is a result found from symbolic computational verification.

Lemma 3.3.

Let for any for , then is represented exactly in the form


for some finite .


It is sufficient to check


is of the correct form for all as is made from weighted linear combinations of these. For there are checks to complete, which is easily computed in a symbolic manipulation software. ∎

Putting all of our work together we can prove the following bound.

Theorem 3.4.

Let be positive functions whose squares are once weakly differentiable and have Fourier co-efficients , then


for the level truncation of .

Currently, we have the strong assumption that our actual constraint set has the form , for a convex set . This is not in general true, but we can note that


where we have the constraint set . An appropriate convex relaxation of this problem is to consider the cost of the form


for in such that and . Such a relaxation still has uniqueness and existence of minimiser, but as soon as we introduce the constraint of equality on we lose uniqueness.

Example 3.5.

Consider the case

with the dummy variables

which obey the relations


We will consider the three dimensional space spanned by . Recall granting that the constraints given correspond to a parametric parabloid, meanwhile the constraint


is equivalently the plane given through the implicit equation


The intersection of these curves corresponds to a one-dimensional subspace given through the parametric equation


for . Now consider the cost function


This is a convex function in the parameters , but our constrained optimisation yields two global minima. Consider the sublevel sets


where it is clear is closed and convex for all as is convex. We may reform our optimisation criteria as


This occurs at showing the minimum is one, and there are two points of intersection corresponding to being either or , granting two minima. This occurs as the constraint set is not convex. For the case, by Bézout’s bound we are guaranteed at most two solutions as there are at most four points of intersection of these quadratic curves, but when counted with multiplicity there are at most two intersection points due to the minimality condition.

Finally, we can state the (finite-dimensional) optimisation to be performed.


subject to


The cost function is now convex and coercive in the variables , granting existence of a minimiser as it positive, but due to the non-convex constraints, we may lose uniqueness. Using standard least-squares optimisation techniques, we may find global minima. The results for operators of orders are presented in the Appendix. These are presented in rational arithmetic.

By a simple bootstrap argument, we may also consider the minimisation


subject to


for a weighted norm that is induced by the inner product


for . This can be done through noting that for we have


which can be formulated with an identical procedure. Of particular interest is the weight function


for some maximum desired frequency . Such an optimisation would yield operators which are more accurate for lower frequencies at the expense of high-frequency accuracy. We can easily target which part of the spectrum we wish to have the most accuracy with this bootstrap. This is not the focus for the current study, and so our -DRP operators will be made with the weight , but we will give a small example.

Example 3.6.

Consider the weight function which weights the importance of the high-frequencies higher than that of the low-frequencies. We will use the approximation space so in this case the value is small and the accuracy of this operator for low frequencies would not be suitable for application. The objective for our optimisation has the form


where are appropriately weighted by the the Fourier cosine representation of . Let be the found dispersion relation, then we can compute the point-wise relative error to be . We plot and compare this relative error with an identical optimisation where the weight was used.

Figure 1: Dispersion errors on with the maximal dispersion error of the high-frequency scheme plotted in red.

As we can see, our weight function penalises the higher-end of the spectrum more and yields a smaller error for the -mode at the expense of worse error bounds for lower frequencies.

Remark 3.7.

A heuristic observed for designing interior operators is that higher accuracy can be found with increasing

and for different parts of the spectrum. We observe:

  1. higher values correspond to higher-order test functions and generally lead to operators with better accuracy in the low-frequency region ;

  2. the value can be chosen sufficiently larger than to allow the free variables in the linearly independent cosine approximating functions to minimise the error in the high-frequency region.

In Example 3.6, there is a larger error in the low-frequency region than the high-frequency region as is small but .

We may also use the weight to be more accurate at the high-frequency end of the spectrum, at the expense of accuracy for low frequencies. This sacrifice could be compensated for by setting the value sufficiently large.

Last, we consider the simplest dispersion relation here, but this can be generalised to more complicated relations in the same framework. For instance, in fluid dynamics dispersion relations of the form are common. In Figure 2 is an approximant we have made for this case. This is a second order operator, common in fluid modelling, with three additional stencil points. The maximal dispersion error is no greater than .

Figure 2: A second order dispersion relation preserving scheme for . This function is difficult to approximate due to the poor Fourier representations of but can be done through Stejiles representations.

3.2 Boundary Stencil

In this section we show how to close the boundary stencil for a given interior stencil so that we keep the dual-pairing SBP framework. As we have used upwind stencils for our interior scheme, we will close the boundaries so that the final dual-pairing SBP scheme is also an upwind scheme. We assume the ansatz


where . For example, in the case