# Optimized Finite Difference Formulas for Accurate High Frequency Components

Authors: Manuel Kindelan; Miguel Moscoso; Pedro Gonzalez-Rodriguez
Journal: Mathematical Problems in Engineering
 (2016)

---

## 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.


---

## Body

## 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 valueE. 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 E.Based on an observation regarding the locations of the maxima and minima errors for two consecutive stencil sizes, we propose acombined method that uses different stencil sizes for odd and even time steps. Thecombined method often results in higher accuracy at the same computational cost.The paper is organized as follows. Section2 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 orderp of the resulting formula. That is, the local truncation error approaches zero as O(hp) where h 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 equation(1)ut+cux=0,with constant velocity c. This equation admits solutions of the form(2)ux,t=eik′x-ωt,where k′ is the wavenumber and ω is the angular frequency. Thus, the dispersion relation that determines the angular frequency in terms of the wavenumber is simply ω(k′)=ck′. The group velocity, which is the rate of change of ω(k′) with respect to k′, that is, Vg(k′)=dω(k′)/dk′=c, and the phase velocity, which is just Vp(k′)=ω(k′)/k′=c, are also constant for (1).The key point is that, in general, the angular frequency is some functionω=ω(k′) 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 functionf on a conventional equispaced grid using a symmetric stencil of N+1 nodes, where N is an even number, with weights dj. Then,(3)dfdx≈∑j=1N/2djfx+jh-fx-jhh.The effective wavenumber of this differentiator can be computed by substitutingf(x)=eik′x in (3) (we consider time integration exact in (1), so the angular frequency is exact too). Thus,(4)ik′eik′x≈2ieik′x∑j=1N/2djsink′jhh≡ik~heik′x,where we have defined the nondimensional effective wavenumber k~ as(5)k~k=2∑j=1N/2djsin⁡kj.Equation (5) gives the numerical approximation k~ to the nondimensional analytic wavenumber k=k′h.If, instead of (3), we use an approximation to the first derivative using a staggered stencil with N+1 nodes(6)dfdx≈∑j=1N/2djfx+2j-1h/2-fx-2j-1h/2h,where N is an even number, then the effective wavenumber of the resulting FD scheme is given by(7)k~=2∑j=1N/2djsin⁡k2j-12.Note that both approximations to the analytical wavenumber (5) and (7) depend on the choice of the weights dj,j=1,…,N/2, and, therefore, these approximations can be optimized in a wavenumber interval of interest. Indeed, for a given FD formula with weights dj,j=1,…,N/2, the dispersion error(8)ϵk=k~k-kdetermines 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 k will propagate with velocity ω(k~(k))/k, with effective wavenumber k~ given by (5) or (7). Accordingly, the relative error in phase velocity Ep(k)=(Vp(k~)-Vp(k))/Vp(k) when using a FD scheme to approximate the spatial derivative in (1) is given by(9)Epk=k~k-1.Similarly, the relative error in group velocity Eg(k)=(Vg(k~)-Vg(k))/Vg(k) is given by(10)Egk=dk~dk-1.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:(11)GridtypeObjectivefunctionconventionalϵ=2∑j=1N/2djsin⁡kj-k≤E,Ep=2k∑j=1N/2djsin⁡kj-1≤E,Eg=2∑j=1N/2djjcos⁡kj-1≤E,staggeredϵ=2∑j=1N/2djsin⁡k2j-12-k≤E,Ep=2k∑j=1N/2djsin⁡k2j-12-1≤E,Eg=∑j=1N/2dj2j-1cos⁡k2j-12-1≤E.Here, E is the maximum error allowed in the wavenumber interval [0,kc]. The problem is to find the weights dj to maximize the length of this interval. kc is the largest wavenumber for which the objective function is less than or equal to E. 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 N+1 stencil in a conventional grid. We impose the conditions that at N/2 extrema |ϵ|=±E, and that at these extrema the derivative of ϵ with respect to k is zero. Thus, we solve the following system of N nonlinear equations:(12)2∑j=1N/2djsin⁡kij-ki+E-1N/2+i+1=0,i=1,2,…,N2,2∑j=1N/2djjcos⁡kij-1=0,i=1,2,…,N2.Solving this system of N equations we find the N/2 sought weights dj, and the N/2 locations of the extrema ki. 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 anN+1 symmetric stencil (N even number) with weights dj,j=1,…,N/2, can be computed as follows. Substituting f(x)=eik′x in(13)d2fdx2≈d0fxh2+∑j=1N/2djfx+jh+djfx-jhh2,and setting k=k′h, we obtain(14)-k2eikx≈eikxd0+2∑j=1N/2djcos⁡kj≡-eikxk~2.From this equation we find that the square of the effective wavenumber is(15)k~2=-d0+2∑j=1N/2djcos⁡kj.Thus, the error in the computation of the second derivative is given by(16)εk=k~2-k2=-d0+2∑j=1N/2djcos⁡kj-k2.Similarly to what was done for the first derivative, we compute the optimal FD weights dj,j=1,…,N/2, for the second derivative by maximizing the length of the wavenumber interval [0,kc] for which the approximation error (16) is bounded by a specified value E. Thus(17)-d0+2∑j=1N/2djcos⁡kj-k2≤E,0≤k≤kc.The nonlinear equations to be solved are(18)-d0-2∑j=1N/2djcos⁡kij-ki2+E-1N/2+i+1=0,i=1,…,N2,(19)2∑j=1N/2djjsin⁡kij-2ki=0,i=1,…,N2,(20)d0+2∑j=1N/2dj=0.This is a system of N+1 equations with N+1 unknowns (N/2+1 weights dj, and N/2 locations of the extrema ki). Equation (20) is used to enforce the approximation error to be zero at k=0. This condition was also used in [6]. Slightly better results are obtained if condition (20) is relaxed and we allow the error at k=0 to be different from zero, though. In this case, equation (20) has to be replaced by(21)-d0-2∑j=1N/2dj+E-1N/2+1=0.


