# How to calculate dispersion relation from a Finite Difference (FD) wave simulation

I have a python code that calculates the solution of the inhomogeneous acoustic wave equation for a 2D medium with any velocity and source configuration. It was implemented using Finite Differences from the following set of equations:

$$\rho( \mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) - \frac{1}{c^2(\mathbf{r})} \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t} = -s(\mathbf{r},t)$$

where,

$$s(\mathbf{r},t) = \rho( \mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t}$$

and with the propagation velocity

$$c(\mathbf{r}) = \sqrt{\frac{\kappa(\mathbf{r})}{\rho(\mathbf{r})}}$$

Where $\mathbf{\kappa}$ is the adiabatic compression modulus of the medium, $i_v$ is the source and $\rho$ is the density of the medium.

The FD schema uses 2nd order discretization in time and 4th in space and is implemented for 2D space.

My question is how do I calculate from my numerical simulation the dispersion relation of my code? In fact i want to calculate the phase velocity $v(k) = \frac{\omega(k)}{k}$ dispersion.

I want to do that to compare with the expected relation from literature. I expect it varying with different grid angles, source frequencies and stability parameters. I know how to input those in my simulation but I dont know how to use the results (2D time panels) to establish the dispersion relation.

• Is there a reason you want to do this from numerical results and not analytically like it is usually done? The dispersion (and dissipation) can be looked at analytically for linear equations. – tpg2114 Dec 11 '15 at 15:40
• Would Computational Science be a better home for this question? – Qmechanic Dec 11 '15 at 15:44
• @tpg2114 due numerical solution approximation of F.D. the dispersion relation is not expected to be linear anymore it depends on the parameters chosen for its stability and consequently dispersion. I don't know if it is that what u mean. I want to know how my code is performing comparing with the analytical approximations for F.D. wave codes. – eusoubrasileiro Dec 11 '15 at 16:17
• Depending on what your solutions look like, you can define a hilbert transform of your dependent variable, to get a phase, and hence (if geometrical optics holds) a frequency and wavenumber. – Nick P Feb 10 '16 at 7:26