bugGNU Octave - Bugs: bug #57603, quadv gives wrong results...

 
 

bug #57603: quadv gives wrong results integrating sin^2 or cos^2 on multiple of 8 periods

Submitter:  None
Submitted:  Wed 15 Jan 2020 10:59:14 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  rik5
Originator Name:  Olivier Amoignon Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 21 Jan 2020 04:56:27 PM UTC, comment #8: 

Thanks Rik for taking care about this!  All examples of this item are now working.

As Software "primarily intended for numerical computations", I feel a little unhappy if Octave is struggling with such.  Even tough there exist functions to handle this case explicitly.

Kai Torben Ohlhus <siko1056>
Group Member
Fri 17 Jan 2020 05:01:54 AM UTC, comment #7: 

Turns out that there was a lot of room for improvement.  I re-wrote a lot of quadv.m in this changeset https://hg.savannah.gnu.org/hgweb/octave/rev/6d35ab1cd544.

The function now correctly handles periodic integrands.  I also implemented the Richardson extrapolation algorithm which reduces the number of required function evaluations by ~3.5X.

This will all be in the new 6.1 release of Octave which is probably in February.

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Thu 16 Jan 2020 03:57:13 PM UTC, comment #6: 

I took a look at this yesterday and discovered quite a few things that were undesirable about the implementation.

1) The algorithm always begins by bisecting the interval, which can cause problems if the integrand is periodic.  I also came up with the same idea as Olivier in comment #5, that is, Octave should divide the initial integral in to two unequal regions.  Ideally, the division should be an irrational number so that no set of integer factors N/D can ever be periodic.  At this point, I ran some experiments and it appears that Matlab divides the initial integral in to three separate regions.  The left and right regions are approximately 27% of the absicssa distance (b - a), and the middle region is 46% of (b - a).  I'm not sure why trisection would be particularly better than bisection, but it is also easy to implement.

2) The Octave algorithm is not using Richardson extrapolation which makes each integral interval just a bit more accurate.  For this integral


quadv(@(t)sin(t).^2,0,2.5*pi, 1e-6, true)


The current method requires 217 function evaluations versus just 61 if extrapolation is used.  That is a tremendous difference and needs to be coded in to Octave.

3) The algorithm is not correctly working to implement the user's requested tolerance.

We want to bound Q to +/- tol, and we calculate Q by dividing the region in two and then calculating q1 and q2 where Q = q1 + q2.  In order to preserve the user's request, the tolerance needs to halve at each interval.  Then


Q = (q1 +/- tol/2) + (q2 +/- tol/2) = q1 + q2 +/- tol.


4) The tracing output display could be better aligned.

I've made some progress on this already so I've assigned the bug to myself.

Rik <rik5>
Group administrator
Thu 16 Jan 2020 01:19:28 PM UTC, comment #5: 

After having a look at the code in scripts/general/quadv.m implementing an adaptive Simpson algorithm the cause of the wrong result for certains intervals of periodic functions being periodically equal to 0 is straight-forward: at nfun=5 the value of q is still q=q0 and thus the calculation stops and returns the last value q which is 0 because of the choice of the new quadrature points (d and e) as middle points. With b-a=8pi and a=0 when integrating sin(x), or a=pi/2 when integrating cos(x), the value of sin or cos at the new points (d and e) is still 0. To avoid the problem, which could occur with the current implementation for any harmonic function on specific integration intervals, a simple cure would be to initialize quadv starting with a division of the interval [a,b] in intervals of different length. It rises the initial cost of quadv from 3 to 5 function "f" calls, if initializing with a two intervals of different length, but it works fine.

Olivier Amoignon <olamo>
Thu 16 Jan 2020 09:59:32 AM UTC, comment #4: 

The singularity reported (the result is wrong on 4 periods=8Pi and not 3 periods=6Pi or 5 periods=10Pi etc...) is a bug because independently of the accuracy of the method the integration of sin^2 or cos^2 cannot be 0 on an entire period, so even less so on 4 periods. Also I should mention that the error appeared when I conducted tests comparing matlab and octave. In matlab quadv behaves as all other integration methods (quad, quadgk etc ..) and gives the righ result.

Olivier Amoignon <olamo>
Wed 15 Jan 2020 06:06:54 PM UTC, comment #3: 

Check the Octave manual for guidance about which integrands are appropriate for which numerical integration function.  Quoting


The best quadrature algorithm to use depends on the integrand. If you have empirical data, rather than a function, the choice is trapz or cumtrapz. If you are uncertain about the characteristics of the integrand, quadcc will be the most robust as it can handle discontinuities, singularities, oscillatory functions, and infinite intervals. When the integrand is smooth quadgk may be the fastest of the algorithms.

Function
 Characteristics
quad
 Low accuracy with nonsmooth integrands
quadv
 Medium accuracy with smooth integrands
quadl
 Medium accuracy with smooth integrands. Slower than quadgk.
quadgk
 Medium accuracy (1e-6 – 1e-9) with smooth integrands.
Handles oscillatory functions and infinite bounds
quadcc
 Low to High accuracy with nonsmooth/smooth integrands
Handles oscillatory functions, singularities, and infinite bounds


For an oscillatory integrand the programmer should be using quadcc, or maybe quadgk.

However, I do agree that quadv can probably be improved.  All of the other integration routines seem to be able to handle the integrand and bounds mentioned.  There must be something quite specific about quadv that is preventing it from reaching the correct result.


Rik <rik5>
Group administrator
Wed 15 Jan 2020 11:43:34 AM UTC, comment #2: 

Thanks for spotting this one.  I think Octave should be able to handle periodic functions with these typical border values.

In the current stable release 5.1.90 and the default development branch this error is present as well.

Kai Torben Ohlhus <siko1056>
Group Member
Wed 15 Jan 2020 11:07:59 AM UTC, comment #1: 

Repost comment #0, as Savannah messed up the format.

original submission:

> Tested with version 4.2.2 under Ubuntu 18.04LTS and version 3.8.2 under Centos 7:
> examples are


> quadv(@(t)sin(t).^2,0,8*pi)


> which gives 8.0412e-30 whereas


> quadv(@(t)sin(t).^2,0,10*pi)


> gives 15.708
> Another example is


> quadv(@(t)cos(t).^2,0.5*pi,8.5*pi)


> which gives 2.2364e-29 instead of 12.566, the result obtained by quad or quadgk.

Kai Torben Ohlhus <siko1056>
Group Member
Wed 15 Jan 2020 10:59:14 AM UTC, original submission:  

Tested with version 4.2.2 under Ubuntu 18.04LTS and version 3.8.2 under Centos 7:
examples are

quadv(@(t)sin(t).^2,0,8*pi)
+verbatim+
which gives 8.0412e-30 whereas
+verbatim+
quadv(@(t)sin(t).^2,0,10*pi)
+verbatim+
gives 15.708
Another example is
+verbatim+
quadv(@(t)cos(t).^2,0.5*pi,8.5*pi)
+verbatim+
which gives 2.2364e-29 instead of 12.566, the result obtained by quad or quadgk.

Anonymous

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by olamo (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by siko1056 (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 group members can vote.

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-01-17 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2020-01-16 rik5 Severity4 - Important 3 - Normal
        StatusConfirmed In Progress
        Assigned toNone rik5
    2020-01-15 siko1056 Severity3 - Normal 4 - Important
        StatusNone Confirmed
        Release4.2.2 dev

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code