multidimensional numerical integration

multidimensional numerical integration

[Kevin Jackman]

Hi all,

A couple of questions followed by some useful information for others (in my
opinion):

Q: I realize "epsabs" and "epsrel" control the desired accuracy of numerical
integration, but how? Mainly, how do they each affect the GSL routines, and
how do they differ from each other? How about on a Quadpack level? I think
the GSL manual is kind of vague on how these parameters individually
control numerical integration convergence. (No hard feelings, ok?)


Some Useful Info:

Looking through the GSL documentation, archives, and source code, I soon
discovered that GSL does not have specific routines for multidimensional
numerical integration, except using Monte Carlo methods.  I saw some
efforts to initially incorporate cubature routines:

http://sourceware.org/ml/gsl-discuss/2007-q2/msg00046.html
http://www.cygwin.com/ml/gsl-discuss/2005-q2/msg00020.html

However, as cubature has not been formally incorporated into GSL, and that
topic is far from my specialty, I moved on. Below is a little tidbit of code
which demonstrates how to recursively integrate using 1D QAGS routines and
thereby yield an effective 3D integral. The code integrates the volume of a
sphere in spherical coordinates. I realize this method can be time and computer
resource consuming for complex integrals, but I have yet to see a better GSL
alternative suited for my purposes. (Monte Carlo is not good for my needs,
long explanation).  Also, the absolute error estimation may be lost by doing it
this way, though I'm not sure.

Aside from telling me which hand I should wish in and which hand I should...,
I'd like to see something like this incorporated into the GSL manual, so others
don't have to reinvent the wheel. You are welcome to use my code, modify it,
etc... to make it up to your manual standards. No props necessary even.


compile and build with something like:
gcc -o simpleQAGS_3D simpleQAGS_3D.c -L/usr/lib -I/usr/include -lgsl
-lgslcblas -lm


simpleQAGS_3D.c:

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

double Phi (double phi, void * params) {

    return 1.0;
}


double Theta (double theta, void * params) {

    gsl_integration_workspace * w3 = gsl_integration_workspace_alloc (1000);

    double result3, error3;
    double expected3 = 2.0*M_PI;
    double alpha3 = 1.0;
    double TwoPi = 2.0*M_PI ;


    gsl_function F3;
    F3.function = &Phi;
    F3.params = &alpha3;

    gsl_integration_qags (&F3, 0, TwoPi, 1e-5, 1e-7, 1000, w3,
&result3, &error3);

    printf ("W3 result          = %.18f\n", result3);
    printf ("W3 expected result = %.18f\n", expected3);
    printf ("W3 estimated error = %.18f\n", error3);
    printf ("W3 actual error    = %.18f\n", result3 - expected3);
    printf ("W3 intervals       = %d\n", w3->size);

    gsl_integration_workspace_free (w3);

    return sin(theta)*result3;
}


double R (double r, void * params) {

    gsl_integration_workspace * w2 = gsl_integration_workspace_alloc (1000);

    double result2, error2;
    double expected2 = M_PI;
    double alpha2 = 1.0;
    double Pi = M_PI;

    gsl_function F2;
    F2.function = &Theta;
    F2.params = &alpha2;

    gsl_integration_qags (&F2, 0, Pi, 1e-5, 1e-7, 1000, w2, &result2, &error2);

    printf ("W2 result          = %.18f\n", result2);
    printf ("W2 expected result = %.18f\n", expected2);
    printf ("W2 estimated error = %.18f\n", error2);
    printf ("W2 actual error    = %.18f\n", result2 - expected2);
    printf ("W2 intervals       = %d\n", w2->size);

    gsl_integration_workspace_free (w2);

    return pow(r,2.0)*result2;
}


int main (void)
{
    gsl_integration_workspace * w1 = gsl_integration_workspace_alloc (1000);

    double result1, error1;
    double alpha1 = 1.0;
    double radius = 2.0;
    double expected1 = (4.0*M_PI*pow(radius,3.0)) / 3.0;

    gsl_function F1;
    F1.function = &R;
    F1.params = &alpha1;

    gsl_integration_qags (&F1, 0, radius, 1e-5, 1e-7, 1000, w1,
&result1, &error1);

    printf ("\nW1 result          = %.18f\n", result1);
    printf ("W1 expected result = %.18f\n", expected1);
    printf ("W1 estimated error = %.18f\n", error1);
    printf ("W1 actual error    = %.18f\n", result1 - expected1);
    printf ("W1 intervals       = %d\n", w1->size);

    gsl_integration_workspace_free (w1);

    return 0;
}

