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

Submitted by:  None
Submitted on:  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:  Rik <rik5>
Originator Name:  Olivier Amoignon Originator Email:  -email is unavailable-
Open/Closed:  Closed Release:  dev
Operating System:  GNU/Linux

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( 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>
Project 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>
Project AdministratorIn charge of this item.
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>
Project AdministratorIn charge of this item.
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>
Project AdministratorIn charge of this item.
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>
Project 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>
Project 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

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)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project 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.5