## 2.1. First Derivative

Consider the nondispersive first-order linear transport equation(1)ut+cux=0,with constant velocity c. This equation admits solutions of the form(2)ux,t=eik′x-ωt,where k′ is the wavenumber and ω is the angular frequency. Thus, the dispersion relation that determines the angular frequency in terms of the wavenumber is simply ω(k′)=ck′. The group velocity, which is the rate of change of ω(k′) with respect to k′, that is, Vg(k′)=dω(k′)/dk′=c, and the phase velocity, which is just Vp(k′)=ω(k′)/k′=c, are also constant for (1).The key point is that, in general, the angular frequency is some functionω=ω(k′) 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 functionf on a conventional equispaced grid using a symmetric stencil of N+1 nodes, where N is an even number, with weights dj. Then,(3)dfdx≈∑j=1N/2djfx+jh-fx-jhh.The effective wavenumber of this differentiator can be computed by substitutingf(x)=eik′x in (3) (we consider time integration exact in (1), so the angular frequency is exact too). Thus,(4)ik′eik′x≈2ieik′x∑j=1N/2djsink′jhh≡ik~heik′x,where we have defined the nondimensional effective wavenumber k~ as(5)k~k=2∑j=1N/2djsin⁡kj.Equation (5) gives the numerical approximation k~ to the nondimensional analytic wavenumber k=k′h.If, instead of (3), we use an approximation to the first derivative using a staggered stencil with N+1 nodes(6)dfdx≈∑j=1N/2djfx+2j-1h/2-fx-2j-1h/2h,where N is an even number, then the effective wavenumber of the resulting FD scheme is given by(7)k~=2∑j=1N/2djsin⁡k2j-12.Note that both approximations to the analytical wavenumber (5) and (7) depend on the choice of the weights dj,j=1,…,N/2, and, therefore, these approximations can be optimized in a wavenumber interval of interest. Indeed, for a given FD formula with weights dj,j=1,…,N/2, the dispersion error(8)ϵk=k~k-kdetermines 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 k will propagate with velocity ω(k~(k))/k, with effective wavenumber k~ given by (5) or (7). Accordingly, the relative error in phase velocity Ep(k)=(Vp(k~)-Vp(k))/Vp(k) when using a FD scheme to approximate the spatial derivative in (1) is given by(9)Epk=k~k-1.Similarly, the relative error in group velocity Eg(k)=(Vg(k~)-Vg(k))/Vg(k) is given by(10)Egk=dk~dk-1.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:(11)GridtypeObjectivefunctionconventionalϵ=2∑j=1N/2djsin⁡kj-k≤E,Ep=2k∑j=1N/2djsin⁡kj-1≤E,Eg=2∑j=1N/2djjcos⁡kj-1≤E,staggeredϵ=2∑j=1N/2djsin⁡k2j-12-k≤E,Ep=2k∑j=1N/2djsin⁡k2j-12-1≤E,Eg=∑j=1N/2dj2j-1cos⁡k2j-12-1≤E.Here, E is the maximum error allowed in the wavenumber interval [0,kc]. The problem is to find the weights dj to maximize the length of this interval. kc is the largest wavenumber for which the objective function is less than or equal to E. 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 N+1 stencil in a conventional grid. We impose the conditions that at N/2 extrema |ϵ|=±E, and that at these extrema the derivative of ϵ with respect to k is zero. Thus, we solve the following system of N nonlinear equations:(12)2∑j=1N/2djsin⁡kij-ki+E-1N/2+i+1=0,i=1,2,…,N2,2∑j=1N/2djjcos⁡kij-1=0,i=1,2,…,N2.Solving this system of N equations we find the N/2 sought weights dj, and the N/2 locations of the extrema ki. 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 anN+1 symmetric stencil (N even number) with weights dj,j=1,…,N/2, can be computed as follows. Substituting f(x)=eik′x in(13)d2fdx2≈d0fxh2+∑j=1N/2djfx+jh+djfx-jhh2,and setting k=k′h, we obtain(14)-k2eikx≈eikxd0+2∑j=1N/2djcos⁡kj≡-eikxk~2.From this equation we find that the square of the effective wavenumber is(15)k~2=-d0+2∑j=1N/2djcos⁡kj.Thus, the error in the computation of the second derivative is given by(16)εk=k~2-k2=-d0+2∑j=1N/2djcos⁡kj-k2.Similarly to what was done for the first derivative, we compute the optimal FD weights dj,j=1,…,N/2, for the second derivative by maximizing the length of the wavenumber interval [0,kc] for which the approximation error (16) is bounded by a specified value E. Thus(17)-d0+2∑j=1N/2djcos⁡kj-k2≤E,0≤k≤kc.The nonlinear equations to be solved are(18)-d0-2∑j=1N/2djcos⁡kij-ki2+E-1N/2+i+1=0,i=1,…,N2,(19)2∑j=1N/2djjsin⁡kij-2ki=0,i=1,…,N2,(20)d0+2∑j=1N/2dj=0.This is a system of N+1 equations with N+1 unknowns (N/2+1 weights dj, and N/2 locations of the extrema ki). Equation (20) is used to enforce the approximation error to be zero at k=0. This condition was also used in [6]. Slightly better results are obtained if condition (20) is relaxed and we allow the error at k=0 to be different from zero, though. In this case, equation (20) has to be replaced by(21)-d0-2∑j=1N/2dj+E-1N/2+1=0.


