Hello,
in 4.8 Multidimensional integrals, how can I easily modify the code in "quad3d" to take functors in the case of three demensional integrals?

...how can I easily...

Well, "easily," like "beauty," very well may be in the eye of the beholder...

The function quad3d has five arguments that are functions. The first one is already defined by a template; you have to work on the others.

Starting from the bottom-up, suppose we want to calculate the integral of the function "r^2" over a sphere with a given radius using quad3d.


This is the problem is illustrated by xquad3d in the NR2 CPP distribution. The "old-fashioned" way of NR2 required us to store the radius in a global variable that was accessed by the individual functions used for calculating y and z limits. (And quad3d.h in the NR3 distribution does it the same way.)

Now, instead of using functions, we can use functors in a more modern approach to C++ and eliminate all requirements for any of those dreaded and despised global variable thingies.

The source for a test program could look like:

//
// Functors used with quad3d
//
// davekw7x
//
// The NR3 Distribution files
#include "../code/nr3.h"
#include "../code/qgaus.h"

// A special templatized version of NR3 quad3d.h
#include "./quad3d_t.h"

// The following is not actually necessary, since it is in nr3.h, but...
using namespace std;

//
// r^2 = x^2 + y^2 + z^2: This is the function being integrated
// over a spherical surface.
//

//
// For this example, it could be defined as a function,
// but I show that it can be a functor
//
struct Ftor
{
Doub operator() (const Doub x, const Doub y, const Doub z)
{
return x*x + y*y + z*z;
}
};

//
// Use functors to pass functions to a templated version of
// quad3d. This way the radius can be incorporated in the
// calculations without making it a global variable
//
// For a given value of x, this calculates
// the lower limit for integration over y
//
struct Fty1
{
Doub radius;

Fty1(Doub r):radius(r) {}

Doub operator() (const Doub x) const
{
return -sqrt(radius*radius - x*x);
}
};

//
//For a given value of x, this calculates
// the upper limit for integration over y
//

struct Fty2
{
Doub radius;

Fty2(Doub r):radius(r) {}

Doub operator() (const Doub x) const
{
return +sqrt(radius*radius - x*x);
}
};

//
//For given values of x and y, this calculates
// the lower limit for integration over z
//
struct Ftz1
{
Doub radius;
Ftz1(Doub r):radius(r) {}
Doub operator() (const Doub x, const Doub y) const
{
return -sqrt(radius*radius - x*x - y*y);
}
};

//
//For given values of x and y, this calculates
// the upper limit for integration over z
//
struct Ftz2
{
Doub radius;
Ftz2(Doub r):radius(r) {}
Doub operator() (const Doub x, const Doub y) const
{
return +sqrt(radius*radius - x*x - y*y);
}
};

// Some <cmath> headers have M_PI defined; some do not
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

//
// The analytical value of this integral is 4*pi*(radius^5)/5
// Calculate it here for comparison purposes.
//
Doub actual(Doub r)
{
return 4.0 * M_PI * pow(r, 5.0) / 5.0;
}

int main(void)
{
//
// Use this many values of radius up to a value of 1.0
//
int num_values = 10;

cout << "Integration of r^2 over a spherical volume using quad3d_t"
<< endl << endl;
cout << setw(11) << "Radius" << setw(15) << "Approx";
cout << setw(15) << "Actual" << endl;
cout << scientific << setprecision(6);
for (int i = 1; i <= num_values; i++) {
Doub radius = Doub(i)/num_values;

// Limits of integration over x
Doub xmin = -radius;
Doub xmax = radius;

// Instantiate functors
Ftor func;
Fty1 y1(radius);
Fty2 y2(radius);
Ftz1 z1(radius);
Ftz2 z2(radius);

// Now perform the 3D integration
Doub result = quad3d_t(func, xmin, xmax, y1, y2, z1, z2);

cout << setw(15) << xmax
<< setw(15) << result
<< setw(15) << actual(radius) << endl;
}
return 0;
}


Now, we have to work on the quad3d stuff. I made a copy of quad3d.h and called it quad3d_t.h so that the original can not be corrupted by any misguided efforts on my part.

Start with the quad3d function itself. Originally it had one argument that was a template and four that were functions. Now it has to have five templated arguments, corresponding to our five functors:


// The original
//template <class T>
//Doub quad3d(T &func, const Doub x1, const Doub x2, Doub y1(Doub), Doub y2(Doub),
// Doub z1(Doub, Doub), Doub z2(Doub, Doub))

//The templated version
template <class T, class Ty1, class Ty2, class Tz1, class Tz2>
Doub quad3d_t (T &fun, const Doub x1, const Doub x2, Ty1 &y1, Ty2 &y2,
Tz1 &z1, Tz2 &z2)


