bugGNU Octave - Bugs: bug #46349, dblquad/triplequad - periodic...


bug #46349: dblquad/triplequad - periodic functions over the whole period

Submitted by:  None
Submitted on:  Sat 31 Oct 2015 06:19:16 PM UTC
Category:  Octave Function Severity:  4 - Important
Priority:  5 - Normal Item Group:  Performance
Status:  Fixed Assigned to:  None
Originator Name:  Robert Grigo Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * GNU/Linux
* Mandatory Fields

Add a New Comment (Rich Markup)

Comment Type & Canned Response:


( Jump to the original submission )

Mon 25 Sep 2017 07:11:55 PM UTC, comment #11: 

I changed quadcc to use both a relative tolerance and an absolute tolerance for the stopping criteria.  The cset is here (http://hg.savannah.gnu.org/hgweb/octave/rev/71dad5be765a).  This fixes this bug and now quadcc is about 12X faster than quadgk.  Closing report.

Rik <rik5>
Project Administrator
Tue 21 Mar 2017 04:10:09 PM UTC, comment #10: 

looking at the integrators and came across this. in 4.2.1 dblquad works quickly for all but quadcc on the test in Comment #1.  It actually does finish with quadcc, but takes quite some time. since quadcc is the default, dblquad without specifying the integrator takes a long time.

>> tic;dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1, [], 'quadgk');toc
Elapsed time is 0.188 seconds.

>> tic;dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1, [], 'quadv');toc
Elapsed time is 0.072 seconds.

>> tic;dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1, [], 'quadl');toc
Elapsed time is 0.0909998 seconds.

>> tic;dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1, [], 'quadcc');toc
Elapsed time is 49.752 seconds.

>> tic;dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1);toc
Elapsed time is 48.5679 seconds.

so it's incorrect to say that it doesn't terminate, but a factor of 200-300x in computation time is surely not ideal.

Nicholas Jankowski <nrjank>
Project Member
Mon 11 Jul 2016 01:25:13 PM UTC, comment #9: 

Marco, the reason that

quadcc(@(x) sin(2*pi*x),0,1)

converges is that  nivals  (the number of intervals being further refined) goes to 0.

The current exit criterion is

  while (nivals > 0 && err > 0.0 && err > fabs (igral) * tol
         && !(err_final > fabs (igral) * tol
              && err - err_final < fabs (igral) * tol))

and on exit the values are:
nivals 0
err 3.06866e-16
fabs(igral) 5.55112e-17
tol 1e-06
err_final 3.06866e-16
err - err_final 0

The criterion for increasing the number of intervals seems to be independent of the tolerance.

Have there been any examples of quadcc giving trouble, or is it just that it uses a relative error termination criterion?  If Rik's patches fix dblquad/triplequad, and no actual failures have been seen in quadcc, I'd be inclined to close this report.


Lachlan Andrew <lachlan>
Tue 10 Nov 2015 07:43:30 AM UTC, comment #8: 

I find quite strange that quadcc succeeds in the computation of

quadcc(@(x) sin(2*pi*x),0,1)

I have to say that I do not understand very much the conditions at lines 1787-1789 (what is err_final?). Anyway, since in this case there is no compatibility problem (quadcc in not in Matlab), I would explicitely introduce a relative and an absolute tolerance, like in quadgk and change the condition to

&& err > fabs(igral) * reltol && err > abstol

Sounds reasonable?

Marco Caliari <caliari>
Project Member
Mon 09 Nov 2015 04:39:16 PM UTC, comment #7: 

