3
$\begingroup$

I came across a possible bug in Mathematica's evaluation of integrals that I thought I'd share. Consider the code:

zmax = 4;
a =  Integrate[Exp[-x^2] Abs[Sin[x]], {x, -z   , z   }] /. z -> zmax
b = NIntegrate[Exp[-x^2] Abs[Sin[x]], {x, -zmax, zmax}]
a - b

If zmax is smaller than $\pi$ it evaluates close to $0$ (to the precision limit) while for zmax greater than $\pi$ it evaluates to $-4.50908\cdot 10^{-6}$ (which happens to be roughly $4\cdot \int_{\pi}^{\infty}e^{-x^2}\sin(x){\rm d}x$). The condition $z< \pi$ is not listed in Matematica's (11.2.0.0) analytical expression for the integral a above:

ConditionalExpression[( Sqrt[\[Pi]] (2 Erfi[1/2] - Abs[Sin[z]] Csc[z] (Erfi[1/2 - I z] + Erfi[1/2 + I z])))/(2 E^(1/4)), Re[z] >= 0 && Im[z] == 0]

Anybody knows why this happens? If I were to guess, based on the result above, it looks like Mathematica uses the symmetry to rewrite the integral as $2$ times the integral over $[0,z]$ and then for some reason removes the absolute value over $\sin(x)$ (which is allowed only if $z<\pi$), but forgets to mention it.

$\endgroup$
  • 2
    $\begingroup$ @ Winther Just for completeness and because you don't mention where you "came across" this problem. Here's my guess: it was my painful answer to a question in MSE: math.stackexchange.com/questions/2832912/… $\endgroup$ – Dr. Wolfgang Hintze Jun 29 '18 at 17:58
  • $\begingroup$ @Dr.WolfgangHintze Yes. I added the link to the question. $\endgroup$ – Winther Jun 29 '18 at 21:30
  • 1
    $\begingroup$ @ Winther: the discontinuity of the antiderivative at z = n pi was first pointed out by Mariusz Iwaniuk in trying to find the error in my answer. A correct formula for the integral taking into account the jumps in the antiderivatives was given there independently by him and myself. $\endgroup$ – Dr. Wolfgang Hintze Jun 30 '18 at 7:03
2
$\begingroup$

There is a discontinuity in the antiderivative:

Plot[(Sqrt[π] (2 Erfi[1/2] - 
    Abs[Sin[z]] Csc[z] (Erfi[1/2 - I z] + Erfi[1/2 + I z]))) / (2 E^(1/4)),
 {z, 3, 3.25}]

Mathematica graphics

If a concrete limit greater than Pi in magnitude is supplied, Integrate is able to account for the discontinuity; otherwise the general formula is used that is not valid:

Integrate[Exp[-x^2] Abs[Sin[x]], {x, -4, 4}] // N
(*  0.848877 + 0. I  *)

Integrate[Exp[-x^2] Abs[Sin[x]], {x, -z, z}, 
   Assumptions -> -z < x < z] /. z -> 4 // N
(*  0.848873 + 0. I  *)

The use of Assumptions merely removes the ConditionalExpression from the result:

Integrate[Exp[-x^2] Abs[Sin[x]], {x, -z, z}]
(*
  ConditionalExpression[(
   Sqrt[π] (2 Erfi[1/2] - 
      Abs[Sin[z]] Csc[z] (Erfi[1/2 - I z] + Erfi[1/2 + I z])))/(2 E^(1/4)), 
   C[1] ∈ Integers && z >= 2 π C[1] && C[1] >= 0]
*)

The expression gives the erroneous result 0.848873 for C[1] -> 0; otherwise it is undefined.

Further, since Integrate can handle the discontinuity in the concrete case, one is left wondering whether the condition that Abs[z] < Pi could be added to the ConditionalExpression.

$\endgroup$
  • 1
    $\begingroup$ @ Michael E2: Not only at z = pi, the antiderivative has jumps at z = k pi (k=1,2,3..). You might wish to see the original discussion at math.stackexchange.com/questions/2832912/…. Here the correct expression for the integral is given by correcting for the jumps. $\endgroup$ – Dr. Wolfgang Hintze Jun 30 '18 at 7:21
  • $\begingroup$ Thanks for the answer. Interesting that it handles it correctly for a concrete input. The fact that it's discontinuous means it's not a correct antiderivative (an antiderivative is continuous) so the conditions $z<\pi$ should (not just could) be added to result IMO. $\endgroup$ – Winther Jun 30 '18 at 7:30
  • $\begingroup$ @ Micheal E2 You will remember that we had this discussion several times before. I'd like to suggest now that Mathematica gives a warning to the user. Something like "Warning: the result of this integral may be incorrect in the case where the antiderivative contains spurious jumps. Please check for jumps and correct them manually (using Limit)." Such a warning should be issued in cases where there are "conspicious" functions in the integrand, like Abs etc. What do you think? $\endgroup$ – Dr. Wolfgang Hintze Jun 30 '18 at 13:33
  • $\begingroup$ @Dr.WolfgangHintze Yes, but is it really that important? We have seen many times that Integrate struggles with even one discontinuity in the integral. $\endgroup$ – Michael E2 Jun 30 '18 at 17:56
