Hypergeometric function integration

Question

I have a problem with an integral: $$\small \int_0^\infty \alpha^{-\frac{5}{4}}\left(\alpha^{ik_2}+\alpha^{-ik_2}\right) {}_{2}F_1\left(\frac{1}{4}-ik_1-ik_2,\frac{1}{4}+ik_1+ik_2,\frac{1}{2},-\alpha\right) {}_{2}F_1\left(\frac{1}{2}-ik_3-ik_2,\frac{1}{2}+ik_3+ik_2,1,-\frac{1}{\alpha}\right) \, d\alpha$$ where $ {}_{2}F_1$ is the hypergeometric function. Wolfram Mathematica gives answer but it's tedious. Could you give some information about how to solve evaluate such integral? $k_i\in \mathbb{R}^+$.

Some information connected to the problem which I know: There is a useful integral for this case: $$\int_0^\infty \alpha^{\gamma-1} {}_{2}F_1\left(\frac{\gamma}{2}-ik_1,\frac{\gamma}{2} + ik_1,\gamma,-\alpha\right) {}_{2}F_1\left(\frac{\gamma}{2}-ik_2,\frac{\gamma}{2} + ik_2,\gamma,-\alpha\right) \, d\alpha\propto \delta(k_1-k_2)f(k_1)$$ where $f$ is a product of $\Gamma$-functions.
It is connected with sturm-liouville operator: $-\frac{d}{d\alpha}\alpha(\alpha+1)\frac{df}{d\alpha}+\frac{(1-\gamma)^2}{4\alpha}f$

Function: $\alpha^{\frac{1}{2}(\gamma-1)} {}_{2}F_1(\frac{\gamma}{2}-ik,\frac{\gamma}{2}+ik,\gamma,-\alpha)$ is an eigenfunction of this operator with eigenvalues $k^2+\frac{1}{4}$ If you have some information about that I'll be glad to know it.

Answer 1

To evaluate something this tough, I'd start with Mathematica. Let's see what it (version 10.4.1.0, Linux x86 (64-bit)) does with this integral.

integrand[ k_, \[Alpha]_ ] := \[Alpha]^(-5/4) *
   ( \[Alpha]^(I k[[2]]) + \[Alpha]^(-I k[[2]]) ) *
   Hypergeometric2F1[
      1/4 - I (k[[1]] + k[[2]]), 
      1/4 + I (k[[1]] + k[[2]]), 
      1/2, -\[Alpha]] *
   Hypergeometric2F1[
      1/2 - I (k[[3]] + k[[2]]), 
      1/2 + I (k[[3]] + k[[2]]), 
      1, -1/\[Alpha]];
result = Assuming[{
      k1 \[Element] Reals, k2 \[Element] Reals, k3 \[Element] Reals,
      k1 > 0, k2 > 0, k3 > 0
   },
   FullSimplify[
      Integrate[
         integrand[{k1, k2, k3}, y], 
         {y, 0, \[Infinity]}
      ]
   ]
]