Cheers,
Kevin

Re: multidimensional numerical integration

[Robert G. Brown]

On Mon, 22 Oct 2007, Kevin Jackman wrote:

> Hi all,
>
> A couple of questions followed by some useful information for others (in my
> opinion):
>
> Q: I realize "epsabs" and "epsrel" control the desired accuracy of numerical
> integration, but how? Mainly, how do they each affect the GSL routines, and
> how do they differ from each other? How about on a Quadpack level? I think
> the GSL manual is kind of vague on how these parameters individually
> control numerical integration convergence. (No hard feelings, ok?)
>
>
> Some Useful Info:
>
> Looking through the GSL documentation, archives, and source code, I soon
> discovered that GSL does not have specific routines for multidimensional
> numerical integration, except using Monte Carlo methods.  I saw some
> efforts to initially incorporate cubature routines:
>
> http://sourceware.org/ml/gsl-discuss/2007-q2/msg00046.html
> http://www.cygwin.com/ml/gsl-discuss/2005-q2/msg00020.html

I'm guessing that this is one of the source packages from those
discussions:

  http://www.phy.duke.edu/~rgb/libcubature-0.1.0.tgz

I actually built rpms (source and binary) for it as well, if that helps.
The algorithm works quite well for certain classes of integrand, but as
always, there are functional patterns that will fool it and/or give
poor performance.  In those cases it may be wise to further subdivide
the volume of integration.  It is, however, reasonably efficient.

    rgb

>
> However, as cubature has not been formally incorporated into GSL, and that
> topic is far from my specialty, I moved on. Below is a little tidbit of code
> which demonstrates how to recursively integrate using 1D QAGS routines and
> thereby yield an effective 3D integral. The code integrates the volume of a
> sphere in spherical coordinates. I realize this method can be time and computer
> resource consuming for complex integrals, but I have yet to see a better GSL
> alternative suited for my purposes. (Monte Carlo is not good for my needs,
> long explanation).  Also, the absolute error estimation may be lost by doing it
> this way, though I'm not sure.
>
> Aside from telling me which hand I should wish in and which hand I should...,
> I'd like to see something like this incorporated into the GSL manual, so others
> don't have to reinvent the wheel. You are welcome to use my code, modify it,
> etc... to make it up to your manual standards. No props necessary even.
>
>
> compile and build with something like:
> gcc -o simpleQAGS_3D simpleQAGS_3D.c -L/usr/lib -I/usr/include -lgsl
> -lgslcblas -lm
>
>
> simpleQAGS_3D.c:
>
> #include <stdio.h>
> #include <math.h>
> #include <gsl/gsl_integration.h>
>
> double Phi (double phi, void * params) {
>
>    return 1.0;
> }
>
>
> double Theta (double theta, void * params) {
>
>    gsl_integration_workspace * w3 = gsl_integration_workspace_alloc (1000);
>
>    double result3, error3;
>    double expected3 = 2.0*M_PI;
>    double alpha3 = 1.0;
>    double TwoPi = 2.0*M_PI ;
>
>
>    gsl_function F3;
>    F3.function = &Phi;
>    F3.params = &alpha3;
>
>    gsl_integration_qags (&F3, 0, TwoPi, 1e-5, 1e-7, 1000, w3,
> &result3, &error3);
>
>    printf ("W3 result          = %.18f\n", result3);
>    printf ("W3 expected result = %.18f\n", expected3);
>    printf ("W3 estimated error = %.18f\n", error3);
>    printf ("W3 actual error    = %.18f\n", result3 - expected3);
>    printf ("W3 intervals       = %d\n", w3->size);
>
>    gsl_integration_workspace_free (w3);
>
>    return sin(theta)*result3;
> }
>
>
> double R (double r, void * params) {
>
>    gsl_integration_workspace * w2 = gsl_integration_workspace_alloc (1000);
>
>    double result2, error2;
>    double expected2 = M_PI;
>    double alpha2 = 1.0;
>    double Pi = M_PI;
>
>    gsl_function F2;
>    F2.function = &Theta;
>    F2.params = &alpha2;
>
>    gsl_integration_qags (&F2, 0, Pi, 1e-5, 1e-7, 1000, w2, &result2, &error2);
>
>    printf ("W2 result          = %.18f\n", result2);
>    printf ("W2 expected result = %.18f\n", expected2);
>    printf ("W2 estimated error = %.18f\n", error2);
>    printf ("W2 actual error    = %.18f\n", result2 - expected2);
>    printf ("W2 intervals       = %d\n", w2->size);
>
>    gsl_integration_workspace_free (w2);
>
>    return pow(r,2.0)*result2;
> }
>
>
> int main (void)
> {
>    gsl_integration_workspace * w1 = gsl_integration_workspace_alloc (1000);
>
>    double result1, error1;
>    double alpha1 = 1.0;
>    double radius = 2.0;
>    double expected1 = (4.0*M_PI*pow(radius,3.0)) / 3.0;
>
>    gsl_function F1;
>    F1.function = &R;
>    F1.params = &alpha1;
>
>    gsl_integration_qags (&F1, 0, radius, 1e-5, 1e-7, 1000, w1,
> &result1, &error1);
>
>    printf ("\nW1 result          = %.18f\n", result1);
>    printf ("W1 expected result = %.18f\n", expected1);
>    printf ("W1 estimated error = %.18f\n", error1);
>    printf ("W1 actual error    = %.18f\n", result1 - expected1);
>    printf ("W1 intervals       = %d\n", w1->size);
>
>    gsl_integration_workspace_free (w1);
>
>    return 0;
> }
>
> Cheers,
> Kevin
>

