Comparison of Finite Difference Schemes for the Wave Equation Based on Dispersion

Finite difference techniques are widely used for the numerical simulation of time-dependent partial differential equations. In order to get better accuracy at low computational cost, researchers have attempted to develop higher order methods by improving other lower order methods. However, these types of methods usually suffer from a high degree of numerical dispersion. In this paper, we review three higher order finite difference methods, higher order compact (HOC), compact Padé based (CPD) and non-compact Padé based (NCPD) schemes for the acoustic wave equation. We present the stability analysis of the three schemes and derive dispersion characteristics for these schemes. The effects of Courant Friedrichs Lewy (CFL) number, propagation angle and number of cells per wavelength on dispersion are studied.


Introduction
Real life situations in many disciplines including engineering, physics, economics, biosciences, etc, can be described through mathematical models.Differential equations play an important role in modelling interactions in physics, engineering, and biological processes, from ground motion (earthquake) to interactions between neurons.Differential equations can either be categorized as ordinary differential equations (ODEs), which are for example, used to model dynamical systems and their applications to science and engineering or partial differential equations (PDEs), which most importantly used to describe a wide variety of phenomena in physics; such as sound, heat, fluid dynamics, elasticity and in most field of engineering in general.However, most differential equations, such used to solve real-life problems, have no analytical solutions or are unrealistic to solve analytically.Therefore, solutions for these problems are attempted by means of numerical approximations.To this end various numerical methods are developed to provide better approximate solutions to these equations.
Numerical techniques are a powerful tool for handling large systems of equations involving complex geometries that are otherwise impossible to handle analytically [1].There are several types of numerical methods, each with their own pros and cons.Finite difference method (FDM) is the most popular numerical technique which is used to approximate solutions to differential equations using finite difference equations [2].These techniques are widely used for the numerical solutions of time-dependent partial differential equations.However, when we approximate derivatives using finite difference methods truncation of terms produces inaccuracies.Due to these truncations an inherent problem arises, undesirable ripples are produced due in the process of discretization.Sometimes, due to truncation of higher order terms, frequency dependent numerical errors called dispersion is produced.A large amount of work has been developed for higher order accurate scheme with low dispersion; for instance, Krum [3] developed time-domain schemes based on multiresolution analysis, Lee [4] introduced a low-dispersive schemes based on Pseudo-spectral approach.
The aim of this research is to explore a higher order compact finite difference scheme with comparably low dispersive characteristics.Therefore, in this paper, firstly, we discussed a conventional higher order finite difference schemes.Then we discussed another higher order finite difference scheme developed by Das [5] and Liao [6].We also discussed the dispersion properties of these methods by comparing with a conventional higher order finite difference scheme.The paper explores comparably low dispersive scheme with among the finite difference schemes.

Higher Order Finite Difference Discretization for the Wave Equation
The two dimensional version of the wave equation with velocity v and acoustic pressure u in homogeneous media can be written as , , 0, , along with the initial conditions and the boundary condition

Finite Difference Operators
To discretize the above PDE we consider a uniform grid equally spaced in all spacial direction with a grid spacing .
A higher order finite difference approximation for the second derivative is given by Liu [8]. , where 0 a and m a are central difference coefficients.Some values of these coefficients are presented in Table 1, and the formula to find all coefficients can be found in Liu [8].
Though this approximation is accurate up to order 2M, its not efficient because when M becomes large, the stencil that need to be computed becomes large.
The standard second-order finite difference approximation for the acoustic wave equation in (1) using the above operators is ( ) , ( )

Implicit Finite Difference Method
A fourth order accurate implicit finite difference scheme for one dimensional wave equation is presented by Smith [9].We extend the idea for two-dimensional case as discussed below.
Consider two dimensional wave equation, using Taylor's series expansion of ( )

, , u t h x y
+ and ( ) , , u t h x y − about the point ( ) If u is a solution of (1), then we have the following Substituting equations ( 13) and ( 14) in to equation (12) and using the notations ( )  , , , Following Strikwerda [10] and using Taylor's series expansion we obtain and We also have 1 .
Substituting ( 16)-( 19) into (15) gives a scheme which is fourth order accurate where v r h τ = as denoted above.
Applying the operator to both sides of (20), we obtain Equation ( 21) can be written as an ADI-scheme and can be referred to as HOC-ADI, which can be presented more precisely as, 12 We have introduced the higher order finite difference alternating direction implicit scheme for two-dimension-al wave equation.In the next section we will derive another higher order scheme based on Padé approximations.