(*
    (Sqrt[\[Pi]] (
       Cosh[(k2 + k3) \[Pi]] Gamma[1/4 - I k1 - I k2] Gamma[1/4 + I k1 + I k2] (
          Csc[\[Pi]/4 + I k2 \[Pi]] * 
          Gamma[1/4 - 2 I k2 - I k3] * 
          Gamma[1/4 + I k3] * 
          HypergeometricPFQRegularized[
             {1/4 - I k1 - I k2, 1/4 + I k1 + I k2, 1/4 - 2 I k2 - I k3, 1/4 + I k3}, 
             {1/2, 3/4 - I k2, 3/4 - I k2}, 1] + 
          Gamma[1/4 - I k3] * 
          Gamma[1/4 + 2 I k2 + I k3] * 
          HypergeometricPFQRegularized[
              {1/4 - I k1 - I k2, 1/4 + I k1 + I k2, 1/4 - I k3, 1/4 + 2 I k2 + I k3}, 
              {1/2, 3/4 + I k2, 3/4 + I k2}, 1] *
          Sec[\[Pi]/4 + I k2 \[Pi]]
       ) + 
       Sqrt[2] \[Pi] Sech[2 k2 \[Pi]] (
          -Gamma[1/2 - I k1] * 
          Gamma[1/2 + I k1 + 2 I k2] * 
          HypergeometricPFQRegularized[
             {1/2 - I k1, 1/2 + I k1 + 2 I k2, 1/2 - I k2 - I k3, 1/2 + I k2 + I k3}, 
             {1, 3/4 + I k2, 5/4 + I k2}, 
             1] * 
          (Cosh[k2 \[Pi]] - I Sinh[k2 \[Pi]]) - 
          Gamma[1/2 + I k1] * 
          Gamma[1/2 - I k1 - 2 I k2] * 
          HypergeometricPFQRegularized[
             {1/2 + I k1, 1/2 - I k1 - 2 I k2, 1/2 - I k2 - I k3, 1/2 + I k2 + I k3}, 
             {1, 3/4 - I k2, 5/4 - I k2}, 
             1] * 
          (Cosh[k2 \[Pi]] + I Sinh[k2 \[Pi]])
       )
    )) / 
    (Gamma[1/4 - I k1 - I k2] Gamma[1/4 + I k1 + I k2])
*)

... which certainly is a mess. Let's simplify that at $k_2 = 1$ and graph it to get a glimpse of how much numerical treachery we are facing.

plot1 = Assuming[{
      k1 \[Element] Reals, k2 \[Element] Reals, k3 \[Element] Reals,
      k1 > 0, k2 > 0, k3 > 0
   },
   FullSimplify[result /. k2 -> 1]
];
(* If you inspect the value of plot1, you will see that the regularization has 
   been split out of the hypergeometrics as extra factors of Gamma[]. *)
Plot3D[Abs[plot1], {k1, 0, 5}, {k3, 0, 7}]

... so, a lot of treachery. Sometimes increasing WorkingPrecision (above the default, $MachinePrecision = 15.9546) helps, especially for functions that have substantial cancellation in the sums that define them.

Plot3D[Abs[plot1], {k1, 0, 5}, {k3, 0, 7}, WorkingPrecision -> 24]

The non-smooth surface normals in the corner of large k1 and k3 suggest we are near the edge of losing the meager precision shown, so we add a few more digits, giving a result that is better, but not perfect.

Plot3D[Abs[plot1], {k1, 0, 5}, {k3, 0, 7}, WorkingPrecision -> 40]

So if we are going to use Mathematica to evaluate this integral, we should be careful about specifying precision and accuracy goals. As an example, with $k_1 = k_2 = k_3 = 100$ ("medium" choices of the $k_i$):

ParallelTable[
   Block[
      {$MaxExtraPrecision = Max[2^(k + 4), 50]}, 
      Timing[
         N[result /. {k1 -> 100, k2 -> 100, k3 -> 100}, 
            {2^(k + 4), 2^(k + 4)}
         ]
      ]
   ], 
   {k, 0, 5}
]

(*  Errors about hitting $MaxExtraPrecision  *)

(*  
    {   {90.044, 0.*10^346 + 0.*10^346 I}, 
        {56.188, 0.*10^330 + 0.*10^330 I}, 
        {88.068, 0.*10^284 + 0.*10^284 I}, 
        {79.136, 0.*10^171 + 0.*10^171 I}, 
        {124.788, -0.0030119028984024280992543112994012888... + 0.*10^-61 I}, 
        {198.756, -0.0030119028984024280992543112994012888... + 0.*10^-348 I}}
*)

So we get no digits of precision until the PrecisionGoal and AccuracyGoal are somewhat greater than the $k_i$s and the computation time will be long.

---

This question and answer over at CS.SE discuss how to evaluate the various (regularized) hypergeometric functions appearing in your integrand and result. References there are to Abramowitz and Stegun's "Handbook of Mathematical Functions" which has continued life as NIST's DLMF (Digital Library of Mathematical Functions). Chapter 15 covers the Gauss hypergeometric function ${}_2F_1(a,b;c;z)$ and its regularized version.

