The EPS method: A new method for constructing pseudospectral derivative operators

Abstract

We develop a new type of derivative matrix for pseudospectral methods. The norm of these matrices grows at the optimal rate O(N²) for N-by-N matrices, in contrast to standard pseudospectral constructions that result in O(N⁴) growth of the norm. The smaller norm has a big advantage when using the derivative matrix for solving time dependent problems such as wave propagation. The construction is based on representing the derivative operator as an integral kernel, and does not rely on the interpolating polynomials. In particular, we construct second derivative matrices that incorporate Dirichlet or Neumann boundary conditions on an interval and on the disk, but the method can be used to construct a wide variety of commonly used operators for solving PDEs and integral equations. The construction can be used with any quadrature, including traditional Gauss–Legendre quadratures, but we have found that by using quadratures based on prolate spheroidal wave functions, we can achieve a near optimal sampling rate close to two points per wavelength, even for non-periodic problems. We provide numerical results for the new construction and demonstrate that the construction achieves similar or better accuracy than traditional pseudospectral derivative matrices, while resulting in a norm that is orders of magnitude smaller than the standard construction. To demonstrate the advantage of the new construction, we apply the method for solving the wave equation in constant and discontinuous media and for solving PDEs on the unit disk. We also present two compression algorithms for applying the derivative matrices in O(N log N) operations.

Introduction

Pseudospectral methods have proved useful in modeling a wide variety of phenomena including wave propagation and fluid dynamics. (See e.g. [12], [6], [13] for an overview.) By differentiating a fixed set of basis functions analytically, pseudospectral methods can provide superior accuracy compared to finite difference methods assuming that the basis is well-suited for the modeled phenomena. For example, using periodic trigonometric basis functions gives optimal performance for periodic problems, but are less suitable for non-periodic problems. For non-periodic problems a polynomial basis, such as Chebyshev polynomials, is often used. Recently, prolate spheroidal wave functions (PSWF) have also been proposed as a basis for pseudospectral methods [26], [7], [5]).

When using pseudospectral methods for solving PDEs, it is critical to construct accurate and efficient numerical representations of derivative operators with respect to the basis of the pseudospectral method. A major problem with such pseudospectral derivative matrices is that the norms of these matrices may become large, effectively restricting the time step used when modeling time dependent phenomena [24]. Controlling the norm of pseudospectral derivative matrices, without adversely affecting accuracy, is therefore of great importance when designing pseudospectral based algorithms.

Our goal in this paper is to construct accurate and efficient derivative matrices for functions restricted to a finite domain while controlling the norm of the derivative matrix. For the second derivative matrix the norm of our derivative matrix D² will scale as ∥D²∥ ∼ O(N²) where N denotes the number of nodes. The size of the norm, O(N²), can be considered optimal in the sense that the eigenvalues for the analytic second derivative operator with zero Dirichlet or Neumann conditions on the interval [−1, 1] are given by $-\frac{m^2{\pi }^2}{4}\text{,}m=1\text{,}2\text{,}\dots$. Since the analytic second derivative operator is symmetric, the norm is given by the spectral radius of the operator, which therefore must grow as O(N²).

In standard pseudospectral methods derivative matrices are typically constructed by analytically differentiating the interpolating basis functions of the method (see e.g. [12]). As an example, consider the second derivative matrix based on the Nth degree Chebyshev polynomial with roots $\{{\theta }_k{\}}_{k=1}^N$. We first expand a function f(x) with respect to the basis of Lagrange’s interpolating polynomials such that $f(x)={\sum }_{l=1}^{N}f({\theta }_l)L_l(x)$ where L_l(x) is an N − 1 degree polynomial with the property L_l(θ_k) = δ_kl. The second derivative matrix is now given by $(D^2{)}_{\mathit{kl}}=L_l^{″}({\theta }_k)$. Throughout this paper, we will refer to this method as the interpolation based method or the standard construction.

The interpolation based second derivative matrix for Chebyshev polynomials suffers from a large norm, scaling as O(N⁴), which can be related to the spectral radius of the matrix. The spectrum of pseudospectral derivative matrices have been studied in, e.g. [22], [24], [8]. Weideman and Trefethen [24] showed that the large norm is related to so-called “spurious” or “out-lying” eigenvalues. The spurious eigenvalues can be related to the highest frequency the derivative matrix can accurately represent. As an example, consider the eigenvalue problem f_xx = λf, f(±1) = 0 with eigenvalues $-\frac{m^2{\pi }^2}{4}\text{,}m=1\text{,}2\text{,}\dots$, and corresponding eigenfunctions $\sin (\frac{m\pi }{2}(x+1))$. The highest frequency the interpolation based Chebyshev derivative matrix can accurately represent is then (asymptotically) given by N/π. The first 2N/π eigenvalues of the interpolation based Chebyshev derivative matrix will match the analytic eigenvalues with high accuracy (growing as O(N²)), whereas the remaining “spurious” eigenvalues will grow as O(N⁴).

This problem can be overcome by mapping the spurious eigenvalues to smaller, non-zero values. This method, spectral filtering, has earlier been investigated by Weideman and Trefethen [24]. However, it was observed that when using this approach for Chebyshev polynomials, the decreased norm comes at the expense of decreased accuracy.

In this paper we introduce a new method, which we call the Eigen-decomposition Pseudo-Spectral (EPS) method. The EPS method is based on writing the derivative operator as an integral operator with a kernel given by the eigen (spectral) decomposition of the derivative operator with boundary conditions imposed. For example, the second derivative operator with zero Dirichlet boundary conditions on the interval [−1, 1] has the eigen decomposition ${\left\{\left(-\frac{m^2{\pi }^2}{4}\text{,}\sin \left(\frac{m\pi (x+1)}{2}\right)\right)\right\}}_{m=1}^{\infty }$. Taking the second derivative of a function f can be expressed via the integral operator

$$f^{″}(x)={\int }_{-1}^{1}K(x\text{,}y)f(y)\mathit{dy}\text{,}$$

where the kernel K(x, y) is given by

$$K(x\text{,}y)=-\sum _{m=1}^{\infty }\frac{m^2{\pi }^2}{4}\sin \left(\frac{m\pi (x+1)}{2}\right)\sin \left(\frac{m\pi (y+1)}{2}\right)\text{.}$$

Expressing the kernel using the first N terms of this sum, Eq. (1) can be approximated by

$$f^{″}(x)\approx -\sum _{m=1}^N\frac{m^2{\pi }^2}{4}\sin \left(\frac{m\pi (x+1)}{2}\right){\int }_{-1}^1\sin \left(\frac{m\pi (y+1)}{2}\right)f(y)\mathit{dy}\text{.}$$

Hence, for a given rank-N approximation of the kernel, constructing the derivative operator has now been reduced to the problem of constructing quadratures for integrals of the type

$${\int }_{-1}^1\sin \left(\frac{m\pi (y+1)}{2}\right)f(y)\mathit{dy}\text{.}$$

In order to select quadratures, we investigate the use of prolate spheroidal wave functions [19], [20] that have long been known to be optimal for representing bandlimited functions, which refers to functions with finite energy and compactly supported Fourier transform. prolate spheroidal wave functions provide an orthonormal basis for such functions while maximizing the energy of the functions within a fixed finite interval in the spatial domain. In a sense they are therefore optimal for representing bandlimited functions restricted to a finite interval.

For the EPS method, quadratures for prolate spheroidal wave functions are appropriate since these are near optimal for the integrals of the type in (2) for smooth functions f(x). However, we stress that the EPS method can be used with any quadrature, and we will also demonstrate the method for Gauss–Legendre and Chebyshev quadratures.

