Re: Using quad() for multidimensional integration

[Ted Harding]

( Re Message From: Vinayak Dutt )
> 
> I am trying to use Octave to do multidimensional integration using the
> quadrature method: quad(). I have a triple integral of type:
> 
> 
> (TeX Notation)
> 
> \int_{x=0}^{R} F(x)  \int_{y=0}^{f(x)} G(y) \int_{0}^{g(y)}  H(z) dz dy dx
> 
> I guess, th exact integral does not matter, except that its triple integral.
> 
> So I implement this by having the integral function itself involve another
> integration using quad() which again involves another integration.
> 
> a = quad('func1',x1,x2);
> 
> % and
> 
> function f2 = func1(x)
>   y1=0;
>   y2 = f(x);
>   f2 = F(x)*quad('func2',y1,y2);
> endfunction
> 
> % and
> 
> function f3 = func2(y)
>   z1 = 0;
>   z2 = g(y);
>   f3 = G(y)*quad('H',z1,z2);
> endfunction
> 
> 
> But when I do triple integral this way, I find that first two quadratures are 
> not
> evaulated at all. The functions func1 and func2 were evaluated only once. Only
> func2 was used to completely evaluate func1 for that single func1 evaluation.
> It seems the recursion in quad() is somehow killing the hiher level quad() 
> evaluation.
> 
> Could this be due to usage of some globals in quad() which make it 
> non-recursive?
> 
> Any comments?
> 
> I am using Octave-1.1.1 on SPARC-Solaris2.4 with gcc2.6.3.
> 
****************************************************************************
There do seem to be problems here. I tried a simpler task:

   \int_0^1 \int_0^x y dy dx   [Ans = 1/6].

First the function definitions (note the trace prints):--
-----------------------------------------------------------------
#script file
1;

function v=V(y)
fprintf(stdout,"V(%4.2f)\n",y);               ## printout on entry
v=y;
endfunction

function u = U(x)
fprintf(stdout,"U(%4.2f)\n",x);               ## printout on entry
u = quad('V',0,x);
fprintf(stdout,"U(0,%4.2f) = %4.2f\n",x,u);   ## printout on exit
endfunction
-----------------------------------------------------------------
Next, execute "quad('U',0,1)". Results:--
U(0.50)          <-----  "U" is entered here for first time
V(0.25)          <-----  "V" here for first time
V(0.01)
V(0.49)
....      [all "V"] 21 V's in all
V(0.39)
V(0.18)
V(0.32)
U(0,0.50) = 0.12  <----- first "U" completes here
V(0.01)
V(0.49)
V(0.03)
....      [all "V"] 86 V's in all
V(0.00)
V(0.00)
V(0.00)
V(0.00)
ans = 0           <----- end of computation

Therefore it appears that "U" is entered once only, for the bottom half of
the bisection of (0,1). This is similar to what Vinayak observed.

Secondly, however, I am seriously concerned that it should take so many
iterations of [what I presume is] Simpson's rule in order to get the
quadrature of "y dy" (in "V"). It should come out almost at once.

I am wondering, given the number of calls to V, if octave is running out
of stack space.

As a general comment, it would not be considered efficient practice to
recursively apply a crude quadrature rule for repeated integration, but
that is another story. The above example (unless I have made some subtle
error) should have worked.

Comments?

Ted.                                    (address@hidden)