• Section 15.1 covers notation.
• Section 15.19 covers computation methodology.

The DLMF recommends splitting your integral and applying identities so that all evaluations of ${}_2F_1(a,b;c;z)$ have $z$ in the interval $\left[ 0,\frac{1}{2} \right]$.

You have \begin{align*} f(k_1, k_2, k_3) &= {}_{2}F_1\left(\frac{1}{4}-ik_1-ik_2,\frac{1}{4}+ik_1+ik_2,\frac{1}{2},-\alpha\right) \text{,} \\ g(k_1, k_2, k_3) &= {}_{2}F_1\left(\frac{1}{2}-ik_3-ik_2,\frac{1}{2}+ik_3+ik_2,1,-\frac{1}{\alpha}\right) \text{, and} \\ I(k_1, k_2, k_3) &= \int_0^\infty \; \alpha^{-\frac{5}{4}}\left(\alpha^{ik_2}+\alpha^{-ik_2}\right) f(k_1, k_2, k_3) g(k_1, k_2, k_3) \, \mathrm{d}\alpha \text{.} \end{align*} Following the suggestion, we write \begin{align*} I(k_1, k_2, k_3) &= J_1(k_1, k_2, k_3) + J_2(k_1, k_2, k_3) + K_1(k_1, k_2, k_3) + K_2(k_1, k_2, k_3) \\ J_1(k_1, k_2, k_3) &= \int_0^{1/2} \; \alpha^{-\frac{5}{4}}\left(\alpha^{ik_2}+\alpha^{-ik_2}\right) f(k_1, k_2, k_3) g(k_1, k_2, k_3) \, \mathrm{d}\alpha \\ J_2(k_1, k_2, k_3) &= \int_{1/2}^1 \; \alpha^{-\frac{5}{4}}\left(\alpha^{ik_2}+\alpha^{-ik_2}\right) f(k_1, k_2, k_3) g(k_1, k_2, k_3) \, \mathrm{d}\alpha \\ K_1(k_1, k_2, k_3) &= \int_1^2 \; \alpha^{-\frac{5}{4}}\left(\alpha^{ik_2}+\alpha^{-ik_2}\right) f(k_1, k_2, k_3) g(k_1, k_2, k_3) \, \mathrm{d}\alpha \\ K_2(k_1, k_2, k_3) &= \int_2^\infty \; \alpha^{-\frac{5}{4}}\left(\alpha^{ik_2}+\alpha^{-ik_2}\right) f(k_1, k_2, k_3) g(k_1, k_2, k_3) \, \mathrm{d}\alpha \end{align*}

• In $J_1$, the $- \alpha$ in $f$ ranges over $[-1/2,0]$, so we apply identity 15.8.1 (remembering to translate between regularized and unregularized functions).
• In $J_1$, the $-1/\alpha$ in $g$ ranges over $(-\infty, -2]$, so we apply identity 15.8.3.
• et c.

Having done this, there are a profusion of new hypergeometric functions, but no apparent simplifications. Unless there is a sneaky identity application, I think you're stuck with evaluating with software: Mathematica, mpmath (for Python, see hyp2f1(a,b,c,z)), or some other tool that is capable.

---

