0

如果我尝试使用和比较来评估integral 0 to 2pi: (1/2pi) f(t) exp(-i*m*t) dt（f 在 t 中是周期性的，周期为 2pi），它们会给出几乎相同的结果（1e-18 数量级的差异）。但是 3 维的类似情况不匹配。我的问题是为什么？scipy.integrate.quadscipy.fft.fft

第一种情况：

def f(theta): return np.exp(-1j * np.cos(theta))

def g_real(theta, m):
    kick = -np.cos(theta)
    kernel = -m * theta
    return np.cos(kick + kernel)

def g_imag(theta, m):
    kick = -np.cos(theta)
    kernel = -m * theta
    return np.sin(kick + kernel)

def getIntegral(m):
    result_real, err_real = quad(g_real, 0, 2*np.pi, args=(m,))
    result_imag, err_imag = quad(g_imag, 0, 2*np.pi, args=(m,))
    return (result_real + 1j * result_imag) / (2*np.pi)

def getFourier():
    theta = np.arange(64) * 2 * np.pi / 64
    x = f(theta)
    return fft.fftshift(fft.fft(x, norm="forward"))

def main():
    m = 10
    fourier = getFourier()[m+32]
    integral = getIntegral(m)
    print(f"Fourier: {fourier}")
    print(f"Integral: {integral}")
    print(f"Bessel: {(1j)**(-m) * jv(-m, -1)}") # I know the analytical answer to this integral
    print(f"Diff: {np.abs(fourier - integral)}")

# Fourier: (-2.093833800239648e-05+1.2583048337422513e-18j)
# Integral: (-2.0938338002223678e-05-2.820993268635476e-17j)
# Bessel: (-2.093833800238929e-05+0j)
# Diff: 1.7529613390482267e-16

3d案例：

integral 0 to 2pi for all t1, t2, t3: (1/2pi)^3 f(t1, t2, t3) exp(-i * (m1*t1 + m2*t2 + m3*t3)) dt1 dt2 dt3（f 在所有三个变量中都是周期性的，周期为 2pi）。但这并不奏效，如下所示。

ALPHA = 0.5

def f_real(t1, t2, t3, m1, m2, m3):
    kick = -np.cos(t1) * (1 + ALPHA * np.cos(t2) * np.cos(t3))
    kernel = -(m1*t1 + m2*t2 + m3*t3)
    return np.cos(kernel + kick)

def f_imag(t1, t2, t3, m1, m2, m3):
    kick = -np.cos(t1) * (1 + ALPHA * np.cos(t2) * np.cos(t3))
    kernel = -(m1*t1 + m2*t2 + m3*t3)
    return np.sin(kernel + kick)

def f_kick_full(t1, t2, t3):
    kick = -np.cos(t1) * (1 + ALPHA * np.cos(t2) * np.cos(t3))
    return np.exp(1j * kick)

def findIntegral(m1, m2, m3):
    params = (m1, m2, m3)
    ranges = [(0,2*np.pi), (0,2*np.pi), (0,2*np.pi)]
    args_real = (f_real, ranges, params)
    args_imag = (f_imag, ranges, params)
    result_real, abserr_real = nquad(*args_real)
    result_imag, abserr_imag = nquad(*args_imag)
    return (result_real + 1j * result_imag) / (2*np.pi)**3

def getFourier():
    theta = np.arange(64) * 2 * np.pi / 64
    t1, t2, t3 = np.meshgrid(theta, theta, theta)
    x = f_kick_full(t1, t2, t3)
    y = fft.fftshift(fft.fftn(x, norm="forward"))
    return y

def main():
    m1, m2, m3 = (1, 2, 1)
    fourier_coeffs = getFourier()

    fourier = fourier_coeffs[m1+32, m2+32, m3+32]
    integral = findIntegral(m1, m2, m3)
    print("From Fourier:", fourier)
    print("From Integral:", integral)
    print("Difference:", np.abs(fourier - integral))

# From Fourier: (-0.025662302376245+1.1891906739905623e-17j)
# From Integral: (8.325321276060769e-17-5.401283469428714e-18j)
# Difference: 0.02566230237624508

由于我在上面的代码中没有看到任何重大错误，我是否误解了 fftn 或 nquad 的工作原理？

4

0 回答 0