# Numerical integration with special function over the positve half line

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.

• Mathematica file is here uocmath.wikispaces.com/file/view/ma_file.nb/615770393/… [![enter image description here][1]][1] Jul 20 '17 at 7:56
• Please post your actual Mathematica code, not an image of it. Without real code no one will be able to work with it to see what you might have done wrong, nor will they be able to experiment with possible repairs. Jul 20 '17 at 12:12
• @m_goldberg: The nb file with the code is available through the link uocmath.wikispaces.com/file/view/ma_file.nb/615770393/… under the question. Jul 20 '17 at 16:54
• For reference, the integrand can be expressed in terms of the $G$ function: 2 Sqrt[π] - 2 MeijerG[{{1/2, 3/4}, {1/4, 1}}, {{0}, {}}, 4/x, 1/4] Nov 1 '17 at 4:08

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

• If you use PrecisionGoal -> 50, WorkingPrecision -> 150 + 3 i, then the integral decreases to ~10^-81, so the errors in the terms computed by NIntegrate in this answer dominate the result. All it shows is that the integral is probably very small. Dec 10 '17 at 23:39

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.

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

• Interestingly, the terms you have written out exactly cancel. Try, eg. N[31 Sqrt[2/3]/9 - 31 Gamma[1/4] Gamma[3/4]/(9 Sqrt[3] \[Pi]), 10]. Jul 20 '17 at 9:22
• (and FullSimplify[(2 2^(5/6))/Sqrt[3] == (2 2^(1/3) Gamma[1/4] Gamma[3/4])/(Sqrt[3] \[Pi])] ) Jul 20 '17 at 9:23
• @evanb: Thank you for your valuable comment. I changed my mind and edited my answer. The plot produced with WorkingPrecision ->15 is exacter. Jul 20 '17 at 9:27
• The AccuracyGoal -> 10 indicates that NIntegrate should accept a result with an error estimate of 10^-10 or less, which is greater than the result returned in this case. Playing with your option settings, I can get positive and negative values down to 10^-35. The integral might be zero. Dec 10 '17 at 14:56