Since we seem to be stuck with software, and you have access to Mathematica, can we build an automated evaluator? (Note: I have not done the numerical analysis to show that the following works or requires the minimal amount of computation to get the desired accuracy/precision. I'm just leaning on prior experience trying to get combinations of hypergeometric functions to evaluate accurately.) Starting with a cut-and-paste of result from above, and specifying the evaluation only occurs when all arguments are numbers,

evalIntegral[{k1in_?NumericQ, k2in_?NumericQ, k3in_?NumericQ}, targetPrecision_?NumericQ] := Module[{
      digits, realTargetPrecision,
      k1, k2, k3
   },
   realTargetPrecision = Max[$MachinePrecision, targetPrecision];
   digits = 2(k1in + k2in + k3in + realTargetPrecision + 40);
   (* The "+40" is based on what we saw in the original Plot3D[]s. *)
   {k1, k2, k3} = SetPrecision[{k1in, k2in, k3in}, digits];
   N[Block[{
         $MaxExtraPrecision = digits
      },
      N[
         (Sqrt[\[Pi]] (
            Cosh[(k2 + k3) \[Pi]] * 
            Gamma[1/4 - I k1 - I k2] * 
            Gamma[1/4 + I k1 + I k2] (
               Csc[\[Pi]/4 + I k2 \[Pi]] * 
               Gamma[1/4 - 2 I k2 - I k3] * 
               Gamma[1/4 + I k3] * 
               HypergeometricPFQRegularized[
                  {1/4 - I k1 - I k2, 
                   1/4 + I k1 + I k2, 
                   1/4 - 2 I k2 - I k3, 
                   1/4 + I k3}, 
                  {1/2, 3/4 - I k2, 3/4 - I k2}, 
                  1] + 
               Gamma[1/4 - I k3] * 
               Gamma[1/4 + 2 I k2 + I k3] * 
               HypergeometricPFQRegularized[
                  {1/4 - I k1 - I k2, 
                   1/4 + I k1 + I k2, 
                   1/4 - I k3, 
                   1/4 + 2 I k2 + I k3}, 
                  {1/2, 3/4 + I k2, 3/4 + I k2}, 
                  1] * 
               Sec[\[Pi]/4 + I k2 \[Pi]]
            ) + 
            Sqrt[2] \[Pi] Sech[2 k2 \[Pi]] (
               -Gamma[1/2 - I k1] * 
               Gamma[1/2 + I k1 + 2 I k2] * 
               HypergeometricPFQRegularized[
                  {1/2 - I k1, 
                   1/2 + I k1 + 2 I k2, 
                   1/2 - I k2 - I k3, 
                   1/2 + I k2 + I k3}, 
                  {1, 3/4 + I k2, 5/4 + I k2}, 
                  1] * 
               (Cosh[k2 \[Pi]] - I Sinh[k2 \[Pi]]) - 
               Gamma[1/2 + I k1] * 
               Gamma[1/2 - I k1 - 2 I k2] * 
               HypergeometricPFQRegularized[
                  {1/2 + I k1, 
                   1/2 - I k1 - 2 I k2, 
                   1/2 - I k2 - I k3, 
                   1/2 + I k2 + I k3}, 
                  {1, 3/4 - I k2, 5/4 - I k2}, 
                  1] * 
               (Cosh[k2 \[Pi]] + I Sinh[k2 \[Pi]])
            )
         )) / (
            Gamma[1/4 - I k1 - I k2] * 
            Gamma[1/4 + I k1 + I k2]
         ),
         {digits, digits}
      ]
   ],
      realTargetPrecision
]]

And we test a few values. We don't get $MaxExtraPrecision errors this time. We could suppress them with Quiet[], if needed.

evalIntegral[{100, 100, 100}, 10]
evalIntegral[{100, 100, 100}, 20]
evalIntegral[{100, 100, 100}, 40]
evalIntegral[{100, 100, 100}, 80]

(* No N::meprec errors now... *)
(*
    -0.003011902898402428 + 0.*10^-19 I
    -0.0030119028984024280993 + 0.*10^-23 I
    -0.003011902898402428099254311299401288849657 + 0.*10^-43 I
    -0.0030119028984024280992543112994012888496566189627235819934233737819167464566341282 + 0.*10^-83 I
*)

Happily, the digits we claimed we wanted to be good are replicated with the higher precision calculations and the last digits of prior results are rounded in the correct direction.

So, do we get a more encouraging graph?

Plot3D[Abs[evalIntegral[{k1, 1, k3}, 10]], {k1, 0, 5}, {k3, 0, 7}]

Yes, but very slowly.