Bug in quadgk on an infinite domain

[David Bateman]

Marco Caliari wrote:
> On Tue, 9 Jun 2009, David Bateman wrote:
>
>> Marco Caliari wrote:
>>> Dear maintainers,
>>>
>>> there is a bug in the 'Standard Infinite to finite integral 
>>> transformation' in quadgk.m (3.2.0). The correct transformation is
>>>
>>> ##   g(t)  = t / (1 - t^2)
>>> ##   g'(t) =  (1 + t^2) / (1 - t^2) ^ 2
>>>
>>> With the enclosed patch, I get
>>>
>>> octave:1> quadgk(@(x)exp(-x.^2),-10,10)
>>> ans =  1.77245385090546e+00
>>> octave:2> quadgk(@(x)exp(-x.^2),-Inf,Inf)
>>> ans =  1.77245385078031e+00
>>>
>>> that is sqrt(pi)
>>>
>>> octave:3> sqrt(pi)
>>> ans =  1.77245385090552e+00
>>>
>>> whereas, with the original quadgk, I get
>>>
>>> octave:4> quadgk(@(x)exp(-x.^2),-Inf,Inf)
>>> ans =  1.15710256465273e+00
>>>
>>> Best regards,
>>>
>>> Marco
>> Oppps missed that when I wrote this.. There is also a bug in the 
>> waypoint transform for the doubly infinite integrands in that the 
>> waypoint at x=0 should be zero but the transform return NaN.
>
> A much more elegant and stable waypoint transform is
>
> trans = @(x) (2 * x) ./ (1 + sqrt(1 + 4 * x .^ 2));
>
> Best regards,
>
> Marco
>
Ok I committed this waypoint transform credit to you to savannah

D.


-- 
David Bateman                                dbateman at dbateman.org
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)