In this article we develop and implement a split region technique that solves the time-dependent acoustic wave equation with greatly increased efficiency. This method uses a sophisticated Chebyshev propagation scheme in areas where there are interfaces and medium variations, and a simple free space propagator where the medium is homogeneous. Mathematica provides a cohesive and interactive environment, where the mathematical functions and visualization tools required for this work are already built in. The interactive interface allows users to modify the code and study specific problems with ease.

The Acoustic Wave Propagator

Theory

In an earlier paper by Pan and Wang [1], an explicit acoustic wave propagator (AWP) was proposed to describe the time-domain evolution of mechanical waves in various media. This method was based on a similar scheme that was originally developed by Tal-Ezer and Kosloff in [2, 3, 4, 5, 6], who studied seismic wave propagation and a variety of gas-phase reactive scattering and related chemical processes. The AWP method has been successfully applied to study both the propagation of a flexural wave in a thin plate by Peng, Pan, and Sum [7] and an acoustic wave in a room by Sun, Wang, and Pan [8] and [9]. However, this method requires significant computer resources since the AWP propagation scheme utilizes a large set of modified Chebyshev polynomials with Bessel functions of the first kind as the expansion coefficients.

In this article we develop and implement a split region technique in Mathematica that uses the sophisticated Chebyshev propagation scheme in areas where there are interfaces and medium variations, but a simple and much more efficient free space propagator where the medium is homogeneous.

First, we give a brief outline of the AWP and then introduce the Chebyshev propagation scheme. This method is then implemented in Mathematica code. In a later section, the free space propagator is defined, along with a method of integration that uses Fourier transforms. Finally, the split region technique is constructed and propagation is demonstrated across a splitting region with a boundary.

The motion of acoustic waves in air and solids can be described by a partial differential equation known as the acoustic wave equation:

(1)

Integrating this with respect to time, yields a formal solution to this equation:

(2)

where denotes the spatial coordinates collectively and stands for time, with being the initial starting time. is a state vector, while is the system Hamiltonian that describes the physical properties of the propagation and the boundary medium. In the case of a one-dimensional duct, describes the sound pressure and the particle velocity :

(3)

and is of the form:

(4)

where is the speed of sound within the medium, and is the density of the medium.

The AWP is defined as:

(5)

However, this exponential operator is impractical in its current form, and thus it must be expanded as a finite polynomial. In this work we use a Chebyshev polynomial expansion. To ensure the convergence of this expansion, the system Hamiltonian must be normalized:

(6)

where is the maximum eigenvalue of the system operator . In the case of sound pressure in a one-dimensional duct, this value is:

(7)

If we let:

(8)

then the AWP may be expressed as:

(9)

The next step is to make a simple, if somewhat nonintuitive, change of variables. We let:

(10)

Then we expand the exponential operator in terms of the Chebyshev polynomials :

(11)

Using the orthogonality relationship for the Chebyshev polynomials, we find that the coefficients are:

(12)

where and for , and is a Bessel function of the first kind. However, there are complex numbers involved in this expansion, while the state vector and operator are real. Hence, we define a new set of modified Chebyshev polynomials, as given by:

(13)

It can be shown that these satisfy the following recursion relation:

(14)

with and . It is now possible for us to write the AWP, , in the form:

(15)

which only involves real-valued operations. We now obtain our state vector with an expanded AWP:

(16)

which appears in full as:

(17)

The benefit of this method is that the Bessel functions decay exponentially with the coefficient index when . These properties are very useful for numerical computation, as they allow expansions of the exponential function to be accurately calculated for arbitrarily large values of (that is, arbitrarily large time steps).

Numerical Implementation

Discretising the Problem

In this section we closely follow the work of Falloon and Wang [10], with the necessary changes for the AWP. In particular, we need to handle two components representing the sound wave, instead of just one as in the electronic wavefunction.

The first step is to create a one-dimensional grid of length that contains points, which stands for our spatial dimension .

We also a need a similar grid in -space, as the evaluation of derivatives involves taking the Fourier transform of the functions:

This grid has a grid spacing of , and consequently can represent a maximum wave-number of :

