Hi,

I am using qromb( ) to evaluate integrals.
I get following error message after execution of my code:

Numerical Recipes run-time error...
Too many steps in routine qromb
...now exiting to system...

Why do I get this message?

Regards

Hello,

Can you provide some more information?

For instance, what is the function that you are trying to integrate?

Thanks
Matt

Hello,

Can you provide some more information?

For instance, what is the function that you are trying to integrate?

Thanks
Matt

I am using first 10 Zernike Polynomials. But for our discussion let us take our function as a linear combination of 4th and 5th Zernike polynomials.
i.e.

w(x,y) = c4*z4(x,y) + c5*z5(x,y)

where c4, c5 are coefficients and
z4 = sqrt(3)*(2*x*x+2*y*y-1)
z5 = 2*sqrt(6)*x*y

And, x- limits of integration be -5.5 to -4.5 &
y limits be -0.5 to +0.5.

Regards

...

w(x,y) = c4*z4(x,y) + c5*z5(x,y)

where c4, c5 are coefficients and
z4 = sqrt(3)*(2*x*x+2*y*y-1)
z5 = 2*sqrt(6)*x*y

And, x- limits of integration be -5.5 to -4.5 &
y limits be -0.5 to +0.5.
...
Hmmm... The limits of integration are finite, the function looks pretty smooth, and there are no singularities.

So...

Do you have an example problem for which you know the result?

What values of c4 and c5 are you using?

What is the expected output?


Oh, yeah, I almost forgot: Show us your code.

Regards,

Dave

//
// Another demonstration of 2-D integration with Numerical Recipes functions.
// davekw7x
//
#include "../code/nr3.h"
//
// For gaussian integration, include qgaus.h
//
//#include "../code/qgaus.h"
// For romberg integration you need the following
// three #includes
#include "../code/interp_1d.h"
#include "../code/quadrature.h"
#include "../code/romberg.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);
}
};

// For romberg integration, use qromb.
// For gaussian integration, use qgaus
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));
return qromb(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);
return qromb(f1, x1, x2);
}
// End of template stuff
//--------------------------------------------------------------------
// The template stuff can be placed in a separate header and #included
// in your main program file.

//
// Test program for 2D integration
//
// 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;
}


Output:

_______________________

The answer is 42
_______________________


Regards,

Dave