2D numerical integration with infinite limit (C++)

Question

In order to integrate a two dimensional function of the form

$$\int_{1}^\infty \int_{-\sqrt{x^2-1}}^{\sqrt{x^2-1}} e^{-x} \rm{d}y \rm{d}x,$$

I have been attempting to use the following code (written in C++) taken mostly from the Numerical Recipes book which calls a gaussian quadrature routine for the integration:

 static float xsav;
 static float (*nrfunc)(float,float);

 float quad2d(float (*func)(float, float), float x1, float x2)
 {
     float qgaus(float (*func)(float), float a, float b); 
     float f1(float x);

     nrfunc=func;
     return qgaus(f1,x1,x2); 
 }

 float f1(float x) 
 {
     float qgaus(float (*func)(float), float a, float b); 
     float f2(float y);
     float yy1(float),yy2(float);

     xsav=x;
     return qgaus(f2,yy1(x),yy2(x)); 
 }

 float f2(float y)  
 {
     return (*nrfunc)(xsav,y);
 }

This code works fine for two dimensional integrals with finite limits, but fails as the outer limit is taken to infinity. To account for this, I have attempted to use a change of variables:

#define FUNC(x) ((*funk)(-log(x))/(x))

float qgaus(float (*funk)(float), float aa, float bb)
{
    int j;
    float xr,xm,dx,s,a,b;

    b=exp(-aa); 
    a=0.0;

    static float x[]={0.0,0.1488743389,0.4333953941,
        0.6794095682,0.8650633666,0.9739065285};
    static float w[]={0.0,0.2955242247,0.2692667193, 
        0.2190863625,0.1494513491,0.0666713443};

    xm=0.5*(b+a); 
    xr=0.5*(b-a);
    s=0;
    for (j=1;j<=5;j++)
    {
        dx=xr*x[j];
        s += w[j]*(FUNC(xm+dx)+FUNC(xm-dx)); 
    }
    return s *= xr; 
}


float f(float x, float y)
{
    float a = exp(-x);
    return a;
}

float yy1(float x)
{
    float y = -sqrt(x*x-1);
    return y;   
}

float yy2(float x)
{
    float y = sqrt(x*x-1);
    return y;
}



static float xsav;
static float (*nrfunc)(float, float);

float quad2d(float (*func)(float, float), float x1, float x2)
{
    float qgaus(float (*func)(float), float aa, float bb);
    float f1(float x);

    nrfunc=func;
    float t = qgaus(f1,x1,x2);
    return t;
}

float f1(float x)
{
    float qgaus(float (*func)(float), float aa, float bb);
    float f2(float y);
    float yy1(float);
    float yy2(float);

    xsav=x;
    float r = qgaus(f2,yy1(x),yy2(x));
    return r;
}

float f2(float y) 
{
    float k = (*nrfunc)(xsav,y);
    return k;
}

int main ()
{
    float z;
    z = quad2d(f, 1.0, 20.0);
    cout << z << endl;
}

but this still doesn't give the right answer. It should be $2 \times \rm{BesselK}[1,1] \approx 1.20381$, but instead gives 2.15501.

Any suggestions on how I could modify this code to account for the infinite limit would be greatly appreciated!

Answer 1

There is a standard way to deal with infinite limits in integrals. This is explained nicely in the README of the cubature package (which you should use for low-but-multi-dimensional integration).

Infinite intervals

Integrals over infinite or semi-infinite intervals is possible by a change of variables. This is best illustrated in one dimension.

To compute an integral over a semi-infinite interval, you can perform the change of variables $x=a+t/(1-t)$:

$$\int_a^\infty f(x) dx = \int_0^1 f\biggl(a+\frac{t}{1-t}\biggr) \frac{1}{(1-t)^2} dt.$$

For an infinite interval, you can perform the change of variables $x=t/(1 - t^2)$:

$$\int_{-\infty}^\infty f(x) dx = \int_{-1}^1 f\biggl(\frac{t}{1-t^2}\biggr) \frac{1+t^2}{(1-t^2)^2} dt.$$

Note the Jacobian factors multiplying $f(\dots)$ in both integrals, and also that the limits of the $t$ integrals are different in the two cases.

In multiple dimensions, one simply performs this change of variables on each dimension separately, as desired, multiplying the integrand by the corresponding Jacobian factor for each dimension being transformed.

---

In your specific case where you can easily work out the $y$ integration by hand, I'd use the GSL package to compute the integral with the existing quadrature package, which can deal with infinite intervals.

Using $$\int_{-\sqrt{x^2-1}}^{\sqrt{x^2-1}} \exp(-x) dy = 2 \exp(-x) \sqrt{x^2-1}$$ one then has

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f(double x, void * params)
{
  return 2*exp(-x)*sqrt(x*x-1);
}

int main()
{
  gsl_integration_workspace * w = gsl_integration_workspace_alloc(100);
  double result, error;

  gsl_function F;
  F.function = &f;

  gsl_integration_qagiu(&F, 1, 1e-7, 1e-7, 100, w, &result, &error);

  printf ("result          = %.18f\n", result);
  printf ("estimated error = %.18f\n", error);
  printf ("intervals       = %zu\n", w->size);

  gsl_integration_workspace_free (w);
}