Bug in quadgk on an infinite domain

[David Bateman]

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. I took your 
patch, fixed the waypoint transform, added a test case and changelog and 
committed this to savannah, but it should probably go into 3.2.1 as well...

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)