How many quadrature points should I use?

I am trying to compute the following integration $$\int_0^\infty e^{-y}y^{a/2}L_c^b(y)L_e^d(y)dy$$ using the generalized Gauss-Laguerre quadrature routine in the GNU Scientific Library. Here the $$L$$'s are the generalized Laguerre functions, and $$a>-1, b>0, d>0$$. I am not sure how to choose the number of quadrature points, different numbers give different values.

Any suggestions?

(Edit: MWE)

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_laguerre.h>

struct data { double b; int c; double d; int e; };

double f(double y, void* userdata) {
struct data *d = (struct data *) userdata;
return gsl_sf_laguerre_n(d->c, d->b, y) * gsl_sf_laguerre_n(d->e, d->d, y);
}

int main() {
const double a = -0.5;
const double b = 0.5;
const int c = 0;
const double d = 0.5;
const int e = 0;
const size_t num_nodes = (c + e + 1) / 2;

const gsl_integration_fixed_type *T = gsl_integration_fixed_laguerre;
gsl_integration_fixed_workspace *w
= gsl_integration_fixed_alloc(T, num_nodes, 0.0, 1.0, 0.5 * a, 0.0);

struct data params = { b, c, d, e };
gsl_function F;
F.function = &f;
F.params = &params;

double result;
gsl_integration_fixed(&F, &result, w);
printf("%12.e\n", result);

gsl_integration_fixed_free(w);

}

• The fractional power of $y$ yields the non-decaying part of the integrand generally non-polynomial. You will not get the right answer with this quadrature. You will want to use a generalized Gauss-Laguerre quadrature, which can integrate functions of the form $e^{-y}y^{\alpha}f(y),\alpha\in\mathbb{R}$. See also: people.sc.fsu.edu/~jburkardt/m_src/gen_laguerre_rule/… – smh Sep 24 '18 at 17:06
• @smh I am using the generalized Gauss-Laguerre quadrature, the GNU Scientific Library implements at the generalized Gauss-Laguerre quadrature. – monstergroup42 Sep 24 '18 at 17:24
• If that's the case, then because it's a Gaussian quadrature and the remaining portion of the integral is polynomial, you need only $(c+e+1)/2$ points to integrate the Laguerre polynomial product. The rest of the integrand is the weight function, which is baked into the integration weights. You should also update your question to reflect the fact that you are using generalized Gauss-Laguerre quadrature. – smh Sep 24 '18 at 18:28
• That is what I am doing, but the answers are not matching with Mathematica. Also these are the generalised Laguerre functions, not the Laguerre polynomials. – monstergroup42 Sep 24 '18 at 19:23
• What Mathematica function are you using to evaluate your $L$ functions? Your original question says you're using "associated" Laguerre functions (which are polynomials), but your last comment says you're using "generalized" Laguerre functions (which are not). Which is it? – smh Sep 24 '18 at 21:10

The way printf formats work, a call like

printf("%12.e\n", result);


instructs printf to print a number with 12 digits/characters/spaces (of width 12), with NO digits after the decimal point, just like in the precision specifier "12.".

Changing it to

printf("%.12e\n", result);


results in printing 1.225416702465e+00, which is probably what was intended, and which matches vibe's answer.

(optional) . followed by integer number or *, or neither that specifies precision of the conversion. In the case when * is used, the precision is specified by an additional argument of type int. If the value of this argument is negative, it is ignored. If neither a number nor * is used, the precision is taken as zero. See the table below for exact effects of precision. [Emphasis added]

• You are right. It is now giving the right answer, even if I use num_nodes = (c + e + 1) / 2 + 1. – monstergroup42 Sep 27 '18 at 15:59

The following program will integrate your integral with GSL. I was unable to reproduce the result you quoted with Mathematica. Note that when $$c = e = 0$$, the Laguerre polynomials are always $$1$$ for any value of $$x$$. The code below does a slightly more interesting case of $$c = 2, e = 3$$. The output is:

result = 1.539250228072e+00


Mathematica's output is $$1.539250228072043$$. When $$c = e = 0$$ as you originally posted, both GSL and Mathematica produce the result $$1.225416702465178$$.

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_specfunc.h>
#include <gsl/gsl_integration.h>

struct data
{
double b;
int c;
double d;
int e;
};

double
f(double x, void * params)
{
struct data *d = (struct data *) params;
double Lcb = gsl_sf_laguerre_n(d->c, d->b, x);
double Led = gsl_sf_laguerre_n(d->e, d->d, x);
return Lcb * Led;
}

int
main()
{
const gsl_integration_fixed_type * T = gsl_integration_fixed_laguerre;
const double a = -0.5;
const double b = 0.5;
const int c = 2;
const double d = 0.5;
const int e = 3;
const size_t n = (c + e + 1) / 2;
gsl_integration_fixed_workspace * w =
gsl_integration_fixed_alloc(T, n, 0.0, 1.0, 0.5 * a, 0.0);
gsl_function F;
double result;
struct data data_params = { b, c, d, e };

F.function = &f;
F.params = &data_params;

gsl_integration_fixed(&F, &result, w);
fprintf(stderr, "result = %.12e\n", result);

gsl_integration_fixed_free(w);

return 0;
}

• I have now added an example to my post. It is very similar to the code that you have posted here. However I am getting $1$ for $c=e=0$, and $2$ for $c=2, e=3$. – monstergroup42 Sep 27 '18 at 12:59