## 3. Optimal Weights

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


### 3.1. First Derivative

Figure1 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 N+1 node stencils with N equal to 4, 6, 8, 10, and 12. The optimal weights have been obtained by solving equations (12) with E=10-4. Figure 2 enlarges the interval around ϵ=0. Note that the spectral improvement when optimal weights are used is more significant for large N. For example, for N=4 the maximum wavenumber for which the dispersion error is smaller than E=10-4 changes from kc=0.31 with standard weights to kc=0.54 with optimal weights, while for N=12 it changes from kc=1.05 to kc=1.78.Figure 1
Dispersion errorϵ(k) for conventional stencils with N+1 nodes (from left to right N=2,4,6,8,10, and 12). Dashed lines: standard FD formulas. Solid lines: optimal FD formulas for E=10-4.Figure 2
Detail of Figure1. Dashed lines: standard FD formulas. Solid lines: optimal FD formulas for E=10-4. From left to right N=2,4,6,8,10, and 12.We observe in Figure2 that, for fixed N, the locations of k where the errors are ±E approximately coincide with the locations of k where the errors are ∓E for N+2. This observation, which happens to be true for N 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 Table1 also show the optimal weights dj for N=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 kc for which the dispersion errors are smaller than E. For instance, for N=8, the weights computed in [6] result in kc=1.3009, while the weights in Table 1 result in kc=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 (N=26 and 28), as we show in Figure 3.Table 1
Two first rows: maximum wavenumber for which the dispersion error is less than or equal toE=10-4. Bottom rows: corresponding optimized weights dj.
N = 4 N = 6 N = 8 N = 10 N = 12 k c standard 0.31309 0.5475 0.74763 0.91369 1.052 k c optimal 0.54129 0.96452 1.3112 1.5812 1.7888  d 1 0.678935 0.779124 0.842632 0.884271 0.912322 d 2 −0.089700 −0.174992 −0.246700 −0.302511 −0.3444207 d 3  0.023744 0.061634 0.102893 0.140765 d 4   −0.008624 −0.026887 −0.050463 d 5    0.0040000 0.013796 d 6     −0.002218Figure 3
Dispersion errorϵ(k) for optimized FD formulas with N+1 conventional nodes (N=26 and 28).Table2 summarizes, for all the cases analyzed in this paper, the values of the largest wavenumbers kc for which the corresponding error is bounded by E=10-4. For each of these cases we compare the values of kc using the standard FD weights (rows with s in the second column) with those obtained using the optimal weights (rows with o in the second column). The different values of kc shown in Table 2 have been obtained by solving the systems of nonlinear equations in Appendix A.Table 2
Largest wavenumberskc for which the corresponding error is bound by E=10-4. These values have been obtained by solving the nonlinear systems shown in Appendix A for the optimal weights.
N = 2 N = 4 N = 6 N = 8 N = 10 N = 12 Conventional nodes  ϵ s 0.0843 0.3136 0.5480 0.7481 0.9141 1.0525 ϵ o 0.1338 0.5415 0.9651 1.3121 1.5814 1.7465 E p s 0.0236 0.2338 0.4944 0.7204 0.9051 1.0571 E p o 0.0332 0.3919 0.8624 1.2582 1.5622 1.7922 E g s  0.1561 0.3565 0.5441 0.7061 0.8445 E g o  0.2621 0.6236 0.9589 1.2350 1.4586  Staggered nodes  ϵ s 0.1339 0.4644 0.7654 1.0022 1.1889 1.3358 ϵ o 0.2125 0.8011 1.3396 1.7337 2.0120 2.2114 E p s 0.0490 0.38301 0.7313 1.0025 1.2112 1.3724 E p o 0.0693 0.6404 1.2648 1.7224 2.0368 2.2551 E g s 0.0283 0.2559 0.5269 0.7563 0.9427 1.0928 E g o 0.0400 0.4282 0.9155 1.3138 1.6136 1.8385Figure4 shows the optimal weights for a stencil with N+1 nodes (N=4,6,8,10, and 12) on a conventional grid as a function of the error bound E.Figure 4
Optimal weightsdj as a function of E, for which |ϵ|<E. The different curves correspond to N=4,6,8,10,12.In each of the plots depicted in this figure, the optimal weightsdj increase in absolute value as N increases. As expected, the optimal weights approach the values corresponding to standard FD formulas when the error E 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 boundE. For example, for N=4, the optimal weights are given to leading order by(22)d1≈23+0.2193E3/8+⋯,d2≈-112-0.1096E3/8+⋯.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 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 bandwidthkc determines the number of grid points per wavelength(23)PPW=2πkc=2πkc′h=λc′hrequired to resolve the shortest wavelength component λc′ for which the corresponding error is bounded by a given threshold E. The values of PPW as a function of E 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 N+1. From top to bottom, the curves depicted in this figure correspond to optimized FD formulas with N=2,4,6,8,10, and 12. For each of these curves, the points per wavelength are greatly reduced as the required E 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
Points per wavelength as a function of the maximum dispersion error allowedE for N=2 (top curve) to N=12 (bottom curve).Figure5 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 E 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 h is given one should perform the following steps:(1)
Choose the maximum admissible errorE.(2)
Compute the Fourier transform of the initial condition.(3)
Pick the maximum significative wavenumberkc to be resolved below the error E.(4)
The corresponding shortest wavelengthλc gives the number of points per wavelength PPW through (23).(5)
Using Figure5, PPW and E define the number of nodes N+1 of the FD formula that fulfills both requirements.For example, let the internodal distance and the maximum error beh=1 and E≤10-4. Consider an initial condition whose maximum significative wavenumber kc 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 N=10. Figure 5 refers to the dispersion errors, but any other error could be used as well following the same procedure.If the discretization lengthh 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 h.


