Re: Integration of periodic functions

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

* Jacques Gélinas [2018-12-11 01:03]:
> Responding to http://pari.math.u-bordeaux.fr/archives/pari-dev-1812/msg00003.html
> 
> For periodic functions integrated over a period, the trapezoidal or
> midpoint methods can give much more accurate results than other
> numerical integration methods. Compare, for example, the number of
> bits LOST in
> 
> exr(x,y) = [exponent( if(!x,y,1-y/x) ), 0] - exponent(0.)*[1,1];
> 
> { localbitprec(64*4);
> 
> rm = intnumromb(X=0,2*Pi,sin(X)/(5+3*sin(X)));
> 
> de = intnum(X=0,2*Pi,sin(X)/(5+3*sin(X)));
> 
> exr(-Pi/6,rm) == [9,256] && exr(-Pi/6,de) == [158,256]
> 
> }

You can improve the result accuracy by choosing a smaller integration
step [ see ??intnum ] :

  de1 = intnum(X=0,2*Pi,sin(X)/(5+3*sin(X)), 1 /* h -> h/2 */);
  de2 = intnum(X=0,2*Pi,sin(X)/(5+3*sin(X)), 2 /* h -> h/4 */);

  (12:32) gp > exr(-Pi/6,de1)
  %6 = [60, 256] \\ better but not perfect
  (12:32) gp > exr(-Pi/6,de2)
  %7 = [1, 256]  \\ perfect

This is due to the pole at asin(-5/3) ~ -Pi/2  + 1.098...*I
which is (relatively) close to the integration interval and
causes larger error terms.

PARI not being a CAS doesn't know anything about your function [ it only
gets a black box able to evaluate to given accuracy at various points ]. 
So we can't compensate for nearby poles on our own, the user must do it...

Cheers,

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