0
$\begingroup$

I'm trying to evaluate the following code

mp = 50;
\[Omega] = 10^-2;
me = 1;
M = mp + me;
\[Gamma] = 1/2*M*\[Omega];
{\[Alpha], \[Beta]} = \
{0.98023965221867725983173613713006488978862762451171875`50., 
   0.00014974969686706935865193324186606105286045931279659271240234375\
`50.};
norm = 0.140101117848786449079471089613455306450;

funE[re_] := 
  1/(8 Sqrt[2]
     mp re \[Gamma] (M^2 \[Beta] + mp^2 \[Gamma]) Sqrt[\[Beta] + (
     mp^2 \[Gamma])/M^2]) E^(-2 re^2 \[Gamma]) M^2 norm^2 \[Pi]^(
   3/2) (E^((M \[Alpha] - 2 mp re \[Gamma])^2/(
      2 (M^2 \[Beta] + 
         mp^2 \[Gamma]))) ((-1 + E^((4 M mp re \[Alpha] \[Gamma])/(
           M^2 \[Beta] + mp^2 \[Gamma]))) M \[Alpha] + 
        2 (1 + E^((4 M mp re \[Alpha] \[Gamma])/(
           M^2 \[Beta] + mp^2 \[Gamma]))) mp re \[Gamma]) + 
     E^((M \[Alpha] - 2 mp re \[Gamma])^2/(
      2 (M^2 \[Beta] + mp^2 \[Gamma]))) (M \[Alpha] - 
        2 mp re \[Gamma]) Erf[(M \[Alpha] - 2 mp re \[Gamma])/(
       Sqrt[2] M Sqrt[\[Beta] + (mp^2 \[Gamma])/M^2])] - 
     E^((M \[Alpha] + 2 mp re \[Gamma])^2/(
      2 (M^2 \[Beta] + mp^2 \[Gamma]))) (M \[Alpha] + 
        2 mp re \[Gamma]) Erf[(M \[Alpha] + 2 mp re \[Gamma])/(
       Sqrt[2] M Sqrt[\[Beta] + (mp^2 \[Gamma])/M^2])]);
ksorbE[re_] := Sqrt[funE[re]];
lapE[re_] := D[ksorbE[re], {re, 2}] + 2/re D[ksorbE[re], re];
res = -(1/(2*me))*2*\[Pi]*norm^2*
  NIntegrate[
   re^2*Sin[\[Theta]]*lapE[re]*ksorbE[re], {\[Theta], 0, \[Pi]}, {re, 
    0, \[Infinity]}]

where I want to obtain the result of res but it doesn't return any result. I see in the error dialog box that Mathematica assumes boundaries of re 0 to 0.35 and I really don't know how it reaches such boundary! Any idea?

Addendum

The way that Akku14 suggested works just for above case but it fails for the following case:

mp = 1836;
\[Omega] = 10^-2;
me = 1;
M = mp + me;
\[Gamma] = 1/2*M*\[Omega];
{\[Alpha], \[Beta]} = {66833261/66879672, 1050947/7018024223};
norm = (16*\[Pi]^2*(
     Sqrt[\[Pi]/
      2] (-2 \[Alpha] Sqrt[\[Beta]] + 
        E^(\[Alpha]^2/(2 \[Beta])) Sqrt[
         2 \[Pi]] (\[Alpha]^2 + \[Beta]) Erfc[\[Alpha]/(
          Sqrt[2] Sqrt[\[Beta]])]))/(
     128 \[Beta]^(5/2) \[Gamma]^(3/2)))^(-1/2) // Rationalize[#, 0] &;
funP[rp_] = (E^(-2 rp^2 \[Gamma]) M^2 norm^2 \[Pi]^(
      3/2) (E^((M \[Alpha] - 2 me rp \[Gamma])^2/(
         2 (M^2 \[Beta] + 
            me^2 \[Gamma]))) ((-1 + E^((4 M me rp \[Alpha] \[Gamma])/(
              M^2 \[Beta] + me^2 \[Gamma]))) M \[Alpha] + 
           2 (1 + E^((4 M me rp \[Alpha] \[Gamma])/(
              M^2 \[Beta] + me^2 \[Gamma]))) me rp \[Gamma]) + 
        E^((M \[Alpha] - 2 me rp \[Gamma])^2/(
         2 (M^2 \[Beta] + me^2 \[Gamma]))) (M \[Alpha] - 
           2 me rp \[Gamma]) Erf[(M \[Alpha] - 2 me rp \[Gamma])/(
          Sqrt[2] M Sqrt[\[Beta] + (me^2 \[Gamma])/M^2])] - 
        E^((M \[Alpha] + 2 me rp \[Gamma])^2/(
         2 (M^2 \[Beta] + me^2 \[Gamma]))) (M \[Alpha] + 
           2 me rp \[Gamma]) Erf[(M \[Alpha] + 2 me rp \[Gamma])/(
          Sqrt[2] M Sqrt[\[Beta] + (me^2 \[Gamma])/M^2])]))/(8 Sqrt[2]
       me rp \[Gamma] (M^2 \[Beta] + me^2 \[Gamma]) Sqrt[\[Beta] + (
       me^2 \[Gamma])/M^2]) // Simplify[#, rp > 0] &;

