9
$\begingroup$

Bug introduced after 8.0.4 and before 9.0.1, and fixed in 11.1.0.


In Mathematica 10.2 I'm trying to integrate this piecewise continuous function, but Integrate and NIntegrate seem to disagree on the branch cut of Sqrt[]:

 Integrate[Sqrt[Exp[I*t]^2 - 1], {t, 0, Pi}]
 (* 2 Sqrt[2] - I Pi *)

 NIntegrate[Sqrt[Exp[I*t]^2 - 1], {t, 0, Pi}]
 (* 1.06568 - 4.51028*10^-16 I *)

It seems that NIntegrate's answer is correct if Sqrt[] has its branch cut on the negative real axis. Is there a way to force Integrate to use the same branch cut?

$\endgroup$
6
  • $\begingroup$ Same result in v9.0.1 and v10.4.1, but 1/2 (4 Sqrt[2] + Log[17 - 12 Sqrt[2]]) in v8.0.4, which is correct. I think it's safe to say it's a bug. $\endgroup$
    – xzczd
    May 27, 2016 at 2:34
  • 1
    $\begingroup$ 5.2 gives consistent results, tho. $\endgroup$ May 27, 2016 at 2:41
  • 1
    $\begingroup$ @xzczd, I mean, 5.2 yields (after simplification) 2 (Sqrt[2] - ArcSinh[1]). So, something went wrong in between 5.2 and 9. $\endgroup$ May 27, 2016 at 2:50
  • 2
    $\begingroup$ @J.M. To be more precise, between 8.0.4 and 9.0.1 :) $\endgroup$
    – xzczd
    May 27, 2016 at 2:53
  • 2
    $\begingroup$ @xzczd, for any other software, it would be weird to keep older versions hanging around your computer... :) $\endgroup$ May 27, 2016 at 2:57

2 Answers 2

6
$\begingroup$

Some insight can be gained by plotting Sqrt[Exp[I*t]^2 - 1] in the complex plane.

Plot3D[Evaluate[ReIm[Sqrt[Exp[I*t]^2 - 1] /. t -> tr + I ti]], {tr, 0, Pi}, 
    {ti, -1, 1}, AxesLabel -> {tr, ti, f}]

enter image description here

Branch points occur at t == n Pi, n an integer, with branch cuts extending from the branch points to t == n Pi + I ∞. Visibly, there also are branch cuts at t an odd half-integer of Pi, even though there are no corresponding branch points. They are, in other words, simply an artifact of Mathematica's definition of Sqrt[]. Nonetheless, they must be taken into account.

The impact of the branch cut at t == Pi/2 is apparent from the following plot of the integrand.

Plot[Evaluate[ReIm[Sqrt[Exp[I*t]^2 - 1]]], {t, 0, Pi}]

enter image description here

The real part of the integrand is symmetric, the imaginary part antisymmetric, which is consistent with the result of NIntegrate in the Question. So, what does Integrate do as it crosses the branch cut?

Integrate[Sqrt[Exp[I*t]^2 - 1], {t, Pi/2, Pi/2 + 10^-10}]
(* Sqrt[2] - I π + I (-Sqrt[-1 - E^(I/5000000000)] + 
   ArcTan[Sqrt[-1 - E^(I/5000000000)]]) + ArcTanh[Sqrt[2]] *)
% // N
(* 1.76275 - 3.14159 I *)

Evidently, integrating an infinitesimal distance produces a finite value. Because the contour cannot be distorted around this branch cut, Integrate must have some built-in rule to add this value as the branch cut is crossed. In any case, Integrate can be made to give the correct value by excluding the region at t == Pi/2.

Integrate[Sqrt[Exp[I*t]^2 - 1], {t, 0, Pi/2 - 10^-10}] + 
    Integrate[Sqrt[Exp[I*t]^2 - 1], {t, Pi/2 + 10^-10, Pi}]
(* I (-Sqrt[-1 - E^(-(I/5000000000))] + ArcTan[Sqrt[-1 - E^-(I/5000000000))]]) + 
   I (Sqrt[-1 - E^(I/5000000000)] - ArcTan[Sqrt[-1 - E^(I/5000000000)]]) *)
% // N
(* 1.06568 + 0. I *)

which is the same value as obtained from NIntegrate, up to roundoff.

Version 11.1.0 update.

Version 11.1.0 returns the correct answer.

Integrate[Sqrt[Exp[I*t]^2 - 1], {t, 0, Pi}]
(* 1/2 (4 Sqrt[2] + Log[17 - 12 Sqrt[2]]) *)
% // N
(* 1.06568 *)
$\endgroup$
8
  • 1
    $\begingroup$ Hmm. Limit[-I (Sqrt[-1 - E^(-2 I ε)] - ArcTan[Sqrt[-1 - E^(-2 I ε)]]) + I (Sqrt[-1 - E^(2 I ε)] - ArcTan[Sqrt[-1 - E^(2 I ε)]]), ε -> 0] yields 2 Sqrt[2] - I Pi in the current version, but (after simplification) 2 (Sqrt[2] - ArcSinh[1]) in 5.2. $\endgroup$ May 27, 2016 at 2:56
  • $\begingroup$ @J.M. Yes, I noticed the first part too. What simplification did you use to obtain second. I can obtain it only be setting ε to a very small number. $\endgroup$
    – bbgodfrey
    May 27, 2016 at 4:05
  • 2
    $\begingroup$ I used FullSimplify[TrigToExp[expr]] in 5.2 . $\endgroup$ May 27, 2016 at 4:10
  • $\begingroup$ @J.M. This does not work for 10.4.1, returning expr. However, Limit[expr // TrigToExp, ϵ-> 0] // FullSimplify reproduces the second result in your first comment. Incidentally, Direction -> 1 flips the sign of both answers. $\endgroup$
    – bbgodfrey
    May 27, 2016 at 4:17
  • 2
    $\begingroup$ @J.M. With regard to the two limits in your first comment, the first almost certainly is another bug. In fact, because this erroneous limit is the same expression as the erroneous value from Integrate, it is possible that Integrate is somehow using this limit in its calculation. $\endgroup$
    – bbgodfrey
    May 27, 2016 at 18:00
1
$\begingroup$

ComplexExpand with TargetFunctions -> Abs help us to solve this integral.

f[t_] := Sqrt[Exp[I*t]^2 - 1];
complex = ComplexExpand[f[t], TargetFunctions -> Abs]

Sqrt[Abs[-1 + E^(2 I t)]] Cos[1/2 (-I Log[-1 + E^(2 I t)] + I Log[Abs[-1 + E^(2 I t)]])] + I Sqrt[Abs[-1 + E^(2 I t)]] Sin[1/2 (-I Log[-1 + E^(2 I t)] + I Log[Abs[-1 + E^(2 I t)]])]

sol = Integrate[complex, {t, 0, Pi}]

2 (Sqrt[2] - ArcSinh[1])

sol//N 

1.06568

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