A norm function is defined for both the position and the -space grids:

The following functions are defined to plot the pressure and velocity distributions in both the position and -spaces.

The position and distribution functions form a Fourier transform pair and are related by the transformations:

(18)

For discretised functions, such as those being used here, the Fast Fourier Transform (FFT) algorithm may be used. The built-in Fourier functions of Mathematica utilise this algorithm. The functions defined next allow us to transform distributions between the position and -spaces. The

command allows us to easily shift our distributions from the interval to , in order to provide a centred transform. These distributions are normalised by the factor .

To test the AWP, we use Gaussian distributions for both the pressure and the velocity components of the wave description. A Gaussian distribution is of the form:

(19)

where its width is determined by . This defines a Gaussian function of this form:

The speed of sound in the medium () and the density of the material () should not change throughout the calculation of the propagation. Other parameters that should not change are the initial pressure pulse spread () and position (0), the duct length (), the propagation time (time), and the numbers of grid points for the entire region (num) and the Chebyshev region (nint). These have been defined using global variables, denoted with the Mathematica $name notation.

Here we turn off some unimportant Mathematica warning messages.

Defining the Propagator

Both components of the system Hamiltonian in the AWP apply a first-order derivative to each component of the state vector . Thus, a good starting place for the implementation of the propagator is to define a discrete differentiation operator. DelGrid takes a given function, labelled as grid, and computes its derivative via Fourier transformations.

The function Propagator implements the Chebyshev propagation scheme. First a number of definitions are made, then a Do loop is constructed to perform and repeat the Chebyshev expansion the number of times determined by the term M. The pressure and velocity distributions are handled simultaneously, due to their related definitions. The whole set of commands is contained within a Module function so that all definitions remain local. Finally, we use Compile on the overall function, which keeps the calculation time down by keeping all values numerical.

Testing the Propagator

Now we carry out a few tests of this propagation scheme. Here the initial pressure and velocity distributions are defined, followed by their respective derivatives. First the pressure:

Then the velocity:

Now we run the Chebyshev expansion propagation scheme for these distributions. This returns the pressure and velocity distributions of the acoustic wave after it has propagated for seconds.

Now we compare the initial and final pressure and velocity distributions and can see that the acoustic wave has indeed propagated to the right.

The plots show a Gaussian sound pulse propagating to the right, travelling a distance of approximately 12 metres in seconds. This corresponds to a wave speed of approximately , the expected speed of sound in the material. This observation agrees perfectly with expectations—for a sound pulse like that defined earlier, the exact solution is a Gaussian pulse travelling to the right at the speed of sound in the medium. A more exact measurement of the speed and distance travelled gives:

This measurement confirms that the propagation speed of the sound wave agrees with the specified speed of sound in air.

The Region Splitting Propagator

The aim of our work has been to develop an algorithm that propagates an acoustic wave through a region that is mainly free space, but also has regions where there are changes in the media and obstructions to the wave’s propagation. In principle, the Chebyshev expansion scheme is sufficient to perform any calculation that is necessary in this situation, but for a large grid, the computational effort that it requires becomes prohibitive. The exact solution that is available for regions of free, homogeneous space involves less computational effort, but is not appropriate for regions where there are changes in the media or obstructions to the wave. Because of this, we construct an algorithm where the space is divided into multiple regions—regions of free space are handled by the exact solution, and for those where there are variations, the Chebyshev scheme is used.

The Free Space Propagator

An exact solution to the problem of the propagation of a sound wave through free space is available and applicable to any situation that is free of medium changes and interfaces.

It can be shown that the pressure () and particle velocity () waves that satisfy the equation:

(20)

with the parameters described earlier, also satisfy the second-order wave equation:

(21)

which has the initial conditions:

(22)

The solution to this partial differential equation is of the form:

(23)

Due to the excellent accuracy and efficiency of the numerical differentiation technique based upon the FFT [11], we seek a similar technique for the integration required in the previous equation. We substitute the function with its transform:

(24)

Thus the integral can be computed numerically by taking the transform of the original function, multiplying it by , and then taking the inverse transform. We implement the Fourier integration method as the function FourierIntegrate:

