Hello everyone,

I'm a novice user to NR3 and need help modifying Quad3d.h. I would like to evaluate a double integral, but the Quad3d evaluates only triple integrals. The following code is the changes I have so far of quad3d_t.h (a templated version), which I call quad2d_t.h. It does compile but there is an error in the result. I don't know what I did wrong:
//--------------------------------------------------------------------

template <class T>
struct NRf3 {
Doub xsav,ysav;
T *func2d;
Doub operator()(const Doub z)
{
return (*func2d)(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>
template <class T, class Y1, class Y2>

struct NRf1 {
Y1 &y1;
Y2 &y2;
//NRf2<T,Z1,Z2> f2;
NRf3<T> f3;


//NRf1(Y1 &yy1, Y2 &yy2, Z1 &z1, Z2 &z2) : y1(yy1),y2(yy2), f2(z1,z2) {}
NRf1(Y1 &yy1, Y2 &yy2) : y1(yy1),y2(yy2) {}

Doub operator()(const Doub x)
{
//f2.f3.xsav=x;
f3.xsav=x;

//return qgaus(f2,y1(x),y2(x));
return qgaus(f3,y1(x),y2(x));

}
};

//template <class T, class Y1, class Y2, class Z1, class Z2>
template <class T, class Y1, class Y2>

//Doub quad3d(T &func, const Doub x1, const Doub x2, Y1 &y1, Y2 &y2, Z1 &z1,
//Z2 &z2)
Doub quad2d(T &func, const Doub x1, const Doub x2, Y1 &y1, Y2 &y2)
{
//NRf1<T,Y1,Y2,Z1,Z2> f1(y1,y2,z1,z2);
NRf1<T,Y1,Y2> f1(y1,y2);

//f1.f2.f3.func2d=&func;
f1.f3.func2d=&func;


return qgaus(f1,x1,x2);
}
//--------------------------------------------------------------------

Thanks for any input,
Bob

...there is an error ...
The top function has two parameters, not three

#include "../code/nr3.h"
#include "../code/qgaus.h"

//--------------------------------------------------------------------
// After debugging/testing, you can put the template stuff in
// its own header

//--------------------------------------------------------------------
// Start of template stuff
// Derived from the quad3d_t in post number 3 of this thread:
//http://www.nr.com/forum/showthread.php?t=1663
//
template < class T >
struct NRf2
{
Doub xsav;
T *func2d;
Doub operator() (const Doub y) {
return (*func2d) (xsav, y);
}
};

template < class T, class Y1, class Y2 >
struct NRf1 {
Y1 & y1;
Y2 & y2;
NRf2 < T > f2;
NRf1(Y1 & yy1, Y2 & yy2):y1(yy1), y2(yy2) { }
Doub operator() (const Doub x) {
f2.xsav = x;
return qgaus(f2, y1(x), y2(x));
}
};

template < class T, class Y1, class Y2 >
Doub quad2d_t(T & func, const Doub x1, const Doub x2, Y1 & y1, Y2 & y2)
{
NRf1 < T, Y1, Y2 > f1(y1, y2);
f1.f2.func2d = &func;
return qgaus(f1, x1, x2);
}
// End of template stuff
//--------------------------------------------------------------------

// Test program for 2D integration
//
// Will integrate r**2 over a circular region.
//
// Here are the user functions for the test case
//
//
// r**2 = x**2 + y**2: This is the function being integrated.
//
// For this example, it could be defined as a function
// But I show how it can be a functor in case
// it needs parameters in addition to x and y
//

struct Ftor
{
Doub operator() (const Doub & x, const Doub & y) const
{
return x*x + y*y;
}
};

//
// Use functors to pass functions to quad2d.
//
// 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(const 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(const Doub & r):radius(r) {}

Doub operator() (const Doub &x) const
{
return +sqrt(radius*radius - x*x);
}
};
//
// Some compiler library math.h headers define M_PI, some do not
//
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif


//
// The integral of r**2 over a circle with a given radius
// is obtained easily with polar coordinates. Here's
// the analytical result.
//
Doub actual(Doub r)
{
return 2.0*M_PI * pow(r,4)/4.0;
}

int main()
{
Doub num_values = 10;
cout << "Integration of r^2 over a circle using quad2d_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;
Doub xmax = radius;
Doub xmin = -xmax;

// Instantiate functors
Ftor func;
Fty1 y1(radius);
Fty2 y2(radius);

// Calculate the integral for this radius
Doub result = quad2d_t(func, xmin, xmax, y1, y2);
cout << setw(15) << xmax
<< setw(15) << result
<< setw(15) << actual(radius) << endl;

}

return 0;
}


Output:

Integration of r^2 over a circle using quad2d_t

Radius Approx Actual
1.000000e-01 1.572260e-04 1.570796e-04
2.000000e-01 2.515617e-03 2.513274e-03
3.000000e-01 1.273531e-02 1.272345e-02
4.000000e-01 4.024987e-02 4.021239e-02
5.000000e-01 9.826628e-02 9.817477e-02
6.000000e-01 2.037649e-01 2.035752e-01
7.000000e-01 3.774997e-01 3.771482e-01
8.000000e-01 6.439979e-01 6.433982e-01
9.000000e-01 1.031560e+00 1.030599e+00
1.000000e+00 1.572260e+00 1.570796e+00


Regards,

Dave

It works! Thanks Dave, your the best.

I would like to pack the template part of your code into a class and then call quad2d_t from another class...

In order to enclose the templates into a class (say class A), I needed paste the qgauss code inside the class and make it static.

But when I try to use members of class A (namely quad2d_t) from a different class (say class B), I get errors while compiling.

Inside class B I am using function pointers and pass them to quad2d_t, e.g.
FuncPtr1=&B::function;
FuncPtr2=&B::yy1;
FuncPtr3=&B::yy2;

A integral;
result=integral.quad2d_t(FuncPtr1,xx,yy,FuncPtr2,F uncPtr3);

Error I get concerns calling qgauss from structs:

error: must use ‘.*’ or ‘->*’ to call pointer-to-member function in ‘((A::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>*)this)->A::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>::y1 (...)’

error: must use ‘.*’ or ‘->*’ to call pointer-to-member function in ‘((A::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>*)this)->A::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>::y2 (...)’

Any clue what should I do?

I would like to...

I can't figure out what you would like to do. I mean, the quad2d_t constructor doesn't take function pointers as arguments. It is takes references to functions or functors.

Here's where I start.

First of all: Here's my header file for my so-called quad2d_t class

//
// quad2d_t.h from davekw7x
//
template < class T >
struct NRf2
{
Doub xsav;
T *func2d;
Doub operator() (const Doub y) {
return (*func2d) (xsav, y);
}
};

template < class T, class Y1, class Y2 >
struct NRf1 {
Y1 & y1;
Y2 & y2;
NRf2 < T > f2;
NRf1(Y1 & yy1, Y2 & yy2):y1(yy1), y2(yy2) { }
Doub operator() (const Doub x) {
f2.xsav = x;
return qgaus(f2, y1(x), y2(x));
}
};

template < class T, class Y1, class Y2 >
Doub quad2d_t(T & func, const Doub x1, const Doub x2, Y1 & y1, Y2 & y2)
{
NRf1 < T, Y1, Y2 > f1(y1, y2);
f1.f2.func2d = &func;
return qgaus(f1, x1, x2);
}


Now, here's the main program file that uses functors from ExplorerBob's example:

//
// Test program for 2D integration
// davekw7x
//

#include "../code/nr3.h"
#include "quad2d_t.h"

// Here are the user functions for the test case
//
// This is the function that will be integrated
//
struct Ftor
{
Doub c4, c5;
Ftor(Doub a, Doub b): c4(a), c5(b) {}

Doub operator() (const Doub & x, const Doub & y) const
{
Doub z4 = sqrt(3.0)*(2.0*x*x+2.0*y*y-1.0);
Doub z5 = 2.0*sqrt(6.0)*x*y;
return c4*z4 + c5*z5;
}
};

//
// Use functors to pass functions to quad2d.
//
// For a given value of x, this calculates
// the lower limit for integration over y.
// Since it's a rectangular region, you just
// give a constructor with ymin and that's
// what you get.
//
struct Fty1
{
Doub ymin;
Fty1(Doub y):ymin(y){};

Doub operator() (const Doub & x) const
{
return ymin;
}
};


//
// For a given value of x, this calculates
// the upper limit for integration over y.
// Since it's a rectangular region, you just
// give a constructor with ymax and that's
// what you get.
//
struct Fty2
{
Doub ymax;
Fty2(Doub y): ymax(y){}
Doub operator() (const Doub &x) const
{
return ymax;
}
};
int main()
{
// Limits of integration over a rectangular region
Doub xmin = -5.5;
Doub xmax = -4.5;
Doub ymin = -0.5;
Doub ymax = +0.5;
//
// Constants used in the function to be integrated
//
Doub c4 = 0.4915279318776544;
Doub c5 = 0.4915279318776544;

// Instantiate the functors
Ftor func(c4, c5);
Fty1 y1(ymin);
Fty2 y2(ymax);

//
// Do the integration.
//
Doub result = quad2d_t(func, xmin, xmax, y1, y2);

cout << "The answer is " << result
// If you know actual, calculate it here
// << endl
// << "The actual value is " << actual(xmin, xmax, ymin, ymax)
<< endl;

return 0;
}


Here's the output:

_________________

The answer is 42
_________________

Now, can you tell us (better yet: show us) exactly what it is that you want to do?


Regards,

Dave

Thank you for your swift reply.

Please bear in mind that I have limited experience with templates, so perhaps what I am doing is completely wrong.

Here is an working example, where quad2d templates are inside a class:

#include <iostream>
#include <iomanip>
#include <cmath>

//http://www.nr.com/forum/showthread.php?t=1920
//Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919

using namespace std;

class integrator {
public:

integrator(){};
~integrator(){};

template <class T>
static double qgaus(T &func, const double a, const double b){

//here I implemented code for 96 point legendre quadrature integration, found in Abramovitz
//removed it so that the post is accepted and not too long
//basically the same as qgauss, just on 96 points instead of 5
}

template < class T >
struct NRf2
{
double xsav;
T *func2d;
double operator() (const double y) {
return (*func2d) (xsav, y);
}
};

template < class T, class Y1, class Y2 >
struct NRf1 {
Y1 & y1;
Y2 & y2;
NRf2 < T > f2;
NRf1(Y1 & yy1, Y2 & yy2):y1(yy1), y2(yy2) { }
double operator() (const double x) {
f2.xsav = x;
return qgaus(f2, y1(x), y2(x));
}
};

template < class T, class Y1, class Y2 >
double quad2d_t(T & func, const double x1, const double x2, Y1 & y1, Y2 & y2)
{
NRf1 < T, Y1, Y2 > f1(y1, y2);
f1.f2.func2d = &func;
return qgaus(f1, x1, x2);
}


};

static double eta;
static double M;

double yy1(double W){
double a=1-W/eta/M;
a=1-sqrt(a);
return M*M*eta*a*a;
}
double yy2(double W){
double a=1-W/eta/M;
a=1+sqrt(a);
return M*M*eta*a*a;
}
double funct(double x, double y){
return x*x+y*y;
}

int main(){

eta=1;
M=1;

double xmin=-1;
double xmax=1;

double (*FuncPtr1)(double, double)=&funct;
double (*FuncPtr2)(double)=&yy1;
double (*FuncPtr3)(double)=&yy2;


integrator integral;
double out=integral.quad2d_t(FuncPtr1,xmin,xmax,FuncPtr2, FuncPtr3);

cout<<"result = "<<out<<endl;

return 0;
}



Compiling (g++ this_code.cc) will not give an error and running the program will succeed.

Now if I try to access the class integrator from inside another class, where functions yy1, yy2 and funct are defined, the compiling fails. The code now is:


#include <iostream>
#include <iomanip>
#include <cmath>

//http://www.nr.com/forum/showthread.php?t=1920
//Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919

using namespace std;

class integrator {
public:

integrator(){};
~integrator(){};

template <class T>
static double qgaus(T &func, const double a, const double b){
// code for qgauss here
}

template < class T >
struct NRf2
{
double xsav;
T *func2d;
double operator() (const double y) {
return (*func2d) (xsav, y);
}
};

template < class T, class Y1, class Y2 >
struct NRf1 {
Y1 & y1;
Y2 & y2;
NRf2 < T > f2;
NRf1(Y1 & yy1, Y2 & yy2):y1(yy1), y2(yy2) { }
double operator() (const double x) {
f2.xsav = x;
return qgaus(f2, y1(x), y2(x));
}
};

template < class T, class Y1, class Y2 >
double quad2d_t(T & func, const double x1, const double x2, Y1 & y1, Y2 & y2)
{
NRf1 < T, Y1, Y2 > f1(y1, y2);
f1.f2.func2d = &func;
return qgaus(f1, x1, x2);
}


};

class B{
public:
B(){
//fictitious numbers at the moment
eta=1;
M=1;

FuncPtr1=&B::funct;
FuncPtr2=&B::yy1;
FuncPtr3=&B::yy2;

};
~B(){};

//fictitious functions
double yy1(double W){
double a=1-W/eta/M;
a=1-sqrt(a);
return M*M*eta*a*a;
}
double yy2(double W){
double a=1-W/eta/M;
a=1+sqrt(a);
return M*M*eta*a*a;
}
double funct(double x, double y){
return x*x+y*y;
}

double calculate(){
double xmin=-1;
double xmax=1;

double out=integral.quad2d_t(FuncPtr1,xmin,xmax,FuncPtr2, FuncPtr3);

return out;
}

private:
static double eta;
static double M;

double (B::*FuncPtr1)(double, double);
double (B::*FuncPtr2)(double);
double (B::*FuncPtr3)(double);

integrator integral;
};


int main(){
B myB;
cout<<"result = "<<myB.calculate()<<endl;

return 0;
}


The output from g++ is:
classy_quad2d.cc: In member function ‘double integrator::NRf1<T, Y1, Y2>::operator()(double) [with T = double (B::*)(double, double), Y1 = double (B::*)(double), Y2 = double (B::*)(double)]’:
classy_quad2d.cc:82: instantiated from ‘static double integrator::qgaus(T&, double, double) [with T = integrator::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>]’
classy_quad2d.cc:115: instantiated from ‘double integrator::quad2d_t(T&, double, double, Y1&, Y2&) [with T = double (B::*)(double, double), Y1 = double (B::*)(double), Y2 = double (B::*)(double)]’
classy_quad2d.cc:159: instantiated from here
classy_quad2d.cc:106: error: must use ‘.*’ or ‘->*’ to call pointer-to-member function in ‘((integrator::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>*)this)->integrator::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>::y1 (...)’
classy_quad2d.cc:106: error: must use ‘.*’ or ‘->*’ to call pointer-to-member function in ‘((integrator::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>*)this)->integrator::NRf1<double (B::*)(double, double), double (B::*)(double), double (B::*)(double)>::y2 (...)’

Since I can't explain myself, I hope that this example will show you, what I would like to do.

Best, etc
Matej