Re: dblquad or similar in Octave?

[David Bateman]

Jens Benecke wrote:
> Hi everybody,
>
> I need to perform a double quadrature in Octave. In Matlab, this is the
> code:
>
> function retval=induk(z1,z2,r1,th)
>     tol  = [];
>     retval = r1* dblquad(@(p1,p2) ...
>         (sin(p1).*sin(p2).*cos(th)+cos(p1).*cos(p2))./ ...
>         sqrt( (z2*sin(th)+cos(p2)*cos(th)-r1*cos(p1)).^2.+...
>               (sin(p2)-r1*sin(p1)).^2. +...
>               (z2*cos(th)-cos(p2)*sin(th)-z1).^2. ), ...
>    0, pi, 0, 2*pi, tol, @quadl)/(2.*pi)
> end
>
>
> I have searched the net but all I could find was several references
> to "missing functions" in Octave (compared to Matlab).
>
> I use Octave 2.1.71 under Linux (SuSE and Ubuntu). How would I perform this
> kind of quadrature with Octave?
>
> Thank you!
>
>
>   
Ideally the quad function should be made re-entrant, so that you could
nest one quad function inside another. Unfortunately, the quad function
is based on the fortran code quadpack and f77 doesn't define recursive
functions.

I suggest you use the quadl function instead as it is a pure m-file and
so can be nested. It is likely to be slow though. As an example consider

function z = fun2 (y)
  z = y
endfunction
function z = fun1 (x, f, s, e)
  z = x * quadl(f, s, e, x);
endfunction
quadl(@(x) fun1(x, @(y) fun2(y), 0, 1), 0, 1)

though you might need 2.9.9...

D.


quadl


-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary