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 .
We observe in Figure 2 that, for fixed , the locations of where the errors are approximately coincide with the locations of where the errors are for . This observation, which happens to be true for less than or equal to , can be exploited in time dependent problems to improve the accuracy of the numerical approximations without additional computational cost (see Section 4 for a detailed discussion on this issue).
The bottom rows of Table 1 also show the optimal weights for , and . These values are in good agreement with those shown in Table of [6] that are computed using a simulated annealing optimization procedure. We mention, however, that we obtain slightly larger wavenumber ranges for which the dispersion errors are smaller than . For instance, for , the weights computed in [6] result in , while the weights in Table 1 result in . We stress, though, that the advantage of the method used here compared to the one used in [6] is that the former is much simpler and computationally less expensive. This allows us to compute optimal weights for FD formulas with a higher number of nodes ( and ), as we show in Figure 3.
Table 2 summarizes, for all the cases analyzed in this paper, the values of the largest wavenumbers for which the corresponding error is bounded by . For each of these cases we compare the values of using the standard FD weights (rows with in the second column) with those obtained using the optimal weights (rows with in the second column). The different values of shown in Table 2 have been obtained by solving the systems of nonlinear equations in Appendix A.
Figure 4 shows the optimal weights for a stencil with nodes (, and ) on a conventional grid as a function of the error bound .
In each of the plots depicted in this figure, the optimal weights increase in absolute value as increases. As expected, the optimal weights approach the values corresponding to standard FD formulas when the error decreases to zero.
To facilitate the use of the optimized FD formulas we provide closed form equations for the optimal weights as a function of the error bound . For example, for , the optimal weights are given to leading order byIn Appendix B we provide closed form equations for the optimal weights in all the cases analyzed in this paper. For instance, in the case of dispersion error for conventional stencils, the values computed using the equations in Table 5 are undistinguishable from those shown in Figure 4.
Optimized FD schemes are designed to require only a few points per wavelength for accurate resolution of PDEs. The (nondimensional) spatial bandwidth determines the number of grid points per wavelengthrequired to resolve the shortest wavelength component for which the corresponding error is bounded by a given threshold . The values of as a function of for the case of weights that optimize the dispersion error on a conventional grid are shown in Figure 5 for stencils with different number of nodes . From top to bottom, the curves depicted in this figure correspond to optimized FD formulas with , and . For each of these curves, the points per wavelength are greatly reduced as the required increases. This fact is more relevant for optimized FD formulas with small support (upper curves). On the other hand, as expected, the FD formulas that use a larger support to approximate the first derivative (lower curves), and therefore require more computational effort, perform better since they require fewer points per wavelength.
Figure 5 or similar figures for the different cases analyzed are very useful to determine which finite difference formula should be used to solve a particular problem. In fact, given a PDE with initial and/or boundary conditions, the problem is to approximate its solution within a given dispersion, phase, or group velocity error for all the significant wavenumbers involved in it. Furthermore, one would like to use the simplest possible FD formula, that is, the one with as less nodes as possible, in order to minimize the computational cost. If the discretization length is given one should perform the following steps:(1)Choose the maximum admissible error .(2)Compute the Fourier transform of the initial condition.(3)Pick the maximum significative wavenumber to be resolved below the error .(4)The corresponding shortest wavelength gives the number of points per wavelength through (23).(5)Using Figure 5, and define the number of nodes of the FD formula that fulfills both requirements.
For example, let the internodal distance and the maximum error be and . Consider an initial condition whose maximum significative wavenumber is ; then, according to (23), the number of points per wavelength is . Finally, 5 shows that optimized FD formula with less nodes that fulfills all the requirements is the one with . Figure 5 refers to the dispersion errors, but any other error could be used as well following the same procedure.
If the discretization length is not fixed, then you choose the stencil length and the maximum admissible error and use Figure 5 to determine the value of PPW. Then use equation (23) to determine the value of .
3.2. Second Derivative
Figure 6 shows the errors in the approximation to the second derivative (16) for standard (dashed lines) and optimized (solid lines) FD formulas as a function of the wavenumber . The latter have been computed by solving (18)–(20) with in order to determine the optimal weights . Compared to the standard FD formulas, there is a significant increase in wavenumbers for which the error is kept below .
The top rows of Table 3 compare the maximum wavenumber of the standard and optimal FD formulas for the second derivative in the case . The second row corresponds to (18)–(20), and the third row to (18)-(19) and (21). In the first case (see (18)–(20)), the ratio between the optimal and the standard is in the range . Notice that there is a significant improvement when using instead (18)-(19) and (21) for small values of . This advantage becomes, however, negligible for large values of .
Table 3 also contains the optimal weights for second derivative using (18)–(20). These are in good agreement with the values shown in Table of [6] using a simulated annealing optimization procedure. As was the case with the first derivative, with our procedure the results are slightly better. For instance, in the case the weights computed in [6] result in , while using the weights in Table 3 results in . Similar improvements of around % are obtained for other values of .
4. Numerical Experiments
To illustrate the performance of the optimized schemes proposed in this paper, we will consider the 1D advection equation (1) with velocity and different initial conditions. In all the cases considered here, we use conventional stencils with step size , and periodic boundary conditions. Time integration is carried out using a fourth-order Runge–Kutta method with time step . Hence, the numerical errors of the approximations presented here are dominated by the spatial discretization. In fact, we have repeated the numerical experiment with obtaining identical results.
Let us start by considering the propagation of a pure wave with wavenumber . Thus, we solve (1) with initial condition
Figure 7 shows the exact solution at (solid line) together with the numerical solution obtained using optimized FD formulas of different support and maximum allowed error . Notice that the exact solution is the initial wave (24) centered at , and therefore . The squares correspond to the solution obtained with the optimal weights corresponding to nodes conventional stencil with . Notice that instead of zero (large square) and, therefore, the propagation velocity of the numerical solution is slightly higher than one. In fact, using (9) the error in phase velocity is and, therefore, at the wave has travelled a distance . This value exactly coincides with the displacement observed in the numerical solution . The circles correspond to the solution obtained with the optimal weights corresponding to nodes conventional stencil with . In this case, (large circle) and, therefore, the propagation velocity is smaller than one. In fact, and , which again coincides with the numerical solution . These results are shown in Table 4.
An interesting possibility is to efficiently combine two stencils of consecutive sizes. In fact, it was pointed out in Section 3 that the errors for two consecutive stencil sizes tend to cancel each other: at locations where the error is using a stencil of size , it is approximately using a stencil of size (see Figure 2). Thus, one can alternate between two stencil sizes during the time integration in order to reduce the errors. That is, one can use a stencil with nodes at time step , a stencil with nodes at , a stencil with nodes at , and so on. This technique has been used to obtain the numerical solution represented by triangles in Figure 7. It has been obtained by integrating odd step times with a conventional stencil, and even time steps with a conventional stencil. Notice that the results obtained with this combined method are much more accurate () and have the same computational cost.
We consider next the initial conditionwhere is the dominant wavelength and is the half-width of the Gaussian function. This is the same initial condition considered in [10]. We consider the case , . Figure 8(a) shows the initial condition (25) (notice that there are about four points per wavelength), and Figure 8(b) the corresponding normalized power spectrum which is peaked at .
(a)
(b)
Figure 9 displays the numerical solution at . Figure 9(a) shows the numerical solution (solid line) obtained using the standard weights corresponding to a stencil with equispaced nodes. The exact solution is shown with dashed lines. Notice that there are significant errors because of both numerical dissipation and dispersion. In fact, for , the error of the group velocity using standard nodes stencil is . Thus, the group velocity is and the location of the maximum at is , which is in very good agreement with the location of the maximum shown in Figure 9(a).
(a)
(b)
Notice in Table 2 that for stencil with optimized weights, the maximum wavenumber for which is . Since the spectrum of the initial condition (25) peaks at , we can not use in this case those optimized weights. Instead, we use the optimized weights corresponding to , for which . These optimized weights can be obtained using the formulas in Table 7. The numerical results at , displayed in Figure 9(b) with a solid line, show a remarkable good agreement with the exact solution (dashed lines).
To quantify the agreement between the exact and the numerical solution we computewhere is the numerical solution and is the exact solution. For the solution obtained with optimized weights (Figure 9(b)), and .
Finally, we consider the initial condition defined in (25) with , . Figure 10(a) displays the corresponding normalized power spectrum. In this case, the normalized power spectrum peaks at , and the largest significant wavenumbers are around . We have integrated (1) with this initial condition until using optimized weights corresponding to and sizes and . In both cases is well above the wavenumbers contained in the initial condition. The errors in the case areand in the case However, if we combine both stencils by integrating odd step times with the stencil, and the even time steps with the stencil, the results are much more accurate: The numerical solution for using the combined method is shown in Figure 10(b) (solid line) together with the exact solution (dashed line).
(a)
(b)
5. Conclusions
In this paper, we have computed optimized weights for FD formulas for the first and the second derivatives. This weights are optimized in the sense that they allow small errors (dispersion, phase, or group velocity errors) for low frequencies but maximize the interval of wavenumbers for which the error is smaller than a given threshold. We also provide explicit formulas that give the values of the optimized weights as functions of the allowed error.
As a byproduct, we have also showed that integrating a time dependent equation combining stencils with two consecutive numbers of nodes gives much more accurate results than using only one of them. The gain in accuracy is achieved without an increase of the computational cost.
Appendix
A. Nonlinear System of Equations
The nonlinear equations to be solved in the other cases are as follows.
Conventional, :
Conventional, :
Conventional, :
Staggered, :
Staggered, :
Staggered, :
B. Closed Form Equations for Optimal Finite Difference Weights
In this appendix we present closed form equations to compute the optimal weights for all the cases analyzed. Figure 4 shows that, as approaches zero, the difference between the optimal weight and the standard finite difference weight also approaches zero. In fact, it is straightforward to check that this difference approaches zero as . We compute the exponent numerically as
We find that for the optimal weights corresponding to phase velocity, , and group velocity, . This result was derived analytically in [9] for the case of group velocity with staggered stencils. For the optimal weights corresponding to dispersion error, , we find that the exponent is slightly smaller ( for ).
Once the exponent is determined we derive a closed form formula for the optimal weights; namely,The weights are given in Tables 5–10 for all the cases analyzed. They have been derived by computing the projection of the numerical values on the functional space spanned by the functions . Thus, they are slightly different from the weights of the Taylor series of in powers of , but they are more precise in the range , if only a small number of terms () are used to compute .
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work has been supported by the Spanish MICINN Grant FIS2013-41802-R.