-- 
Robert G. Brown
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone(cell): 1-919-280-8443
Web: http://www.phy.duke.edu/~rgb
Lulu Bookstore: http://stores.lulu.com/store.php?fAcctID=877977

Re: multidimensional numerical integration

[Brian Gough]

At Mon, 22 Oct 2007 12:59:16 -0600,
Kevin Jackman wrote:
> Q: I realize "epsabs" and "epsrel" control the desired accuracy of numerical
> integration, but how? Mainly, how do they each affect the GSL routines, and
> how do they differ from each other? How about on a Quadpack level? I think
> the GSL manual is kind of vague on how these parameters individually
> control numerical integration convergence. (No hard feelings, ok?)

Hello, 

There is a definition of epsabs and epsrel in the introductory section
of the Integration chapter but maybe it is not clear enough.  

The original book on quadpack has a complete explanation of them. It
is worth getting for serious use, the GSL routines are a faithful
implementation of the algorithms in the book. 

I think the book also has a section on error bounds for computing
multidimensional integrals recursively.

-- 
Brian Gough

GNU Scientific Library -
http://www.gnu.org/software/gsl/

Re: multidimensional numerical integration

[Kevin Jackman]

Thanks for the info about the quadpack book. I should definitely get that.
Any suggestions on where? You did hit the nail right on the head Brian.
The definition of epsabs and epsrel in the introductory section of the
Integration chapter just seemed a little vague to me.

Kevin

On 10/23/07, Brian Gough <bjg@gnu.org> wrote:
> At Mon, 22 Oct 2007 12:59:16 -0600,
> Kevin Jackman wrote:
> > Q: I realize "epsabs" and "epsrel" control the desired accuracy of numerical
> > integration, but how? Mainly, how do they each affect the GSL routines, and
> > how do they differ from each other? How about on a Quadpack level? I think
> > the GSL manual is kind of vague on how these parameters individually
> > control numerical integration convergence. (No hard feelings, ok?)
>
> Hello,
>
> There is a definition of epsabs and epsrel in the introductory section
> of the Integration chapter but maybe it is not clear enough.
>
> The original book on quadpack has a complete explanation of them. It
> is worth getting for serious use, the GSL routines are a faithful
> implementation of the algorithms in the book.
>
> I think the book also has a section on error bounds for computing
> multidimensional integrals recursively.
>
> --
> Brian Gough
>
> GNU Scientific Library -
> http://www.gnu.org/software/gsl/
>

Re: multidimensional numerical integration

[Brian Gough]

At Tue, 23 Oct 2007 17:05:48 -0600,
Kevin Jackman wrote:
> 
> Thanks for the info about the quadpack book. I should definitely get that.
> Any suggestions on where? 

It's out of print so you'd have to get it via a library or abe.com

-- 
Brian Gough

Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/