1

I want to write a Python (with SymPy lib) code that it can make calculations back and forth between Fourier (FT) and inverse Fourier (iFT) Transformation of a function. Let's say there is a function f(r), FT of it g(q) = FT[f(r)] and iFT of new function f(r) = iFT[g(q)].

My Python code using sympy could not achieve this but a simple Mathematica code can do this easily.

  • am I missing something in my Python code?
  • or Sympy is not capable of what I am trying to do?

The Python Code: 

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)))

Output: 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 Code:(To try, copy&paste to 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]

Output: 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)

hesap
  • 11
  • 1
  • SymPy already has `jn`, the spherical Bessel function of the first kind. – asmeurer Oct 25 '16 at 16:10
  • How did you define `invFT` for SymPy? Also, SymPy already has `fourier_transform` and `inverse_fourier_transform`. – asmeurer Oct 25 '16 at 16:15
  • Thanks for the response: First of all before returning my problem: `jn()` of sympy gives for `jn(0,0) = nan` where it should `jn(0,0)=1`. And second, I did come accros to `jn` until I see your comment. I know `fourier_transform` and inverse of it are exist but what I am trying to achive is FT and invFT in Spherical Coordinates. – hesap Oct 25 '16 at 21:01
  • I've opened [an issue](https://github.com/sympy/sympy/issues/11776) for the `jn(0, 0)` problem. – asmeurer Oct 26 '16 at 16:59

0 Answers0