Now, we have to follow through with the template stuff all the way down through the lower functions.

So the beginning of quad3d_t can look like:


template <class T, class Ty1, class Ty2, class Tz1, class Tz2>
Doub quad3d_t (T &fun, const Doub x1, const Doub x2, Ty1 &y1, Ty2 &y2,
Tz1 &z1, Tz2 &z2)
{
T func(fun);
NRf1 <T,Ty1,Ty2,Tz1,Tz2> f1(y1,y2,z1,z2);


See how it goes? The arguments to quad3d_t that are defined by templates have to propagate to the lower ones as templates.

So the NRf1 becomes a templated class, and the beginning of NRf1 can now look like:


template <class T, class Ty1, class Ty2, class Tz1, class Tz2>
struct NRf1 {
Ty1 y1;
Ty2 y2;
NRf2 <T, Tz1,Tz2> f2;
NRf1(Ty1 yy1, Ty2 yy2, Tz1 z1, Tz2 z2):y1(yy1),y2(yy2),f2(z1,z2){}


Get the pattern? Now NRf2 has to be templated, and its beginning can look like


template <class T, class Tz1, class Tz2>
struct NRf2 {
NRf3 <T> f3;
Tz1 z1;
Tz2 z2;
NRf2(Tz1 zz1, Tz2 zz2) : z1(zz1), z2(zz2) {}


Now are you getting the pattern? Here's a possible start of NRf3 (and this is the last change):

template <class T>
struct NRf3 {
T func3d;


With these changes applied to quad3d_t, after compiling with GNU g++, the output of the test program looks like:


Integration of r^2 over a spherical volume using quad3d_t

Radius Approx Actual
1.000000e-01 2.515219e-05 2.513274e-05
2.000000e-01 8.048699e-04 8.042477e-04
3.000000e-01 6.111981e-03 6.107256e-03
4.000000e-01 2.575584e-02 2.573593e-02
5.000000e-01 7.860058e-02 7.853982e-02
6.000000e-01 1.955834e-01 1.954322e-01
7.000000e-01 4.227328e-01 4.224060e-01
8.000000e-01 8.241868e-01 8.235497e-01
9.000000e-01 1.485211e+00 1.484063e+00
1.000000e+00 2.515219e+00 2.513274e+00


Regards,

Dave

First, let me say that I think this qualifies as a bug report. The code for quad3d.h in the book is templated on func and the intention was that the user could supply func either as a function or a functor. However, the code works only for a function and not for a functor.

The code in the previous post works for a functor (and templates the functions that supply the limits so that they can be functors as well.) However, it doesn't work for a function. Writing a single piece of code that works for either functions or functors (or any combination, if we consider all 5 function inputs) is a little tricky. As far as I can tell, the only way to do it is to make func3d in NRf3 be a pointer. Here is the complete code (normally we don't post copyrighted code in the forum, but this seems like a good time to make an exception):

template <class T>
struct NRf3 {
Doub xsav,ysav;
T *func3d;
Doub operator()(const Doub z)
{
return (*func3d)(xsav,ysav,z);
}
};
template <class T, class Z1, class Z2>
struct NRf2 {
NRf3<T> f3;
Z1 &z1;
Z2 &z2;
NRf2(Z1 &zz1, Z2 &zz2) : z1(zz1), z2(zz2) {}
Doub operator()(const Doub y)
{
f3.ysav=y;
return qgaus(f3,z1(f3.xsav,y),z2(f3.xsav,y));
}
};
template <class T, class Y1, class Y2, class Z1, class Z2>
struct NRf1 {
Y1 &y1;
Y2 &y2;
NRf2<T,Z1,Z2> f2;
NRf1(Y1 &yy1, Y2 &yy2, Z1 &z1, Z2 &z2) : y1(yy1),y2(yy2), f2(z1,z2) {}
Doub operator()(const Doub x)
{
f2.f3.xsav=x;
return qgaus(f2,y1(x),y2(x));
}
};
template <class T, class Y1, class Y2, class Z1, class Z2>
Doub quad3d(T &func, const Doub x1, const Doub x2, Y1 &y1, Y2 &y2, Z1 &z1,
Z2 &z2)
{
NRf1<T,Y1,Y2,Z1,Z2> f1(y1,y2,z1,z2);
f1.f2.f3.func3d=&func;
return qgaus(f1,x1,x2);
}

The code in the previous post...doesn't work for a function...

Saul: Thank you for the corrections. I had so much fun functorizing the limits functions it that I didn't think to test with unfunctorized arguments.

Regards,

Dave

thank you very much for your helps, i have benefit a lot

regards

worrygemini