Bug in quadgk on an infinite domain

[Marco Caliari]

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