4
$\begingroup$

(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 *)
$\endgroup$
  • $\begingroup$ its (mathematically) ok as z->-3/k because Erf[Infinity] is 1. This is part of the difficulty numerically integrating the thing though. $\endgroup$ – george2079 May 7 '18 at 18:04
3
$\begingroup$

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

|improve this answer|||||
$\endgroup$
  • $\begingroup$ 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 . $\endgroup$ – nutan May 7 '18 at 17:37
  • 1
    $\begingroup$ It "seems" the Exp[18/k^2]*Erf.. term integrates to zero. No luck trying to show that in a convincing way however. $\endgroup$ – george2079 May 7 '18 at 18:34
  • $\begingroup$ Is it possible to increase the working precision? I guess that should resolve the issue upto some extent. $\endgroup$ – nutan May 7 '18 at 19:11
  • $\begingroup$ 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. $\endgroup$ – george2079 May 7 '18 at 21:04
  • $\begingroup$ The "exact" 0.6 is actually 0.5610991269147620893`0.9908795242543437, which suggests the input has reached the boundary of catastrophic precision loss. $\endgroup$ – Michael E2 May 7 '18 at 22:10
2
$\begingroup$

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}

|improve this answer|||||
$\endgroup$
  • $\begingroup$ (+1) Please comment on the methodological similarity with the approach your took here. $\endgroup$ – Anton Antonov May 27 '19 at 22:33
  • $\begingroup$ @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. $\endgroup$ – Roman May 28 '19 at 4:46
1
$\begingroup$

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.297088667477530482343599956196447156134025802833492431136165342000934397235532650950740214993203116`100.*^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} *)
|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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