Karim Belabas on Tue, 11 Dec 2018 12:46:45 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: Integration of periodic functions |
* 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] `