1

我想编写一个 Python（带有 SymPy lib）代码，它可以在函数的傅里叶 (FT) 和逆傅里叶 (iFT) 转换之间来回进行计算。假设有一个函数f(r)，它的 FT g(q) = FT[f(r)]和新函数f(r) = iFT[g(q)]的 iFT 。

我使用 sympy 的 Python 代码无法做到这一点，但一个简单的 Mathematica 代码可以轻松做到这一点。

  • 我在我的 Python 代码中遗漏了什么吗？
  • 或者 Sympy 无法完成我想要做的事情？

Python代码： 

import sympy as sym

q, r, alpha = sym.symbols('q r alpha', positive=True)
def sph_bessel(n, r): 
    if r<>0:
        sph = sym.sqrt((sym.pi)/(2.*r))*sym.besselj(n + 1./2., r)
    elif n==0 and r==0:
        sph = 1
    elif n<>0 and r==0:
        sph = 0
    return sph

FT = lambda f, n, q, r: 4*sym.pi*sym.integrate(f*sph_bessel(n, q*r) * r**2, (r, 0, sym.oo))
invfunc = lambda q, alpha:  alpha**4/(alpha**2 + q**2)**2.

print "              f(r) = %s"%sym.N(func(r, alpha))
print "F(q) = FT(f(r))    = %s"%sym.simplify(sym.N(FT(func(r, alpha), 0, q, r)))
print "f(r) = invFT(F(q)) = %s"%sym.simplify(sym.N(invFT(invfunc(q, alpha), 0, q, r)))

输出： f(r) = 0.0397887357729738*alpha**3*exp(-alpha*r) F(q) = FT(f(r)) = 1.0*alpha**4*(alpha**2 + 1.0*q**2)**(-2.0) f(r) = invFT(F(q)) = 0.0498677850501791*sqrt(2)*I*alpha**3.0*(alpha*r*jn(-1.0, alpha*r*exp_polar(I*pi/2)) + I*sinh(alpha*r))/sqrt(pi) Mathematica 代码：（要尝试，复制并粘贴到 Mathematica）

FourierSphrTF =   Integrate[#1 SphericalBesselJ[#2, #4 #3] #3^2, {#3, 0, \[Infinity]},  Assumptions -> #5] &;
FourierSphericalT = 4 \[Pi] FourierSphrTF[#1, #2, #3, #4, #5] &;  InvFourierSphericalT =   1/(2 \[Pi]^2) FourierSphrTF[#1, #2, #3, #4, #5] &; 
FourierSphericalT[\[Alpha]^3 E^(-r \[Alpha])/(8 \[Pi] ) , 0, r, q,   q > 0 && \[Alpha] > 0]
InvFourierSphericalT[\[Alpha]^4/(q^2 + \[Alpha]^2)^2, 0, q, r,  r > 0 && \[Alpha] > 0]

输出： f(r) = 0.0397887357729738*alpha**3*exp(-alpha*r) F(q) = FT(f(r)) = 1.0*alpha**4*(alpha**2 + 1.0*q**2)**(-2.0) f(r) = invFT(F(q)) = 0.0397887357729738*alpha**3*exp(-alpha*r)

4

0 回答 0