Numerical integration with special function over the positve half line

Question

Is there a way to compute the following integral numerically?,

$\int_0^{\infty}\left(-\frac{2 x \Gamma \left(\frac{1}{4}\right) \, _1F_3\left(\frac{1}{4};\frac{1}{2},\frac{3}{4},\frac{5}{4};\frac{x^4}{256}\right)}{\sqrt{\pi }}+x^2 \, _1F_3\left(\frac{1}{2};\frac{3}{4},\frac{5}{4},\frac{3}{2};\frac{x^4}{256}\right)+\sqrt{\frac{2}{\pi }} \Gamma \left(\frac{1}{4}\right) \Gamma \left(\frac{3}{4}\right)\right)dx$

I know that $-\frac{2 x \Gamma \left(\frac{1}{4}\right) \, _1F_3\left(\frac{1}{4};\frac{1}{2},\frac{3}{4},\frac{5}{4};\frac{x^4}{256}\right)}{\sqrt{\pi }}+x^2 \, _1F_3\left(\frac{1}{2};\frac{3}{4},\frac{5}{4},\frac{3}{2};\frac{x^4}{256}\right) \rightarrow -\sqrt{\frac{2}{\pi }} \Gamma \left(\frac{1}{4}\right) \Gamma \left(\frac{3}{4}\right)$ as $x \rightarrow \infty.$

Mathematica file is here.

Answer 1

The integral over the half-oscillations between roots of the integrand seems to be decreasing rapidly, so one can estimate the whole integral as an alternating sum:

(* Integrand from @user64494 *)
g = x^2*HypergeometricPFQ[{1/2}, {3/4, 5/4, 3/2}, x^4/256] - 
   2*Gamma[1/4]/(Sqrt[Pi])*x*
    HypergeometricPFQ[{1/4}, {1/2, 3/4, 5/4}, x^4/256] + 
   Sqrt[2/π] Gamma[1/4] Gamma[3/4];

(* Find roots of the integrand *)
(roots = With[{f = g}, 
     Reap@NDSolve[{t'[x] == 1, t[0] == 0, 
        WhenEvent[f == 0, Sow[x]]}, {}, {x, 0, 200}, 
       PrecisionGoal -> 7, WorkingPrecision -> 500, 
       MaxStepSize -> 1/2]][[2, 1]];
 roots2 = FindRoot[g, {x, #}, WorkingPrecision -> 500] & /@ roots;
) // AbsoluteTiming
(*  {64.5725, Null}  *)

(* terms to sum 0..Infinity *)
ClearAll[term];
mem : term[i_?NumberQ] := 
  mem = Module[{iRes}, 
    Check[iRes = 
      NIntegrate[g, {x, x /. roots2[[i]], x /. roots2[[1 + i]]}, 
       PrecisionGoal -> 20, AccuracyGoal -> Infinity, 
       WorkingPrecision -> 100 + 2 i, 
       Method -> {"GaussKronrodRule", "Points" -> 15}], 
     Print["i = ", i]];
    iRes];
term[0] = NIntegrate[g, {x, 0, x /. roots2[[1]]}, PrecisionGoal -> 20, 
   AccuracyGoal -> Infinity, WorkingPrecision -> 200, 
   Method -> {"GaussKronrodRule", "Points" -> 15}];

NSum[term[Round@i], {i, 0, ∞}, Method -> "AlternatingSigns",
   "VerifyConvergence" -> False, WorkingPrecision -> 100] // AbsoluteTiming
(*
  {61.4593, 
   -9.41471026473132271164542801174142339911608047122339101*10^-58}
*)

Checking the last term calculated by NSum gives an upper bound on the error, assuming the terms continue to be alternating and monotonically decreasing:

DownValues[term][[-2]]
(*  HoldPattern[term[131]] :> -1.76932...88*10^-104  *)

Answer 2

Instead of summing an alternating series, for which NIntegrate seems to have trouble doing accurately enough, we can integrate out to a finite value. The integrand seems to be oscillating and decreasing rapidly, so the error of trunction will be less than the area between the next hump in its graph and the x axis. It turns out that as I pushed the interval of integration further and further, the value of the integral would decrease and always be around the error estimate. It makes me think the integral might actually be zero. If so, one would rather want a proof than an approximation.

We can use the Fourier DCT to approximate the Chebyshev series of g (see this answer by J.M.). The series is easily integrated. The biggest trouble is evaluating the function accurately. I used N[g /. x -> exactValue, precision] to hopefully ensure that accurately values of g are computed.

(* Integrand from @user64494 *)
g = x^2*HypergeometricPFQ[{1/2}, {3/4, 5/4, 3/2}, x^4/256] - 
   2*Gamma[1/4]/(Sqrt[Pi])*x*
    HypergeometricPFQ[{1/4}, {1/2, 3/4, 5/4}, x^4/256] + 
   Sqrt[2/π] Gamma[1/4] Gamma[3/4];

Block[{$MaxExtraPrecision = 500},
  nn = 1024;
  xmax = 220;
  gvals = N[g /. x -> #, 150] & /@
    Rescale[Sin[Pi/2*Range[nn, -nn, -2]/nn], {-1, 1}, {0, xmax}];
  ] // AbsoluteTiming
(*  {45.0455, Null}  *)

(* approximate Chebyshev series *)
cc = foo = Sqrt[2/nn] FourierDCT[gvals, 1];
cc[[{1, -1}]] /= 2;
cc = Drop[cc, -LengthWhile[Reverse@cc, Precision[#] < 2 &]];
Length@cc
(*  610  *)

(* Integrate a Chebyshev series over its complete domain *)
i0Cheb[c_, {a_, b_}] := Total[(b - a)/(1 - Range[0, Length@c - 1, 2]^2) c[[1 ;; ;; 2]]];
i0Cheb[cc, {0, xmax}]
(*  -1.508526852*10^-138  *)

Notes: (1) N[g /. x -> 220, 150] needs a $MaxExtraPrecision of at least 400 and takes nearly a second to evaluate on my machine. The integrand is numerically very hard to evaluate. For N[g /. x -> 1000, 150], we need a $MaxExtraPrecision of at least 3079 and nearly 7 seconds. NIntegrate is going to take a long time, as will any accurate method. (2) The precision 150 allows the integral to be computed up to about 150 decimal places after the decimal point, since Max@Abs@cc is about 0.196. The answer shows 147 decimal places of Accuracy.

Answer 3

This works for me.

g = x^2*HypergeometricPFQ[{1/2}, {3/4, 5/4, 3/2}, x^4/256] -  
2*Gamma[1/4]/(Sqrt[Pi])*x*HypergeometricPFQ[{1/4}, {1/2, 3/4, 5/4}, x^4/256] + 
Sqrt[2/\[Pi]] Gamma[1/4] Gamma[3/4];
NIntegrate[g, {x, 0, Infinity}, Method -> "ExtrapolatingOscillatory", 
 AccuracyGoal -> 10, WorkingPrecision -> 30]

2.53237432183970183831603063456*10^-11