# Why this numerical integration does not a result?

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.9802396522186772598317361371300648897886276245117187550.,
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?

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]

• Look at function funE[re], Plot[Evaluate[funE[re] ], {re, 0, 15}, PlotRange -> All] shows numerical instability Commented Feb 17, 2022 at 9:10
• 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. Commented Feb 17, 2022 at 9:11

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.9802396522186772598317361371300648897886276245117187550.,
0.0001497496968670693586519332418660610528604593127965927124023437550.} // 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   *)

• 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? Commented Feb 17, 2022 at 16:15
• 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. Commented Feb 17, 2022 at 16:23
• Oh yes, you're right, I got your mean. Commented Feb 17, 2022 at 16:25
• 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. Commented Feb 17, 2022 at 16:40
• 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..... Commented Feb 19, 2022 at 7:47

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]}}]


a useful convergence 0.0038417.

• 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. Commented Feb 17, 2022 at 16:11
• @Wisdom It seems only a factor 2 difference between the two solutions, I will check it later Commented Feb 17, 2022 at 16:44
• 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. Commented Feb 17, 2022 at 17:11

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}]):

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?
• 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. Commented Feb 17, 2022 at 9:31