There were a number of issues with the quadl function.  The biggest was that it was using a relative tolerance, when it is clear from the Matlab documentation that an absolute tolerance should be used.  I overhauled the function here (http://hg.savannah.gnu.org/hgweb/octave/rev/655816377845) and now it completes for the single and double integration examples in this bug report.

The only remaining integrator which fails is quadcc, which, unfortunately, is the default and will be encountered frequently.  Like quadl, it uses a relative tolerance.  It may be that this needs to be augmented like quadgk which uses a combination of relative and absolute tolerances to determine convergence.  Or maybe the inner integrator can be quadcc, but the outer needs to to be something else like quadgk.

Rik <rik5>
Project Administrator
Sun 08 Nov 2015 06:22:00 PM UTC, comment #6: 

The stop condition of qudl is:

is+(i1-i2) == is

What about adding more condition? For example:

is+(i1-i2) == is || abs(i1-i2) < tol

see the attached patch.

(file #35416)

Avinoam Kalma <avinoam>
Project Member
Sat 07 Nov 2015 11:16:02 AM UTC, comment #5: 

About quadl:

In the documentation: "quadl will be removed in a future release.".
It also has limits on step size and on number of times that the
integrand is evaluated, to avoid endless loop.

I did not find quadcc in matlab documentation.

Avinoam Kalma <avinoam>
Project Member
Sat 07 Nov 2015 12:44:50 AM UTC, comment #4: 

@Marco: Confirmed for quadl.  Matlab uses an absolute error tolerance and we should too for compatibility, and to fix this bug.

For quadcc, I will quote from the documentation.

The optional argument TOL defines the relative tolerance used to
stop the integration procedure.  The default value is 1e^-6.

So quadcc is going to have problems and require more work to fix.

Rik <rik5>
Project Administrator
Thu 05 Nov 2015 09:31:16 AM UTC, comment #3: 

I think it is related to the exit criteria used by quadl and quadcc, which, from the documentation, are based on the relative error. In fact,

quadl(@(x) sin(2*pi*x),0,1)

does not terminate, as well (but quadcc does terminate).

Long time ago I submitted a bug report (or I wrote to the mainling list, don't remember) about the different exit criteria used by the different quadrature functions. Moreover, quadl in Matlab uses an absolute tolerance. Probably for quadl the best solution is to move to an absolute tolerance. I did not analyze quadcc, yet.


Marco Caliari <caliari>
Project Member
Sun 01 Nov 2015 07:19:02 PM UTC, comment #2: 

This also applies to triplequad.  Even though just one of the integrals is over the whole period, the entire integration fails.

z = triplequad(@(x,y,z) sin(2 * pi * x ) * sin(2 * pi * y) * sin (2*pi*z),0,1,0.2,0.4,0.2,0.4, 1e-1, 'quadcc');

Rik <rik5>
Project Administrator
Sun 01 Nov 2015 02:29:49 PM UTC, comment #1: 

Confirmed.  There appears to be a problem only with certain quad solvers.  The extended syntax

dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1, [], SOLVER);

can be used to check which ones work.

Working: quadv, quadgk
Failing: quadl, quadcc

I checked back to version 3.2.4 and this has always been a problem, apparently, but not discoverd until now.

Rik <rik5>
Project Administrator
Sat 31 Oct 2015 06:19:16 PM UTC, original submission:  

maybe there is a problem on integrating
periodic functions over the whole period with

dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1);

dblquad don't terminate.
Whoever it works e.g. with the 'quadv' option.

dblquad(@(x,y) sin(2 * pi * x * y),0,1,0,1);

will work.




(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #35416:  quadl.diff added by avinoam (2KiB - application/octet-stream)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by avinoam (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only project members can vote.





    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-09-25 rik5 Item GroupSegfault, Bus Error, etc. => Performance
        StatusConfirmed => Fixed
        Open/ClosedOpen => Closed
    2016-07-01 rik5 Severity3 - Normal => 4 - Important
    2015-11-08 avinoam Attached File- => Added quadl.diff, #35416
    2015-11-01 rik5 Summarydblquad - periodic functions over the whole period => dblquad/triplequad - periodic functions over the whole period
    2015-11-01 rik5 StatusNone => Confirmed

    Back to the top

    Powered by Savane 3.9