Mathematical Problems in Engineering

Volume 2016, Article ID 7860618, 15 pages

http://dx.doi.org/10.1155/2016/7860618

## Optimized Finite Difference Formulas for Accurate High Frequency Components

Universidad Carlos III de Madrid, Madrid, Spain

Received 8 August 2016; Accepted 9 November 2016

Academic Editor: Xesús Nogueira

Copyright © 2016 Manuel Kindelan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We present a method to obtain optimal finite difference formulas which maximize their frequency range of validity. The optimization is based on the idea of keeping an error of interest (dispersion, phase, or group velocities errors) below a given threshold for a wavenumber interval as large as possible. To find the weights of these optimal finite difference formulas we solve a system of nonlinear equations. Furthermore, we give compact formulas for the optimal weights as function of the error bound. Several numerical experiments illustrate the performance of the obtained finite difference formulas compared to the standard ones.

#### 1. Introduction

Standard finite difference (FD) formulas approximate differential operators by a weighted sum of the values of the function at a set of neighboring nodes (stencil) so as to maximize the numerical accuracy order (the approximation is exact for all polynomials of as high degree as possible). Classical methods to compute the coefficients of the weighted sum are based either on Lagrange’s interpolation polynomial, or on Taylor expansions, or on monomial test functions [1]. However, these methods are computationally very expensive and therefore restricted to relatively small stencil sizes. For equispaced nodes, a short Padé-based algorithm using symbolic algebra was proposed in [2]. The algorithm has the capability of computing the weights for derivatives of any order, approximated to any level of accuracy. For nonequispaced nodes an efficient, recursive algorithm was also proposed in [2]. We refer to [1] and references therein for other efficient algorithms to compute standard FD weights.

A different approach to analyze FD formulas is based on their frequency response. Standard formulas are exact or very accurate for low frequencies but they generate large errors at high frequencies. In wave propagation simulations, for example, this leads to dispersive errors, that is, a progressive distortion of the waveforms as the time and propagation distance increase. To overcome this problem, Holberg [3] proposed in a pioneering work to maximize the range of frequencies for which the group velocity error was bounded within a prescribed value. Since then, many optimization-based FD formulas can be found in the literature (see [4–6] for a review of proposed methods). They differ in the derivative approximated, stencil type and size, objective function, norm minimized, range of frequencies, and optimization method. For example, Holberg [3] uses the group velocity as objective function while Tam and Webb [7] use the phase velocity. Other works include weights into the objective function to enhance the performance of some relevant frequencies [8]. The objective functions are also minimized using different strategies. For example, Kindelan et al. use Newton’s method [9], while Zhang and Yao [6] use a simulated annealing algorithm. On the other hand, Liu [5] uses least squares to derive the optimal FD weights.

In this work, we consider both conventional and staggered equispaced stencils for first and second derivatives. In the case of first derivatives we consider the dispersion error, the error in phase velocity, and the error in group velocity as objective functions. In the case of second derivatives we consider the error in the computation of the second derivative as the objective function. In all cases, we use the infinity norm to determine the weights that maximize the spectral frequency band (or wavenumber range) for which the objective function is bounded by a specified value . To find the optimal weights we use a simple and fast method which requires the solution of a nonlinear system of equations. Furthermore, we provide compact formulas for the weights of the optimal finite difference formulas as a function of .

Based on an observation regarding the locations of the maxima and minima errors for two consecutive stencil sizes, we propose a* combined* method that uses different stencil sizes for odd and even time steps. The* combined* method often results in higher accuracy at the same computational cost.

The paper is organized as follows. Section 2 describes the method used to compute the optimal weights (Appendix A contains the nonlinear system of equations that are solved in each case). Section 3 presents the optimal weights for the different objective functions used to derive optimal weights (Appendix B contains closed form equations for these weights). Section 4 describes several numerical experiments to illustrate the capabilities of the proposed method. Finally, Section 5 contains the main conclusions of the paper.

#### 2. Methodology

The weights of a FD scheme can be computed by using Taylor series expansions in order to maximize the order of the resulting formula. That is, the local truncation error approaches zero as where is the internodal distance. Using this approach, one might reasonably expect that one FD scheme is superior to another if the local truncation error is of higher order. However, a different analysis can be carried out if one is interested in optimizing the FD scheme focusing in application-specific objectives. For instance, in the case of propagating waves the “goodness” of a numerical scheme might not be well quantified by the local truncation error but by the accuracy of the resulting propagating velocity. Indeed, because the numerical velocity of high frequency waves in the grid is different from the true velocity, numerical dispersion related to grid spacing (and size of the time step) occurs, which has a detrimental effect on the accuracy of the propagating velocity of the waves.

In the following, we describe the method used to optimize the frequency response of a FD formula for the first and second derivatives.

