Re: Using quad() for multidimensional integration

[Ted Harding]

( Re Message From: Vinayak Dutt )
> 
> I just found that even that simple double integral does not work OK.
> 
> But I still wonder why would quad run out of stack space on simple
> double integrals (forget triple integral). I don't see that many levels
> of recursion here?
> 
Nor do I. But they happen. Or at least there is a totally unnecessary
number of calls to the function being integrated. Take the simple example

function u=U(x)
fprintf(stdout,"U(%4.2f)\n",x);
u=0.5*x*x;
endfunction

and do " quad('U',0,1) " and you get (packed onto three lines):

U(0.50) U(0.01) U(0.99) U(0.07) U(0.93) U(0.16) U(0.84) U(0.28) U(0.72)
U(0.43) U(0.57) U(0.00) U(1.00) U(0.03) U(0.97) U(0.11) U(0.89) U(0.22)
U(0.78) U(0.35) U(0.65) ans = 0.16667

i.e. 21 calls in order to integrate a quadratic! Now (if octave used, like
matlab, Simpson's rule) Simpson's rule is theoretically EXACT IN ONE PASS
for a quadratic. That is to say that if

   f(x) = a + b*x + c*x^2

then   \int_0^1 f(x) dx = [f(0) + 4*f(0.5) + f(1)]/6   exactly.

Of course "quad" wouldn't know it was quadratic, but on splitting the
range into 2 halves (which requires 2 further calls to U) and combining
the integrals (also exact) over the two halves, and comparing the second
result with the first, convergence would be exact. So just 5 calls are
required.

[NOTE: when I run this in matlab, U is evaluated at 4 points only, namely
       0, 0.25, 0.125 and 0.625. Hmm]

Octave "help -i quad" says the integration is done using "Quadpack".
I have a feeling I would like to know more about this "Quadpack".

I don't have the octave source (I downloaded the binary only).

Does anyone out there know how octave does "quad"?

Ted.                                    (address@hidden)