1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – b.gates.you.know.what May 23 '18 at 18:38
  • $\begingroup$ @b.gatessucks Why is this happening? $\endgroup$ – user85798 May 23 '18 at 21:43
4
$\begingroup$

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}]

enter image description here

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