2
$\begingroup$

This is not an answer but an extended comment.

Summary

The results of the study below can be summarized as follows

1) the root cause of problems (i.e. discontinuities in the antiderivative) is the apperance of f[x]Abs[] in the integrand when f[x]!=const. (The EDIT in the end done afterwards points to the necessity of Abs[Sin[x]] or Abs[Cos[x]]) to get in trouble.

2) The antiderivative is returned only if the integral is two-sided AND symmetric.

3) As the one sided integral does not have the numerical problems pointed out in this OP we conclude tentatively (and boldly) that the numerical results are correct if Mathematica can NOT find an anitidrivative. Sounds a little paradoxial ...

4) Replacing Sin[x] by Cos[x] gives an antiderivative with clearly visible jumps

Study

Let us consider some related cases and see if there are problems or not.

Case (a1) replacing the Exp[-x^2] factor by unity: no problems

zmax = 4; 
a = Integrate[ Abs[Sin[x]], {x, -z, z}] /. z -> zmax ; 
b =  NIntegrate[ Abs[Sin[x]], {x, -zmax, zmax}] ; 
a - b

-2.44249*10^-15

Hence the problem is not due to Abs[] alone.

The antiderivative in this case

Integrate[ Abs[Sin[x]], {x, -z, z}]

Out[22]= ConditionalExpression[
 4 IntegerPart[z/\[Pi]] - 
  2 (-1 + Cos[\[Pi] FractionalPart[z/\[Pi]]]) Sign[
    FractionalPart[z/\[Pi]]], z \[Element] Reals]

despite its apperance seems to have no jumps.

Case (a2) Replacing Exp[-x^2] by x^2 : problems

zmax = 4; 
a =  Integrate[x^2 Abs[Sin[x]], {x, -z, z}] /. z -> zmax ; 
b =  NIntegrate[x^2 Abs[Sin[x]], {x, -zmax, zmax}] ; 
a - b

-31.4784

The reason is already known: the antideriative

Integrate[x^2 Abs[Sin[x]], {x, -z, z}]

ConditionalExpression[
 2 (-2 + Abs[Sin[z]] (2 z - (-2 + z^2) Cot[z])), 
 Re[z] >= 0 && Im[z] == 0]

has jumps at z = k pi.

This is perhaps the most elementary example exhibiting the problems here.

Case (a3) Replacing Exp[-x^2] by Abs[x]: problems

zmax = 4; a = 
 Integrate[Abs[x] Abs[Sin[x]], {x, -z, z}] /. z -> zmax ; b = 
 NIntegrate[Abs[x] Abs[Sin[x]], {x, -zmax, zmax}] ; a - b

Out[21]= -12.5664

The antiderivative

Integrate[Abs[x] Abs[Sin[x]], {x, -z, z}]

ConditionalExpression[-2 Abs[Sin[z]] (-1 + z Cot[z]), 
 Re[z] >= 0 && Im[z] == 0]

has jumps at z = k pi.

Case (b) the one-sided integral: no problems

zmax = 4; 
a =  Integrate[Exp[-x^2] Abs[Sin[x]], {x, 0, z}] /. z -> zmax ; 
b =  NIntegrate[Exp[-x^2] Abs[Sin[x]], {x, 0, zmax}] ; a - b

6.05072*10^-15 + 0. I

This is surprising (at least for me).

It is intesting that neither for the one sided integral

Integrate[Exp[-x^2] Abs[Sin[x]], {x, 0, z}]

nor for the indefinite integral

Integrate[Exp[-x^2] Abs[Sin[x]], x]

nor for the unsymmtric two sided integral, e.g.

Integrate[E^-x^2 Abs[Sin[x]], {x, -z, 2 z}]

Mathematica provides an antiderivative but returns the input unchanged.

For comparison we repeat here the antidrivative of the two sided integral:

Integrate[Exp[-x^2] Abs[Sin[x]], {x, -z, z}]

Out[30]= ConditionalExpression[(
 Sqrt[\[Pi]] (2 Erfi[1/2] - 
    Abs[Sin[z]] Csc[z] (Erfi[1/2 - I z] + Erfi[1/2 + I z])))/(
 2 E^(1/4)), Re[z] >= 0 && Im[z] == 0]

EDIT 30.08.18 15:15

case "truncated sin"

This case makes me wonder if my conclusions are correct.

Taking a "truncated sine" x(1 - (x/\[Pi])^2) instead of Sin[x] gives an antideriative without jumps:

Integrate[E^-x^2 Abs[x (1 - (x/\[Pi])^2)], {x, -z, z}]

Out[53]= ConditionalExpression[(
 E^-z^2 (1 - \[Pi]^2 + E^z^2 (-1 + \[Pi]^2) + z^2))/\[Pi]^2, 
 0 < Re[z] < \[Pi] && Im[z] == 0]