ksorbP[rp_] = Sqrt[funP[rp]];
lapP[rp_] = D[ksorbP[rp], {rp, 2}] + 2/rp D[ksorbP[rp], rp];

facP = -(1/(2*mp))*2*\[Pi]*norm^2*
   Integrate[Sin[\[Theta]], {\[Theta], 0, \[Pi]}];

fP[re_] = rp^2*lapP[rp]*ksorbP[rp] // Simplify[#, rp > 0] &;
resP = facP* 
  NIntegrate[fP[rp], {rp, 0, \[Infinity]}, WorkingPrecision -> 100, 
   MaxRecursion -> 70, AccuracyGoal -> 20]
$\endgroup$
2
  • 1
    $\begingroup$ Look at function funE[re], Plot[Evaluate[funE[re] ], {re, 0, 15}, PlotRange -> All] shows numerical instability $\endgroup$ Feb 17, 2022 at 9:10
  • 1
    $\begingroup$ You should really at least test your functions you are using inside the integrand on some numerical values to make sure they give something meaningful. For example see what lapE[1] and lapE[0] give,. Also see what ksorbE[0] give. $\endgroup$
    – Nasser
    Feb 17, 2022 at 9:11

3 Answers 3

2
$\begingroup$

This is a problem of numerical precision. Rationalize parameters.

mp = 50;
\[Omega] = 10^-2;
me = 1;
M = mp + me;
\[Gamma] = 1/2*M*\[Omega];
{\[Alpha], \[Beta]} = {0.98023965221867725983173613713006488978862762451171875`50., 
0.00014974969686706935865193324186606105286045931279659271240234375`50.} // Rationalize[#, 0] &;
norm = 0.140101117848786449079471089613455306450 //Rationalize[#, 0] &;

funE[re_] = 
1/(8 Sqrt[
    2] mp re \[Gamma] (M^2 \[Beta] + 
     mp^2 \[Gamma]) Sqrt[\[Beta] + (mp^2 \[Gamma])/
      M^2]) E^(-2 re^2 \[Gamma]) M^2 norm^2 \[Pi]^(3/
   2) (E^((M \[Alpha] - 
         2 mp re \[Gamma])^2/(2 (M^2 \[Beta] + 
          mp^2 \[Gamma]))) ((-1 + 
        E^((4 M mp re \[Alpha] \[Gamma])/(M^2 \[Beta] + 
             mp^2 \[Gamma]))) M \[Alpha] + 
     2 (1 + E^((4 M mp re \[Alpha] \[Gamma])/(M^2 \[Beta] + 
             mp^2 \[Gamma]))) mp re \[Gamma]) + 
  E^((M \[Alpha] - 
         2 mp re \[Gamma])^2/(2 (M^2 \[Beta] + 
          mp^2 \[Gamma]))) (M \[Alpha] - 
     2 mp re \[Gamma]) Erf[(M \[Alpha] - 
       2 mp re \[Gamma])/(Sqrt[
        2] M Sqrt[\[Beta] + (mp^2 \[Gamma])/M^2])] - 
  E^((M \[Alpha] + 
         2 mp re \[Gamma])^2/(2 (M^2 \[Beta] + 
          mp^2 \[Gamma]))) (M \[Alpha] + 
     2 mp re \[Gamma]) Erf[(M \[Alpha] + 
       2 mp re \[Gamma])/(Sqrt[
        2] M Sqrt[\[Beta] + (mp^2 \[Gamma])/M^2])]) // 
