Hello,

I am having trouble setting up the needed function/functor and passing the arguments from my main{} into the template Trapzd with the struct Quadrature. Could anyone nudge me in the right direction with how to correctly pass functions/functors and the needed arguments from the main{} into quadrature.h?

Thanks
Matt D.

Hello,

I am having trouble setting up the needed function/functor...

Here's how to use qsimp and qtrap with a double precision template function:

//
// Demo of qsimp and qtrap
//
// davekw7x
//
#include "../code/nr3.h"
#include "../code/quadrature.h"
//
// A convenient function for x-square
//
inline Doub sqr(Doub x)
{
return x * x;
}

//
// Test function to integrate
//
// x^4 * sin(x) - 2 * x^2 * sin(x)
//
//
// The template for qsimp and qtrap requires a function with a single argument
//
Doub f(Doub x)
{
Doub temp = sqr(x);
return temp * (temp - 2.0) * sin(x);
}

//
// For error calculations:
// The analytic anti-derivative of f(x)
//
Doub fint(Doub x)
{
return 4.0*x*(sqr(x)-7.0)*sin(x)-(sqr(x)*sqr(x)-14.0*sqr(x)+28.0)*cos(x);
}

int main()
{
const Doub pio2 = 1.570796326794896619;
Doub value;
Doub err;

//
// Will integrate from a to b
//
Doub a = 0.0, b = pio2;

Doub actual = fint(b) - fint(a);

//
// Create a Trapzd object with our function.
//
// The template is a function that has a Doub return value.
// The function has a single argument that is a Doub.
//
Trapzd <Doub(Doub)> trapzd(f, a, b);

value = qtrap(f, a, b);

err = value - actual;
cout << scientific << setprecision(15);
cout << "From qtrap, value = " << value << ", error = " << err << endl;

value = qsimp(f, a, b);

err = value - actual;
cout << "From qsimp, value = " << value << ", error = " << err << endl;

return 0;
}


Output

From qtrap, value = -4.791588101032077e-01, error = 3.986089236462931e-12
From qsimp, value = -4.791588101065034e-01, error = 6.903366767119223e-13


If you really need a Trapzd object, you can set it up the same way with the same template.


If your function has parameters other than the argument (the value of x), you can use a functor.


Regards,

Dave
Footnote: Since you asked, I showed how to create a Trapzd object here just in case you need one some day, but I didn't actually use it in this program.

The functions qtrap and qsimp just need your function, and they "do their own thing" with a Trapzd object of their own creation.

So, simply to use qtrap and/or qsimp you can delete the line that creates trapzd. I'm sorry if the extra stuff created confusion.

Thank you very much. The fog has lifted.

Matt D.

No confusion with the extra stuff. Actually, seeing how the Trapzd line was used gave me great insight into how to pass the arguments when using structs and templates.

Thanks again for you help.

Matt

This is a quite old thread but I do have an actual question to the integration.

My problem is that I do not declare the function in a functor but as a class method and that does not work! :confused:


Doub calculatingCrossSection::f(Doub x)
{
double t = dParameters[0][0];
double F1 = dParameters[0][1];
double F2 = dParameters[0][2];
double S = dParameters[0][3];
double I = dIonizationEnergy[0];

return S * F1 *
( 1/pow(1+x/I,3) + 1/pow(t-x/I,3) - 1/pow((1+x/I)*(t-x/I),1.5) )
+ S * F2 *
( 1/pow(1+x/I,2) + 1/pow(t-x/I,2) - 1/( (1+x/I) * (t-x/I)) );

}


And later in the program I simply try to insert the function f as a class method:


int main()
{

calculatingCrossSection CrossSection(100)

Doub a = 0;
Doub b = 50;
Doub value = 0;

value = qsimp(CrossSection.f, a, b)

return 0;
}


But that does not work. I get the error:
no matching function for call to ‘qsimp(<unresolved overloaded function type>, Doub&, Doub&)’

The function is okay, I checked it with a "normal" functor but it is not possible to use it as a class method. Or do I do something wrong?

