## Abstract

In this paper, a symplectic method for structure-preserving modelling of the damped acoustic wave equation is introduced. The equation is traditionally solved using non-symplectic schemes. However, these schemes corrupt some intrinsic properties of the equation such as the conservation of both precision and the damping property in long-term calculations. In the method presented, an explicit second-order symplectic scheme is used for the time discretization, whereas physical space is discretized by the discrete singular convolution differentiator. The performance of the proposed scheme has been tested and verified using numerical simulations of the attenuating scalar seismic-wave equation. Scalar seismic wave-field modelling experiments on a heterogeneous medium with both damping and high-parameter contrasts demonstrate the superior performance of the approach presented for suppression of numerical dispersion. Long-term computational experiments display the remarkable capability of the approach presented for long-time simulations of damped acoustic wave equations. Promising numerical results suggest that the approach is suitable for high-precision and long-time numerical simulations of wave equations with damping terms, as it has a structure-preserving property for the damping term.

## 1. Introduction

The study of the long-term modelling of wave motion forms the basis of the solutions to many physical problems in the realm of wave propagation, such as noise analysis, modelling of the Earth's free oscillations, seismic noise propagation modelling and long-time modelling of seismic waves. These studies usually require damped wave equations with variable coefficients to be solved; these equations need structure-preserving numerical procedures to avoid accumulated errors in precise or long-time numerical simulations. In the following, taking high-precision and long-term modelling of scalar seismic waves in attenuating media, for example, we discuss the question of structure-preserving schemes for long-term modelling of damped acoustic wave equations with variable coefficients.

Like modelling of undamped wave equations in the time domain using direct methods, modelling of damped wave equations also involves discretization of both spatial and temporal derivatives. In this paper, the emphasis is placed on the structure-preserving time discretization. In recent decades, classical finite difference (FD) methods for temporal discretizations have been widely applied. Because the classical FD approaches for time discretizations are not structure-preserving schemes, it is quite hard to avoid accumulated errors in precise or long-time numerical simulations for wave equations. As for some numerical solutions to differential equations, the corresponding structures (such as symmetry, positive definiteness, conservation, etc.) can be preserved by using some numerical methods; this is called the structure-preserving property of a numerical scheme. In theory, symplectic schemes have a structure-preserving property. A numerical algorithm for Hamiltonian dynamical systems can be termed a symplectic scheme if the resulting solution is also a symplectic mapping. Some symplectic schemes for partial differential equations have been proposed and applied, such as Lax–Wendroff algorithms [1,2] and Nyström algorithms [3–9]. Chen [10] minutely discussed the structure-preserving property of these algorithms. These symplectic algorithms are only suitable for modelling of undamped wave equations in the case of long-term computation. For the long-time modelling of damped wave equations, therefore, new structure-preserving algorithms are required. This sort of method should have a structure-preserving property for damping terms of damped wave equations.

In this paper, we extend the structure-preserving calculation into the problem of damped waves and present a structure-preserving method for accurately and efficiently modelling damped acoustic wave equations with variable coefficients. In the method presented, a truncated and optimized discrete singular convolution differentiator method (DSC) [11,12] is used for space discretizations. In order to improve the performance of wave-field modelling methods for long-time simulations, we substitute the second-order symplectic scheme for the second-order FD scheme in time discretizations. For temporal discretization, the scheme presented is a second-order operator, which requires the same amount of computational time as the second-order FD time discretization.

As an example, we apply the method presented to scalar seismic wave-field modelling in attenuating heterogeneous media. Our numerical results indicate that the method presented is suitable for large-scale numerical modelling because it effectively suppresses numerical dispersion by discretizing the damped acoustic wave equation when coarse grids are used. The numerical results also confirm that the method presented in this paper is superior in solving long-time simulation problems of damped acoustic wave equations.

## 2. Theoretical method and model

