Numerical integration of a hypergeometric function

Question

The Task

Let $z_1, z_2, z_3$ be positive real numbers and define

$$ r(\mathbf{z}):= \int_0 ^\infty (t+z_1)^{-3/2}(t+z_2)^{-3/2}(t+z_3)^{-1/2}\text{d}t. $$

The task is to compute $r$ numerically in an efficient way. Absolute error on the order of $10^{-8}$ is acceptable.

Attempts at a Solution

In the notation of DLMF 19.16, $r$ can be related to the general multivariate hypergeometric function $$ R_{-a}(\mathbf{b};\mathbf{z} ) := \frac{1}{B(a,a')}\int_0^\infty t^{a'-1} \prod_{j=1}^n (t+z_j)^{-b_j}\ \text{d}t, $$ where $B$ is the beta function and $a'=-a+\sum_{j=1}^nb_j$, under the restrictions that $a,a' > 0$ and $z_j \in \mathbb{C}\setminus (-\infty,0]$, by $$ r(\mathbf{z})= \frac{2}{5}R_{-\frac{5}{2}}\left( \frac{3}{2}, \frac{3}{2}, \frac{1}{2} ; \mathbf{z} \right). $$ There is a generalized series representation of $R_{-a}$, but I would guess that it's more expensive to compute than using a canned routine on the indefinite integral above (which is currently too expensive). The DLMF lists several special cases of $R_{-a}$, many of which permit inexpensive recursive evaluations through addition theorems, but the integral here is not one of them. I can, however, rearrange $r$ into two forms that are special cases, when $z_1 = z_2$ and when $z_1 \ne z_2$, as follows.

Case 1: $z_1 = z_2$

Assuming we can pass the limit through the integral, $$ \lim_{z_2 \to z_1} r(\mathbf{z})= \int_0 ^\infty (t+z_1)^{-3}(t+z_3)^{-1/2}\text{d}t $$ which we can rewrite in DLMF 19 notation as $R_{-\frac{5}{2}}\left( 3,\frac{1}{2};z_1,z_3\right)$ and identify with the Gauss series $\,_2F_1(a,b;c;z)$: $$ R_{-a}\left(b_1,b_3;z_1,z_3\right) = z_3^{-a}\ _2F_1\left(a,b_1;b_1+b_3;1-(z_1/z_3)\right) $$ for which there exist recurrence relations.

Case 2: $z_1 \ne z_2$

When $z_1 \ne z_2$ we can express $r$ as $$ r=\frac{\chi_1-\chi_2}{z_2-z_1} $$ where \begin{eqnarray} \chi_1&:=&\int_0^\infty \frac{1}{(t+z_1)\sqrt{(t+z_2)(t+z_3)}}dt \\ \chi_2&:=&\int_0^\infty \frac{1}{(t+z_2)\sqrt{(t+z_1)(t+z_3)}}dt \end{eqnarray} are elliptic integrals corresponding to the special case $R_D(x,y,z)$ in DLMF 19 notation (up to a constant), for which there is another recurrence relation. (See this question for details).

Splitting the Cases Is Unstable

In principle one could use the first method when $z_1 = z_2$ and the second when $z_1 \ne z_2$; however, in practice I need to compute the integral over a discretization of $z_2$ in which $z_2 \to z_1$, and the method above becomes unstable when $z_2 \simeq z_1$. Is there a way to reduce the integral to a form for which there exists a relatively inexpensive recurrence, as there is for each of the two cases?

Answer 1

First of all, from the first paragraph of your attempts at a solution, I assume that the $z_j$ are non-negative? In that case, the integrand has no real problematic points (it's monotonous, decreasing). The only difficulty is the infinite integration range. And then it all depends on the accuracy you want... Do you want machine epsilon precision (or there about) or is it ok for you to aim for $10^{-8}$ or ... ?

I'm impressed by the work you did on finding hypergeometric representations of the integral and trying to find a recurrence for cheap evaluation. However, in this case with such a simple integrand, I would first check if brute force quadrature isn't fast enough.

I would suggest you to try the so-called double-exponential quadrature for an infinite integral. You can find an overview paper of the double-exponential method and its applications in the paper The double-exponential transformation in numerical analysis . This paper contains references to the original papers in which the method was developed. A very nice and free (as in free beer) implementation can be found on professor Ooura's webpage. He has "didactic" implementations (so you can map the implementation to the algorithm described in the papers) and a fast implementation (a bit more obscure code).

I quickly (and dirtily) implemented the function above and linked to the C version of the fast implementation, using the intdei function for infinite range integration. For a result with a relative error tolerance of $10^{-15}$, the code needs typically 100 function evaluations (if $z_j < 1000$, more if $z_j$ are larger) , which in your case are extremely cheap.

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

#include "intde2.h"

int nfunc;
double z1, z2, z3;

double f(double x)
{
    extern int nfunc;
    extern double z1,z2,z3;

    nfunc++;

    const double p1 = -3.0/2;
    const double p2 = -1.0/2;
    const double f1 = pow(x+z1,p1);
    const double f2 = pow(x+z2,p1);
    const double f3 = pow(x+z3,p2);

    return f1*f2*f3;
}

void main()
{
    extern double z1,z2,z3;
    extern int nfunc;
    double tiny, aw[8000], i, err;

    const double a = 10000.0;

    srand(time(NULL));

    z1 = (double)rand()/(double)(RAND_MAX/a);
    z2 = (double)rand()/(double)(RAND_MAX/a);
    z3 = (double)rand()/(double)(RAND_MAX/a);
    lenaw = 8000;
    tiny = 1.0e-307;

    intdeiini(lenaw, tiny, 1.0e-15, aw);
    nfunc = 0;
    intdei(f, 0.0, aw, &i, &err);
    printf("z1 = %20.12e\nz2 = %20.12e\nz3 = %20.12e\n",z1,z2,z3);
    printf(" I =  %20.12e\nerr= %20.12e\nN= %d\n", i, err, nfunc);
}