### 3.2. Second Derivative

Figure6 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 k. The latter have been computed by solving (18)–(20) with E=10-4 in order to determine the optimal weights dj. Compared to the standard FD formulas, there is a significant increase in wavenumbers for which the error is kept below E=10-4.Figure 6
Errors in the approximation to the second derivative (16) using conventional stencils with N+1 nodes as a function of the wavenumber k. Dashed lines: standard FD formulas (from left to right N=2,4,6,8,10, and 12). Solid lines: optimal FD formulas for E=10-4 (from left to right N=2,4,6,8,10, and 12).The top rows of Table3 compare the maximum wavenumber kc of the standard and optimal FD formulas for the second derivative in the case E=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 kc 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 N. This advantage becomes, however, negligible for large values of N.Table 3
Optimized weights for second derivative andE=10-4 using (18)–(20).
N = 2 N = 4 N = 6 N = 8 N = 10 N = 12 k c standard 0.1855 0.4573 0.7039 0.9073 1.0736 1.2110 k c opt. (see (20)) 0.1855 0.7801 1.2294 1.5813 1.8449 2.0442 k c opt. (see (21)) 0.3120 0.8057 1.2500 1.5954 1.8550 2.0508  d 0 −2 −2.556451 −2.824639 −2.978478 −3.071991 −3.131331 d 1 1 1.371592 1.579082 1.708933 1.792168 1.846830 d 2  −0.093366 −0.184639 −0.261022 −0.319217 −0.361654 d 3   0.017876 0.046789 0.077604 0.104965 d 4    −0.005460 −0.016873 −0.030957 d 5     0.002313 0.007709 d 6      −0.001228Table3 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 N=12 the weights computed in [6] result in kc=2.017, while using the weights in Table 3 results in kc=2.0442. Similar improvements of around 1% are obtained for other values of N.


