我试图将标量（或简单声学）波动方程的有限差分解与解析解进行比较。

为此，我使用旧论文Accuracy of the Finite-Difference Modeling of the Accuracy - Geophysics 1974 - RM Alford 等旧论文中提出的以下解析解。人_

该论文中等式 (9) 的解，对于等速介质，在柱坐标中由下式给出：

$$\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}$$ （一种）

$H_0^{(2)}$是零阶的第二类汉克尔函数。

$| \sigma - \sigma_s|$是从源到观察点的距离。

还$k = \frac{\omega}{C_0}$来自$C_0$转换为傅立叶域的介质速度，在原始标量波动方程中定义为：

$$\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)$定义源位置$\phi$X$\rho$飞机。

$F(\omega)$来自源函数$f(t)$

问题：

• 当以数字方式实现A时，我没有得到有意义的结果。使用A获得的速度与$C_0$用作输入

我的代码示例如下所示，我正在使用 python 评估A。solution0我在零处和距离 D=8000 处计算解solutionD。然后我使用全局最大值来计算速度的近似值$C_0$

可能这是我正在做的一个非常简单的错误。如果有人可以为我指出。我得到的速度表明方程（9）可能是错误的，或者我错误地执行它（我认为是这种情况）。

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)