Optimized Finite Difference Formulas for Accurate High Frequency Components

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.


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 optimizationbased FD formulas can be found in the literature (see [4][5][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 2 Mathematical Problems in Engineering 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.

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.

First Derivative. Consider the nondispersive first-order linear transport equation
with constant velocity .This equation admits solutions of the form where   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  + 1 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 k as Equation (5) gives the numerical approximation k to the nondimensional analytic wavenumber  =   ℎ.
If, instead of (3), we use an approximation to the first derivative using a staggered stencil with  + 1 nodes where  is an even number, then the effective wavenumber of the resulting FD scheme is given by Note that both approximations to the analytical wavenumber ( 5) and (7) depend on the choice of the weights   ,  = 1, . . ., /2, and, therefore, these approximations can be optimized in a wavenumber interval of interest.Indeed, for a given FD formula with weights   ,  = 1, . . ., /2, the dispersion error determines 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 ( k())/, with effective wavenumber k given by ( 5) or (7).Accordingly, the relative error in phase velocity   () = (  ( k) −   ())/  () when using a FD scheme to approximate the spatial derivative in (1) is given by Similarly, the relative error in group velocity   () = (  ( k)−   ())/  () is given by The 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 [0,   ].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 +1 stencil in a conventional grid.We impose the conditions that at /2 extrema || = ±, and that at these extrema the derivative of  with respect to  is zero.Thus, we solve the following system of  nonlinear equations: (12) Solving this system of  equations we find the /2 sought weights   , and the /2 locations of the extrema   .The nonlinear systems which are solved for each of the cases analyzed in this paper are detailed in Appendix A.
Mathematical Problems in Engineering 2.2.Second Derivative.The effective wavenumber of a FD formula used to approximate a second derivative on a conventional equispaced grid using an  + 1 symmetric stencil ( even number) with weights   ,  = 1, . . ., /2, can be computed as follows.Substituting () =     in and setting  =   ℎ, we obtain From this equation we find that the square of the effective wavenumber is Thus, the error in the computation of the second derivative is given by Similarly to what was done for the first derivative, we compute the optimal FD weights   ,  = 1, . . ., /2, for the second derivative by maximizing the length of the wavenumber interval [0,   ] for which the approximation error ( 16) is bounded by a specified value .Thus The nonlinear equations to be solved are This is a system of  + 1 equations with  + 1 unknowns (/2 + 1 weights   , and /2 locations of the extrema   ).Equation ( 20) is used to enforce the approximation error to be zero at  = 0.This condition was also used in [6].Slightly better results are obtained if condition ( 20) is relaxed and we allow the error at  = 0 to be different from zero, though.In this case, equation (20) has to be replaced by

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  + 1 node stencils with  equal to 4, 6, 8, 10, and 12.The optimal weights have been obtained by solving equations ( 12) with  = 10 −4 .Figure 2 enlarges the interval around  = 0. Note that the spectral improvement when optimal weights are used is more significant for large .For example, for  = 4 the maximum wavenumber for which the dispersion error is smaller than  = 10 −4 changes from   = 0.31 with standard weights to   = 0.54 with optimal weights, while for  = 12 it changes from   = 1.05 to   = 1.78.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  + 2. This observation, which happens to be true for  less than or equal to 14, 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  = 4, 6, 8, 10, and 12.These values are in good agreement with those shown in Table 1 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  = 8, the weights computed in [6] result in   = 1.3009, while the weights in Table 1 result in   = 1.3112.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 ( = 26 and 28), 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  = 10 −4 .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 +1 nodes ( = 4, 6, 8, 10, and 12) 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  = 4, the optimal weights are given to leading order by In 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 for which the corresponding error is bounded by a given threshold .The values of PPW 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  + 1.From top to bottom, the curves depicted in this figure correspond to optimized FD formulas with  = 2, 4, 6, 8, 10, and 12.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 PPW through (23).
(5) Using Figure 5, PPW and  define the number of nodes  + 1 of the FD formula that fulfills both requirements.
For example, let the internodal distance and the maximum error be ℎ = 1 and  ≤ 10 −4 .Consider an initial condition whose maximum significative wavenumber   is 1.5; then, according to (23), the number of points per wavelength is PPW = 4.1888.Finally, 5 shows that optimized FD formula with less nodes that fulfills all the requirements is the one with  = 10. 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 ℎ.

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  = 10 −4 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  = 10 −4 .
The top rows of Table 3 compare the maximum wavenumber   of the standard and optimal FD formulas for the second  derivative in the case  = 10 −4 .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 1.69-1.74.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 2 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  = 12 the weights computed in [6] result in   = 2.017, while using the weights in Table 3 results in   = 2.0442.Similar improvements of around 1% are obtained for other values of .

Numerical Experiments
To illustrate the performance of the optimized schemes proposed in this paper, we will consider the 1D advection equation ( 1) with velocity  = 1 and different initial conditions.In all the cases considered here, we use conventional stencils with step size ℎ = 1, and periodic boundary conditions.Time integration is carried out using a fourthorder Runge-Kutta method with time step Δ = 0.02.Hence, the numerical errors of the approximations presented here are dominated by the spatial discretization.In fact, we have repeated the numerical experiment with Δ = 0.01 obtaining identical results.Let us start by considering the propagation of a pure wave with wavenumber  = /2.Thus, we solve (1) with initial condition Figure 7 shows the exact solution at  = 20 (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  = 20, and therefore (20) = 0.The squares correspond to the solution obtained with the optimal weights corresponding to  = 12 nodes conventional stencil with |  | <  ≡ 10 −2 .Notice that (20) = −0.2771instead 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   (/2) = 0.008938   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 stencil of size  + 2 (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  + 2 nodes at  +1 , a stencil with  nodes at  +2 , 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  = 10 conventional stencil, and even time steps with a  = 12 conventional stencil.Notice that the results obtained with this combined method are much more accurate ((20) = −0.03677)and have the same computational cost.
We consider next the initial condition where  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  = 4,  = 9.
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  = /2.Figure 9 displays the numerical solution at  = 200.Figure 9(a) shows the numerical solution (solid line) obtained using the standard weights corresponding to a stencil with  = 12 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  = /2, the error of the group velocity using standard  = 12 nodes stencil is   = −0.0693.Thus, the group velocity is V  = 1 − 0.0693 = 0.9307 and the location of the maximum at  = 200 is  = 50 + 200V  = 236.1,which is in very good agreement with the location of the maximum shown in Figure 9(a).
Notice in Table 2 that for  = 12 stencil with optimized weights, the maximum wavenumber for which |  | ≤ 10 −4 is   = 1.4586.Since the spectrum of the initial condition (25) peaks at  = /2, we can not use in this case those optimized weights.Instead, we use the optimized weights corresponding to |  | ≤ 10 −3 , for which   = 1.7229 > /2.These optimized weights can be obtained using the formulas in Table 7.The numerical results at  = 200, 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 compute where  is the numerical solution and   is the exact solution.For the solution obtained with optimized weights (Figure 9(b)),  num = 0.0216 and ‖ −   ‖ ∞ = 0.0166.Finally, we consider the initial condition defined in (25) with  = 8,  = 9. Figure 10(a) displays the corresponding normalized power spectrum.In this case, the normalized power spectrum peaks at  = /4, and the largest significant wavenumbers are around  = 1.25.We have integrated (1) with this initial condition until  = 200 using optimized weights corresponding to |  | ≤ 10 −3 and sizes  = 10 and  = 12.In both cases   is well above the wavenumbers Mathematical Problems in Engineering The numerical solution for  = 200 using the combined method is shown in Figure 10(b) (solid line) together with the exact solution (dashed line).

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. (A.4) Staggered,   : (A.5) 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 We find that  = 2/ 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 ( = 0.4 for  = 4).
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   ,  = 0, . . ., 5. Thus, they are slightly different from the weights of the Taylor series of   in powers of   , but they are more precise in the range 10 −5 ≤  ≤ 10 −1 ,

Figure 5 :
Figure 5: Points per wavelength as a function of the maximum dispersion error allowed  for  = 2 (top curve) to  = 12 (bottom curve).

Table 1 :
Two first rows: maximum wavenumber for which the dispersion error is less than or equal to  = 10 −4 .Bottom rows: corresponding optimized weights   .

Table 2 :
Largest wavenumbers   for which the corresponding error is bound by  = 10 −4 .These values have been obtained by solving the nonlinear systems shown in Appendix A for the optimal weights., the values computed using the equations in Table5are undistinguishable from those shown in Figure4.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 wavelength stencils