Re: Problème avec intnum

[Damian Rossler <damian.rossler@maths.ox.ac.uk>]

Dear Bill, thank you for this. What do you mean by « an asymptotic singularity at I »? 
The integral is along the real line, so it never meets any singularity of the function (in particular, log(x+I) is always well-defined). 
Best wishes, Damian

> On 12 Jun 2023, at 7:43 pm, Bill Allombert <Bill.Allombert@math.u-bordeaux.fr> wrote:
> 
> 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