Higher Order Finite Difference Scheme Based on Padé Approximation from Taylor's Series Expansion We Have
( ) 2 1 12 1 .
The scheme presented in (30) is a 4 th -order accurate both in time and space for the 2-dimensional acoustic wave equation based on Padé approximation.

Higher Order Compact Finite-Difference Method for the Wave Equation
A compact finite difference scheme comprises of adjacent point stencils of which differences are taken at the middle node, therefore typically 3, 9 and 27 nodes are used for compact finite difference descretization in one, two and three dimensions, respectively.Whereas a non-compact stencil comprises of any number of nodes which are near the vicinity of the middle node.Usually, when the order of accuracy of a non compact scheme is greater than two, then the scheme is said to be higher order.Figure 1 represents a compact and non-compact stencil in two-dimensions.

Higher Order Compact Alternating Direction Scheme for the Wave Equation
The standard alternating direction implicit scheme which is fourth order both in time and space is given in equation ( 22)-(23).Recently, Das [5] developed a low dispersive finite difference method for two-dimensional a coustic wave equation using Padé approximation, we will discuss this method in detail.
Multiplying both sides of (30) by the operator which upon simplification leads to . 12 12 144 6 Further simplification give Neglecting the term This scheme needs to be solved in two directions simultaneously, however it is easier to solve it in one direction at a time, therefore we can write it as x y n . 12 The above scheme is based on Pad'e approximation and which can be referred to as CPD-ADI scheme.Substituting ( 16) into (36) we obtain ( ) Equation ( 37) is a tridiagonal system and therefore it can be solved easily by using any tri-diagonal solver.
In the next section we will derive a numerical scheme which is a combination of the compact and non-compact operators that we introduced in chapter two.

Hybrid Scheme
The higher order compact scheme (35)-(36) derived using Padé approximation is efficient and accurate up to order 4.However, the scheme suffers from a moderate numerical dispersion [5].Therefore, we modify the scheme by substituting one of the spacial grids (in x-direction) in Equation (30) by non-compact one in Equation ( 8) to obtain a hybrid scheme.
Therefore, for the two-dimensional wave equation we have 1 .
1 1 12 12 Following the same procedure as we did for (32)-(34) we obtain ( ) Therefore, for implementation, this can further be written as an ADI scheme x y n Equations ( 40)-( 41) is an ADI scheme obtained by mixing compact and non-compact difference operators and therefore it can be referred to as non-compact padé based scheme (NCPD-ADI).
For the implementation, we will apply Equation ( 40)-(41) for the interior nodes, whereas at the nodes near the horizontal boundary we will use (35)-(36).Therefore, the combination of the two schemes can be called compact Padé based interlinked alternating direction implicit scheme (IPD-ADI).
At this stage, we would like to mention that a typical numerical method may produce spurious solution, a phenomena which is greatly attributed to dispersion.Thus, in the next section we will study this phenomena in detail.

Dispersion Relation and Stability Analysis
Finite difference schemes are widely used in solving time dependent partial differential equations numerically, however in approximation of derivatives using Taylor's series, truncation of terms produces inaccuracies in finding solutions for such differential equations.Due to these inaccuracies undesirable ripples are produced in the process of discretization.Eventhough such equations are not frequency dependent, due to truncation of higher order terms frequency dependent numerical errors are produced.This numerical phenomena is called dispersion.Thus, even for non-dispersive partial differential equations their discrete approximations are dispersive.Moreover in descritizing a differential equation, depending on our interest in accuracy, we ignore higher order terms, due to these inaccuracies errors are introduced which results in instability of the finite difference scheme.