## 3.1. First Derivative

Figure1 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 N+1 node stencils with N equal to 4, 6, 8, 10, and 12. The optimal weights have been obtained by solving equations (12) with E=10-4. Figure 2 enlarges the interval around ϵ=0. Note that the spectral improvement when optimal weights are used is more significant for large N. For example, for N=4 the maximum wavenumber for which the dispersion error is smaller than E=10-4 changes from kc=0.31 with standard weights to kc=0.54 with optimal weights, while for N=12 it changes from kc=1.05 to kc=1.78.Figure 1
Dispersion errorϵ(k) for conventional stencils with N+1 nodes (from left to right N=2,4,6,8,10, and 12). Dashed lines: standard FD formulas. Solid lines: optimal FD formulas for E=10-4.Figure 2
Detail of Figure1. Dashed lines: standard FD formulas. Solid lines: optimal FD formulas for E=10-4. From left to right N=2,4,6,8,10, and 12.We observe in Figure2 that, for fixed N, the locations of k where the errors are ±E approximately coincide with the locations of k where the errors are ∓E for N+2. This observation, which happens to be true for N 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 Table1 also show the optimal weights dj for N=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 kc for which the dispersion errors are smaller than E. For instance, for N=8, the weights computed in [6] result in kc=1.3009, while the weights in Table 1 result in kc=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 (N=26 and 28), as we show in Figure 3.Table 1
Two first rows: maximum wavenumber for which the dispersion error is less than or equal toE=10-4. Bottom rows: corresponding optimized weights dj.
N = 4 N = 6 N = 8 N = 10 N = 12 k c standard 0.31309 0.5475 0.74763 0.91369 1.052 k c optimal 0.54129 0.96452 1.3112 1.5812 1.7888  d 1 0.678935 0.779124 0.842632 0.884271 0.912322 d 2 −0.089700 −0.174992 −0.246700 −0.302511 −0.3444207 d 3  0.023744 0.061634 0.102893 0.140765 d 4   −0.008624 −0.026887 −0.050463 d 5    0.0040000 0.013796 d 6     −0.002218Figure 3
Dispersion errorϵ(k) for optimized FD formulas with N+1 conventional nodes (N=26 and 28).Table2 summarizes, for all the cases analyzed in this paper, the values of the largest wavenumbers kc for which the corresponding error is bounded by E=10-4. For each of these cases we compare the values of kc using the standard FD weights (rows with s in the second column) with those obtained using the optimal weights (rows with o in the second column). The different values of kc shown in Table 2 have been obtained by solving the systems of nonlinear equations in Appendix A.Table 2
Largest wavenumberskc for which the corresponding error is bound by E=10-4. These values have been obtained by solving the nonlinear systems shown in Appendix A for the optimal weights.
N = 2 N = 4 N = 6 N = 8 N = 10 N = 12 Conventional nodes  ϵ s 0.0843 0.3136 0.5480 0.7481 0.9141 1.0525 ϵ o 0.1338 0.5415 0.9651 1.3121 1.5814 1.7465 E p s 0.0236 0.2338 0.4944 0.7204 0.9051 1.0571 E p o 0.0332 0.3919 0.8624 1.2582 1.5622 1.7922 E g s  0.1561 0.3565 0.5441 0.7061 0.8445 E g o  0.2621 0.6236 0.9589 1.2350 1.4586  Staggered nodes  ϵ s 0.1339 0.4644 0.7654 1.0022 1.1889 1.3358 ϵ o 0.2125 0.8011 1.3396 1.7337 2.0120 2.2114 E p s 0.0490 0.38301 0.7313 1.0025 1.2112 1.3724 E p o 0.0693 0.6404 1.2648 1.7224 2.0368 2.2551 E g s 0.0283 0.2559 0.5269 0.7563 0.9427 1.0928 E g o 0.0400 0.4282 0.9155 1.3138 1.6136 1.8385Figure4 shows the optimal weights for a stencil with N+1 nodes (N=4,6,8,10, and 12) on a conventional grid as a function of the error bound E.Figure 4
Optimal weightsdj as a function of E, for which |ϵ|<E. The different curves correspond to N=4,6,8,10,12.In each of the plots depicted in this figure, the optimal weightsdj increase in absolute value as N increases. As expected, the optimal weights approach the values corresponding to standard FD formulas when the error E 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 boundE. For example, for N=4, the optimal weights are given to leading order by(22)d1≈23+0.2193E3/8+⋯,d2≈-112-0.1096E3/8+⋯.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 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 bandwidthkc determines the number of grid points per wavelength(23)PPW=2πkc=2πkc′h=λc′hrequired to resolve the shortest wavelength component λc′ for which the corresponding error is bounded by a given threshold E. The values of PPW as a function of E 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 N+1. From top to bottom, the curves depicted in this figure correspond to optimized FD formulas with N=2,4,6,8,10, and 12. For each of these curves, the points per wavelength are greatly reduced as the required E 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
Points per wavelength as a function of the maximum dispersion error allowedE for N=2 (top curve) to N=12 (bottom curve).Figure5 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 E 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 h is given one should perform the following steps:(1)
Choose the maximum admissible errorE.(2)
Compute the Fourier transform of the initial condition.(3)
Pick the maximum significative wavenumberkc to be resolved below the error E.(4)
The corresponding shortest wavelengthλc gives the number of points per wavelength PPW through (23).(5)
Using Figure5, PPW and E define the number of nodes N+1 of the FD formula that fulfills both requirements.For example, let the internodal distance and the maximum error beh=1 and E≤10-4. Consider an initial condition whose maximum significative wavenumber kc 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 N=10. Figure 5 refers to the dispersion errors, but any other error could be used as well following the same procedure.If the discretization lengthh 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 h.


