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

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;
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!

• afaik, there are two standard methods of dealing with such integrals: (1) use a variable transformation that maps your infinite interval to a finite one, or (2) use a special quadrature scheme that can deal with (semi-) infinite integration intervals, such as Gauss-Laguerre. Commented Aug 9, 2017 at 19:14
• You may want to look at one of the following posts: this one or this one
– Paul
Commented Aug 9, 2017 at 20:30
• That change of variables looks slightly dubious to me, can you post the whole code please? Moreover, can you specify how it fails, what the actual output is? (It can fail in different ways, and you're not being specific enough there.) Whatever the reason, try your code first on a simpler testcase, such as $\int_1^\infty x^{-2}\,\mathrm{d}x$—it will make your problem easier to diagnose. Commented Aug 9, 2017 at 23:46
• Sorry, that should be $\int_1^\infty\int_0^1 x^{-2}\,\mathrm{d}y\,\mathrm{d}x$. Commented Aug 9, 2017 at 23:53
• I tested the code I posted above with the function you suggested and it outputs 5.03506 instead of 1. Any idea what might be going on? Commented Aug 10, 2017 at 16:10

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);
}