# The same formula performed on two different versions of Mathematica yields different results. [Life is so tough]

Let me show the code first.

 Gama[kx_, ky_, u_] := NIntegrate[96*Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(-7/2)*1/(
6!*Sqrt[Pi])*(((Sqrt[2*u]/(
Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^6 +
15/2*(Sqrt[2*u]/(
Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^4 +
45/4*(Sqrt[2*u]/(
Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^2 + 15/
8)*Sqrt[Pi]*
Exp[(Sqrt[2*u]/(
Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^2]*
Erfc[Sqrt[2*u]/(
Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2))] - (Sqrt[
2*u]/(Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^5 -
7*(Sqrt[2*u]/(
Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^3 -
33/4*(Sqrt[2*u]/(
Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))), {p, -Pi,
Pi}];


Then you can use

 Plot[Gama[0, 0, u], {u, 0, 1}, PlotRange -> All]


You may find that the values are very large.

Then I try another computer, I get

These days, I am really lucky. Life is tough!

• Those two graphs look the same to me. What are we looking for? Jul 12, 2020 at 13:39
• see Fig1 and Fig3 @ChrisK Jul 12, 2020 at 13:48
• Here I upload three figures. The last two figures are the same. I do not know how to delete one of them. Jul 12, 2020 at 13:49
• @flinty version:11.3.0.0 11.2.0.0 Jul 12, 2020 at 14:43
• @Blueka The same result I got with versions 12.0.0 (your first picture) and 12.1.1 (your second picture). It means that result is unstable. We can check it just call Table[{u, Gama[0., 0., u]}, {u, .0, 1., .1}]. There are pole at p=Pi/2 and p=-Pi/2. So we should use some strategy to compute this integral. Jul 12, 2020 at 16:56

It's a question of numerical precision.

Set ky and ky to zero and plot the integrand ii, enlarge working precision and rationalize variable in Plot.

ii[p_, u_] =
96*Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(-7/2)*1/(6!*
Sqrt[Pi])*(((Sqrt[
2*u]/(Cos[
p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^6 +
15/2*(Sqrt[
2*u]/(Cos[
p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^4 +
45/4*(Sqrt[
2*u]/(Cos[
p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^2 +
15/8)*Sqrt[Pi]*
Exp[(Sqrt[
2*u]/(Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^2]*
Erfc[Sqrt[
2*u]/(Cos[
p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2))] - (Sqrt[
2*u]/(Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^5 -
7*(Sqrt[2*
u]/(Cos[p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))^3 -
33/4*(Sqrt[
2*u]/(Cos[
p/2]*(1 + I*kx*Cos[p] + I*ky*Sin[p])^(1/2)))) /. {kx -> 0,
ky -> 0} // Simplify[#, -Pi < p < Pi && 0 < u < 1] &

Plot3D[ii[p, u], {p, -Pi, Pi}, {u, 0, 1}, PlotRange -> {0, .1},
PlotPoints -> 50, WorkingPrecision -> 30]

Ga00[u_] := NIntegrate[ii[p, u], {p, -Pi, Pi}, WorkingPrecision -> 30]

Plot[Ga00[Rationalize[u, 0]], {u, 0, 1}, PlotRange -> All]


• Cool! Thanks! I get the point. Jul 13, 2020 at 8:47