Re: Slow integration with dblquad

[twinclouds]

Rik, Dan and Alec:
Thank you all for your reply.  I think the situation is clear to me now, even though I still have some questions.
The difference between 3.2.4 and 3.6.2/3 is that the default integrator of dbldquad in 3.2.4 is quadgk and it is not (probably it is quadcc) in 3.6.x.  What my observation is that quadgk is much faster in my case even for the relatively well behaved cases.  However, when the function is not well behaved, quadgk will simply give up and it tied to get a solution in the dblquad in 3.6.x.  Sometime, it would take forever to make not very useful.  BTW, if I set the default integrator to be quadgk, dblquad in 3.6.x will have the same behavior as in 3.2.4.
As for error of integrand over range, I think the reason is the function that I used in integration is f(x,y)*log2(f(x,y)), where f(x,y) takes the form of exp(-((x-x0)^2+(y-y0)^2)/2s).  When f(x,y) is close or equal to 0, f(x,y)*log2(f(x,y)) should be approach zero but log2(f(x,y)) goes to infinity.  This situation could happen in a number of cases, i.e., x ~=x0, y ~= y0 or x and y become very large at the same time.  After put a statement to detect f(x,y) is approximately equal to zero, everything seems works.
The question I still have is that why dblquad in 3.6.x does not work as well as the counterpart in 3.2.4, even for relatively well behaved case.   This seems contradict to what Rik's observation?  I can provide my latest scripts if anybody is interested.
Thanks again for everybody's help.  I am new to Octave but this forum seems to be a very responsive one.  I am certainly would like to use it more.

--- On Sun, 1/20/13, Rik <address@hidden> wrote:

From: Rik <address@hidden>
Subject: Re: Slow integration with dblquad
To: "twinclouds" <address@hidden>
Cc: address@hidden
Date: Sunday, January 20, 2013, 11:31 AM

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