A discontinuity in an integral

Question

I have a function which I'm trying to integrate, which is perfectly smooth and Mathematica sees it as well defined everywhere in the region over which I'm integrating it. The function is:

 Sqrt[((1/2)^2 - Csc[x]^2)/(-1 + (1/2)^2 Cos[x]^2)^3]Sin[x]^2

and the graph of this over 0 to pi looks like

If I calculate the indefinite integral of this function it gives me a discontinuous result, as it seems to think that there is a problem at x=pi/2. The integral of the function is:

Indeed this is an antiderivative of the original function other than at pi/2 where the integral is not defined. I have tried to add an integration constant on the piecewise part for x>pi/2, but I can't take the limit of the function as x->pi/2 from below. Mathematica tells me that it's infinite where clearly the limit from below seems to be well defined.

Is there any way I can put a constraint on the indefinite integral to give me the natural, smooth antiderivative which is defined everywhere?

Answer 1

$Version

(* "12.0.0 for Mac OS X x86 (64-bit) (April 7, 2019)" *)

f[x_] = Sqrt[((1/2)^2 - Csc[x]^2)/
     (-1 + (1/2)^2 Cos[x]^2)^3] Sin[x]^2;

The indefinite integral is

if[x_] = Integrate[f[x], x] // Simplify[#, Element[x, Reals]] &

(* 1/2 Abs[Csc[x]] (Sqrt[49 - Cos[2 x]^2]/(-7 + Cos[2 x]) - (
   2 EllipticE[ArcSin[Sqrt[7 - Cos[2 x]]/Sqrt[14]], 7/4])/Abs[Cos[x]]) Sin[
  2 x] *)

The limits are

(lim0 = Limit[if[x], x -> 0, Direction -> "FromAbove"]) // N // Chop

(* -2.35813 *)

(lim = Limit[if[x], x -> Pi/2, Direction -> #] & /@
     {"FromBelow", "FromAbove"} // Quiet) // N

(* {-1.29897, 1.29897} *)

Adding a piecewise constant to the indefinite integral to form a continuous function

if2[x_] = Piecewise[{
     {if[x] - lim0, 0 <= x < Pi/2},
     {lim[[1]] - lim0, x == Pi/2},
     {if[x] - lim[[2]] + lim[[1]] - lim0, Pi/2 < x <= Pi}}] // Simplify;

Plot[{f[x], if2[x]}, {x, 0, Pi},
 Epilog -> {Red, AbsolutePointSize[4],
   Point[{Pi/2, if2[Pi/2]}]},
 WorkingPrecision -> 15,
 PlotLegends -> Placed["Expressions", {0.25, 0.7}]]

Answer 2

This is a known issue in Mathematica (and indeed in most computer-algebra software). Here's a blog post about the issue (from 2008, but still relevant):

Most calculus students might think that if one could compute indefinite integrals, it would always be easy to compute definite ones. After all, they might think, the Fundamental Theorem of Calculus says that one just has to subtract the values of the indefinite integral at the end points to get the definite integral.

So how come inside Mathematica there are thousands of pages of code devoted to working out definite integrals–beyond just subtracting indefinite ones? ...

The article goes through the issues that arise when doing integration, and comes to the following conclusions:

The problem is that there’s actually no way to produce an antiderivative that has this property (i.e., one without discontinuities) for all definite integrals one might want to compute.

Your function has some pretty hairy properties viewed as a function of the complex numbers (poles, branch cuts). For such functions, it's not actually possible to define an antiderivative that doesn't have discontinuities somewhere.

If all you're concerned about is numerical values, you can define a function in terms of an explicit numerical integral over a fixed range:

f[x_] = Sqrt[((1/2)^2 - Csc[x]^2)/(-1 + (1/2)^2 Cos[x]^2)^3] Sin[x]^2
antideriv[x_] := NIntegrate[f[t], {t, 0, x}];
Plot[antideriv[x], {x, 0, \[Pi]}]

However, if you replace the NIntegrate with Integrate above, forcing Mathematica to look for a symbolic solution, it takes much much longer for Mathematica to calculate each step. Calculating antideriv[1] analytically, for example, takes a little over a minute on my computer and yields the following hot mess:

(1/(6 Sqrt[-1 + 194 E^(4 I) - E^(
  8 I)]))E^-I (3 + 45 E^(2 I) + 45 E^(4 I) + 3 E^(6 I) + 
   4 E^I (Sqrt[-3 + 582 E^(4 I) - 3 E^(8 I)] - 
      3 I (Sqrt[-1 + 194 E^(4 I) - E^(
           8 I)] (EllipticE[I ArcSinh[1/Sqrt[3]], -(3/4)] - 
            EllipticF[I ArcSinh[1/Sqrt[3]], -(3/4)]) + 
         Sqrt[-1 + 14 E^(2 I) - E^(4 I)] Sqrt[
          1 + 14 E^(2 I) + E^(
           4 I)] (EllipticE[
             I ArcSinh[1/2 Sqrt[1/3 (2 + E^(-2 I) + E^(2 I))]], -(3/
              4)] - EllipticF[
             I ArcSinh[1/2 Sqrt[1/3 (2 + E^(-2 I) + E^(2 I))]], -(3/
              4)]))))

Moreover, Mathematica is unable to analytically evaluate antideriv[2] at all; on my machine, it churns for about 45 seconds before returning the integral unevaluated.

Answer 3

f = Sqrt[((1/2)^2 - Csc[x]^2)/(-1 + (1/2)^2 Cos[x]^2)^3] Sin[x]^2;

F = Cos[x] Sqrt[49 - Cos[2 x]^2]/(-7 + Cos[2 x]) - 
   2 (Sign[Pi/2 - x] (EllipticE[ArcSin[Sqrt[7 - Cos[2 x]]/Sqrt[14]], 7/4] - 
        EllipticE[ArcSin[Sqrt[7 - Cos[Pi]]/Sqrt[14]], 7/4]));

D[F, x] - f /. Sign' -> (0 &) // 
 FullSimplify[#, 0 < x < Pi && Sin[x] > 0] &
(*  0  *)

Plot[F, {x, 0, Pi}]