Re: References for Numerical Integration

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Sun, Oct 21, 2018 at 08:35:04PM +0300, kevin lucas wrote:
> PARI obviously has some very powerful tools for numerical integration. I
> have, however, keenly felt the absence of a text treating this area in more
> detail than the manual. 

> I was recently trying to compute
> intnum(x=0,oo, (sin(x)^4)/(x)^2)
> I only get about 3 digits (the integral evaluates to log(2)). 

Hello Kevin,

Do you mean Pi/4 ?

> There's probably a neat way to get more digits in this particular instance.

The documentation ??intnum explain:
<http://pari.math.u-bordeaux.fr/dochtml/html-stable/Sums__products__integrals_and_similar_functions.html#intnum>

    The  last  two  codes  are  reserved for oscillating functions.
 Let k > 0 real,  and g(x) a non-oscillating function tending slowly to 0 
 (e.g. like a negative power of x), then
 * alpha = k * I assumes that the function behaves like cos(kx)g(x).
 * alpha = -k* I assumes that the function behaves like sin(kx)g(x).
 
 Here it is critical to give the exact value of k.  If the oscillating
 part is not a pure sine or cosine,  one must expand it into a Fourier
 series,  use the above codings, and sum the resulting contributions.
 Otherwise you will get nonsense. Note that cos(kx), and similarly
 sin(kx), means that very function, and not a translated version such as
 cos(kx+a).
and later see 'Oscillating functions.'

Here we have g(x)=1/x^2 and 
sin(x)^4 = (3-4*cos(2*x)+cos(4*x))/8

a formula you can find by doing
? lindep([sin(x)^4,1,cos(2*x),cos(4*x)])
%9 = [-8,3,-4,1]~

To avoid the division by 0, we split the integral at 1 and do:

A=intnum(x=0,1,sin(x)^4/x^2);
B=intnum(x=1,oo,1/x^2);
C=intnum(x=1,[oo,2*I],cos(2*x)/x^2);
D=intnum(x=1,[oo,4*I],cos(4*x)/x^2);
A+(3*B-4*C+D)/8
%5 = 0.78539816339744830961566084581987572104
? Pi/4
%6 = 0.78539816339744830961566084581987572105

Cheers,
Bill