# Why NIntegrate does not give correct answer for smaller values of input

(In the following code NIntegration does not give correct result when input falls below 0.9.)

mu22[k_] := NIntegrate[
z  *1/Sqrt[
2 π] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(
2 (3 + k*z))]*(Erf[z/Sqrt[2 (1 + (k*z)/3)]] +
Exp[18/k^2]*Erf[-(z + 6/k)/Sqrt[2 (1 + (k*z)/3)]]),
{z, -3/k, ∞}];


Fails:

mu22[.5]


NIntegrate::ncvb: NIntegrate failed to converge...near {z} = {1.39454}....

(* Out: 8.433254185*10^14 *)


Succeeds:

mu22[.9]
(* Out: 0.5493575446 *)

• its (mathematically) ok as z->-3/k because Erf[Infinity] is 1. This is part of the difficulty numerically integrating the thing though. May 7, 2018 at 18:04

The integral over negative z has a large positive value and the integral over the positive domain has a large nearly equal negative value. I managed to get high precision results by breaking the integral into two parts:

mu22[k_] :=
NIntegrate[
z*1/Sqrt[2 \[Pi]] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(2 (3 + k*z))]*
(Erf[z/Sqrt[2 (1 + (k*z)/3)]] +
Exp[18/k^2]*Erf[-(z + 6/k)/Sqrt[2 (1 + (k*z)/3)]]),
{z, -3/k, 0}, WorkingPrecision -> 50] +
NIntegrate[
z*1/Sqrt[2 \[Pi]] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(2 (3 + k*z))]*
(Erf[z/Sqrt[2 (1 + (k*z)/3)]] +
Exp[18/k^2]*Erf[-(z + 6/k)/Sqrt[2 (1 + (k*z)/3)]]),
{z, 0, Infinity}, WorkingPrecision -> 50]
Table[mu22[i/10], {i, 1, 9}]


{err, err, err, 0.6, 0.55939805, 0.557353,0.55498491455072,0.55231,0.54935762525154}

I'm not convinced this is a correct result, but it seems reasonable. The "exact" .6 is interesting..

• It is correct for 5/10 to 9/10. For 0.1, it should be 0.56399, and in fact for k=0 it is 0.56419. Thanks . May 7, 2018 at 17:37
• It "seems" the Exp[18/k^2]*Erf..  term integrates to zero. No luck trying to show that in a convincing way however. May 7, 2018 at 18:34
• Is it possible to increase the working precision? I guess that should resolve the issue upto some extent. May 7, 2018 at 19:11
• playing with working precision is how I arrived at "It seems". I guess it depends on your application and how confident you need to be if you want to believe that kind of numeric result. May 7, 2018 at 21:04
• The "exact" 0.6 is actually 0.56109912691476208930.9908795242543437, which suggests the input has reached the boundary of catastrophic precision loss. May 7, 2018 at 22:10

As @george2079 says, the second term in the parenthesis integrates to zero:

Assuming[k > 0,
Integrate[z*1/Sqrt[2π] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(2 (3 + k*z))]*
(Exp[18/k^2]*Erf[-(z + 6/k)/Sqrt[2 (1 + (k*z)/3)]]), {z, -3/k, ∞}]]


0

So we only need to integrate the first term in the parenthesis, which is easy:

mu22[k_] := NIntegrate[z*1/Sqrt[2π] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(2 (3 + k*z))]*
(Erf[z/Sqrt[2 (1 + (k*z)/3)]]), {z, -3/k, ∞}]

Table[mu22[k], {k, 0.1, 1, 0.1}]


{0.563994, 0.563409, 0.562441, 0.561099, 0.559398, 0.557354, 0.554985, 0.552312, 0.549358, 0.546144}

• (+1) Please comment on the methodological similarity with the approach your took here. May 27, 2019 at 22:33
• @AntonAntonov I strongly believe that before we explore numerical options (WorkingPrecision, Method, etc.) we must exhaust analytic transformations. In other words, instead of using a faster car on the exact route proposed by the OP, first check if we can find a shortcut. May 28, 2019 at 4:46

Sort of an application of the approaches listed in this answer of "NIntegrate doesn't evaluate".

We increase the working precision, specify modest precision goal (that corresponds to the machine precision computation), and we specify large enough MinRecursion.

Clear[mu22]
mu22[k_, opts : OptionsPattern[]] :=
NIntegrate[z*(1/Sqrt[2*Pi])*(3/(3 + z*k))^(3/2)*
Exp[(-3*z^2)/(2*(3 + k*z))]*(Erf[z/Sqrt[2*(1 + (k*z)/3)]] +
Exp[18/k^2]*Erf[-(z + 6/k)/Sqrt[2*(1 + (k*z)/3)]]),
{z, -3/k, Infinity},
opts, PrecisionGoal -> 6, MinRecursion -> 4, WorkingPrecision -> 30];

mu22[4/10, MinRecursion -> 6, WorkingPrecision -> 50]
(* 0.56109912691473474278093871991934804555057780817151 *)

mu22[5/10]
(* 0.559398049241961146450206736125 *)

mu22[6/10]
(* 0.557353711744127129697499796898 *)

mu22[9/10]
(* 0.549357625251542853028279744451 *)


The comparison computation from the previous answer:

Table[TimeConstrained[mu22[i/10, MinRecursion -> 6, WorkingPrecision -> 50], 200], {i, 1, 9}]

(* During evaluation of In[19]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. *)

(* During evaluation of In[19]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {1.095980305962580454464980064514005730364537970159427154551846177428701345757689625003204536254789809}. NIntegrate obtained -5.055602990319483773686220288847878705656214979951156205424982827525690691913603370306862497743044896*10^25 and 9.297088667477530482343599956196447156134025802833492431136165342000934397235532650950740214993203116100.*^27 for the integral and error estimates. *)

(* {$$Aborted,$$Aborted, \
-5.0556029903194837736862202888478787056562149799512*10^25, \
0.56109912691473474278093871991934804555057780817151, \
0.55939805009602787878559785757854780600727159330663, \
0.55735371174412712980920126629089419941485487591229, \
0.55498491455072038849610463055656444825264621013418, \
0.55231226836469861118654749178597901122820612288778, \
0.54935762525154285303487507916956834797845874625737} *)