## 3.2. Second Derivative

Figure6 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 k. The latter have been computed by solving (18)–(20) with E=10-4 in order to determine the optimal weights dj. Compared to the standard FD formulas, there is a significant increase in wavenumbers for which the error is kept below E=10-4.Figure 6
Errors in the approximation to the second derivative (16) using conventional stencils with N+1 nodes as a function of the wavenumber k. Dashed lines: standard FD formulas (from left to right N=2,4,6,8,10, and 12). Solid lines: optimal FD formulas for E=10-4 (from left to right N=2,4,6,8,10, and 12).The top rows of Table3 compare the maximum wavenumber kc of the standard and optimal FD formulas for the second derivative in the case E=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 kc 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 N. This advantage becomes, however, negligible for large values of N.Table 3
Optimized weights for second derivative andE=10-4 using (18)–(20).
N = 2 N = 4 N = 6 N = 8 N = 10 N = 12 k c standard 0.1855 0.4573 0.7039 0.9073 1.0736 1.2110 k c opt. (see (20)) 0.1855 0.7801 1.2294 1.5813 1.8449 2.0442 k c opt. (see (21)) 0.3120 0.8057 1.2500 1.5954 1.8550 2.0508  d 0 −2 −2.556451 −2.824639 −2.978478 −3.071991 −3.131331 d 1 1 1.371592 1.579082 1.708933 1.792168 1.846830 d 2  −0.093366 −0.184639 −0.261022 −0.319217 −0.361654 d 3   0.017876 0.046789 0.077604 0.104965 d 4    −0.005460 −0.016873 −0.030957 d 5     0.002313 0.007709 d 6      −0.001228Table3 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 N=12 the weights computed in [6] result in kc=2.017, while using the weights in Table 3 results in kc=2.0442. Similar improvements of around 1% are obtained for other values of N.


