[Octave-bug-tracker] [bug #52074] Missing integration functions integral2, integral3

[Rik]

Follow-up Comment #10, bug #52074 (project octave):

I reviewed the patches, made some changes, and checked it all in here
(http://hg.savannah.gnu.org/hgweb/octave/rev/c723faa56ab4).

We can close this report now, or leave it open for other issues on integral2
and integral3.

One remaining issue is the fact that integral2,3 can accept function handles
for the limits of integration.  I'm not sure this is the optimal solution in
terms of computational efficiency, but one solution would be to apply a
masking function to the integrand itself.  As an example, see the help
documentation for integral2
(http://www.mathworks.com/help/matlab/ref/integral2.html).

There they perform a double integral of 


fun = @(x,y) 1./( sqrt(x + y) .* (1 + x + y).^2 )


over the region 0 <= x <= 1, 0 <= y <= 1 - x.

The answer they get for this integral is 0.2854.

If I do the same thing in Octave, which integrates over the rectangular
region, I obviously get a larger answer.


fun = @(x,y) 1./( sqrt(x + y) .* (1 + x + y).^2 );
integral2 (fun, 0, 1, 0, 1)
ans =  0.36953


But now let me define a masking function fmask which is 0 outside of the area
of integration, and one inside.


fmask = @(x,y) y <= (1 - x);


The masking function is implicitly vectorized so if x is a series of points,
fmask will be a series of mask values.

The new integrand with mask applied is


fun2 = @(x,y) fun (x,y) .* fmask (x,y)


The dot operator is used so I am always multiplying by either 0 or 1.

Finally,


integral2 (fun2, 0, 1, 0, 1)
ans =  0.28540


This might be a quick way to add support for limits which are function
handles.  Although, my guess is we would still like to find a way to just
integrate over the relevant area beause using the masking function is about 6X
slower than integrating over the rectangular region.  Although, it is only 4X
slower if I incorporate the fmask definition directly into the definition of a
new function.  That suggests that Octave's composition of functions is slow.

To be a completely generalizable, the mask function probably needs to be


fmask = @(x,y) y >= ymin(x,y) & y <= ymax(x,y)


Anyways, this is just an idea.


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?52074>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/