##### 2.1. First Derivative

Consider the nondispersive first-order linear transport equationwith constant velocity . This equation admits solutions of the formwhere is the wavenumber and is the angular frequency. Thus, the dispersion relation that determines the angular frequency in terms of the wavenumber is simply . The group velocity, which is the rate of change of with respect to , that is, , and the phase velocity, which is just , are also constant for (1).

The key point is that, in general, the angular frequency is some function of the wavenumber, and, therefore, also the group and phase velocities are functions of the wavenumber. Hence, the dispersion (and the dissipation) properties of FD schemes depend on their numerical or effective wavenumber and angular frequency.

Consider, for example, a numerical approximation to the first derivative of a function on a conventional equispaced grid using a symmetric stencil of nodes, where is an even number, with weights . Then,

The effective wavenumber of this differentiator can be computed by substituting in (3) (we consider time integration exact in (1), so the angular frequency is exact too). Thus,where we have defined the nondimensional effective wavenumber asEquation (5) gives the numerical approximation to the nondimensional analytic wavenumber .

If, instead of (3), we use an approximation to the first derivative using a staggered stencil with nodeswhere is an even number, then the effective wavenumber of the resulting FD scheme is given byNote that both approximations to the analytical wavenumber (5) and (7) depend on the choice of the weights , and, therefore, these approximations can be optimized in a wavenumber interval of interest. Indeed, for a given FD formula with weights , the dispersion errordetermines the difference between the effective and the exact wavenumbers, whose absolute value can be kept below a given threshold using appropriate weights. Moreover, a wave with wavenumber will propagate with velocity , with effective wavenumber given by (5) or (7). Accordingly, the relative error in phase velocity when using a FD scheme to approximate the spatial derivative in (1) is given bySimilarly, the relative error in group velocity is given byThe absolute values of (9) and (10) can also be bounded below a specified value using optimized weights.

Almost all previous optimized schemes use the 2-norm to construct the objective function [6]. In this paper, we use the infinity norm of (8), (9), and (10) for all the cases analyzed. In this way we maximize the wavenumber range for which a peak error is bounded by a specified value. This is the same norm used by Holberg [3] in his seminal paper. In particular, we consider the following problems:Here, is the maximum error allowed in the wavenumber interval . The problem is to find the weights to maximize the length of this interval. is the largest wavenumber for which the objective function is less than or equal to . To solve this optimization problem Holberg [3] used a constrained optimization procedure (although, as pointed out in [6], in practice he considers the 4-norm to use least squares to determine the optimal weights). In [6], the authors use the infinity norm and they recur to a simulated annealing algorithm to solve the corresponding optimization problem. In this paper, we use instead a simple procedure which reduces the problem to the solution of a small system of nonlinear equations. This approach was first proposed in [9] using staggered nodes and the relative error in group velocity as the objective function.

To illustrate the proposed method, let us consider the optimization problem for the dispersion error using an stencil in a conventional grid. We impose the conditions that at extrema , and that at these extrema the derivative of with respect to is zero. Thus, we solve the following system of nonlinear equations:Solving this system of equations we find the sought weights , and the locations of the extrema . The nonlinear systems which are solved for each of the cases analyzed in this paper are detailed in Appendix A.

##### 2.2. Second Derivative

The effective wavenumber of a FD formula used to approximate a second derivative on a conventional equispaced grid using an symmetric stencil ( even number) with weights , can be computed as follows. Substituting inand setting , we obtainFrom this equation we find that the square of the effective wavenumber isThus, the error in the computation of the second derivative is given bySimilarly to what was done for the first derivative, we compute the optimal FD weights , for the second derivative by maximizing the length of the wavenumber interval for which the approximation error (16) is bounded by a specified value . ThusThe nonlinear equations to be solved areThis is a system of equations with unknowns ( weights , and locations of the extrema ). Equation (20) is used to enforce the approximation error to be zero at . This condition was also used in [6]. Slightly better results are obtained if condition (20) is relaxed and we allow the error at to be different from zero, though. In this case, equation (20) has to be replaced by

#### 3. Optimal Weights

In this section we present the results of our method. Appendix B contains closed form equations for the optimal weights.

##### 3.1. First Derivative

Figure 1 compares the dispersion errors (8) for FD formulas in conventional grids using standard (dashed lines) and optimized (solid lines) weights. From left to right, the curves correspond to node stencils with equal to 4, 6, 8, 10, and . The optimal weights have been obtained by solving equations (12) with . Figure 2 enlarges the interval around . Note that the spectral improvement when optimal weights are used is more significant for large . For example, for the maximum wavenumber for which the dispersion error is smaller than changes from with standard weights to with optimal weights, while for it changes from to .