The function FreePropagator implements the free space propagation of an acoustic wave in a homogeneous medium. First, the constants and parameters of the system are defined, then the correct shift of grid spacings in the first two terms of equation (23), nshift, is calculated, and the initial pressure and velocity distributions are converted to the left and right travelling forms. Finally, the resulting pressure and velocity distributions are calculated, as well as the proper time of propagation, actualt. All of these commands are gathered into a single Module function, and Compile is used to ensure that they run as efficiently as possible.

The FreePropagator function is used to calculate the propagation of the initial pressure and velocity distributions.

Here is the final pressure distribution.

Here is the final velocity distribution.

Here we compare the initial and final pressure distributions and the initial and final velocity distributions.

The Split Region Propagator

The splitting technique that we use involves multiplying the total acoustic wave distribution by what is essentially a step function, where the step occurs at the intended boundary between regions. This method divides the wave distribution into two sections that can be propagated separately, by the most appropriate technique for that region. Due to the linear nature of the splitting, these two sections can be easily recombined by adding the two resulting distributions together. The whole process can then be repeated if more than one time step is desired.

Ideally, the splitting function would be a step function:

(25)

where one region lies between and , and the other is everything outside that range. This is often called a “top hat function”. However, the discontinuities that such a splitting function would introduce into the wave distribution would create large inaccuracies in the numerical techniques that are to be used. To avoid this consequence, a more gentle splitting function is used, where the step function is convolved with a Gaussian curve, so that the splitting actually occurs over some width , centred on . This eliminates the discontinuities of the pure step function and the associated difficulties and errors. We first define the step function:

Now the function SplittingGrid is defined to create a grid of values from the convolution of the step function and the Gaussian. All values of this lie in the range from 0 to 1. The wave distribution is multiplied by this grid to achieve the required splitting. This calculation makes use of the fact that the transform of the convolution of two functions is the product of the transforms of the individual functions. Later, the new smooth top hat function is plotted.

For the subsequent calculations, we use the following split region barrier. It is important to note that there are as yet no variations in the medium—the speed of sound and the density remain the same across all three regions.

Now we define the function SplittingPropagator, which implements the split region propagation technique. The SplittingGrid function (here it is labelled as fgrid) is used to split the distributions, but it is the quantities nlower, npad, and nupper that divide the propagator into separate regions via a Take selection. In the outer regions, where there is only free space, the FreePropagator function is used, while in the central region, the Chebyshev method is used in the form of the Propagator function. A Do command is used to repeat this propagation a number of times determined by steps.

And now we propagate the acoustic wave through both free space and the splitting region, using SplittingPropagator. Note that the grid should be large enough that a shift of in either direction does not move any significant (that is, nonzero) data outside the boundaries of the grid. We arbitrarily choose the propagation time to be 0.5 second and split the propagation into five time steps with 0.1 second each.

The splitting propagator scheme works as expected. The sound wave propagates through the splitting and Chebyshev interaction regions without introducing any discontinuity and the norm is conserved throughout. Again, the following measure confirms that the propagation speed of the sound wave agrees well with the specified speed of sound in air.

Note that it is vital to ensure that the sizes of the various regions are appropriate—the interaction region must be larger than the splitting region, as the interaction wave packet must have enough space to get clear of the splitting zone (inside the barrier region defined by fgrid) before leaving the interaction zone (defined by nint). Otherwise the wave packet is wrapped around inside the interaction zone due to the Fourier methods used in the calculation, thus never leaving that area.

It is clear that the sizes of the splitting region, the interaction zone, and the time step need careful determination. The buffer zone in the interaction region surrounding the split region needs to be large enough that portions of the wave packet that are at the edge of the split region at the start of the time step do not travel further than the edge of the buffer region, otherwise they are artificially wrapped, thus destroying the accuracy of the scheme. The size of the buffer region needs to be greater than the distance travelled by the wave in the length of the time step:

This length is taken from the point where the splitting grid effectively falls to zero.

Conclusion