As mentioned above, the major goal of this work is to produce a structure-preserving scheme to suppress accumulated errors for long-time modelling of acoustic waves with damping (or to extend the structure-preserving calculation into modelling of wave equations with dissipation terms (damping terms)) and not just to study the problem of the damped acoustic wave itself. Taking seismology for example, here, we present a symplectic scheme for high-precision and long-term modelling of damped seismic waves. Sometimes, the scalar wave equation (acoustic wave equation) can be used to model the seismic-wave field. For mathematical convenience, the scalar wave equation with damping terms (the scalar wave equation with the D'Alembert damping model) for two-dimensional heterogeneous media in the time domain can be written as
*p* is the pressure field, *v* is the velocity of the wave, *D* is the damping coefficient, *f* is the body force, *x* and *y* are Cartesian coordinates and *t* is the time. Here *D*=−*ωQ*^{−1}≤0.0 s^{−1}, where *ω* is the dominant frequency of the wave field and *Q* is the quality factor of the medium.

As in previous work [12,13], we can express equation (2.1) as Hamiltonian PDEs. Writing *v*^{2}((∂^{2}/∂*x*^{2})+(∂^{2}/∂*z*^{2}))=*L* and d*p*/d*t*=*q*, then equation (2.1) can be written as
*D*(d*p*/d*t*)), it is not a typical Hamiltonian system.

For the time variable *t*, the solution of equation (2.3) can be written as
*c*_{i},*d*_{i},*i*=1,⋯ ,*k* are real numbers. Although equation (2.2) is not a typical Hamiltonian system, *k*th-stage symplectic scheme with *m*th-order accuracy of equation (2.4) can be obtained,
*c*_{i} and *d*_{i} are called symplecticity coefficients.

By means of Taylor expansion and the time discretization, the discrete solution of equation (2.6) can be expressed as

where *n* is the index of the discrete time point, Δ*t* is the discrete time step and *i* is the stage index of scheme (2.6). Because the time error of equation (2.5) is *O*(*t*^{m+1}), the time accuracy of equation (2.6) is of *m*th-order. When *k*=*m*=2, equation (2.6) is a second-order operator, that is, it is

When calculating spatial derivatives of equation (2.1), the following discrete singular convolutional differentiator [11,12] is adopted:
*x* is the mesh size, 2*W*+1 is the computational bandwidth and *w*(*i*) is a Hanning window function for avoiding the Gibbs phenomenon:
*α* (0.5≤*α*≤1) and *β* allow a family of different windows to be considered. To obtain an optimal balance between computational efficiency and accuracy of the discrete singular convolutional approach, we chose nine-point explicit operators on regular grids via the discrete Fourier analysis. The nine-point explicit operator used in this paper is accurate to eighth order in space.

In the convolutional differentiator method, the spatial derivatives of *p* in equation (2.1) can be written as
*x* and *L*(*p*^{(n)}) in equation (2.9) can be written as
*m* and *n* are indices along the discrete *x* and *z* axes, Δ*x*, Δ*z* and Δ*t* are sampling rates along the *x*-, *z*- and *t*-axes, and *mx* and *nz* are the half differentiator lengths in the sampling number along the *x*- and *z*-axes. Similarly, equation (2.9) can be expressed as

In subsequent sections, the performance of scheme (2.15) will be evaluated through a set of numerical experiments which compare the accuracy, numerical dispersion, stability and computation cost of scheme (2.15) with a high-order FD scheme. The high-order FD scheme can be written as
** L** is an operator of the second-order partial derivate for two-dimensional spatial discretization. In this paper,

**is an eighth-order FD operator. Scheme (2.16) is a non-symplectic scheme which does not have a structure-preserving property.**

*L*## 3. Stability condition

