Numerical Integral from 0 to infinity

Question

My aim is to calculate the numerical integral of a probability distribution function (PDF) of the distance of an electron from the nucleus of the hydrogen atom in C programming language. I have written a sample code however it fails to find the numerical value correctly due to the fact that I cannot increase the limit as much as its necessary in my opinion. I have also included the library but I cannot use the values stated in the following post as integral boundaries: min and max value of data type in C . What is the remedy in this case? Should switch to another programming language maybe? Any help and suggestion is appreciated, thanks in advance.

Edit: After some value I get the error segmentation fault. I have checked the actual result of the integral to be 0.0372193 with Wolframalpha. In addition to this if I increment k in smaller amounts I get zero as a result that is why I defined r[k]=k, I know it should be smaller for increased precision.

#include <stdio.h>
#include <math.h>
#include <limits.h>
#define a0 0.53 
int N = 200000;
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[], long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {

    long double P[N], r[N], a;
    // Declare and initialize the loop variable
    int k = 0;
    for (k = 0; k <  N; k++)
    {
        r[k] = k ;
        P[k] = r[k] * r[k] * exp( -2*r[k] / a0);
        //printf("%.20Lf \n", r[k]);
        //printf("%.20Lf \n", P[k]);
    }
    a = trapezoid(r, P);    
    printf("%.20Lf \n", a);
}

Last Code:

#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#define a0 0.53 
#define N LLONG_MAX
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[],long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {
    printf("%Ld", LLONG_MAX);
    long double * P = malloc(N * sizeof(long double));
    long double * r = malloc(N * sizeof(long double));
    // Declare and initialize the loop variable
    int k = 0;
    long double integral;
    for (k = 1; k <  N; k++)
        {
        P[k] = r[k] * r[k] * expl( -2*r[k] / a0);
        }
    integral = trapezoid(r, P);
    printf("%Lf", integral);

}

Edit last code working:

#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#define a0 0.53 
#define N LONG_MAX/100
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[],long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {
    printf("%Ld \n", LLONG_MAX);
    long double * P = malloc(N * sizeof(long double));
    long double * r = malloc(N * sizeof(long double));
    // Declare and initialize the loop variable
    int k = 0;
    long double integral;
    for (k = 1; k <  N; k++)
        {
        r[k] = k / 100000.0;
        P[k] = r[k] * r[k] * expl( -2*r[k] / a0);
        }
    integral = trapezoid(r, P);
    printf("%.15Lf \n", integral);
    free((void *)P);
    free((void *)r);
}

In particular I have changed the definition for r[k] by using a floating point number in the division operation to get a long double as a result and also as I have stated in my last comment I cannot go for Ns larger than LONG_MAX/100 and I think I should investigate the code and malloc further to get the issue. I have found the exact value that is obtained analytically by taking the limits; I have confirmed the result with TI-89 Titanium and Wolframalpha (both numerically and analytically) apart from doing it myself. The trapezoid rule worked out pretty well when the interval size has been decreased. Many thanks for all the posters here for their ideas. Having a value of 2147483647 LONG_MAX is not that particularly large as I expected by the way, should the limit not be around ten to power 308?

Answer 1

Numerical point of view

The usual trapezoid method doesn't work with improper integrals. As such, Gaussian quadrature rules are much better, since they not only provide 2n-1 exactness (that is, for a polynomial of degree 2n-1 they will return the correct solution), but also manage improper integrals by using the right weight function.

If your integral is improper in both sides, you should try the Gauss-Hermite quadrature, otherwise use the Gauss-Laguerre quadrature.

The "overflow" error

long double P[N], r[N], a;

P has a size of roughly 3MB, and so does r. That's too much memory. Allocate the memory instead:

long double * P = malloc(N * sizeof(long double));
long double * r = malloc(N * sizeof(long double));

Don't forget to include <stdlib.h> and use free on both P and r if you don't need them any longer. Also, you may not access the N-th entry, so f[N] is wrong.

Using Gauss-Laguerre quadrature

Now Gauss-Laguerre uses exp(-x) as weight function. If you're not familiar with Gaussian quadrature: the result of E(f) is the integral of w * f, where w is the weight function.

Your f looks like this, and:

f x = x^2 * exp (-2 * x / a)

Wait a minute. f already contains exp(-term), so we can substitute x with t = x * a /2 and get

f' x = (t * a/2)^2 * exp(-t) * a/2

Since exp(-t) is already part of our weight function, your function fits now perfectly into the Gauss-Laguerre quadrature. The resulting code is

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

/* x[] and a[] taken from 
 * https://de.wikipedia.org/wiki/Gau%C3%9F-Quadratur#Gau.C3.9F-Laguerre-Integration
 * Calculating them by hand is a little bit cumbersome
*/
const int gauss_rule_length = 3;
const double gauss_x[] = {0.415774556783, 2.29428036028, 6.28994508294};
const double gauss_a[] = {0.711093009929, 0.278517733569, 0.0103892565016};

double f(double x){
    return x *.53/2 * x *.53/2 * .53/2;
}

int main(){
    int i;
    double sum = 0;
    for(i = 0; i < gauss_rule_length; ++i){
        sum += gauss_a[i] * f(gauss_x[i]);
    }
    printf("%.10lf\n",sum); /* 0.0372192500 */
    return 0;
}