> I need to solve multiple integrals (two-dimensional at the moment). I have
> read that quadpack does not work for multiple integrals in the present
> version. Is this still true? I have also read that one should be able to
> download quadpack source code, but I can not find the correct
> location/address. Perhaps I could fix the problem myself, if I had the
> source code and I knew what the problem was. Otherwise, are there any
> advise where I could find some other code that is working?
The source code for quadpack is in netlib.
The reason that QuadPack is difficult to use for multiple integrals
is that it is compiled Fortran, which is non-reentrant
(i.e. a function cannot call itself). The solution is to
use another quadpack function to compute the inner integral.
This, however, did not work with the version of quadpack
which came with the previous version of Ox.
Below is an example program which I recently used to
evaluate a bivariate integral.
Jurgen.
////////////////////////////////////////////////////////////////////////////////
//////////
#include <oxstd.h>
#include <quadpack.h>
output(const sFunc, const result, const abserr, const code)
{
println(sFunc, result, " abserr=", abserr, " code=", code);
}
bias1(const x)
{
return x / sqrt(cosh(x));
}
bias2point(const x1, const x2)
{
return (sqr(x1) - sqr(x2)) * (tanh(x1)/x1 + tanh(x2)/x2) /
sqrt(cosh(x1) * cosh(x2));
}
static decl s_x1;
bias2given1(const x2)
{
return bias2point(s_x1, x2);
}
bias2(const x1)
{
decl result, abserr, retval;
s_x1 = x1;
retval = QAGS(bias2given1, 0, x1, &result, &abserr);
return result;
}
main()
{
decl result, abserr, retval;
// univariate integral
retval = QAGI(bias1, 0, 1, &result, &abserr);
output("QAGI:", 1 - result / 2, abserr, retval);
// bivariate integral
retval = QAGI(bias2, 0, 1, &result, &abserr);
output("QAGI:", 1 - result / 8, abserr, retval);
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|