Bug in analytical expression of integral containing Abs function

Question

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.

Answer 1

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

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.

Answer 2

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

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}

Answer 3

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