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$
  • $\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 '16 at 2:34
  • 1
    $\begingroup$ 5.2 gives consistent results, tho. $\endgroup$ – J. M. will be back soon May 27 '16 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$ – J. M. will be back soon May 27 '16 at 2:50
  • 2
    $\begingroup$ @J.M. To be more precise, between 8.0.4 and 9.0.1 :) $\endgroup$ – xzczd May 27 '16 at 2:53
  • 2
    $\begingroup$ @xzczd, for any other software, it would be weird to keep older versions hanging around your computer... :) $\endgroup$ – J. M. will be back soon May 27 '16 at 2:57
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$
  • 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$ – J. M. will be back soon May 27 '16 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 '16 at 4:05
  • 2
    $\begingroup$ I used FullSimplify[TrigToExp[expr]] in 5.2 . $\endgroup$ – J. M. will be back soon May 27 '16 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 '16 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 '16 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.