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

• Have you done a literature search? It's not difficult to find papers that show in great detail how to compute the dispersion relation for finite difference discretizations... – Wolfgang Bangerth Dec 11 '15 at 19:34
• @WolfgangBangerth Yes I have but you found something I haven't please share with me. – eusoubrasileiro Dec 14 '15 at 14:17
• It's been many years since I dealt with the subject, so I don't have anything at hand any more. But my question stands: there are so many papers that answer how "generally to approach this problem" that I find your question overly broad. What concretely is it that you have trouble with? – Wolfgang Bangerth Dec 14 '15 at 22:19
• @WolfgangBangerth thanks for the effort. The concrete question is that I don't know any of the many methodologies to use in my 2D space time panels from my simulation to create a graph of phase velocity in function of frequency $v(f)$ . I have found how to do it analytically from the F.D. discretization of the wave equation but what I need is to compare my code with it. So I need to find another way to calculate that curve. – eusoubrasileiro Dec 16 '15 at 15:57
• But that's not a good argument. You're saying "I don't know these methods", but what I'm reading is "I'm too lazy to read the papers, can someone do it for me?" Presumably your code uses a numerical method that you can write down. What stops you from using the existing and well-described approaches to what you write down as the method used in your code? – Wolfgang Bangerth Dec 18 '15 at 6:15

I assume you'll get $p(\mathbf{r},t)$ as solution. Take only one dimension for example i.e. writing as $p(x,t)=\Sigma_n \Sigma_m a_{mn} Sin(k_mx-\omega_nt)$. If you Fourier decompose it in space first, you'll $k_m$s. Then for each $k_m$, do Fourier transform in time to get $\omega_n$. Now plot of $\omega_n$ vs $a_{mn}$ for one n at a time. The peak will correspond to the leading $k_m$ for that $\omega_n$ (Here I'll argue that only leading order matters; this can be refined to find second and later order dependence on $k_m$). Think about this way-- if some $a_{mn}$ is numerically zero for some m (call $m_0$) for one particular n (call $n_o$), this means that the component of wave, $Sin(k_{m_0}x-\omega_{n_{0}}t)$ doen't exist in the full wave. Similarly, for other n. This way you'll end up with the pairs $(n,m)$. In general, you'll get different m for each n. Finally, find a fit to this curve or use the plot to define your dispersion relation.
Please look into care and considerations of discrete fourier transform in space and time which has some limits on your choice of wavelengths ($k_m$) and time period ($\omega_n$) by cell size and time step.
Sweep long $x, y$ and $z$ separately to get $(k_x,k_y,k_z)$ in each cell to expand this to three dimensions.