Re: Problème avec intnum

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

On Mon, Jun 12, 2023 at 06:39:37PM +0100, Damian Rössler wrote:
> Bonjour,
> 
> ce message pour vous signaler une problème que j’ai découvert dans Pari/GP lorsque j’ai essayé de faire un calcul avec intnum.
> Voici : lorsque j’évalue

[[Hello Damian, Sorry to reply in English for the benefit of the list]]
 
> intnum(x=-1000,1000,log(x+I)/(x^2+1)) 
> 
> sur mon ordinateur, j’obtiens
> 
> 0.22992031686600422539427385649876229873 + 99.246582527150702120676228264996637697*I

This is a known issue, there is a asymptotic singularity at x=I, so it is better to 
cut the integral at 0

intnum(x=-1000,0,log(x+I)/(x^2+1)) + intnum(x=0,1000,log(x+I)/(x^2+1))
%1 = 2.1617705842398946206845289156864102853+4.9316606089380497070728417671647276968*I

Also, to get the result at the full precision you need to do
? intnum(x=-1000,0,log(x+I)/(x^2+1),2) + intnum(x=0,1000,log(x+I)/(x^2+1),2)
%3 = 2.1617705842396943881732021681813425165+4.9316606089382864390572986833518100961*I

(basically, try ,0) ,1) ,2) until it stabilizes)

> En revanche, la commande 
> 
> NIntegrate[Log[x+I]/(x^2+1),{x,-1000,1000}]
> 
> sur Mathematica donne
> 
> 2.16177 + 4.93166 I
> 
> qui est effectivement proche de 
> 
> \pi\log(2)+i\pi^2/2 =(approx.)  2.1775860903036021305006888982376139473 + 4.9348022005446793094172454999380755677*i
> 
> (que l’on obtiendrait en intégrant de -\infty à \infty).

Indeed! By splitting the integral again:

? intnum(x=-oo,0,log(x+I)/(x^2+1)) + intnum(x=0,oo,log(x+I)/(x^2+1))
%8 = 2.1775860903036021305006888982376139473+4.9348022005446793094172454999380755677*I
? Pi*log(2)+I*Pi^2/2
%9 = 2.1775860903036021305006888982376139473+4.9348022005446793094172454999380755677*I

> Une autre bizarrerie est que Pari/GP en ligne ( https://pari.math.u-bordeaux.fr/gp.html <https://pari.math.u-bordeaux.fr/gp.html> ) donne en revanche
> 
> 0.1951079189482813284692074920 + 120.3479637305210234105619723*I
> 
> lorsqu’on effectue la même commande que plus haut, ce qui ne coïncide pas avec ce que la version de Pari/GP installée sur mon ordinateur calcule.

Yes, the online version is 32bit due to limitation of javascript, which currently implies
that the default precision is 28 digits instead of 38 digits.

> J’ai aussi essayé de calculer l’intégrale avec des bornes plus élevées et on obtient des résultats encore plus éloignés de la vraie valeur. 
> Par ex., sur mon ordinateur, la commande 
> 
> intnum(x=-100000,100000,log(x+I)/(x^2+1)) 
> 
> donne
> 
> 0.0046079812092805937682847726807162324991 + 9916.7863769592496702503585536475594909*I
> 
> Je sais bien que les intégrales oscillantes sont difficiles à approcher numériquement mais l’intégrant n’est pas oscillant ici 
> (mais il converge assez lentement vers 0).

PARI uses the double exponential method that does a change of variable. Unfortunately this causes the singularity
at I to get closer and closer to the integration path when N goes to infinity.

Thanks for using PARI/GP!
Bill