# why this analytical formula can't yields the answer of numerical integration?

I have an initial function as follows

intrp[r_, re_] := (
E^(-2 r (\[Alpha] +
r \[Beta])) (-E^(-((
2 (me re^2 + mp (r + re) (r + re + 2 me re)) \[Gamma])/M^2)) +
E^(-((2 \[Gamma] (me (-mp r^2 + re^2 + mp re^2) + (1 +
me) mp RealAbs[r - re]^2))/M^2))) M^2)/(
4 (1 + me) mp \[Gamma]);


where the constants are

me = 1;
mp = 1;
M = mp + me;
\[Omega] = 0.0001;
{\[Alpha], \[Beta]}={0.50000000206214434417262282295268960297107696533203,
1.1754907046612932983601436560949782261786822346039*10^-8};
\[Gamma] = SetPrecision[0.5*M*\[Omega], 50];
norm=0.00014216419373200363465950006063;


I go to the next step using both Integrate and NIntegrate and I expect to get the same result, but I don't! For the next step I have

int[re_] :=(2*\[Pi])/re norm^2*Integrate[r intrp[r, re], {r, 0, \[Infinity]}];
nint[re_] := (2*\[Pi])/re norm^2*
NIntegrate[r intrp[r, re], {r, 0, \[Infinity]}];


Now if you evaluate the Integrate command in int[re] and then call its result as int2[re] which has the following form

int2[re_] := (E^(-((2 (me + mp + 2 me mp) re^2 \[Gamma])/M^2))
M^2 norm^2 \[Pi]^(
3/2) (-E^(((2 \[Alpha] - (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2)))) M^2 \[Alpha] +
E^((2 \[Alpha] + (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2))) M^2 \[Alpha] +
2 E^((2 \[Alpha] - (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2))) mp re \[Gamma] +
2 E^((2 \[Alpha] + (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2))) mp re \[Gamma] +
2 E^((2 \[Alpha] - (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2))) me mp re \[Gamma] +
2 E^((2 \[Alpha] + (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2))) me mp re \[Gamma] +
E^((2 \[Alpha] - (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2)))
M^2 \[Alpha] Erf[(M^2 \[Alpha] - 2 (1 + me) mp re \[Gamma])/(
Sqrt[2] M^2 Sqrt[\[Beta] + (mp \[Gamma])/M^2])] -
2 E^((2 \[Alpha] - (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2)))
mp re \[Gamma] Erf[(
M^2 \[Alpha] - 2 (1 + me) mp re \[Gamma])/(
Sqrt[2] M^2 Sqrt[\[Beta] + (mp \[Gamma])/M^2])] -
2 E^((2 \[Alpha] - (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2)))
me mp re \[Gamma] Erf[(
M^2 \[Alpha] - 2 (1 + me) mp re \[Gamma])/(
Sqrt[2] M^2 Sqrt[\[Beta] + (mp \[Gamma])/M^2])] -
E^((2 \[Alpha] + (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2)))
M^2 \[Alpha] Erf[(M^2 \[Alpha] + 2 (1 + me) mp re \[Gamma])/(
Sqrt[2] M^2 Sqrt[\[Beta] + (mp \[Gamma])/M^2])] -
2 E^((2 \[Alpha] + (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2)))
mp re \[Gamma] Erf[(
M^2 \[Alpha] + 2 (1 + me) mp re \[Gamma])/(
Sqrt[2] M^2 Sqrt[\[Beta] + (mp \[Gamma])/M^2])] -
2 E^((2 \[Alpha] + (4 (1 + me) mp re \[Gamma])/M^2)^2/(
8 (\[Beta] + (mp \[Gamma])/M^2)))
me mp re \[Gamma] Erf[(
M^2 \[Alpha] + 2 (1 + me) mp re \[Gamma])/(
Sqrt[2] M^2 Sqrt[\[Beta] + (mp \[Gamma])/M^2])]))/(8 Sqrt[
2] (1 + me) mp re \[Gamma] (M^2 \[Beta] +
mp \[Gamma]) Sqrt[\[Beta] + (mp \[Gamma])/M^2]);


when I tried nint[10] ,int[10]and int2[10] for example I get two different results.nemely:

int[10]
int2[10]
nint[10]
(* 4.97596488191011311257756217339320398823879335502331227.\
85176023066818*^-7
0-2135.7355398006894
4.975964881910046*^-7 *)


I know the answer of NIntegrate is the real answer. It seems so odd why the analytic form can't yield the right answer?

I would be very grateful if somebody can help me, because I need to use the analytic form.

• "you get different answers!": for what further calculations? Can you show the exact code in which you used the results of Integrate? Commented Oct 9, 2021 at 14:05
• I updated the question, please read it again Commented Oct 9, 2021 at 14:10

Clear["Global*"]