In this article we developed a highly efficient technique that solves the time-dependent acoustic wave equation using split regions, where a free space solution is used in free space, and the more sophisticated Chebyshev expansion solution is used in the interaction regions. The next stage is to incorporate variations and boundaries in the medium in order to simulate a “real-life” wave passing through air, liquids, and solids. This should only require minor adjustments to the split region technique as it is presented in this article. Preliminary work indicates the reflection and transmission of a wave passing through a barrier, as well as a change in speed as the wave passes through a denser medium.

Acknowledgment

The authors would like to acknowledge support from the Australian Research Council and the National Science Foundation of China (Project No. 60340420325).

For this work we used Mathematica 5.0 and a Pentium 4 computer running Windows XP with a 2.8 GHz CPU and 512 MB of RAM.

References

[1]	J. Pan and J. B. Wang, “Acoustic Wave Propagator,” Journal of the Acoustical Society of America, 108(2), 2000 pp. 481-487.
[2]	H. Tal-Ezer and R. Kosloff, “An Accurate and Efficient Scheme for Propagating the Time Dependent Schrödinger Equation,” Journal of Chemical Physics, 81(9), 1984 pp. 3967-3971.
[3]	H. Tal-Ezer, “A Pseudospectral Legendre Method for Hyperbolic Equations with an Improved Stability Condition,” Journal of Computational Physics, 67(1), 1986 pp. 145-172.
[4]	H. Tal-Ezer, “Spectral Methods in Time for Hyperbolic Equations,” SIAM Journal on Numerical Analysis, 23(1), 1986 pp. 11-26.
[5]	H. Tal-Ezer, “An Accurate Scheme for Seismic Forward Modeling,” Geophysical Prospecting, 35, 1987 pp. 479-490.
[6]	H. Tal-Ezer, “Polynomial Approximation of Functions of Matrices and Applications,” Journal of Scientific Computing, 4(1), 1989 pp. 25-60.
[7]	S. Z. Peng, J. Pan, and K. S. Sum, “Acoustical Wave Propagator Method for Time Domain Analysis of Flexural Wave Scattering and Dynamic Stress Concentration in a Heterogeneous Plate with Multiple Cylindrical Patches,” Tenth International Congress on Sound and Vibration (ICSV10), Stockholm, Sweden, 2003.
[8]	H. M. Sun, J. B. Wang, and J. Pan, “An Effective Algorithm for Simulating Acoustical Wave Propagation,” Computer Physics Communications, 151, 2003 pp. 241-249.
[9]	H. M. Sun, J. B. Wang, and J. Pan, “Acoustical Wave Propagator Method for Two-dimensional Sound Propagation,” Third International Conference on Acoustics, University of Cadiz, Spain, 2003.
[10]	P. Falloon and J. B. Wang, “Electronic Wave Propagation with Mathematica,” Computer Physics Communications, 134(2), 2001 pp. 167-182.
[11]	J. B. Wang, “Numerical Differentiation Using Fourier,” The Mathematica Journal, 8(3), 2002 pp. 383-388.
N. Riste, B. McGrath, J. Wang, and J. Pan, “Acoustic Wave Propagator—A Split Region Implementation,” The Mathematica Journal, 2012. dx.doi.org/10.3888/tmj.10.3-7.

About the Authors

Neil Riste and Bradley McGrath are Ph.D. students at The University of Western Australia, studying and theoretically modelling the dynamical response of both quantum and classical systems to external disturbance.

Associate Professor Jingbo Wang leads the Quantum Dynamics Theory Group at The University of Western Australia and has a wide range of interests, from atomic physics, molecular and chemical physics, spectroscopy, acoustics, chaos, nanostructured electronic devices, and mesoscopic physics to quantum information and computation.

Professor Jie Pan is the director of the Center for Acoustics, Dynamics and Vibration at The University of Western Australia. His research areas are room acoustics, structural acoustics, and active noise and vibration control.

Neil Riste
Bradley McGrath
Jingbo Wang
Jie Pan
School of Physics
School of Mechanical Engineering
The University of Western Australia
35 Stirling Highway
Crawley WA 6009, Australia
wang@physics.uwa.edu.au