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) $$


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

  • 2
    $\begingroup$ 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... $\endgroup$ Commented Dec 11, 2015 at 19:34
  • $\begingroup$ @WolfgangBangerth Yes I have but you found something I haven't please share with me. $\endgroup$
    – imbr
    Commented Dec 14, 2015 at 14:17
  • $\begingroup$ 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? $\endgroup$ Commented Dec 14, 2015 at 22:19
  • $\begingroup$ @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. $\endgroup$
    – imbr
    Commented Dec 16, 2015 at 15:57
  • $\begingroup$ 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? $\endgroup$ Commented Dec 18, 2015 at 6:15

1 Answer 1


I'm trying to answer it with physical intuition and arguments.

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.