傅里叶空间中三维网格的数值积分

问：

我正在尝试实现本文中概述的模型：

一般静磁形状-形状相互作用

具体来说，我正在评估在z方向上被磁化并在x方向上位移的球体的等式（4）。这个方程对任意形状的对象采用“形状函数”，它只是一个二进制的三维数组，在对象存在的地方值为 1，在对象不存在的地方值为 0。然后对这些张量进行傅里叶变换，并可用于计算静磁自能和相互作用能。

我试图对方程（4）进行数值积分，并将其进行无量纲化，并将其与论文中给出的在z轴上磁化并在x轴上分离的球体的简单无量纲化解析结果进行比较，称为方程（27）。

使用粒子半径的长度单位，我编写了一个 Python 脚本来执行积分。

法典

import numpy as np

#Create the shape function of a circle

gridSize = 128 #Speed up FFT algo with a power of 2 size
D_r = np.zeros([gridSize,gridSize,gridSize]) #Preallocate grid for the shape function

radius = 64 #Define sphere size 

for xind in range(gridSize):
    x = xind - (gridSize-1)/2
    
    for yind in range(gridSize):
        y = yind - (gridSize-1)/2
        
        for zind in range(gridSize):
            z = zind - (gridSize-1)/2
            dist =  np.sqrt((x)**2+(y)**2+(z)**2)
        
            if dist < radius:
                    D_r[xind,yind,zind] = 1

"Define the frequency vectors"
unitsPerPixel = 1/radius #Rescale into units of particle radii

freqs = np.fft.fftshift(np.fft.fftfreq(gridSize,d=unitsPerPixel))
#Stretch these out into grids to define the x, y, and z frequencies at every point of the grid D(k)
kx, ky, kz = np.meshgrid(freqs, freqs, freqs, indexing='ij') 


"Define the shape amplitudes"
D_k = np.fft.fftshift(np.fft.fftn(D_r)) #Particle 1
D_k2 = D_k #Particle 2, identical in this case

"Define the magnetization and separation vectors"
m1 = np.array([0,0,1]) #particle 1 magnetization
m2 = np.array([0,0,1]) #particle 2 magnetization
rho = np.array([11,0,0]) #paricle separaion, must be >= 2

"Calculate the analytical solution for refeence"
R= 1 #Sphere radius
V = 4/3*np.pi*R**3 #Sphere volume
Ea = 2/(4*np.pi*rho[0]**3) #Dimensionless interaction energy

"Computer the numerical solution"
Etest = 0

for xind in range(gridSize):
    for yind in range(gridSize):
        for zind in range(gridSize):
            #Build the frequency vector
            k = np.array([kx[xind,yind,zind],ky[xind,yind,zind],kz[xind,yind,zind]])
            
            #Check if this is the zero point 
            if k[0] == 0 and k[1]== 0 and k[2]== 0:
                continue
            else:

                Etest += 1/V**2*D_k[xind,yind,zind]*np.conjugate(D_k2[xind,yind,zind]) \
                          *k[2]*k[2]/(k[0]**2+k[1]**2+k[2]**2) \
                             *np.exp(1j*np.dot(rho,k))


L = gridSize*unitsPerPixel
print('The numerical solution is:')
print(np.real(Etest * 1 / L**3))
print('The analytical solution is:')
print(Ea)

我的输出是巨大的，并且不像物理学所暗示的那样对粒子分离具有反立方依赖性。显然有些事情出了大问题。两个问题是我的输出太大了，而且我的结果实际上不是对象分离的强函数。与论文相比，我无法判断我做了什么不同或错误的事情。

问题

1. 对于数量级问题，我认为这可能源于我对 L 的定义。但我直接使用了作者的定义。我对 L 的定义正确吗？

2. 对于位移依赖性问题，我担心这源于我使用 NumPy 的 fftfreq 函数对频率矢量的定义，因为位移仅显示在指数傅里叶位移因子中。这听起来对吗？我的频率网格实现是否正确？

感谢您的阅读，非常感谢任何帮助。