Dispersion Relation
We will consider a scalar wave Equation ( 1 where k v ω = ⋅ is the analytical dispersion relation and we will derive the numerical dispersion relation for two-dimensional wave equation for the mentioned schemes.
In numerical analysis the relation between the numerical angular frequency ω and k, v and numerical parameters τ and h is called numerical dispersion relation.
Dispersion error therefore can be considered as the difference between the numerical and analytical frequencies or the ratio of these frequencies.Consequently, it can be measured as the difference between, or the ratio of numerical and analytical velocities.Since the latter method is the most convenient to work with, we will use this method for our dispersion measurement.Thus, mathematically, where v is numerical velocity and δ is the normalized phase error.Therefore 1 δ = imply the numerical scheme that we are evaluating is non-dispersive, where as for values largely deviating from 1 indicates the scheme is dispersive.
Applying finite difference operators to (43) gives as , 4sin , , 2 , 4sin , . 2 For two-dimensional HOC-ADI scheme the phase relation is derived as follows.Substituting (43) into (21) we obtain Similarly the normalized phase relation for CPD and NCPD schemes in two dimensions respectively are given by where ( ) The plots for these normalized dispersion relations are given in the subsequent sections.

Stability Analysis
A numerical scheme is said to be stable if the errors produced, due to descritization, at one time-step of the calculation do not propagate as computation continues.For a stable scheme, the errors decay and eventually damp out, whereas for the unstable scheme the errors grow with time.
In this section we will analyse stability for HOC, CPD and NCPD schemes using Von Neumann stability analysis.
Therefore, the CFL condition for the HOC scheme is derived as follows Equation ( 43) can be written as We consider n z ζ = , then (57) can be written as  To examine the effect of number of cells per wavelength, we plot the dispersion errors versus θ for different values of number of cells per wavelength for the three schemes. In

Comparison of Phase Error
In this section we presented plots of the 2 L -norm and relative phase error for these three schemes.Presented in Figure 11 is a comparison of the plots for 2 L -norm for HOC-ADI, CPD-ADI and NCPD-ADI for kh ranges from 0 to π.It is clearly observed in the plot that all the schemes are periodic in θ with a period of π.It is seen in the plot that CPD and NCPD schemes have least error than HOC scheme and it is also observed that NCPD scheme has slightly smaller error than CPD scheme.
The relative phase error is calculated with the following formula Relative error =  the reference value and meas P is the numerically calculated value.In Figure 12 we presented the relative phase errors of the three schemes.As can be seen on the plot, NCPD-ADI scheme has a lower dispersion errors than CPD scheme, and CPD scheme has significant lower dispersion than HOC scheme.
We have discussed the relative phase errors for different schemes in the individual plots above.As we have discussed in the earlier section, phase error is dependent of three factors.Therefore, it is better to see the relative phase error for the combined effect of these factors.The plots in Figures 13-15 show two-dimensional relative  phase errors, where kh represents polar radius, θ is the polar angle and the plot shows for four values r.
From Figures 13-15 we clearly observe that smaller r values show lower dispersion.It is also shown that waves whose propagation direction makes an angle of 45 degree with the axis are the most accurate for CPD and NCPD schemes, whereas for HOC scheme waves propagating in the axial direction shows relatively low error  for some values of r.Moreover, as clearly indicated, CPD and NCPD schemes show large low dispersion region than HOC scheme.

Conclusion
In this thesis, three ADI schemes are presented for 2D and 3D acoustic wave equation with constant velocity.We derived HOC-ADI, CPD-ADI and NCPD-ADI schemes.The former scheme is based on Taylor approximation and the latter schemes are based on Padé approximation.We perform stability analysis for the mentioned   schemes.Dispersion relations for these schemes have been also obtained followed by a thorough discussion of these relations by plotting them for different parameter sets with respect to different values of kh and θ.From this analysis, it is observed that an increase in these parameters increases the dispersion error and vice-versa.It is observed that numerical dispersion is dependent of these parameters in general for the mentioned schemes.Comparison of phase errors for these schemes is also carried out with the given parameters.Phase error comparison of the three schemes showed that NCPD scheme has relatively lower dispersion error than CPD and HOC schemes.However, in this work we could only compare dispersion properties for three higher order numerical schemes, therefore in the future, we will analyse dispersion properties of a large group of schemes so that we will be able to explore higher order accurate low dispersive schemes.
) which has solution of the form grid spacing, τ is time-step and discrete directional wave numbers are given by

Figures 8 -
10 we presented the normalized phase error versus propagation angle θ for different values of kh.It is observed that the phase error increases as the propagation angle increases for all the schemes.Though all the schemes show lower dispersion error for π 8 kh = ; moreover, the CPD and NCPD schemes are relatively low dispersive for π 4 kh = than HOC scheme.

Figure 8 .
Figure 8. Normalized phase error for different values of kh for HOC-ADI scheme, r = 0.6.

Figure 9 .
Figure 9. Normalized phase error for different values of kh for CPD-ADI scheme, r = 0.6.

Figure 10 .
Figure 10.Normalized phase error for different values of kh for NCPD-ADI scheme, r = 0.6.

Figure 12 .
Figure 12.Relative phase error for HOC, CPD and NCPD schemes.

Figure 13 .
Figure 13.Relative phase error for HOC-ADI scheme for different values of r.

Figure 14 .
Figure 14.Relative phase error for CPD-ADI scheme for different values of r.

Figure 15 .
Figure 15.Relative phase error for NCPD-ADI scheme for different values of r.

Table 1 .
Coefficients of central difference formula.