$Version (* "12.3.1 for Mac OS X x86 (64-bit) (June 19, 2021)" *)  Use exact numbers whenever possible and the precision of norm should be consistent with the other values. me = 1; mp = 1; M = mp + me; ω = 10^-4; {α, β} = {0.50000000206214434417262282295268960297107696533203, 1.1754907046612932983601436560949782261786822346039*10^-8}; γ = M*ω/2; norm = SetPrecision[0.00014216419373200363465950006063, 50]; intrp[r_, re_] := (E^(-2 r (α + r β)) (-E^(-((2 (me re^2 + mp (r + re) (r + re + 2 me re)) γ)/M^2)) + E^(-((2 γ (me (-mp r^2 + re^2 + mp re^2) + (1 + me) mp RealAbs[ r - re]^2))/M^2))) M^2)/(4 (1 + me) mp γ) int[re_] := (2*π)/re norm^2*Integrate[r intrp[r, re], {r, 0, ∞}]; nint[re_?NumericQ, wp_ : MachinePrecision] := (2*π)/re norm^2* NIntegrate[r intrp[r, re], {r, 0, ∞}, WorkingPrecision -> wp]; int[10] (* 4.975964881910113117485834958341383768420770382996*10^-7 *) nint[10] (* 4.97596*10^-7 *) int[10] - nint[10] (* 6.67038*10^-21 *)  Using arbitrary precision rather than MachinePrecision for nint int[10] - nint[10, 49] (* 0.*10^-55 *)  The form of int2 requires extreme precision to obtain an accurate result. Clear["Global*"] M = mp + me; γ = M*ω/2; intrp[r_, re_] := (E^(-2 r (α + r β)) (-E^(-((2 (me re^2 + mp (r + re) (r + re + 2 me re)) γ)/M^2)) + E^(-((2 γ (me (-mp r^2 + re^2 + mp re^2) + (1 + me) mp RealAbs[ r - re]^2))/M^2))) M^2)/(4 (1 + me) mp γ) int2[re_] = (2*π)/re norm^2* Integrate[r intrp[r, re], {r, 0, ∞}] // Simplify; Block[{$MaxExtraPrecision = 2000},
int2[10] /.
Rationalize[{me -> 1,
mp -> 1, α ->
0.50000000206214434417262282295268960297107696533203, β ->
1.1754907046612932983601436560949782261786822346039*10^-8,
norm -> SetPrecision[0.00014216419373200363465950006063, 50], ω ->
10^-4}, 0] // N[#, 250] &] // NumberForm[#, {55, 49}] &

(* N::meprec: Internal precision limit $MaxExtraPrecision = 2000. reached while evaluating -(<<1>>/<<1>>). *)  EDIT: using your expression for int2 provides the same results when using appropriate precision. int2[re_] := (E^(-((2 (me + mp + 2 me mp) re^2 γ)/ M^2)) M^2 norm^2 π^(3/ 2) (-E^(((2 α - (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/M^2)))) M^2 α + E^((2 α + (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/M^2))) M^2 α + 2 E^((2 α - (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/M^2))) mp re γ + 2 E^((2 α + (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/M^2))) mp re γ + 2 E^((2 α - (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/M^2))) me mp re γ + 2 E^((2 α + (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/M^2))) me mp re γ + E^((2 α - (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/ M^2))) M^2 α Erf[(M^2 α - 2 (1 + me) mp re γ)/(Sqrt[ 2] M^2 Sqrt[β + (mp γ)/M^2])] - 2 E^((2 α - (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/ M^2))) mp re γ Erf[(M^2 α - 2 (1 + me) mp re γ)/(Sqrt[ 2] M^2 Sqrt[β + (mp γ)/M^2])] - 2 E^((2 α - (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/ M^2))) me mp re γ Erf[(M^2 α - 2 (1 + me) mp re γ)/(Sqrt[ 2] M^2 Sqrt[β + (mp γ)/M^2])] - E^((2 α + (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/ M^2))) M^2 α Erf[(M^2 α + 2 (1 + me) mp re γ)/(Sqrt[ 2] M^2 Sqrt[β + (mp γ)/M^2])] - 2 E^((2 α + (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/ M^2))) mp re γ Erf[(M^2 α + 2 (1 + me) mp re γ)/(Sqrt[ 2] M^2 Sqrt[β + (mp γ)/M^2])] - 2 E^((2 α + (4 (1 + me) mp re γ)/ M^2)^2/(8 (β + (mp γ)/ M^2))) me mp re γ Erf[(M^2 α + 2 (1 + me) mp re γ)/(Sqrt[ 2] M^2 Sqrt[β + (mp γ)/M^2])]))/(8 Sqrt[ 2] (1 + me) mp re γ (M^2 β + mp γ) Sqrt[β + (mp γ)/M^2]); Block[{$MaxExtraPrecision = 2000},
int2[10] /. γ -> M*ω/2 /. M -> mp + me /.
Rationalize[{me -> 1,
mp -> 1, α ->
0.50000000206214434417262282295268960297107696533203, β ->
1.1754907046612932983601436560949782261786822346039*10^-8,
norm -> SetPrecision[0.00014216419373200363465950006063, 50], ω ->
10^-4}, 0] // N[#, 250] &] // NumberForm[#, {55, 49}] &

(* N::meprec: Internal precision limit \$MaxExtraPrecision = 2000. reached while evaluating <<1>>/<<1>>. *)


• Thanks but I said in the question that int[re] gives the right answer, but when I use the analytic form of its integration (namely int2[re]) it gives the wrong result. You have again used the integrating form for int2[re] which is not my desire one. Now when I replace int2[re] in your answer by the analytic form I get the wrong answer again! Commented Oct 9, 2021 at 17:04
• Please use the analytic form (int2[re] in the question)for trying the solution. Commented Oct 9, 2021 at 17:05
• Your form provides the correct result if appropriate precision is used. Presumably you got your form for int2 by doing the symbolic integration which is exactly what I did for my int2 Commented Oct 9, 2021 at 17:22
• Thanks. Is there a way to specify the desire precision at the first of the file? because I should specify the mp, me and omega and then the appropriate alpha, beta, norm and other constants are produced, so I can't specify their values for each calculation. Also I need to plot this function. How about plot I can produce the right answer? Commented Oct 9, 2021 at 17:36
• Define all the constants with exact numbers when possible and high precision otherwise. To plot just use nint and include the working precision (second argument in my definition above) if needed. Commented Oct 9, 2021 at 18:31