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

You are not logged in

Discussion

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

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

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>
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
Low accuracy with nonsmooth integrands
Medium accuracy with smooth integrands
Medium accuracy with smooth integrands. Slower than quadgk.
Medium accuracy (1e-6 – 1e-9) with smooth integrands.
Handles oscillatory functions and infinite bounds
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>
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>
Wed 15 Jan 2020 11:07:59 AM UTC, comment #1:

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

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

> which gives 8.0412e-30 whereas

> gives 15.708
> Another example is

Kai Torben Ohlhus <siko1056>
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

Attached Files

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

Attach Files:

Comment:

No files currently attached

Dependencies

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?
This task has 0 encouragements so far.

Only project members can vote.