[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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
> 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)
- Using quad() for multidimensional integration, Vinayak Dutt, 1995/08/16
- Re: Using quad() for multidimensional integration,
Ted Harding <=
- Re: Using quad() for multidimensional integration, Ted Harding, 1995/08/17
- Re: Using quad() for multidimensional integration, Vinayak Dutt, 1995/08/17
- Re: Using quad() for multidimensional integration, Vinayak Dutt, 1995/08/17
- Re: Using quad() for multidimensional integration, Ted Harding, 1995/08/17
- Re: Using quad() for multidimensional integration, Vinayak Dutt, 1995/08/17
- Re: Using quad() for multidimensional integration, Ted Harding, 1995/08/17
- Re: Using quad() for multidimensional integration, Jim Van Zandt, 1995/08/18
- Re: Using quad() for multidimensional integration, Ted Harding, 1995/08/18