# Contour integration with Mathematica

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.

• A bit tricky; compare Integrate[Sqrt[16 E^(4 I*t) - 2 E^(I*t)]*2 I*E^(I*t), {t, 0, 2 Pi}] //N to Integrate[Sqrt[16 E^(4 I*t) - 2 E^(I*t)]*2 I*E^(I*t), {t, -Pi, Pi}] //N. May 23, 2018 at 18:38
• @b.gatessucks Why is this happening? May 23, 2018 at 21:43

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   *)