Xiao et al. [26] and Beylkin and Monzon [4] constructed numerical representations of prolate spheroidal wave functions. In this paper, we will use the construction in [4], which we will refer to as approximate prolate spheroidal wave functions (APSWF). In [26], the authors observed that the norm of second derivative matrices for prolate spheroidal wave functions appeared to scale as O(N³) compared to O(N⁴) for Chebyshev derivative matrices. (However, see also [7] for a discussion of this result.) By using APSWFs, one can achieve lower sampling rates even for non-periodic functions (asymptotically two points per wavelengths as N → ∞). Despite the seemingly modest improvement in sampling rate (about 50%) compared to Chebyshev based pseudospectral methods, this improvement leads to significant gains in multidimensional problems. For example, an oversampling factor of 1.5 in each dimension in a three-dimensional problem, results in 1.5³ ≈ 3.4 times larger memory requirements, and typically (at least) 1.5⁴ ≈ 5 times more arithmetic operations for solving time dependent problems. Another advantage with APSWFs is its lower node clustering rate compared to polynomial based quadratures [5].

In [5] the authors constructed APSWF derivative matrices by analytically differentiating the interpolating basis functions derived by Beylkin and Monzon [4]. The boundary conditions were incorporated by using an integration by parts method that had earlier been used for multiwavelets [1].

The construction in [5] was supplemented with a spectral projection step that limits the norm of the derivative matrix, similar to the spectral filtering method used in [24]. This is done by projecting the matrix onto eigenfunctions associated with frequencies within the bandlimit of the basis, which effectively maps spurious eigenvalues to zero. In contrast to the spectral filtering method applied to Chebyshev-based pseudospectral derivative matrices in [24], the spectral filtering applied to APSWF derivative matrices leads to a slight increase in accuracy (see Fig. 7 in [5]).

The construction in [5] is somewhat complicated, and potentially unstable for large problems since the method requires multiplication of highly oscillatory matrices that leads to severe cancellation. The construction therefore requires extended precision when constructing large derivative matrices for high accuracy.

The method introduced in this paper, is a more direct and stable approach than the one in [5], and can also be used with other types of quadratures, such as Gauss–Legendre.

The approach in this paper has several important properties:

• Controlled norm. The norm grows at the rate of the analytic operator, which means that the construction allows larger time steps when solving time dependent problems.

• Easily generalized. The method can be generalized to other operators on the form ${\int }_{-1}^{1}K(x\text{,}y)f(y)\mathit{dy}$ with appropriate kernels.

• No need to include endpoints to enforce boundary conditions. Whereas in traditional pseudospectral methods it is sometimes convenient to use quadrature nodes that include the endpoints (such as Chebyshev–Lobatto nodes) to enforce boundary conditions, the construction in this paper does not require the quadrature nodes to include the endpoints to enforce boundary conditions.

• Fast application algorithms. There are (at least) two algorithms for applying the derivative matrix in O(N log N) operations.

Although the method can be used for a wide range of operators, we focus on first and second derivative operators with Dirichlet or Neumann boundary conditions.

The paper is organized as follows. In Section 2 we describe the general construction of discrete representations of operators of the type in Eq. (1) and provide two error estimates in Section 3.

In Section 4 we provide several examples. In our first application we describe the construction of first and second derivatives on an interval, using both APSWF and Gauss–Legendre quadratures. This section also includes a brief review and comparison of node clustering properties for APSWF and Gauss–Legendre quadratures.

We demonstrate the properties of this construction by solving a wave propagation problem on an interval in both constant and discontinuous velocity media.

Our second application describes the construction of the radial part of the Laplacian operator on a disk in 2D using Chebyshev nodes, which we use to solve Poisson’s equation on the unit disk. We also demonstrate the advantage of using the new construction for time dependent problems in the disk. For each example we provide a detailed description of the constructions and provide several numerical examples illustrating the accuracy and norm of the derivative matrices.

In Section 5 we compare memory requirements, time step, and computational costs for solving wave propagation problems using the EPS method, a Chebyshev pseudospectral method, and an 8th order finite difference scheme.

In Section 6 we outline two fast algorithms that can be used for applying the derivative operator on an interval in O(N logN) operations, and conclude the paper with a discussion in Section 7 and a conclusion in Section 8.