Because it is difficult to obtain an explicit treatment for the stability of scheme (2.15), we here evaluate the stability using the Courant number. The following is a numerical calculation for the Courant number. In the numerical calculation, the model parameters are the wave velocity of *V* =3000 m s^{−1} and the damping coefficient of *D*=−0.25. The number of grid points is 301×301; the model size is 3000×3000 m. The mesh size is 10 m. An explosive source is located at (*x*_{s}, *z*_{s})=(2550 m,2550 m), which is a Ricker wavelet and can be written as
*f*_{0} is the dominant frequency and *t*_{0} is the starting time. This is a medium model without any type of absorbing or transmitted boundary conditions. Generally, the Courant number can be defined as
*V*_{max} is the maximum value of the wave velocity, Δ*h*=Δ*x*=Δ*z*. For the above homogeneous medium, the stability condition or the stability limit of equation (3.3) is Δ*t*≤0.6039Δ*h*/*V* or the maximum Courant number *r*_{max}=Δ*tV*/Δ*h*≤0.6039.

The maximum Courant number is the most common numerical stability condition. In the above-mentioned calculation of the maximum Courant number, the model parameters are the wave velocity of *V* =3000 m s^{−1}, a mesh size of Δ*h*=10 m and a maximum time increment of Δ*t*=2.013 ms. Theoretically, scheme (2.15) is equivalent to an optimized FD method.

For scheme (2.16), similarly, the stability condition can be obtained in the same numerical model, as follows.

The maximum time increment of Δ*t*=1.845 ms, Δ*t*≤0.5539Δ*h*/*V* or the maximum Courant number *r*_{max}=Δ*tV*/Δ*h*≤0.5539.

From comparison of the above two stability conditions, it can be found that the stability of the scheme is superior.

## 4. Numerical experiments, results and discussion

In general, the performance of numerical schemes is evaluated by considering the numerical dispersion as a function of the number of grid points per wavelength. However, it is difficult to evaluate the performance of scheme (2.15) using this method because the damping term makes it complicated. Even though scheme (2.15) in complicated cases is usually not known analytically, the overall performance can still be judged qualitatively. In the following, two numerical examples of scalar seismic-wave modelling are given for evaluating the performance of scheme (2.15).

The high-order FD scheme is one of the most widely used methods, which is an efficient scheme for treating heterogeneous media. We compared the numerical results found using scheme (2.15) with those from the conventional high-order FD (scheme (2.16)) for a two-layered medium with a high-velocity contrast. The model consists of two different wave velocity regions separated by a rough inclined interface (figure 1). The model parameters were a velocity of *V*_{1}=1750 m s^{−1} for the upper layer with the source, a velocity of *V*_{2}=3000 m s^{−1} for the lower layer and a damping coefficient of *D*=−3.1416 (*Q*≈50) for both of the two layers. The number of grid points was 301×301, the model size was 3000×3000 m, and the wave source was located at (*x*_{s}, *z*_{s})=(1500 m,1300 m). The receiver was located at (*x*_{r}, *z*_{r})=(1500 m,1000 m). The mesh size was 10 m and the time increment was 1 ms. The rough interface can be considered a velocity discontinuity because the physical parameter contrast is quite high. The seismic source is approximated as a point source (e.g. equation (3.1)) that includes a band-limited Ricker wavelet; it is located in the upper layer and has an amplitude spectrum peak at 25 Hz and a high-frequency cut at 36 Hz. The Ricker wavelet can be written as equation (3.2).

Figure 2*a* is an acoustic wave-field (scalar seismic wave field) snapshot at 600 ms generated by the symplectic scheme (the nine-point operator of scheme (2.15)). The snapshot in figure 2*a* clearly displays the wavefront of the direct wave and other phases (e.g. the reflected, transmitted and scattered waves from the rough interface). In the snapshot generated by scheme (2.15), it is difficult to find any grid dispersion and distortion despite the fact that there are only 4.86 grids per shortest wavelength.

