[Octave-bug-tracker] [bug #37613] Octave precision/accuracy is much lower for quadgk

[anonymous]

URL:
  <http://savannah.gnu.org/bugs/?37613>

                 Summary: Octave precision/accuracy is much lower  for quadgk
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Mon 22 Oct 2012 01:59:36 PM UTC
                Category: Libraries
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Inaccurate Result
                  Status: None
             Assigned to: None
         Originator Name: Pramod Kumbhar
        Originator Email: address@hidden
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 3.4.3
        Operating System: Any

    _______________________________________________________

Details:

Hello All,

I am new to Octave and running Matlab code on Octave 3.4.3. I am getting
incorrect results for my simulations, basically because of lower
precision/accuracy of Octave quadrature routine quadgk.

For example, say we have 

f = @(x)x.^5.*exp(-x).*sin(x);

In Matlab I get: 

K>> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)

q      = -1.499999999999836e+01
errbnd =  9.438576028474561e-09

But in the Octave, I get

debug> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)
warning: quadgk: Error tolerance not met. Estimated error 0.0329476
q      = -1.49999941897813e+01
errbnd =  3.29475856720323e-02

>From above simple test, I see Octave accuracy is very low. 

To find the cause of this, I loooked into the implementation of quadgk routine
in Octave (here:http://fossies.org/dox/octave-3.6.3/quadgk_8m_source.html). 

The problem is, when one of the limit of integral is Inf, current
implementation don't divide subintervals further due to below condition (line
no 322):

if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps))
   .....
   break
endif 

When a or b is Inf, h0 is Inf and hence code returns results with very low
accuracy.

I am not a mathematician and not sure if this is expected behaviour. But
simply removing above check we get more accurate results similar to Matlab.

Can someone please help in this issue? Suggestions for getting more accurate
results?

Regards,
Pramod.




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?37613>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/