Contour integration with Mathematica

Question

I seem to have difficulty getting Mathematica to evaluate contour integrals correctly. I am just parametrising them and inputting them in as usual but I am having problems.

For example, consider $I=\oint_C \sqrt{z^4-z} \,\,\mathrm d z$ where $C$ is the circle of radius $2$ (anticlockwise). I believe the correct answer is $-i \pi$ (or $+i \pi$ depending on which square root is chosen), but Mathematica returns:

N@Integrate[Sqrt[16 E^(4 I*t) - 2 E^(I*t)]*2 I*E^(I*t), {t, 0, 2 Pi}]
6.79746 - 6.40204 I

NIntegrate[Sqrt[16 E^(4 I*t) - 2 E^(I*t)]*2 I*E^(I*t), {t, 0, 2 Pi}]
7.60503*10^-15 + 1.37562 I

and when I try the closed-form 'Integrate' I just get something ridiculous that will never simplify.

I would be most grateful if somebody could show me how to get Mathematica to give me correct results for such integrals as I currently have no reliable way of checking if my work is correct.

Answer 1

Issues with brach cut have to be taken into account.

f[z_] = Sqrt[z^4 - z]

integrand[t_] = 
   f[2 E^(I t)]*2 I*E^(I*t) // 
        ComplexExpand[#, TargetFunctions -> {Re, Im}] & // 
        FullSimplify[#, 0 < t < 2 Pi] &

(*   2 I Sqrt[2] E^(I t) Sqrt[E^(I t) (-1 + 8 E^(3 I t))]   

or the same

   2 I Sqrt[2] (65 - 16 Cos[3 t])^(1/4) (Cos[
   t + 1/2 ArcTan[-Cos[t] + 8 Cos[4 t], -Sin[t] + 8 Sin[4 t]]] + 
   I Sin[t + 1/2 ArcTan[-Cos[t] + 8 Cos[4 t], -Sin[t] + 8 Sin[4 t]]])   *)

There are branch cuts where y of ArcTan[x,y] is zero and x less zero.

{sol = Solve[-Sin[t] + 8 Sin[4 t] == 0 && -Cos[t] + 8 Cos[4 t] < 0 && 
    0 < t < 2 Pi, t], t /. sol // N}

(*   {{{t -> 2 \[Pi] - 
 2 ArcTan[Sqrt[
   Root[-31 + 227 #1 - 221 #1^2 + 33 #1^3 &, 1]]]}, {t -> 
2 ArcTan[Sqrt[
  Root[-31 + 227 #1 - 221 #1^2 + 33 #1^3 &, 1]]]}, {t -> 
2 \[Pi] - 
 2 ArcTan[Sqrt[
   Root[-31 + 227 #1 - 221 #1^2 + 33 #1^3 &, 3]]]}, {t -> 
2 ArcTan[Sqrt[
  Root[-31 + 227 #1 - 221 #1^2 + 33 #1^3 &, 3]]]}}, {5.51943, 
0.763757, 3.94961, 2.33357}}   *)

Therefore construct a piecewise integrand rotating complex numbers by an angle Pi

integrand2[t_] = 
  Piecewise[{{integrand[t], 
     0 <= t <= (t /. sol[[2]])}, {-integrand[t], (t /. sol[[2]]) < 
     t <= (t /. sol[[4]])}, {integrand[t], (t /. sol[[4]]) < 
     t <= (t /. sol[[3]])}, {-integrand[t], (t /. sol[[3]]) < 
     t <= (t /. sol[[1]])}, {integrand[t], (t /. sol[[1]]) < t <= 
     2 Pi}}, 0]

Now NIntegrate yields the right result at once

NIntegrate[integrand2[t], {t, 0, 2 Pi}]

(*   1.77636*10^-15 - 3.14159 I   *)

Do anylytical integration by hand with the indefinite integral. I therefore used the RuleBasedIntegrtor Rubi of Albert Rich, http://www.apmaths.uwo.ca/~arich/ because Integrate of Mathematica produced more jumps.

l1 = Int[integrand2[t][[1, All, 1]], t]

(*    {2/3 Sqrt[2] E^(I t) Sqrt[-E^(I t) + 8 E^(4 I t)] - 
       1/3 ArcTanh[(2 Sqrt[2] E^(2 I t))/Sqrt[-E^(I t) + 8 E^(4 I t)]],

      -(2/3) Sqrt[2] E^(I t) Sqrt[-E^(I t) + 8 E^(4 I t)] + 
       1/3 ArcTanh[(2 Sqrt[2] E^(2 I t))/Sqrt[-E^(I t) + 8 E^(4 I t)]],

      2/3 Sqrt[2] E^(I t) Sqrt[-E^(I t) + 8 E^(4 I t)] - 
      1/3 ArcTanh[(2 Sqrt[2] E^(2 I t))/Sqrt[-E^(I t) + 8 E^(4 I t)]],

      -(2/3) Sqrt[2] E^(I t) Sqrt[-E^(I t) + 8 E^(4 I t)] + 
       1/3 ArcTanh[(2 Sqrt[2] E^(2 I t))/Sqrt[-E^(I t) + 8 E^(4 I t)]], 

      2/3 Sqrt[2] E^(I t) Sqrt[-E^(I t) + 8 E^(4 I t)] - 
      1/3 ArcTanh[(2 Sqrt[2] E^(2 I t))/Sqrt[-E^(I t) + 8 E^(4 I t)]]}    *)

l2 = integrand2[t][[1, All, 2]];

Since the Int-command of Rubi is not compatible with Piecewise, construct the integrated function.

tp = Transpose[{l1, l2}];

rint[t_] = Piecewise[tp, 0]

Plot shows still two jumps of imaginary part at 2/3 Pi and 4/3 Pi, therefore limits have to be taken.

Plot[{rint[t] // Re, rint[t] // Im}, {t, 0, 2 Pi}, PlotPoints -> 200, 
 PlotStyle -> {Blue, Red}, GridLines -> {{2/3 Pi, 4/3 Pi}, Automatic},
 Ticks -> {{2/3 Pi, 4/3 Pi}, Automatic}]

Limit gives the result

lim = Plus @@ {Limit[rint[t], t -> 2 Pi, 
 Direction -> 1], -Limit[rint[t], t -> 4/3 Pi, Direction -> -1], 
Limit[rint[t], t -> 4/3 Pi, 
 Direction -> 1], -Limit[rint[t], t -> 2/3 Pi, Direction -> -1], 
Limit[rint[t], t -> 2/3 Pi, 
 Direction -> 1], -Limit[rint[t], t -> 0, Direction -> -1]} // 
 FullSimplify

(*  - I Pi   *)