help-octave
[Top][All Lists]

## Re: Using quad() for multidimensional integration

 From: Ted Harding Subject: Re: Using quad() for multidimensional integration Date: Thu, 17 Aug 1995 10:09:47 +0200 (BST)

( Re Message From: Vinayak Dutt )
>
> I am trying to use Octave to do multidimensional integration using the
>
>
> (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.
>
>
> % and
>
> function f2 = func1(x)
>   y1=0;
>   y2 = f(x);
> endfunction
>
> % and
>
> function f3 = func2(y)
>   z1 = 0;
>   z2 = g(y);
> 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?
>
>
> 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
fprintf(stdout,"U(0,%4.2f) = %4.2f\n",x,u);   ## printout on exit
endfunction
-----------------------------------------------------------------
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.