Message Boards Message Boards

Issue with Inverse Fourier transforms of Hypergeometric functions?

Fab Nicol
Posted 2 years ago

I was trying to compute the Hilbert Transform (HT) of Exp[-x^4]
using the well-known formula:

f[w_] := InverseFourierTransform[ I Sgn[x] FourierTransform[g[y], 
           y, x,  FourierParameters -> {1, -1}], 
           x, w,  FourierParameters -> {1, -1}]

Oddly, this gave wrong results.
The first direct Fourier transform step was neatly solved:

FourierTransform[Exp[-x^4], x, w, FourierParameters -> {1, -1}]

yielding:

2 Gamma[5/4] HypergeometricPFQ[{}, {1/2, 3/4}, w^4/256] - 
 1/4 w^2 Gamma[3/4] HypergeometricPFQ[{}, {5/4, 3/2}, w^4/256]

Then the inverse Fourier transform step:

I InverseFourierTransform[ % Sign[w], w, x, 
  Assumptions -> Element[w,  Reals], FourierParameters -> {1, -1}]

gave results with unexpected complex imaginary parts:

(1/(2 Pi Abs[x]))(-1)^(
 3/4) E^x^4 x (Gamma[1/4] Gamma[3/4, x^4] - 
   I Sqrt[2] Pi ((-1 - I) + GammaRegularized[1/4, x^4]))

as testifies the shocking:

% /. x -> 1 // N
-0.261089 - 2.27256 I

Hard as I could try to work this out, I am at a loss what to think. My tentative guess is that there may be a glitch in Inverse Fourier Transform algos. Attached is the notebook. Any idea out there for a workaround?

Notes:

  1. For Exp[-x^2] the process outlined above worked and gave the right result (- 2 DawsonF[x]/Sqrt[Pi]), identical to principal value integration of the convolution with 1/x on the real axis.
  2. For Exp[-x ^4], the HT was not fount by direct symbolic integration
Attachments:
invfthypergeo.nb
POSTED BY: Fab Nicol
6 Replies
Sort By:
Fab Nicol
Fab Nicol
Posted 2 years ago

Thanks Markus, this at least gives a possible idea of the Mathematica bug (I am running Mathematica 13.0).
I am grateful for this workaround. I could replicate it step by step. Running: 

FullSimplify[InverseLaplaceTransform[(t (Sqrt[s] + t^2))/(Sqrt[2] s^(3/4) (s + t^4)), 
             s, a,  Assumptions -> t \[Element] Reals && t > 0] /. a -> 1 , 
                    Assumptions -> t \[Element] Reals && t > 0]

I arrived at: 

 HT =  (E^-t^4 (-4 I \[Pi] - 2 t^3 ExpIntegralE[1/4, -t^4] Gamma[1/4] - 
  2 t ExpIntegralE[3/4, -t^4] Gamma[3/4]))/(4 \[Pi])

and: 

  % /. t -> 1 // N
  (* 0.81944 + 2.60013*10^-17 I *)

Which established replication. The fact that the imaginary part is not numerically significant is established by: 

Mean[Im[Table[HT]^2, {t, 1, 5, 0.1}]]]
 (* 2.66728*10^-18 *)

There still remains to understand where this imaginary part comes from.
POSTED BY: Fab Nicol

Another workaround:

 ClearAll["`*"]

 ML = MellinTransform[1/Pi*Exp[-\[Tau]^4]/(t - \[Tau]), t, s]
 INT = Integrate[ML, {\[Tau], -Infinity, Infinity}, 
    Assumptions -> {s > 0}] // ExpandAll
  InverseMellinTransform[INT, s, t](*Hilbert Transfrom Solution:*)

 (*(-MeijerG[{{0, 1/4, 1/2, 3/4}, {}}, {{0, 0, 1/4, 1/2, 3/4}, {}}, -t, 
    1/4] + MeijerG[{{0, 1/4, 1/2, 3/4}, {}}, {{0, 0, 1/4, 1/2, 3/
     4}, {}}, t, 1/4])/(8 \[Pi]^4)*)

  % /. t -> 1 // N(*For t=1,Only Real part is True*)

   (*0.81944 + 0.367879 I*)
POSTED BY: Mariusz Iwaniuk
Fab Nicol
Fab Nicol
Posted 2 years ago

I consider the Mellin transform result * not * a workaround though, since it gives a numerically significant erroneous imaginary part... Rather, this looks like an instantiation of the same Wolfram algo bug. In itself, it is not more (or less) reliable than the initial InverseFourier-based post.
POSTED BY: Fab Nicol

I wrote: "Real part is True".

I fond another solution expressed by MeijerG function:

 MeijerG[{{0, 1/4, 1/2, 3/4}, {}}, {{0, 0, 1/4, 1/2, 3/4}, {}}, t^4]/(
    8 \[Pi]^4) - 
    MeijerG[{{0, 1/4, 1/2, 3/4}, {1/8, 3/8, 5/8, 7/8}}, {{0, 0, 1/4, 1/
       2, 3/4}, {1/8, 3/8, 5/8, 7/8}}, t^4] /. t -> 1 // N

(*0.81944 + 1.27204*10^-17 I*)
POSTED BY: Mariusz Iwaniuk
Fab Nicol
Fab Nicol
Posted 2 years ago

Yes, but as long as the imaginary part is there, it is impossible to rule out a doubt (unless you performed other converging methods independently, which you did :-) As far as the Mellin/InverseMellin transform method is concerned, the convincing workaround would rather be:

InverseMellinTransform[-Tan[ Pi \[Rho]/2 ] MellinTransform[
   Exp[- t^4], t, \[Rho]], \[Rho], x] // FullSimplify

