Analytic solution 2D scalar wave equation in cylindrical coordinates numerical implementation

Question

I am trying to compare my finite difference's solution of the scalar (or simple acoustic) wave equation with an analytic solution.

For that purpose I am using the following analytic solution presented in the old paper Accuracy of the finite-difference modeling of the acoustic wave equation - Geophysics 1974 - R. M. Alford et. al.

The solution, equation (9) in that paper, for a constant velocity medium is given in cylindrical coordinates by:

\begin{eqnarray} u_s(\rho, \phi, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \{ -i \pi H_0^{(2)}\left(k | \sigma - \sigma_s| \right) F(\omega) e^{i\omega t} d\omega\} \end{eqnarray} (A)

$H_0^{(2)}$ is the Hankel function of the second kind with order zero.

$| \sigma - \sigma_s| $ is the distance from the source to the observation point.

Also $ k = \frac{\omega}{C_0}$ come from $C_0$ the medium velocity converted to the Fourier domain, defined in the original scalar wave equation as:

\begin{eqnarray} \left[ \nabla^2 - \frac{1}{C_0^2}\frac{\partial ^2}{\partial t^2}\right] u(\rho, \phi, t) = -4 \pi \frac{\delta(\rho-\rho_s)\delta(\phi-\phi_s)f(t)}{\rho} \end{eqnarray}

$ (\rho_s, \phi_s) $ defines the source location in the $\phi$ x $\rho$ plane.

$ F(\omega)$ comes from the source function $f(t)$

The problem:

• When implementing numerically A I'm not getting meaningful results. The velocity obtained using A is totally different from the $C_0$ used as input

My code example is bellow where I am evaluating A, using python. I calculate the solution at zero solution0 and at a distance D=8000 solutionD. Then I use the global maximum to calculate an approximation to velocity $C_0$

Probably it's a quite simple error that I am doing. If someone can point it out for me. The velocity I am getting shows that maybe the equation (9) is wrong or I am implementing it wrongly (what I believe to be the case).

from scipy.special import hankel2
import numpy as np

# a simple source function f(t) gauss derivative normalized
def gaussource(time, wlength, delay=None):
        if delay == None: # enough delay time
            delay = 3*wlength
        t = time - delay
        return ((2*np.sqrt(np.e)/(wlength))
               *t*np.exp(-2*(t**2)/(wlength**2)))

# analytic solution equation (A) paper
def waves(source, distance, c):
    n = len(source)
    sourcew = np.fft.fft(source) # source in the frequency domain
    hnk_factor = distance/c # distance to source divided by velocity
    if not distance == 0.0:
        hankelshift = np.complex(0,-np.pi)*np.array([hankel2(0,hnk_factor*omega) for omega in xrange(n)])
        hankelshift[0] = 0. # i infinity in limit
        return np.real(np.fft.ifft(hankelshift*sourcew))
    return source

c = 4000.
fc = 5.0
dt = 0.01
t = np.arange(0., 2., dt) # 2 seconds to evaluate the solution
source = gaussource(t, 1./fc) # source
solution0 = waves(source, 0., c) # solution at zero source (source itself)
solutionD = waves(source, 8000., c) #solution at certain distance D=8000 metres e.g.   
t0 = solution0.argmax(); # get max index time at zero
tD = solutionD.argmax(); # get max index time at distance
print "phase velocity aproximation ", 8000./((tD-t0)*dt) 

Answer 1

After some weeks and no one found the answer to my problem. I found myself the errors and the solution. To help the community I'am sharing it (meta advice).

It was indeed a very simple error.

Really basic sampling theorem, maximum frequency based on sample rate is:

$$ f_{max} = 1/\Delta t $$

or using $ \omega $

$$ \omega_{max} = 2\pi/\Delta t $$

I changed the def waves(source, distance, c): function definition. I used $\Delta t$ to calculate $ \Delta \omega $ to sample the solution in frequency. Also removed the calculation at zero distance since the source is the solution itself.

# analytic solution equation (A) paper
def waves(source, distance, c, dt):
    n = len(source)
    sourcew = np.fft.fft(source) # source in the frequency domain
    hnk_factor = distance/c # distance to source divided by velocity
    dw = np.pi*2.0*(1./dt)*1/n
    hankelshift = np.complex(0,-np.pi)*np.array([hankel2(0.0,hnk_factor*k*dw) for k in xrange(n)])
    hankelshift[0] = 0. # i infinity in limit
    return np.real(np.fft.ifft(hankelshift*sourcew))

t = np.arange(0., 3.5, dt) # 3.5 seconds to evaluate the solution
source = gaussource(t, 1./fc) # source and solution at 0
solutionD = waves(source, 8000., c, dt) # solution at certain distance     D=8000 metres e.g.   
t0 = source.argmax(); # get max index time at zero
tD = solutionD.argmax(); # get max index time at distance

I also increased the solution time spam to 3.5 seconds.

And bellow is how it did work, it even can give closer results changing the sample rate.

phase velocity approximation 3883.49514563