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.

  • $\begingroup$ 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. $\endgroup$ – tpg2114 Dec 11 '15 at 15:40
  • 1
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$ – Qmechanic Dec 11 '15 at 15:44
  • $\begingroup$ @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. $\endgroup$ – eusoubrasileiro Dec 11 '15 at 16:17
  • $\begingroup$ 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. $\endgroup$ – Nick P Feb 10 '16 at 7:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.