1

関数のフーリエ (FT) 変換と逆フーリエ (iFT) 変換の間で計算を行ったり来たりできる Python (SymPy lib を使用) コードを作成したいと考えています。関数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