## 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 c=1 and different initial conditions. In all the cases considered here, we use conventional stencils with step size h=1, and periodic boundary conditions. Time integration is carried out using a fourth-order Runge–Kutta method with time step Δt=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 Δt=0.01 obtaining identical results.Let us start by considering the propagation of a pure wave with wavenumberk=π/2. Thus, we solve (1) with initial condition(24)ux,0=sin⁡π2x.Figure7 shows the exact solution at t=20 (solid line) together with the numerical solution obtained using optimized FD formulas of different support and maximum allowed error E. Notice that the exact solution is the initial wave (24) centered at x=20, and therefore u(20)=0. The squares correspond to the solution obtained with the optimal weights corresponding to N=12 nodes conventional stencil with |Ep|<E≡10-2. Notice that u(20)=-0.2771 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 Ep(π/2)=0.008938 and, therefore, at t=20 the wave has travelled a distance D=20(1+0.008938)=20.17877. This value exactly coincides with the displacement observed in the numerical solution D=20+2arcsin⁡(0.2771)/π=20.17877. The circles correspond to the solution obtained with the optimal weights corresponding to N=10 nodes conventional stencil with |Ep|<E≡10-2. In this case, u(20)=0.2058 (large circle) and, therefore, the propagation velocity is smaller than one. In fact, Ep(π/2)=-0.006597 and D=19.86806, which again coincides with the numerical solution D=20+2arcsin⁡(-0.2058)/π=19.8680. These results are shown in Table 4.Table 4 Weights u ( 20 ) D N=12 −0.2771 20.17877 N=10 0.2058 19.8680 Combined −0.03677 20.0234Table 5
Conventional nodes,ϵ.
N d j E 0 E 0.4 E 0.8 E 1.2 E 1.6 E 2 4 d 1 2/3 0.4924 −0.1558 −0.0718 0.0661 −0.0119 d 2 −1/12 −0.2462 −0.2857 −0.1559 −0.0347 0.0093 N d j E 0 E 0.3 E 0.6 E 0.9 E 1.2 E 1.5  6 d 1 3/4 0.5048 −0.7967 2.0934 −3.7087 2.5988 d 2 −3/20 −0.4059 0.2154 −1.1857 2.5889 −1.8717 d 3 1/60 0.1022 0.1298 0.4946 −0.3993 0.4132 N d j E 0 E 0.2 E 0.4 E 0.6 E 0.8 E 1  8 d 1 4/5 0.2468 0.2738 −1.0570 1.3929 −0.7130 d 2 −1/5 −0.2442 −0.4511 1.0109 −1.1471 0.5844 d 3 4/105 0.1037 0.2999 −0.1453 0.2508 −0.2520 d 4 −1/280 −0.0176 −0.0642 −0.1516 −0.0018 −0.3805 N d j E 0 E 0.2 E 0.4 E 0.6 E 0.8 E 1  10 d 1 5/6 0.4139 −0.7854 1.5663 −1.9306 0.9708 d 2 −5/21 −0.4790 0.6220 −1.2993 1.7253 −0.8936 d 3 5/84 0.2738 −0.0557 0.4723 −0.8206 0.4281 d 4 −5/504 −0.0831 −0.1096 −0.3179 0.3435 −0.0613 d 5 1/1260 0.0111 0.0270 0.1898 −0.0605 0.4262Figure 7
Exact solution att=20 of the 1D advection equation (solid line) along with the solution obtained with optimal weights (|Ep|<E≡10-2) for N=12 conventional nodes (squares), optimal weights (|Ep|<E≡10-2) for N=10 conventional nodes (circles), and alternating N=10 and N=12 optimal weights (triangles).An interesting possibility is to efficiently combine two stencils of consecutive sizes. In fact, it was pointed out in Section3 that the errors for two consecutive stencil sizes tend to cancel each other: at locations where the error is +E using a stencil of size N, it is approximately -E using a stencil of size N+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 N nodes at time step tn, a stencil with N+2 nodes at tn+1, a stencil with N nodes at tn+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 N=10 conventional stencil, and even time steps with a N=12 conventional stencil. Notice that the results obtained with thiscombined method are much more accurate (u(20)=-0.03677) and have the same computational cost.We consider next the initial condition(25)ux,0=sin⁡2πx-50aexp⁡-log⁡2x-50b2,where a is the dominant wavelength and b is the half-width of the Gaussian function. This is the same initial condition considered in [10]. We consider the case a=4, b=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 k=π/2.Figure 8
(a) Initial condition (25) for a=4, b=9. (b) Normalized spatial power spectrum.
(a) (b)Figure9 displays the numerical solution at t=200. Figure 9(a) shows the numerical solution (solid line) obtained using the standard weights corresponding to a stencil with N=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 k=π/2, the error of the group velocity using standard N=12 nodes stencil is Eg=-0.0693. Thus, the group velocity is vg=1-0.0693=0.9307 and the location of the maximum at t=200 is x=50+200vg=236.1, which is in very good agreement with the location of the maximum shown in Figure 9(a).Figure 9
(a)u(x,200) using conventional N=12 nodes stencil weights. (b) u(x,200) using optimal N=12 nodes stencil weights corresponding to Eg and E=10-3. Dashed line: exact solution.
(a) (b)Notice in Table2 that for N=12 stencil with optimized weights, the maximum wavenumber for which |Eg|≤10-4 is kc=1.4586. Since the spectrum of the initial condition (25) peaks at k=π/2, we can not use in this case those optimized weights. Instead, we use the optimized weights corresponding to |Eg|≤10-3, for which kc=1.7229>π/2. These optimized weights can be obtained using the formulas in Table 7. The numerical results at t=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(26)enum=∑u-ue2∑ue21/2,where u is the numerical solution and ue is the exact solution. For the solution obtained with optimized weights (Figure 9(b)), enum=0.0216 and u-ue∞=0.0166.Finally, we consider the initial condition defined in (25) with a=8, b=9. Figure 10(a) displays the corresponding normalized power spectrum. In this case, the normalized power spectrum peaks at k=π/4, and the largest significant wavenumbers are around k=1.25. We have integrated (1) with this initial condition until t=200 using optimized weights corresponding to |Eg|≤10-3 and sizes N=10 and N=12. In both cases kc is well above the wavenumbers contained in the initial condition. The errors in the case N=10 are(27)enum=0.01586,u-ue∞=0.01124,and in the case N=12(28)enum=0.01403,u-ue∞=0.01092.However, if we combine both stencils by integrating odd step times with the n=10 stencil, and the even time steps with the n=12 stencil, the results are much more accurate: (29)enum=0.00463,u-ue∞=0.000335.The numerical solution for t=200 using thecombined method is shown in Figure 10(b) (solid line) together with the exact solution (dashed line).Figure 10
(a) Normalized spatial power spectrum for (25) with a=8, b=9. (b) u(x,200) usingcombined method of N=10 and N=12 stencils corresponding to Eg and E=10-3. Dashed line: exact solution.
(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.



---

*Source: 7860618-2016-12-13.xml*