quad fails on low amplitude oscillatory functions

[David Bateman]

Blair Sutton wrote:
> Bug report for Octave 3.1.50 configured for i686-pc-msdosmsvc
> 
> Description:
> -----------
> 
> When numerically integrating a low amplitude oscillatory function quad 
> outputs "ABNORMAL RETURN FROM DQAGP" to stderr and gives an incorrect 
> answer.
> 
> Repeat-By:
> ---------
> 
> Run the following code in the interpreter. The first quad fails whereas 
> the second is fine.
> 
> function y = quadtest(x, k)
> y = cos(2*x) * exp(-k^2);
> endfunction
> quad(inline(sprintf('quadtest(x,%f)', 0.001)),-300,300)
> quad(inline(sprintf('quadtest(x,%f)', 10)),-300,300)
> 
> Fix:
> ---
> 
> Since an integral is linear, one can scale the integrand to reasonable 
> range of the underlying cpu's floating point capabilities before using 
> the quadrature method. Then scale back the final result. Perhaps even 
> have an adaptive scaling method.
>

The answer is in fact not all that bad. Consider

octave:22> k= 0.001; quad(@(x) cos(2*x) * exp(-k^2), -300, 300)
  ABNORMAL RETURN FROM DQAGP
ans =  0.044182
octave:23> k= 0.001; quadgk(@(x) cos(2*x) * exp(-k^2), -300, 300)
warning: quadgk: Error tolerance not met. Estimated error 8.73317e-06
ans =  0.044182

Yes quadrature of oscillatory integrals is always difficult. But quad 
and quadgk both attempt to be intelligent about them and use adaptive 
techniques. The issue is that both quad and quadk use an absolute 
tolerance for their convergence test. Try using a purely relative error 
instead like

octave:13> k= 0.001; quadgk(@(x) cos(2*x) * exp(-k^2), -300, 300, 
'AbsTol', 0, 'RelTol', 1e-6)
ans =  0.044182
octave:14> k= 0.001; quad(@(x) cos(2*x) * exp(-k^2), -300, 300,  [0, 1e-2])
ans =  0.044182

I presonally prefer quadgk, as its a bit faster in real cases, seems to 
have better convergence properties, and uses much the same algorithm as 
in matlab for the same function.

Cheers
David

-- 
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)