Best regards,
Quantum.

...
My problem is that I do not declare the function in a functor but as a class method and that does not work!...
Here's the thing: A non-static member function of a class has an extra (implicit) parameter, this. The result is that the signature of a class member function is different from a stand-alone function defined with the same (explicit) parameters and same return type.

That is, the qsimp constructor requires a function with return type Doub and a single parameter with type Doub. The class function does not match this signature even though it looks the same to use mere humans.


You can try something like the following: Make a functor that acts as a wrapper for your class function:


//
// Your class definition for calculatingCrossSection goes here
//

class cCSWrapper
{

private:
calculatingCrossSection cCS;

public:
cCSWrapper(const calculatingCrossSection & cr): cCS(cr){}

//The functor that we will use in the qsimp constructor
Doub operator ()(Doub & x) { return cCS.f(x); }

};

int main()
{

// Your constructor
calculatingCrossSection CrossSection(100);

//
// Whatever you need to do with the CrossSection object
//

// Instantiate the wrapper functor object
cCSWrapper csWrapper(CrossSection);

Doub a = 0;
Doub b = 50;

Doub value = qsimp(csWrapper, a, b);
.
.
.



Regards,

Dave

Every time I see someone having to deal with stuff like this, it makes me happy all over again to have moved on to managed languages.

Every time I see someone having to deal with stuff like this, it makes me happy all over again to have moved on to managed languages.

:confused: :confused:

Uh, sorry Dave but I get a lot of errors:


/quadrature.h:1: error: redefinition of ‘struct Quadrature’
/quadrature.h:1: error: previous definition of ‘struct Quadrature’
/quadrature.h:6: error: redefinition of ‘struct Trapzd<T>’
/quadrature.h:6: error: previous definition of ‘struct Trapzd<T>’
/quadrature.h:30: error: redefinition of ‘template<class T> Doub qtrap(T&, Doub, Doub, Doub)’

/quadrature.h:30: error: ‘template<class T> Doub qtrap(T&, Doub, Doub, Doub)’ previously declared here

/quadrature.h:44: error: redefinition of ‘template<class T> Doub qsimp(T&, Doub, Doub, Doub)’

/quadrature.h:44: error: ‘template<class T> Doub qsimp(T&, Doub, Doub, Doub)’ previously declared here

/quadrature.h:60: error: redefinition of ‘struct Midpnt<T>’
/quadrature.h:60: error: previous definition of ‘struct Midpnt<T>’
/quadrature.h:91: error: redefinition of ‘struct Midinf<T>’
/quadrature.h:91: error: previous definition of ‘struct Midinf<T>’
/quadrature.h:102: error: redefinition of ‘struct Midsql<T>’
/quadrature.h:102: error: previous definition of ‘struct Midsql<T>’
/quadrature.h:114: error: redefinition of ‘struct Midsqu<T>’
/quadrature.h:114: error: previous definition of ‘struct Midsqu<T>’
/quadrature.h:126: error: redefinition of ‘struct Midexp<T>’
/quadrature.h:126: error: previous definition of ‘struct Midexp<T>’


What did I do wrong?

Here my Code:


using namespace std;

int main()
{

calculatingCrossSection CrossSection(200);

cCSWrapper cCSWrapper(CrossSection);

Doub a = 0;
Doub b = 100;

Doub value = qsimp(cCSWrapper, a, b);


return 0;

}


And the wrapper

class cCSWrapper
{
public:
calculatingCrossSection cCS;

cCSWrapper(const calculatingCrossSection & cr) : cCS(cr){}

Doub operator()(Doub & x) { return cCS.f(x); }
};


And the part of my class functor:

double calculatingCrossSection::f(double x)
{
double t = dParameters[i][0];
double F1 = dParameters[i][1];
double F2 = dParameters[i][2];
double S = dParameters[i][3]; // Einheit: Angstrom
double I = dIonizationEnergy[i];

return S * F1 *
( 1/pow(1+x/I,3) + 1/pow(t-x/I,3) - 1/pow((1+x/I)*(t-x/I),1.5) )
+ S * F2 *
( 1/pow(1+x/I,2) + 1/pow(t-x/I,2) - 1/( (1+x/I) * (t-x/I)) );
}

...lot of errors:


/quadrature.h:1: error: redefinition of ‘struct Quadrature’
/quadrature.h:1: error: previous definition of ‘struct Quadrature’
.
.
.


The Numerical Recipes distribution has been written in such a way that all source code must be compiled at the same time (in practical terms that means that your source code must be all together in the same file).

In other words, it is not possible (with unmodified NR files) to include quadrature.h from separate C++ source files in any given project.


From author Bill Press in this thread: http://www.nr.com/forum/showthread.php?t=1260

"2. The majority of the NR3 .h files contain (global) function definitions, not just declarations. These are definitely intended to be included only once and in a single compilation unit."


Regards,

Dave

Hi Dave,

one of the main problems is that my function (and thus the functor) has several variables.
To make an example:
Let's integrate the function f(x) = x^2 with the functor
Doub f (Doub x)
{
return x * x;
}

But in my program the function also depens on several variables, which are determined previously.

Something like: f(x) = A * x^2
First the program determines A.
And then you can declare the functor
Doub f (Doub x)
{
return A * x * x;
}

So it is not possible for me to simply declare a functor ahead of the main function, becaus my functor contains variables which are to be determined later in the program. :(

Best regards,
Quantum.

...my function (and thus the functor) has several variables....
Something like: f(x) = A * x^2

A big reason to use a functor instead of a simple function is that you can create a constructor that accepts parameters. The functor will have the proper signature to satisfy qsimp or qtrap and will let the user supply the information required to calculate the function values.

Using your example:


//
// Demo of a functor with a user-supplied parameter for
// use with qsimp and qtrap
//
// davekw7x
//
#include "../code/nr3.h" // Or whatever path you have to the nr3 stuff
#include "../code/quadrature.h"


//
// Test function to integrate: alpha * x^2, where the parameter alpha
// is supplied by the user
//
struct Ftor {
Doub alpha; // The parameter

// The function that will be integrated
Doub operator () (const Doub & x) const
{
return alpha * x * x;
}

//
// The constructor takes an argument that lets the user
// set the value of alpha
//
Ftor(Doub a):alpha(a) { }
};

// Analytic anti-derivative of test function
Doub fint(Doub alpha, Doub x)
{
return alpha * x * x * x / 3.0;
}

int main()
{
Doub alpha, a, b;
Doub approx;
Doub actual;
Doub err;
cout << "Will integrate alpha * x^2 over an interval." << endl << endl
<< "Enter values for interval end points and alpha: " << endl;
cin >> a >> b >> alpha;
cout << endl;

if (!cin) {
cout << "Invalid input." << endl;
return 0;
}

cout << "Will compute the approximate integral from "
<< a << " to " << b << " with alpha = " << alpha << endl;

// Create the functor object with user's value of alpha
Ftor f(alpha);

cout << scientific << setprecision(14);

actual = fint(alpha, b) - fint(alpha, a);
cout << "Actual value of integral is "
<< setw(10) << actual << endl;

approx = qsimp(f, a, b);
err = approx - actual;
cout << "From qsimp, value = " << approx
<< ", absolute error = " << err << endl;

approx = qtrap(f, a, b);
err = approx - actual;
cout << "From qtrap, value = " << approx
<< ", absolute error = " << err << endl;


return 0;
}


Here's a run:

Will integrate alpha * x^2 over an interval.

Enter values for interval end points and alpha:
0 100 1

Will compute the approximate integral from 0 to 100 with alpha = 1
Actual value of integral is 3.33333333333333e+05
From qsimp, value = 3.33333333333333e+05, absolute error = 0.00000000000000e+00
From qtrap, value = 3.33333333343035e+05, absolute error = 9.70129622146487e-06


Regards,

Dave

Ah,
thank you so much Dave.

Now it works so well!

Thanks!

Best regards,
Quantum.