Cos instaed of Sin

If we replace Sin[x] by Cos[x] there is no hiding of tiny jumps anymore

Integrate[E^-x^2 Abs[Cos[x]], {x, -z, z}]

Out[66]= ConditionalExpression[(
 I Sqrt[\[Pi]]
   Abs[Cos[z]] (Erfi[1/2 - I z] - Erfi[1/2 + I z]) Sec[z])/(
 2 E^(1/4)), \[Pi] + z + Conjugate[z] <= 0]

The graph is

enter image description here

The jumps are at the zeroes of Cos[], i.e. at z= (2 k + 1) \[Pi]/2.

And the jump sizes are

j[k_] := -((
  I Sqrt[\[Pi]] (Erfi[1/2 - ((2 k + 1) I \[Pi])/2] - 
     Erfi[1/2 + ((2 k + 1) I \[Pi])/2]))/(2 E^(1/4)))

The sequnce of j[k] is extremely slowing decreasing

Table[N[j[k], 30], {k, 0, 4}] // Chop

{
-1.39137915049102862583968181080, 
-1.380388447038392527457349263471, 
-1.380388447043142974773415246738, 
-1.380388447043142974773415246726, 
-1.380388447043142974773415246726
}

Case explizit numbers for integration limits

Explizit numbers instead of a variable z leads to numerically correct results.

Example is the comparison of a sum with its integral representation

Let

sN = 
 NSum[Exp[-n^2]/(1 - 4 n^2), {n, -\[Infinity], \[Infinity]}, 
  WorkingPrecision -> 20]

Out[210]= 0.7522978984722431448

Then the sum is increasingly well approximated by the integral

Table[N[Sqrt[\[Pi]]/2 Integrate[E^-x^2 Abs[Sin[x]], {x, -n, n}] - sN, 
   20], {n, 1, 6}] 

{-0.22995897048993183416 + 0.*10^-21 I,
 -0.0058037201019733100 +  0.*10^-23 I,
 -3.67451669256173*10^-6 + 0.*10^-26 I,
 -1.99280197442*10^-8 + 0.*10^-28 I,
 -2.2294568*10^-12 + 0.*10^-32 I,
 -6.9*10^-18 + 0.*10^-37 I}
$\endgroup$
  • $\begingroup$ You can try for free MMA 11.3 in wolfram cloud: develop.wolframcloud.com only you must sign-in. I use because they have a better computer than mine->28 Cores 32 GB RAM :). Limitation is time(or another things, well I don't know) to use,about 5 minutes max. $\endgroup$ – Mariusz Iwaniuk Jun 30 '18 at 9:31
  • 1
    $\begingroup$ @ Mariusz Iwaniuk Thanks for the hint. It works fine. $\endgroup$ – Dr. Wolfgang Hintze Jun 30 '18 at 12:32
2
$\begingroup$

OK when the right syntax (RealAbs instead of Abs ) is used:

zmax = 4;a = Integrate[Exp[-x^2] RealAbs[Sin[x]], {x, -z, z}] /. z -> zmax
(-I Sqrt[π] Erf[4 + I/2] + 2 Sqrt[π] Erfi[1/2] + 
 Sqrt[π] Erfi[1/2 + 4 I] - 2 Sqrt[π] Erfi[1/2 - I π] - 
 2 Sqrt[π] Erfi[1/2 + I π])/(2 E^(1/4))
b = NIntegrate[Exp[-x^2] RealAbs[Sin[x]], {x, -zmax, zmax}]

0.84887725359944

a - b

1.21014*10^-14 + 0. I

$\endgroup$
  • $\begingroup$ Thanks for the reply. I would not say RealAbs is the right syntax and Abs is wrong, both should give correct results! But based on this it seems RealAbs is safer than Abs and should be the recommended choice for real integrals which is something I'll keep in mind $\endgroup$ – Winther Jun 30 '18 at 7:36
  • $\begingroup$ @Winter: Your claim "I would not say RealAbs is the right syntax and Abs is wrong, both should give correct results! " does not correspond to reality. For example, compare the results of Integrate[RealAbs[x],x] and Integrate[Abs[x],x]. $\endgroup$ – user64494 Jun 30 '18 at 8:37
  • $\begingroup$ My comment was not on what Mathematica outputs, but what is mathematically correct. Abs[x] and RealAbs[x] are equivalent when x is real. $\endgroup$ – Winther Jun 30 '18 at 8:38
  • $\begingroup$ It’s perfectly fine that using one function does not give an answer while the other one does (as in your example). It’s not fine that both choices gives an answer which disagree with each other. That is a bug. $\endgroup$ – Winther Jun 30 '18 at 8:52
  • $\begingroup$ @ Winther Sorry for intervening, but - as I have no Version 11 - could you please just repeat your cases with RealAbs instead of Abs, and also calculate the antiderivative with RealAbs? $\endgroup$ – Dr. Wolfgang Hintze Jun 30 '18 at 9:16

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.