Re: Slow integration with dblquad

[Rik]

On 01/20/2013 10:00 AM, address@hidden wrote:
> Message: 5
> Date: Sun, 20 Jan 2013 00:00:02 -0600
> From: Daniel J Sebald <address@hidden>
> To: twinclouds <address@hidden>
> Cc: address@hidden
> Subject: Re: Switching default integrator for multiple variable to
>       quadcc
> Message-ID: <address@hidden>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> On 01/19/2013 05:14 PM, twinclouds wrote:
>> > @Rik
>> > I don't know if this is related to the problem I am having right now.  When
>> > I use Octave 3.2.4, the integration is pretty fast.  The only problem is
>> > when the area gets too large, i.e., most area are zeros, it shows "warning:
>> > quadgk: non-finite integrand encountered."  When I use Octave 3.6.3, I 
>> > don't
>> > have this problem but it gets very slow, from a few seconds to a few
>> > minutes.  I would like to change the default integrator but don't know how
>> > to.  Add the name of the integrator as the last argument of dblquad yields
>> > all sorts of errors.  Any suggestions?
> Twinclouds,
>
> The change that Rik made is in the scripts portion of the source tree. 
> The latest source code can be found here:
>
> http://hg.savannah.gnu.org/hgweb/octave/file/d56dd6794a20/scripts/general
>
> There are a small number of related changesets associated with Rik's mod 
> in this log (search CNTRL-F for "dbl"):
>
> http://hg.savannah.gnu.org/hgweb/octave/shortlog/68a59630798d
>
> The most relevant changeset is here:
>
> http://hg.savannah.gnu.org/hgweb/octave/rev/eb4afb6a1a51
>
> which should give a good idea of the changes Rik made and what you can 
> do if you want to make changes.  (Start by getting the most recent 
> version and go from there.)
For reference, the test code inside dblquad.m shows how to use the function
with different integrators besides the default and it is also mentioned in
the documentation.  The relevant test code is:

%% Nasty integrand to show quadcc off
%!assert (dblquad (@(x,y) 1 ./ (x+y), 0, 1, 0, 1), 2*log (2), 1e-6)

%!assert (dblquad (@(x,y) exp (-x.^2 - y.^2) , -1, 1, -1, 1, 1e-6,
@quadgk), pi * erf (1).^2, 1e-6)
%!assert (dblquad (@(x,y) exp (-x.^2 - y.^2) , -1, 1, -1, 1, 1e-6, @quadl),
pi * erf (1).^2, 1e-6)
%!assert (dblquad (@(x,y) exp (-x.^2 - y.^2) , -1, 1, -1, 1, 1e-6, @quadv),
pi * erf (1).^2, 1e-6)

I used function handles here but you could also just use the string name of
the integrator you want such as "quadgk".

> ...
>
> I'm curious how the error you describe is coming about.  "Non-finite 
> integrand" suggests your function might have a singularity somewhere. 
> If you are integrating through a singularity, that might be where the 
> slowness is coming from.  But you described an integrand with a lot of 
> negligible area--that sort of thing is usually not a problem for 
> adaptive integrators and runs pretty fast.  Discontinuities and 
> singularities are typically where numerical integration slows down and 
> often breaks with poor accuracy.  Some tricks you can try is that if you 
> know where there is a singularity or discontinuity, the problem can be 
> broken up into multiple sections.
>
> Dan
What Dan wrote sounds correct to me.  Usually the integrator is slowed by a
singularity, not by large stretches of zero values.  You might try a quick
plot of the function to see if there is anything odd about it.

--Rik