(See Frederick W. King, Hilbert Transforms, Cambridge U.P, p. 274, which works for even functions).

This yields the elegant :

    (1/2 - I/2) E^-x^4 ((-1 + I) + GammaRegularized[1/4, -x^4] - 
         I GammaRegularized[3/4, -x^4])

Now for the reduction to reals, I have not completed the algebra yet, but the reasoning goes as follows: incomplete Gammas with a strictly negative argument are rotated in the complex plane by Pi theta / 2 with respect to the value for positive arguments (see NIST handbook, 8.5.5, Confluent Hypergeometric Representations), where theta is the form index.
Here the increase in the form index is 1/2 for an argument that can only be strictly negative or null, hence a rotation by Pi /4 from the first Gamma to the second. The null case is obvious, and otherwise after some basic algebra it follows that the big bracket is proportional to Exp[ i Pi / 4], which multiplied by the (1 - I) initial factor yields a real.
POSTED BY: Fab Nicol

I only have a workaround:

 ClearAll["`*"]
 LT = LaplaceTransform[Exp[-a x^4], a, s]
 (*1/(s + x^4)*)

Then:

  f[x_] := 1/(s + x^4)
  HT = -(1/Pi)*
    Integrate[(
     f[t + \[Tau]] - f[t - \[Tau]])/\[Tau], {\[Tau], eps, Infinity}, 
     Assumptions -> {t \[Element] Reals, s > 0, 
       eps > 0}](*Hilbert transform*)

  (*-(1/(8 \[Pi] s^(
    3/4) (s + t^4)))(-4 Sqrt[2] \[Pi] Sqrt[s] t - 4 Sqrt[2] \[Pi] t^3 + 
     2 t (Sqrt[2] Sqrt[s] + 2 s^(1/4) t + Sqrt[2] t^2) ArcTan[
       1 + (Sqrt[2] (eps - t))/s^(1/4)] - 
     2 t (Sqrt[2] Sqrt[s] - 2 s^(1/4) t + Sqrt[2] t^2) ArcTan[
       1 + (Sqrt[2] (-eps + t))/s^(1/4)] - 
     2 t (Sqrt[2] Sqrt[s] + 2 s^(1/4) t + Sqrt[2] t^2) ArcTan[
       1 - (Sqrt[2] (eps + t))/s^(1/4)] + 
     2 t (Sqrt[2] Sqrt[s] - 2 s^(1/4) t + Sqrt[2] t^2) ArcTan[
       1 + (Sqrt[2] (eps + t))/s^(1/4)] - 
     Sqrt[2] Sqrt[s]
       t Log[eps^2 + Sqrt[2] eps s^(1/4) + Sqrt[s] - 2 eps t - 
        Sqrt[2] s^(1/4) t + t^2] + 
     Sqrt[2] t^3 Log[
       eps^2 + Sqrt[2] eps s^(1/4) + Sqrt[s] - 2 eps t - 
        Sqrt[2] s^(1/4) t + t^2] + 
     Sqrt[2] Sqrt[s]
       t Log[eps^2 - Sqrt[2] eps s^(1/4) + Sqrt[s] + 2 eps t - 
        Sqrt[2] s^(1/4) t + t^2] - 
     Sqrt[2] t^3 Log[
       eps^2 - Sqrt[2] eps s^(1/4) + Sqrt[s] + 2 eps t - 
        Sqrt[2] s^(1/4) t + t^2] + 
     Sqrt[2] Sqrt[s]
       t Log[eps^2 - Sqrt[2] eps s^(1/4) + Sqrt[s] - 2 eps t + 
        Sqrt[2] s^(1/4) t + t^2] - 
     Sqrt[2] t^3 Log[
       eps^2 - Sqrt[2] eps s^(1/4) + Sqrt[s] - 2 eps t + 
        Sqrt[2] s^(1/4) t + t^2] - 
     Sqrt[2] Sqrt[s]
       t Log[eps^2 + Sqrt[2] eps s^(1/4) + Sqrt[s] + 2 eps t + 
        Sqrt[2] s^(1/4) t + t^2] + 
     Sqrt[2] t^3 Log[
       eps^2 + Sqrt[2] eps s^(1/4) + Sqrt[s] + 2 eps t + 
        Sqrt[2] s^(1/4) t + t^2] - 
     2 s^(3/4)
       Log[eps^4 + s - 4 eps^3 t + 6 eps^2 t^2 - 4 eps t^3 + t^4] + 
     2 s^(3/4)
       Log[eps^4 + s + 4 eps^3 t + 6 eps^2 t^2 + 4 eps t^3 + t^4])*)
  L = Limit[HT, eps -> 0, Assumptions -> t \[Element] Reals]
  InverseLaplaceTransform[L[[1]], s, a, 
     Assumptions -> t \[Element] Reals] /. a -> 1 // Expand(*Hilbert Transform Solution:*)

  (*-(((-1)^(1/4) E^-t^4 t^3)/(Sqrt[2] Abs[t]^3)) - ((-1)^(3/4)
     E^-t^4 t)/(
   Sqrt[2] Abs[t]) + ((-1)^(3/4) E^-t^4 t Gamma[1/4, -t^4])/(
   Sqrt[2] Abs[t] Gamma[1/4]) + ((-1)^(1/4)
     E^-t^4 t^3 Gamma[3/4, -t^4])/(Sqrt[2] Abs[t]^3 Gamma[3/4])*)

    % /. t -> 1 // N(*For t=1,Only Real part exist*)

   (*0.81944 - 5.55112*10^-17 I*)
POSTED BY: Mariusz Iwaniuk
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract