bugGNU Octave - Bugs: bug #37896, quadcc(f,-Inf,Inf) incorrectly...

 
 

bug #37896: quadcc(f,-Inf,Inf) incorrectly returns 0 (quadgk too)

Submitter:  Drew Parsons <dparsons>
Submitted:  Tue 11 Dec 2012 01:04:34 AM UTC
   
 
Category:  None Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Invalid / Not an Octave Bug Assigned to:  None
Originator Name:  Drew Parsons Open/Closed:  * Closed
Release:  * 3.6.2 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Thu 14 Mar 2013 08:40:45 PM UTC, comment #3: 

Computers are great aids, but they can't replace humans.  Try plotting the function to get some intuition about it.  Because of the Heaviside function, the value is 0 for all times less than 0.  And while it has a peak at early times, it quickly decays off to zero.  When you ask quadcc or quadgk to integrate a function over an infinite range it will check a couple of points.  In this case the value at -Inf,0,+Inf which are all zero.  The integration routine simply can't find anything to integrate and returns 0.

If I use an integration range that includes the non-zero portions, say t = [0, 5000], I get


octave:15> quadcc (f, 0, 5000)
ans =  1.2703e+09
octave:16> ans - quadcc (f, 0, Inf)
ans = -1.0156


And notice that +Inf is acceptable as an endpoint as long as it can find something to integrate.

Rik <rik5>
Group administrator
Wed 12 Dec 2012 02:19:05 AM UTC, comment #2: 

Thanks Jordi.  But after applying the #37414 patch, quadcc stills behaves the same for this function, returning 0 when the bounds are bigger than 1e6.

Drew Parsons <dparsons>
Tue 11 Dec 2012 02:16:13 PM UTC, comment #1: 

This looks like a duplicate of bug #37414. Can you apply that fix and confirm?

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Tue 11 Dec 2012 01:04:34 AM UTC, original submission:  

I'm using quadcc (and comparing with quadgk) for an integral which in general has infinite bounds.

The function is
f = @(x) exp (-x .^ 2 / (2 500 ^ 2)) / (sqrt (2 pi) 500) . (150e5 / 2) .* (22 - x) .^ 2 ./ (50 - (22 - x)) .* heaviside (-(22 - x), 1)
(the real coefficients and the Heaviside function relate to the physics that I'm modelling)

It has a large, finite value, which I can estimate with

> lim=1e5; quadcc ( f, -lim, lim, 1e-5 )

ans =  1.2703e+09

But with any larger limit (including Inf), quadcc returns 0:

> lim=1e6; quadcc ( f, -lim, lim, 1e-5 )

ans = 0

> lim=Inf; quadcc ( f, -lim, lim, 1e-5 )

ans = 0



quadgk is also misbehaving at large limits. I get the same solution (with some small error) at a limiting finite bound,

> lim=1e5; quadgk ( f, -lim, lim, 1e-5 )

warning: quadgk: maximum interval count (650) met
warning: quadgk: Error tolerance not met. Estimated error 0.000776351
ans =  1.2703e+09

> lim=1e6; quadgk ( f, -lim, lim, 1e-5 )

warning: quadgk: maximum interval count (650) met
warning: quadgk: Error tolerance not met. Estimated error 0.171285
ans =  1.2703e+09

but "mostly zero" at larger finite limits,

> lim=1e7; quadgk ( f, -lim, lim, 1e-5 )

ans =  5.4000e-131

But then at Inf, quadgk gives a finite value, meeting the correct value within gross error bars,

> lim=Inf; quadgk ( f, -lim, lim, 1e-5 )

warning: quadgk: Error tolerance not met. Estimated error 2.41746e+09
ans =  2.4847e+09


This quadgk behaviour is curious, but I understand quadgk is expected to come with greater error.  I'm mainly filing this bug to address the misbehaviour in quadcc.

Drew Parsons <dparsons>

 

(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 rik5 (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by dparsons (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 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2013-03-14 rik5 StatusNeed Info Invalid / Not an Octave Bug
        Open/ClosedOpen Closed
    2012-12-11 jordigh StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code