Simplify[#, re > 0] &;

ksorbE[re_] = Sqrt[funE[re]];

lapE[re_] = D[ksorbE[re], {re, 2}] + 2/re D[ksorbE[re], re];

fac = -(1/(2*me))*2*\[Pi]*norm^2*
  Integrate[Sin[\[Theta]], {\[Theta], 0, \[Pi]}]

f[re_] = re^2*lapE[re]*ksorbE[re] // Simplify[#, re > 0] &;

Plot[fac f[re], {re, 0, 15}, WorkingPrecision -> 30, PlotRange -> All]

res = fac NIntegrate[f[re], {re, 0, \[Infinity]}, 
     WorkingPrecision -> 50, MaxRecursion -> 50, AccuracyGoal -> 10]

(*   0.0038050253906429147971294149529417266291524768497926   *)
$\endgroup$
5
  • $\begingroup$ Many thanks, your result is a bit less than when I set 0<re<6 in code of the question (0.003820). Regarding highly oscillations in re>12 do you think these oscillations are responsible for this decrease? $\endgroup$
    – Wisdom
    Feb 17, 2022 at 16:15
  • $\begingroup$ There are no oscillations ever for re>12, they were an artefact due to low workingPrecision. See my Plot with high WorkingPrecision. You loose a little negative part, if you stop at re=06. $\endgroup$
    – Akku14
    Feb 17, 2022 at 16:23
  • $\begingroup$ Oh yes, you're right, I got your mean. $\endgroup$
    – Wisdom
    Feb 17, 2022 at 16:25
  • $\begingroup$ unfortunately your method fails for other cases and I encounter the same problem again! Do you have another suggestion? I add it to my question as an addendum. $\endgroup$
    – Wisdom
    Feb 17, 2022 at 16:40
  • $\begingroup$ Curiously, when I run your code on a fresh kernel I get a slightly different result starting with the 9th significant digit: 0.00380502535... instead of 0.00380502539..... $\endgroup$
    – theorist
    Feb 19, 2022 at 7:47
1
$\begingroup$

modified

With the suggestions (2 times //Rationalize) from @Akku14 very good answer this approach (WorkingPrecison&AccuracyGoal not necessary) with decreased integration range gives similar results!

    E^((4 M mp re \[Alpha] \[Gamma])/(M^2 \[Beta] + 
             mp^2 \[Gamma]))) M \[Alpha] + 
     2 (1 + E^((4 M mp re \[Alpha] \[Gamma])/(M^2 \[Beta] + 
             mp^2 \[Gamma]))) mp re \[Gamma]) + 
  E^((M \[Alpha] - 
         2 mp re \[Gamma])^2/(2 (M^2 \[Beta] + 
          mp^2 \[Gamma]))) (M \[Alpha] - 
     2 mp re \[Gamma]) Erf[(M \[Alpha] - 
       2 mp re \[Gamma])/(Sqrt[
        2] M Sqrt[\[Beta] + (mp^2 \[Gamma])/M^2])] - 
  E^((M \[Alpha] + 
         2 mp re \[Gamma])^2/(2 (M^2 \[Beta] + 
          mp^2 \[Gamma]))) (M \[Alpha] + 
     2 mp re \[Gamma]) Erf[(M \[Alpha] + 
       2 mp re \[Gamma])/(Sqrt[
        2] M Sqrt[\[Beta] + (mp^2 \[Gamma])/M^2])]) // 
Simplify[#, re > 0] &;

ksorbE[re_] = Sqrt[funE[re]];

lapE[re_] = D[ksorbE[re], {re, 2}] + 2/re D[ksorbE[re], re];

First you might simplify the integration by integrating the part Sin[\[Theta]], which might be separated, first

Integrate[Sin[\[Theta]], {\[Theta], 0, Pi}]
(* 2 *)

Remains

res = Function[inf, -(1/(2*me))*2*\[Pi]*   norm^2*2 NIntegrate[re^2*lapE[re]*ksorbE[re], {re, 0, inf} , Method -> "LocalAdaptive" ]   ]

which calculates the integral res[inf] in the range 0<re<inf

ListPlot indicates

ListPlot[Table[{ inf, res[inf]}, {inf, 1, 10}], AxesLabel -> {"inf", "res[inf]"}, GridLines -> {None, {res[10]}}]

enter image description here

a useful convergence 0.0038417.

$\endgroup$
3
  • $\begingroup$ Thanks a lot however it seems this method can't recover the correct result. It fails for inf>10 and when I set inf=10 for example it returns 0.003786 which is less than results of @Akku14 's answer and setting 0<re<6 in the question. $\endgroup$
    – Wisdom
    Feb 17, 2022 at 16:11
  • $\begingroup$ @Wisdom It seems only a factor 2 difference between the two solutions, I will check it later $\endgroup$ Feb 17, 2022 at 16:44
  • $\begingroup$ Thanks, Can I replace funE[re] by its numerical integration instead of its analytical form? I wrote its analytical form because I wasn't sure that D[] command is applied to numerical integration. $\endgroup$
    – Wisdom
    Feb 17, 2022 at 17:11
0
$\begingroup$

this is a nasty function you want to integrate.. If this is really the function you'd expect

verify with Plot[re^2*Sin[0.01]*lapE[re]*ksorbE[re], {re, 25, 25.01}]): enter image description here

I think you'd need some smart regularization of this stochastic part of the integrand. Numerical integration works well when the integrand is bounded and of order 1. Your function oscillates between 10^18 and -10^18 in an interval of 10^-8. Do you know if this integral is supposed to converge?

I fear you can't just rely on the automatic numerical integration of NIntegrate but need to something "smart" about your function.

  • Do you know any statistical properties of the stochastic/oscillatory region?
  • Do you know analytical expressions of the integral?
$\endgroup$
1
  • $\begingroup$ Thanks, as @UlrichNeumann said funE[re] oscillates about re>12 and when I limit the boundary to 6 for example it returns a reasonable result. Let me clarify: funE[re] is a density function and must be have all properties of a density probability function. Also by numerical integration I know that funE[re] decays about re=5. $\endgroup$
    – Wisdom
    Feb 17, 2022 at 9: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.