I am experimenting a little bit with the error function and understanding its behavior. For example, by the asymptotic behavior of the error function, the following integral should converge (and its absolute square be constant!). However, mathematica gives me, for example for $t=0.2$, that the integral diverges:

int2[t_] :=  NIntegrate[Erf[Sqrt[(-1/2 - I/2*Cot[t])]*x/Sqrt[2]]*Exp[-x^2/4], {x, 0,Infinity}, WorkingPrecision -> 1000]

For $t=0.19$ however, it converges (as it does for most $t$). Is this due to WorkingPrecision?

  • $\begingroup$ Have you seen the result of Assuming[u > 0, Integrate[Erf[Sqrt[(-1/2 - I/2*u)]*x/Sqrt[2]]*Exp[-x^2/4], {x, 0, Infinity}]]? $\endgroup$ – J. M.'s technical difficulties May 18 '17 at 20:52
  • $\begingroup$ Thanks for the answer! I have, but I wonder where this behaviour in NIntegrate comes from, as this integral is not the actual one I am interested in. The one which I am interested in cannot be done analytically, I think. I saw a similar behaviour in NIntegrate, though, so I cooked it down to this simple integral... Is this caused by 'WorkingPrecision'? $\endgroup$ – Gesbesgue May 18 '17 at 21:07
  • $\begingroup$ WorkingPrecision -> 1000 causes the precision goal to be half that or 500 digits. Do you really want such a high-precision result? It takes forever to run, even if I set PrecisionGoal -> 2 and it discourages me from exploring the situation. $\endgroup$ – Michael E2 May 19 '17 at 10:38
  • $\begingroup$ If I put it lower or leave it out, the integral diverges, which cannot be. $\endgroup$ – Gesbesgue May 19 '17 at 14:41

First we should distinguish the convergence of the integral from the convergence of the numerical procedure used to compute the integral. If the integral diverges, then the numerical procedure will diverge. However, if the integral converges, it is still possible for the numerical procedure to fail to converge to the value of the integral. NIntegrate complains about the failure of the numerical procedure to converge.

The problem with convergence has to do with both precision and the integration method selected. The principal problem is with the method selection. When t is 0.2, the Gauss-Kronrod rule is used and the integration fails. If t is 20/100, then the Levin rule is used and the integration succeeds.

The error function grows quite large quite quickly for t around 0.2, and although it is multiplied by a much smaller gaussian, there seems to be significant precision loss. This may be showing up in the error estimation in the Gauss-Kronrod rule, causing NIntegrate to subdivide the interval {0, Infinity} more, further exacerbating the problem. If the recursion is allowed to go far enough (MaxRecursion -> 10 usually does it), you get overflow. It may be hopeless to try to use Gauss-Kronrod.

On the other hand, the Levin rule seems to work reliably. So one should force NIntegrate to use it. We can actually do the symbolic preprocessing for all t using LevinIntegrandReduce, which may be found in the documentation on the Levin rule. This saves some time, too. Here is the code:

With[{li = NIntegrate`LevinIntegrandReduce[
    Erf[Sqrt[(-1/2 - I/2*Cot[t])]*x/Sqrt[2]] Exp[-x^2/4], x]},
 int3[t0_] := Block[{t = t0},
     Erf[Sqrt[(-1/2 - I/2*Cot[t])]*x/Sqrt[2]] Exp[-x^2/4],
     {x, 0, Infinity},
     Method -> {"LevinRule", "AdditiveTerm" -> li@"AdditiveTerm", 
       "Amplitude" -> li@"Amplitude", "Kernel" -> li@"Kernel", 
       "DifferentialMatrix" -> First@li@"DifferentialMatrices"}]


int2[2/10] // RepeatedTiming  (* OP's *)
int3[0.2] // RepeatedTiming
  {0.066, 1.42631 - 0.370151 I}
  {0.0233, 1.42631 - 0.370151 I}

How to see what rule is used:

int2[t_] := NIntegrate[Erf[Sqrt[(-1/2 - I/2*Cot[t])]*x/Sqrt[2]]*Exp[-x^2/4],
  {x, 0, Infinity},
  IntegrationMonitor :> print]

 print[x_] := (Print[Head@First[x]@"GetRule"]; Clear[print]);
 print[x_] := (Print[Head@First[x]@"GetRule"]; Clear[print]);


(*  1.42631 - 0.370151 I  *)

NIntegrate`GeneralRule (Gauss-Kronrod is inferred)

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections....

(*  1.83887 - 0.382448 I  *)

If you doubt which rule is being used with NIntegrate`GeneralRule, change print to just Print[x] and compare the abscissae and weights at the end with those returned by NIntegrate`GaussKronrodRuleData[5, MachinePrecision].

|improve this answer|||||

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.