Comparing the wave-field snapshot generated by scheme (2.15) (figure 2*a*) with that generated by the conventional high-order (eighth-order) FD (figure 2*b*), one can see that there is hardly any evidence of numerical dispersion and distortion in the symplectic approach, whereas numerical dispersion and distortion are very obvious when using the eighth-order FD method. We can find a similar phenomenon when comparing the synthetic seismograms (figure 3*a* generated by scheme (2.9) and figure 3*b* generated by the eighth-order FD).

By comparing the synthetic seismograms (figure 3*a* generated by scheme (2.15) and figure 3*b* generated by the conventional high-order FD) with the synthetic seismogram of the undamped wave (*D*=0 or

To examine the long-time performance of the symplectic scheme presented, the numerical results computed by scheme (2.15) and those from the eighth-order FD scheme for a two-dimensional attenuating homogeneous medium model were compared. The model parameters are a damping coefficient of *D*=−0.25 s^{−1} (*Q*=500) and a velocity of *V* =3000 m s^{−1}. The number of grid points was 601×601, the model size was 6000×6000 m, and the point source was located at (*x*_{s}, *z*_{s})=(3000 m,3000 m). The receiver was located at (*x*_{r}, *z*_{r})=(3000 m,2000 m). The spatial step sizes were 10 m and the time step was 1 ms. The point source, a band-limited Ricker wavelet, has an amplitude spectrum peak at 20 Hz ( *f*_{0}=20 Hz).

Figure 5*a* displays a synthetic seismogram generated by the eighth-order FD scheme for the case of an undamped wave (*D*=0 or *b* and figure 5*c* show synthetic seismograms generated by scheme (2.15) and the eighth-order FD scheme for the case of an damped wave (*D*=−2.5 s^{−1} or *Q*=50), respectively. In figure 5, what we find is the same as that seen in figure 3. In comparison with the synthetic record of an undamped wave, it can be seen that the synthetic seismogram generated by scheme (2.15) does not have any distortion in either the phase or the wave form and the attenuation behaviour appears clearly, whereas serious distortion of these appears in the synthetic seismogram generated by the conventional high-order FD. Figure 6 shows the frequency spectra of the above-mentioned synthetic seismograms in the case of long-term simulation. In comparison with the frequency spectrum of undamped waves, figure 6 shows that the amplitude attenuation, frequency band narrowing and dominant frequency reducing are quite obvious in the frequency spectrum for scheme (2.15), while the result obtained by using the eighth-order FD scheme has obvious frequency band broadening and dominant frequency increasing. It is clear that the results obtained by the presented symplectic scheme are reasonable in physics.

Figure 7*a* and figure 7*b* display acoustic wave-field snapshots computed by scheme (2.15) after 500 and 10 000 time steps, respectively. Figure 7*c* and figure 7*d* show acoustic wave-field snapshots obtained by the eighth-order FD scheme after 500 and 10 000 time steps, respectively. Figure 7*a* and figure 7*c* show results for weak damping (*D*=−0.25 s^{−1} or *Q*=500). Similarly, figure 7*b* and figure 7*d* show results for strong damping (*D*=−2.5 s^{−1} or *Q*=50). In the case of weak damping, we can see that the wavefront curves generated by the two schemes after 500 time steps are very clear from figure 7*a*,*c*. This means that the two schemes have similar performance for short-time numerical simulations of the acoustic wave equation with weak damping. For long-time numerical simulations of acoustic waves with strong damping, however, the two schemes are quite different in performance and have different error growth. After 10 000 time steps, the wavefront curves generated by scheme (2.15) are still quite clear. In the result computed by the eighth-order FD scheme, however, it can be seen that the wavefront curves have blurred seriously. This comparison shows that the two schemes perform quite differently for long-time modelling of acoustic waves with strong damping, and scheme (2.15) is very appropriate for long-time simulation of wave equations with damping terms.

In order to test the performance of the non-symplectic DSC method for long-time simulations in the case of strong damping, we substitute the second-order FD scheme (scheme (2.16)) for the second-order symplectic scheme (scheme (2.15)) in time discretizations. The model parameters of the numerical model adopted here are the same as those of the previous long-time modelling in this paper. Figure 8 shows a damped acoustic wave-field snapshot computed by the non-symplectic DSC method after 10 000 time steps. This result is for strong damping (*D*=−2.5 s^{−1} or *Q*=50). From the result, it can be found that the wavefront curves have blurred seriously, and numerical dispersion and distortion are very obvious. This result looks almost identical to that of scheme (2.16), which is shown in figure 7*d*. The result also proves that the nice long-time behaviour of the symplectic method is not caused by the DSC spatial discretization.

Theoretically, a spatial discretization operator is not structure-preserving for the time discretization. Therefore, the DSC method with non-symplectic time discretization (such as the second-order FD time discretization or other high-order FD time discretization) does not have any structure-preserving property for the time discretization, although it is able to suppress numerical dispersion efficiently in short-time modelling.

Some robust methods with non-symplectic time discretization, such as RNADM [14] and the DSC method [11], can effectively suppress the numerical dispersion for short-time modelling of wave propagation. Especially RNADM can simulate high-frequency wave propagation for a given grid spacing and automatically suppresses the numerical dispersion for short-time modelling of the acoustic wave equation without the damping term. However, these methods with non-symplectic time discretization are unable to suppress the numerical dispersion for long-time modelling, especially for the case of strong damping. Existing symplectic algorithms are only suitable for modelling of undamped wave equations in the case of long-term computation. So far, no studies into long-time modelling of wave equations with damping terms have been found in the literature. Because the symplectic method presented (scheme (2.15)) has a structure-preserving property for the time discretization of acoustic wave equations with damping terms, it is able to effectively suppress the numerical dispersion for long-time modelling of acoustic wave equations with strong damping. This means that only the symplectic method presented is suitable for long-time modelling of acoustic wave equations with strong damping at present. This is an advance on existing work.

From the above-mentioned numerical results and discussion, it can also be seen that this study extends the structure-preserving calculation into the modelling of an attenuating wave in nature. This is the significant difference between this study and previous work [12].

## 5. Conclusion

Generally, the ultimate goal of a symplectic scheme for long-time modelling of attenuating wave propagation is to avoid accumulated errors in long-time numerical simulations for partial differential equations or to preserve some intrinsic properties of these equations (such as the structure-preserving property for the attenuating term of wave equations). To this end, we have extended the structure-preserving calculation into modelling of wave equations with dissipation terms (damping terms) and have presented a sympletic approach for modelling of the damped acoustic wave equation with variable coefficients. In the method presented, an explicit second-order symplectic scheme is used for the time discretization and physical space is discretized by the DSC method. The approach presented (scheme (2.15)) is accurate to second order in time and to eighth order in space. The performance of the proposed scheme has been tested and verified using numerical simulations of the attenuating acoustic wave equation (attenuating scalar seismic-wave equation). These numerical experiments demonstrate the superior capability of the approach presented for the suppression of numerical dispersion and for long-time simulations of the acoustic wave equation with strong damping terms. Our numerical results indicate that the sympletic approach presented has the structure-preserving property for the damping terms of wave equations despite the fact that the damped wave equation could not be transformed into a typical Hamiltonian system. In this study, the structure-preserving calculation is extended into damped acoustic wave-field modelling. The results here are a positive development not only for future acoustics studies but also for any physical research that requires a structure-preserving numerical solution of wave equations with attenuating terms.

## Authors' contributions

X.L. conceived and designed the mathematical models and the computational method, carried out results analysis, drafted and wrote the manuscript; M.L. and S.L. carried out numerical experiments and computational code writing; S.C. and H.Z. drafted the manuscript, carried out results analysis and participated in the numerical experiments; M.Z. participated in the numerical experiments and paper writing. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This work was supported by the National Science Foundation of China (grant nos. 41174047, 41574053).

- Received February 16, 2015.
- Accepted September 29, 2015.

- © 2015 The Author(s)