Re: Strange behaviour of intnum

[Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>]

* Bill Allombert [2010-02-26 19:37]:
> On Fri, Feb 26, 2010 at 10:30:05AM +0100, Olivier Ramare wrote:
> > Hello,
> >
> >  (1)----------------------------------------------------------------------
> >
> > Can anyone tell me what wrong with
> >
> > intnum(t = [0,0], [1], (sin(t)/t)^2)
> 
> Your function is oscillating at oo: the documentation says:
[...]
> so here you must replace sin(x)^2 by (1-cos(2*x))/2 and sum:
> You get:
> 
> ? intnum(t = 0, 1,if(t,(1-cos(2*t))/2/t^2,1))+intnum(t = 1, [[1],2*I],-cos(2*t)/2/t^2)+intnum(t=1,[1],1/2/t^2)
> %1 = 1.5707963267948966192313216916397514421

Hum, on the other hand :

(01:18) gp > \p100
  realprecision = 105 significant digits (100 digits displayed)
(01:18) gp >  intnum(t = 0, [1], (sin(t)/t)^2)
%1 = 1.571166537031496604360105950346597224302063277903999870302529278008130571434265738562554253368755171

which is (almost) entirely wrong, but in an «expected» way (same answer
at various other accuracies). The error you get simply comes from

(01:21) gp > \p28
   realprecision = 28 significant digits
(01:21) gp > sin(1e38)
  ***   at top-level: sin(1e38)
  ***                 ^---------
  *** sin: precision too low in mpsc1.

We try to reduce the argument mod 2*Pi, but since it allows an absolute
error of about 10^10, we do not go very far...

> > intnum(t = [0,0], 1000, 9*((sin(t)-t*cos(t))/t^3)^2)
> 
> ? intnum(t = 0, 1, 9*((sin(t)-t*cos(t))/t^3)^2)
> %1 = 0.E143
> 
> This one is a bug, I do not see any reason why it fails.

A simpler variant

(01:25) gp > intnum(t = 0, 1, (sin(t)-t*cos(t))/t^3)
%1 = 0.E32

This comes from the apparent singularity at 0 :

(01:27) gp > f(t) = (sin(t)-t*cos(t))/t^3;
(01:27) gp > f(1e-40)
%3 = 0.E52  \\ function evaluated with an enormous loss of accuracy

??intnum suggests replacing the integrand by truncated expansion in this case.

(01:30) gp > F = truncate(f(t) + O(t^7))
%4 = -1/45360*t^6 + 1/840*t^4 - 1/30*t^2 + 1/3

(01:30) gp > intnum(t = 0, 1, if (t < 1e-18, eval(F), f(t)))
%5 = 0.3224571957137131128479832569

Cheers,

    K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite Bordeaux 1          Fax: (+33) (0)5 40 00 69 50
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~belabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`