Hi!

I have some problem with using an adaptive 3D integration method. In the book "Numerical Recipes in C++" there is an example how to make a 3D integral from 1D integrals (page 169) but the 1D integral is not adaptive (just a simple 10 point gaus). When I am trying to use the example with an adaptive 1D quadrature (e.g. romberg, simpson, trapezoid) it just doesn't work! Can anybody tell me what is the problem?

Zoli

Hi Zoli,

The problem is that qromb and hence trapzd are called recursively. The routine trapzd contains the static variable s, so it can't be called properly recursively. (A good compiler should have warned you.) A simple way to fix this is to make s an argument of trapzd (non-const). It's easiest just to make your own version of trapzd with a new name. You'll have to change the call to the function in qromb. This is actually a little tricky to get right - you might want to look at the Fortran 90 version (online), where this is done.

Good luck.
Saul Teukolsky