0
$\begingroup$

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.975964881910113112577562173393203988238793355023312`27.\
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.

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

1 Answer 1

1
$\begingroup$
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>>). *)

enter image description here

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>>. *)

enter image description here

$\endgroup$
5
  • $\begingroup$ 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! $\endgroup$
    – Wisdom
    Commented Oct 9, 2021 at 17:04
  • $\begingroup$ Please use the analytic form (int2[re] in the question)for trying the solution. $\endgroup$
    – Wisdom
    Commented Oct 9, 2021 at 17:05
  • $\begingroup$ 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 $\endgroup$
    – Bob Hanlon
    Commented Oct 9, 2021 at 17:22
  • $\begingroup$ 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? $\endgroup$
    – Wisdom
    Commented Oct 9, 2021 at 17:36
  • $\begingroup$ 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. $\endgroup$
    – Bob Hanlon
    Commented Oct 9, 2021 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.