Physics 230 -- Lab 11 (Complex Analysis)

Fall 2015 - Transtrum

Labs developed by John S. Colton, BYU Physics & Astronomy, based in part on past labs by Branton Cambell and Steve Turley

Complex numbers are essential to the mathematical analysis of physical systems. Because the field of complex numbers is “algebraically closed” while the field of real numbers is not, we can't help running into them from time to time.  Loosely speaking, algebraically closed means that polynomial equations with complex coefficients will have complex solutions--whereas there’s no guarantee that polynomial equations with real coefficients will have only real solutions. In addition, complex numbers are frequently used in physics courses to represent real quantities such as amplitude and phase of a sinusoidal wave in order to simplify the mathematical operations. Thus the study of complex number is an important one to physicists.

Clear["`*" ]

Complex Numbers

Basic complex number calculations

As mentioned in a previous lab Mathematica uses a capital I to represent the imaginary number (Sqrt[-1])in standard input. Output cells use a funny looking lower case i:  i. 

In the cell below, we define the complex number, z0 = 2 + i , and show you several useful complex number functions: the real part (Re), imaginary part (Im), absolute value (Abs), phase in radians (Arg), normalized value (Sign, = z/Abs[z] ) and complex conjugate (Conjugate).

z0 = 2+I

Re[z0]
Im[z0]
Abs[z0]
Sign[z0]
Arg[z0]
Conjugate[z0]

The complex plane

Because real and imaginary numbers are in some sense orthogonal to each other, we can plot the real and imaginary parts of a complex number z as coordinates in a plane, the “complex plane”. The x and y axes represent the real and imaginary components of the number, respectively. To make graphing in the complex plane easier, here are two functions which convert a complex number z into the corresponding {x, y} point, and vice versa.

complexToXY[z_] := {Re[z],Im[z]} ;
complexToXY[2+3I]  (* testing *)

XYToComplex[{x_,y_}] :=x + I y
XYToComplex[{2,3}] (* testing *)

All of the functions just mentioned, Re, Im, Abs, Arg, Sign, and Conjugate, have graphical interpretations. 

Consider the following point, (Subscript[z, 0])=3+4 i, plotted on the complex plane. Because points in a plane can be expressed as rectangular (x and y) or polar coordinates (magnitude and direction), complex numbers can be thought of as either points or vectors. I’ve used the following box to plot (Subscript[z, 0]) as both a point and a vector.

z0 = 3+4 I
z0coords = z0 // complexToXY
plot1=ListPlot[{z0coords},PlotStyle->PointSize[0.03]];
plot2 = Graphics[Arrow[{{0,0},z0coords}] ];
Show[plot1,plot2]

Assignment 1. Graphical interpretation of Re, Im, Abs, Arg, and Conjugate

(a) What are Re[z] and Im[z] for this particular point (Subscript[z, 0]), and what do they represent on the above graph?
(b) What are Abs[z] and Arg[z] for this particular point (Subscript[z, 0]), and what do they represent on the above graph? 
(c) Make a similar combination point & arrow plot for Conjugate[z0]. What does the complex conjugate represent, graphically?

Assignment 2. Random complex numbers and interpretation of Sign.

(a) Use RandomComplex to generate 100 random complex numbers within the rectangle defined by corners at -1-i  and 1+i .  Convert them to {x, y} points (the Map function,  /@ , may be helpful) and ListPlot them using AspectRatio -> 1. 
(b) Now apply the Sign function to your list of complex numbers and plot again. Use your plot to figure out the graphical interpretation of the Sign function.

Complex numbers as vectors: addition and subtraction

As mentioned, when viewed as polar coordinates complex numbers are essentially the same as 2D vectors! This means that adding and subtracting complex numbers is just like adding and subtracting vectors.

Addition: Run the next cell a few times. It plots two random complex numbers (Subscript[z, 1]) and (Subscript[z, 2]) (as vectors), along with the complex sum: (Subscript[z, 1])+(Subscript[z, 2]). The sum is in blue. Notice how the sum of the two complex numbers is exactly the same as the vector sum of the two vectors done through the “tip to tail” method that you undoubtedly learned in Physics 121.

z1 =RandomComplex[{-1 -I , 1 + I}]; 
z2 =RandomComplex[{-1 -I , 1 + I}];
z3 = z1+z2;
objects := {Arrow[{{0,0},z1 //complexToXY}], Arrow[{{0,0},z2//complexToXY }], Blue, Arrow[{{0,0},z3 //complexToXY }] };
 Graphics[objects ,Axes->True,PlotRange-> {{-2,2},{-2,2}},AspectRatio->1,AxesLabel->{"x","y"}]

Subtraction: Run the next cell a few times. It plots two random complex numbers (Subscript[z, 1]) and (Subscript[z, 2]) (as vectors), along with the complex difference: (Subscript[z, 1])-(Subscript[z, 2]). The difference is in red. Because order matters I have plotted (Subscript[z, 1]) in solid black and (Subscript[z, 2])in dashed black, so red arrow = solid - dashed.  Notice how the difference of two complex numbers (red arrow) is exactly the same as the vector subtraction of the two black arrows. (The vector subtraction of (Subscript[z, 1])-(Subscript[z, 2]) is the same as adding (Subscript[z, 1]) plus the opposite of (Subscript[z, 2]).)

z1 =RandomComplex[{-1 -I , 1 + I}]; 
z2 =RandomComplex[{-1 -I , 1 + I}];
z4=z1-z2;
objects := {Arrow[{{0,0},z1 //complexToXY}], Red, Arrow[{{0,0},z4 //complexToXY }] , Black,  Dashed, Arrow[{{0,0},z2//complexToXY }]};
Graphics[objects ,Axes->True,PlotRange-> {{-2,2},{-2,2}},AspectRatio->1,AxesLabel->{"x","y"}]

Exponential to Trig, and vice versa

In a previous lab we discussed Euler’s formula:  ((e )^(iϕ))=cosϕ + i sinϕ. To emphasize the geometrical qualities of this equation, consider a complex number  r (e^(iϕ)) where r and ϕ are real numbers. By Euler’s formula it equals  r (cosϕ + i sinϕ) = r cosϕ + i r sinϕ. Viewed in the complex plane, the x coordinate is then  r cosϕ and the y coordinate is  r sinϕ. But those are exactly the x- and y-coordinates of a point (x,y) written in polar form as (r, ϕ)! This is the mathematical justification for representing complex numbers as vectors.

There is a deep connection between complex exponentials and sines/cosines which goes beyond polar/rectangular conversion, and in many physics classes you will need to switch back and forth between the two equivant representations. The two built-in functions,  TrigToExp and  ExpToTrig aid in this process.

Assignment 3. Exponential and sin/cos formulas

(a) Use ExpToTrig to convert (e^ix) into trigonometric form, to verify that Mathematica “knows” Euler’s formula.
(b) Use ExpToTrig to convert (e^x) into trigonometric form. This is the analog of Euler’s formula for hyperbolic functions sinh(x) and cosh(x).
(c) Use TrigToExp to convert cos(x) and sin(x) into exponential form.
(d) Use TrigToExp to convert cosh(x) and sinh(x) into exponential form. Compare to the exponential forms of cosine and sine.
(e) All six of these results are used so frequently in physics classes that they are worth memorizing, in case you don’t already know them. Be prepared to pass them off with a TA.  (Note: the result for sin(x) is usually written with the “i” in the denominator, using 1/i = -i.)

Circular vs. simple harmonic motion

You may have learned in Phys 121 that there is a connection between circular motion and simple harmonic motion. This connection can be visualized using complex numbers. Simple harmonic motion of the form cosωt (or sinωt) is the real part (or imaginary part) of the complex exponential (e^(iωt)). Many mathematical operations such as derivatives and integrals are easier to do with exponentials than with cosines or sines, so (e^(iωt)) is frequently used to analyze harmonic motion mathematically--though it is implicitly understood that only the real part of the expression is physical (in the case of cosωt).

Evaluate the cell below, which animates the arrow pointing towards z = (e^(i θ)) in the complex plane as θ is varied from 0 to 2π. The angle  θ is changing constantly in time, as in: θ = ωt.  The real and imaginary components of z are also animated as large red and blue dots, respectively. The connection between the circular motion of the arrow (i.e. the complex number z) and simple harmonic motion: the red and blue dots move in simple harmonic motion as the arrow rotates in circular motion.

Because of this connection, Arg of the complex number is called either the angle of the complex number (thinking of the vector), or the phase of the complex number (thinking of the simple harmonic motion)--both terms are used equivalently.

euler[th_] =Exp[I  th] ;
eulerobjects[th_] := {Circle[],Thick,Arrowheads[0.1],Arrow[{{0,0},  euler[th]//complexToXY}], PointSize[0.05],Red,Point[{Cos[th],0}], Blue,Point[{0,Sin[th]}]}    (* see below for what this means *)
Animate[ Graphics[eulerobjects[th] ,Axes->True,AspectRatio->1,AxesLabel->{"x","y"}, Background->White],{th,0,2Pi}]

Things to notice: When used inside a Graphics command, the eulerobjects stuff specifies four graphics objects to be plotted. It can be read as: “Plot a circle of default radius, then using a thick line and a larger-than-normal arrowhead plot an arrow centered at the origin, then using a larger-than-normal dot plot a red point, then using the same size dot plot a blue point.”

Multiplying and dividing complex numbers

Euler's formula makes it possible to express any complex number in exponential (polar) form, in terms of a magnitude and a phase angle. The laws of exponential multiplication can be used to obtain short-cuts for multiplying and dividing complex numbers in polar form. If  (Subscript[z, 1])=(Subscript[A, 1])(e^(Subscript[iϕ, 1]))and (Subscript[z, 2])=(Subscript[A, 2])(e^(Subscript[iϕ, 2])) then:

(Subscript[z, 1]) (Subscript[z, 2])=(Subscript[A, 1])(Subscript[A, 2])(e^(Subscript[iϕ, 1]))(e^(Subscript[iϕ, 2]))
          =((Subscript[A, 1])(Subscript[A, 2])) (e^(i((Subscript[ϕ, 1])+(Subscript[ϕ, 2])))) 

In other words, when two complex numbers multiply, the magnitude of the product is the product of the magnitudes, and the angle of the product is the sum of the angles.  This is an incredibly useful rule!

Similarly, when they divide, the magnitude of the result is (Subscript[A, 1])/(Subscript[A, 2]) and the angle of the result is (Subscript[ϕ, 1])-(Subscript[ϕ, 2]).

Complex numbers and rotations

In a previous lab you learned about rotating vectors via rotation matrices. While that is quite useful, if you are dealing with 2D vectors, Euler’s formula gives you another mathematical tool for rotating vectors. Specifically, consider the product (Subscript[z, 2])=(Subscript[z, 1])×R when R is a complex number having magnitude = 1,   R=(e^(iϕ))  with no amplitude multiplying it.

The multiplication rule you just learned (magnitudes multiply, angles add) means that (Subscript[z, 2]) will have the same magnitude as (Subscript[z, 1]), and it will have an angle equal to angle((Subscript[z, 1])) + angle(R). In other words, (Subscript[z, 2]) when viewed as a vector looks just like (Subscript[z, 1]), but rotated by the angle of R! The complex number R is therefore an analog of the rotation matrices studied in a previous lab.

Assignment 4. A “rotate by angle” function

Create a function, rotate[number_, angle_], which “rotates” a complex number by multiplying it by a unit vector having the desired angle. Test out your function by rotating the vector (Subscript[z, 1])=2+i  by 35 degrees and plotting the results (initial vector and rotated vector) with the complexToXY function and commands modified from the “Complex numbers as vectors” section above. Use PlotRange-> {{-2,2},{-2,2}}  and  AspectRatio->1 as options in your Graphics command so that the rotation looks right.

Application: Reflection/Transmission of light

Fresnel coefficients

As a light wave strikes a boundary between two media, some of the wave reflects and some transmits. The equations which govern how much reflects and how much transmits are called the four “Fresnel equations”. You can find the equations in the Wikipedia article at http://en.wikipedia.org/wiki/Fresnel_equations.  (See the “Formulas” given for the “Amplitude equations”, about half-way through the article, here.)  There are two sets of equations, one for “s-polarized” light and one for “p-polarized” light. Each set contains two equations, describing the amplitude ratios for reflection (written as r) and for transmission (written as t). The equations are given in terms of the index of refraction of each material, (Subscript[n, 1]) and (Subscript[n, 2]), and the incident and transmitted angles, (Subscript[θ, i]) and (Subscript[θ, t]). Although you learned about the index of refraction in Physics 123, you likely did not learn that to represent possible phase shifts between the waves, the index of refraction is actually a complex number. The imaginary part of the index is often ignored for insulators, but becomes especially important for conductors as well as for insulators at wavengths of light that are strongly absorbed.

As mentioned, the s- and p-polarizations behave differently from each other. In case you can’t recall (or didn’t learn) the difference between those two polarizations, they are defined relative to the “plane of incidence”. That is the plane in which all three light rays can be drawn (incident, reflected, and transmitted), as in the second figure in the Wikipedia article. P-polarization is when the electric field in the light wave oscillates in this plane; s-polarization is when the electric field oscillates perpendicular to this plane. Other polarizations can always be broken up into their s and p components.

As you should have learned in Physics 123, there is a special incident angle at which there is no reflected beam for p-polarized light. In other words, if you shine unpolarized light at a surface at that special angle, the reflected light will be entirely s-polarized. That angle is called Brewster’s angle. For example, here is the reflected power for a light wave striking a typical air-glass interface (reflected s-polarization plotted in red, p-polarization plotted in blue); you can see Brewster’s angle for this situation occurs around (Subscript[θ, i])= 1 rad.


Although your Physics 123 textbook likely said (as I just did) that there is no p-polarized reflection at Brewster’s angle, it’s actually not quite true. When complex indices of refraction are involved, p-polarized light will partially reflect. It will still have a minimum when plotted vs. incident angle, however, and the angle which gives the minimum reflection is still called Brewster’s angle.

Assignment 5. Fresnel coefficients

(a) According to this website database, for 532 nm light (green), aluminum has an index of refraction of 0.939 + 6.42i. If you are shining a green laser pointer (532 nm, 1 mW of power) at an aluminum mirror, how much laser power will reflect? The answer will depend on the incident angle (Subscript[θ, i]), so use Mathematica to make a plot of the reflected power as a function of (Subscript[θ, i]) (in radians). Note that the reflecting power is described by R, not r; the two are related to each other by R=|r((|)^2). Plot the reflected powers for both s and p polarizations on the same graph as I’ve done in the graph above. Hint: Since the equations are given in terms of (Subscript[θ, i]) and (Subscript[θ, t]), you’ll have to use Snell’s law to relate (Subscript[θ, t]) to the incident angle (Subscript[θ, i]):  (Subscript[n, 1])(Subscript[sinθ, i])= (Subscript[n, 2]) (Subscript[sinθ, t]),  or (Subscript[θ, t])=(sin^(-1))((Subscript[n, 1])(Subscript[sinθ, i])/(Subscript[n, 2])). The complex index will give you a complex (Subscript[θ, t]), but don’t worry about that--the equations still hold for a complex angle of transmission even if a complex angle can’t be visualized. Mathematica can take sines, cosines, and tangents of complex angles just fine.
(b) Find the angle that corresponds to the minimum value of the p-polarization graph, and how much p-polarized light will still reflect at that “Brewster’s angle”.

Complex Algebra

ComplexExpand

At times you may have a complicated algebraic expression that you need to separate into real and imaginary parts. If you had numbers, it would be no problem--you could just use the Re and Im commands.

complicatedz = ((3 + 4 I)/(2 - 5 I) + (1 + 3 I))/(2 + 2 I)
Re[%] //N
Im[%%] //N

But look what happens when you try to use Re and Im with an algebraic expression containing unknowns:

complicatedz= ((a + b I)/(c - d I) + (e + f I))/(g + h I)
Re[%]
Im[%%]

Mathematica doesn’t know what to do, because it doesn’t know whether the numbers a ... h themselves are real or complex.  The ComplexExpand function comes to our rescue. It will algebraically separate z into real and imaginary parts, assuming that all of the unknown parameters are purely real numbers.

ComplexExpand[complicatedz] (* Notice how in the complicated expression, the real stuff comes first and everything that's multiplied by i comes second. *)

Complex Functions

Functions in the complex plane, part 1

Consider the function, f(x)=(e^x), also written f(x)=exp(x). The function is quite easy to visualize. Among other things, the function is never negative:

Plot[Exp[x],{x,-2,2}]

However, what if x can be imaginary, such as a number like πi ? Then what does f(x) equal? You can do it directly in Mathematica, or you can use Euler’s formula.

f(πi) = (e^(πi))
      = cos(π) + i sin(π)
      = -1

Exp[Pi I] //N

So (e^(πi))=-1!  Kind of cool. Hopefully you’ve seen that before. But notice how that doesn’t seem to fit with the graph above, which indicates the function (e^x) is never negative. So a regular graph of a function such as we’re used to doesn’t tell us everything when it comes to complex numbers.

How about x = 3 + 4i?  Again, we can use Mathematica to find the result, and we can understand the result via Euler’s formula:

f(3 + 4i) = (e^(3+4i))
             = (e^3)(e^(4i))
             = (e^3)(cos(4) + i sin(4)
             = 20.0855 (-0.653644 - 0.756802 i)
             = -13.1288 - 15.2008 i

Exp[3 + 4 I] //N
E^3  * (Cos[4] +  I Sin[4] ) //N   (* same thing *)

The moral of the story is: functions can in general take complex arguments, and they can in general spit out complex values of the function. And you can get non-intuitive results, such as a negative value for a function that doesn’t “look” like it ever goes negative.

Assignment 6. Plotting (e^x) as a complex function

What does it mean to plot a complex function? The seemingly simple matter of plotting a function is made more “complex” (ha ha) in two different ways: 

• Since the function can take a complex value as an input, you technically need to have the whole complex plane as the domain of the function. There are two typical ways of dealing with this: (1) do a Plot3D as a function of x and y, where x and y are the real and imaginary parts of the complex number input, or (2) restrict your range to a linear subset of the complex plane, such as only the real numbers, or only the purely imaginary numbers, or only the numbers that have equal real and imaginary components, and do a regular plot with that linear subset as the “x-axis”.

• Since the function can return a complex value as an output, you need to have the ability to plot two independent numbers as the output. The typical way to deal with this is to plot the real and imaginary parts of the output, separately.

Now for the assignment:

(a) Plot the real and imaginary components of f(x) =(e^x) on two separate graphs, using the real numbers from -2 to 2 as the x-axis. (I basically did this above.)
(b) Plot the real and imaginary components of f(x) =(e^x) on two separate graphs, using the purely imaginary numbers from -2i to 2i as the x-axis.  Be careful in how you do this as you cannot use complex numbers as the limits for the plot variable. So, for example, you cannot simply do this: Plot[Re[Exp[x]],{x,-2 I,2 I}]
(c) Plot the real and imaginary components of f(x) =(e^x) on two separate graphs, using the numbers that have equal real and imaginary components from -2-2i to 2+2i as the x-axis. See the “Be careful” comment of part (b).
(d) Plot the real and imaginary components of f(x) =(e^x) on two separate graphs, using Plot3D to incorporate the entire complex plane as the domain, i.e. a square that incorporates the previous three linear subset domains. Show that the results of (a), (b), and (c) are visible in the two plots.
(e) Use the ComplexExpand command to separate the real and imaginary parts of (e^(x+iy)). From the results, you can explicitly see the two separate functions that make up our complex function f(x)=(e^x). How does the output of ComplexExpand relate to your plots?

Functions in the complex plane, part 2

To reiterate: a complex function can be viewed as two separate functions, namely the real part of the result and the imaginary part of the result. The relationship between the two functions is not arbitrary, though--it’s set by the definition of how the function acts on real numbers, and by the properties of i.

Example 1: f=(x^2). Functions don’t get much more basic than this: f(x)=(x^2).  Look what happens when that function acts on a complex number x + i y.

f(x+i y)=(((x+i y))^2)
              = (x^2)+2 x y i -(y^2)
              = ((x^2)-(y^2)) + i (2 x y)

So the real part of the function is (x^2)-(y^2)  and the imaginary part is 2 x y.  And, as must be the case, if y=0 the function turns back into the regular real function, f(x)=(x^2).

We can use ComplexExpand to get that result more quickly:

(x + I y)^2 //ComplexExpand

Example 2: f=sin(x). Here's another example: f(x) = sin(x), applied again to a complex number.

Sin[x+I y]//ComplexExpand

Notice again that if y = 0, then the function returns to the regular real function, sin(x), because cosh(0) = 1 and sinh(0) = 0.

Let’s graph the real part of that:

Plot3D[ Re[Sin[x+I y]] , {x,-2 Pi, 2 Pi},{y,-2 Pi,2 Pi},PlotRange->All,AxesLabel->{"x","y"}]    (* We're plotting f(x,y) = Cosh[y] Sin[x]  *)

Make sure you can see that for a constant y, it’s just a sine function (with an amplitude that depends on what y you picked). And for a constant x, it’s just a cosh function (with an amplitude that depends on the x you picked). If you’ve forgotten what cosh looks like, now would be a good time to make a quick plot. Still looking at this plot along lines of constant x, why are they sometimes cosh functions that open up and sometimes cosh functions that open down?

Complex Series Expansions

Power series expansions

As you know, functions can be expanded in powers of x, called Taylor series. For  example, here’s the Series expansion for f= ln(x) to 10th order, about the point x=2:

fabout2=Series[Log[x],{x,2,10}]  //Normal   (* "//Normal" makes it cut off the "order x^11" part at the end *)
fabout2 /. x-> 2.1  (* the series, evaluated at x=2.1 *)
Log[2.1]  (* to verify *)

fabout2 /. x-> 5.1  (* the series, evaluated at x=5.1. The expansion is quite accurate. *)
Log[5.1]  (* not even close, because 5.1 isn't close to the point about which we expanded (x=2) *)

Things to notice: When we expand about a point, the expansion is valid for points nearby. In this case, x = 2.1 is near x = 2 where the expansion was done, so the Taylor series about x=2 produces a good approximation for the function. But it’s not a good approximation at x = 5.1

Now here’s the same function, expanded in a series about the point  z = 2i. This time the coefficients are obviously complex.

series = Series[Log[z],{z,2 I,10}] //Normal
% /. z -> 2.13 I - 0.09  (* the series, evaluated at a point close to z = 2i *)
Log[2.13 I - 0.09 ]  (* to verify *)

Things to notice: When we expand about a point in the complex plane, the expansion is valid for points nearby. But “nearby” now takes on a geometric meaning--you can go in any direction in the plane as long as you stay close to the point where the expansion was done. In this case, I picked a point close to z = 2 i, but off in a random direction.

Laurent series and “poles”

Let’s consider the function, f(x)=sin(πx)/(x^3).  Here’s the (10^th)-order power series expansion around the point x = 0.

Series[ Sin[π x]/x^3,{x,0,10}] //Normal
Sin[π x]/x^3 /. x->0 //N   (* the function goes infinite at x = 0 *)

What to notice: This is a real function (or at least not an explicitly complex one). Yet the expansion contains a term that has a negative power of x. Apparently not all functions can be expanded in Taylor series such as  (Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])(x^2)+... !  This is something your calculus teacher likely did not dwell on, but it’s true of many functions. As another example, consider something as simple as y = 1/x. You can’t expand that one in a regular Taylor’s series about x = 0, either.

Negative powers of x in series expansions like this often occur with functions of complex numbers. Such series expansions are given a special name: instead of being called “Taylor series”, which implies only positive powers of x, they are called "Laurent series". Negative powers of x will exist in the series expansions whereever the function “blows up”. Those points are called singularities. Singularities do not only exist at x = 0--depending on the function they can potentially exist at any value of x. Focusing now on functions of complex variables, a singularity can potentially exist at any point (Subscript[z, 0]) on the complex plane. And if the function is expanded in a Laurent series about that point, negative powers of z (or more specifically, z-(Subscript[z, 0])) will appear.

A singularity that diverges like a negative power of z (i.e. (z^(-n))) in the series expansion has some special properties and gets another special name--we call it a "pole". Poles only exist at singularities. That means that a pole that exists in the series expansion about a particular point won’t in general exist when we expand the function about a different point. For example, here’s a 4th order series expansion of the same function around the point z = i. See how there is no pole in the series expansion, nor is there a singularity in the function.

Series[Sin[π z]/z^3,{z,I,4}] //Normal  (* doesn't contain any negative powers of (z-i) *)
Sin[π z]/z^3 /. z->I //N   (* the function has a finite value at z = i *)

The pole at z = 0 is clearly visible if we plot the real and imaginary components of the function:

f=Sin[π z]/z^3   
list={Re[f],Im[f]}   /. z-> x+I y
Plot3D[# , {x,-2 , 2 },{y,-2 ,2 },AxesLabel->{"x","y",None},PlotRange->{-30,30}, PlotLabel->#]&  /@ list

Notice how the function goes infinite at the origin, but at no other points.

Assignment 7. Laurent series and residue of 1/(1-(e^(z-0.5 )))

(a)  Plot the real and imaginary parts of this function:  f=1/(1-(e^(z-0.5))). Feel free to copy & paste the code just above as a starting point. Where is the pole, (Subscript[z, 0])?
(b) Compute the (10^th)-order Laurent series expansion around that pole, i.e. the series expansion of f in powers of (z-(Subscript[z, 0])).
(c) Now use Mathematica’s built-in SeriesCoefficient function to determine the (Subscript[a, -1]) and (Subscript[a, 9]) coefficients of the series expansion of f in powers of (z-(Subscript[z, 0])) and verify that they match what you found from part (a). If you just care about one or two of the coefficients, this is a more efficient method than using Series.
(d) The (Subscript[a, -1]) coefficient of a Laurent series expansion has some special significance as we’ll see below, so there’s a special command just for it: Residue. Verify that the Residue command gives the same result as SeriesCoefficient for the (Subscript[a, -1]) term in the series.

Complex Integration

Path integrals

Because the complex plane is very similar to the regular x-y plane, let’s begin by talking about integration of a function of x and y. You are used to thinking of the integral of a 1D function (function of x only) from the point A to the point B as being the “area under the curve” of the function from A to B.  An integral of a function of x and y can be thought of as an area in the exact same fashion. (Note: in a previous lab we considered double integrals of the form ∫∫f(x,y) dx dy, and thought of the result as a volume. Here we are not considering double integrals, but single integrals only.)

The main difference between an integral in the 2D plane from A to B and a regular 1D integral is that when you go from point A = ((Subscript[x, 0]),(Subscript[y, 0])) to point B = ((Subscript[x, f]),(Subscript[y, f])) in 2D space, there’s more than one way to go from A to B. The relevant area that is the result of the integral depends on the path you take. Thus, such integrals are called “path integrals” or “line integrals”. 

Let me help you visualize the situation with some graphs.  Consider a simple function, f(x,y)=(x^2)+(y^2). We want to integrate it from the point A = (0,0) to the point B = (1,1). To represent the function I’ll use Plot3D (the third dimension is f). The start and end points are also plotted here; rotate the graph around so you can see them.

f = x ^2 + y^2;
plot1 = Plot3D[f,{x,-1.5,1.5},{y,-1.5,1.5} ,BoxRatios->{1,1,1}];
plot2= ListPointPlot3D[{{ 0,0,0},{1,1,2}}, PlotStyle->PointSize[.05], BoxRatios->{1,1,1}];
Show[plot1,plot2]

There are an infinite number of ways to travel from A to B. Here are two of them:

(* first path from A to B*)
plot0=ParametricPlot3D[{x,y,f} /.{ x -> t, y-> t},{t,0,1}, BoxRatios->{1,1,1}];
plot00=ParametricPlot3D[{x,y,0} /.{ x -> t, y-> t},{t,0,1}, BoxRatios->{1,1,1}, PlotStyle->Dashed]; (* projection in the x-y plane *)
Show[plot0,plot00, PlotRange->All]

(* second path from A to B *)
plot0=ParametricPlot3D[{x,y,f} /.{ x -> t, y-> t^6},{t,0,1}, BoxRatios->{1,1,1}];
plot00=ParametricPlot3D[{x,y,0} /.{ x -> t, y-> t^6},{t,0,1}, BoxRatios->{1,1,1}, PlotStyle->Dashed]; (* projection in the x-y plane *)
Show[plot0,plot00, PlotRange->All]

If we integrate the function over the first path, we will get the area represented by this figure:

plot3=ParametricPlot3D[{x,y,z}  /.{ x -> t, y-> t},{t,0,1},{z,0,f /.{ x -> t, y-> t}}, BoxRatios->{1,1,1}]

And if we integrate the function over the second path, we will get this area:

plot4=ParametricPlot3D[{x,y,z}  /.{ x -> t, y-> t^6},{t,0,1},{z,0,f /.{ x -> t, y-> t^6}}, BoxRatios->{1,1,1}]

Those two areas are not the same! The value of the integral of f  depends on the path you take from A to B.

Just for kicks, here are the two areas plotted on the same graph as the function and the points A and B. Look at this from many angles until you fully understand everything that’s being shown.

Show[plot1,plot2, plot3, plot4]

To actually calculate the two areas, one has to understand how I’ve defined the paths. I’ve used a parameter “t” to indicate how x and y are related, just as we did when we learned about ParametricPlot back in the Plotting lab.

Path 1:   x = t;  y = t;  t goes from 0 to 1   (resulting in a straight line for the path in the x-y plane)

Path 2:   x = t;  y = (t^6);  t goes from 0 to 1   (resulting in a y=(x^6) equation for the path in the x-y plane)


To integrate, we substitute the appropriate parametrization information into the integrand, then do the integral with respect to t.  However, there’s a critical extra factor that must be included when you do this, related to “the length of an arc” section we discussed back in the Integration lab. You’ll learn more about this in multivariable calculus, so I’ll just tell you the extra factor. It is: (Sqrt[((((dx/dt)))^2)+((((dy/dt)))^2)]). In other words, for a path parametrized by the variable t, 

area under path = (∫_path)f(x,y)ds =∫f(x(t),y(t))(Sqrt[((((dx/dt)))^2)+((((dy/dt)))^2)])dt 

(I’ve used “ds” on the left hand side just to indicate that it’s a 1D integral along a certain path. The s doesn’t mean anything.)

There doesn’t seem to be a built-in function in Mathematica to calculate a parametrized path integral like this, so I wrote a program to do so:

pathintegral[f_,xoft_,yoft_,t_,tstart_,tend_] :=  Module[{extrafactor,foft},   
	(*  f = f(x,y);  
		xoft = x as func of t for the particular path, 
		yoft = y as func of t for the particular path,  
		t is just the letter t, 
		tstart and tend are limits of t  *)
extrafactor=Sqrt[D[xoft,t]^2 + D[yoft,t]^2];
foft = f /. {x->xoft, y-> yoft} ;
NIntegrate[foft * extrafactor,{t,tstart,tend}]
]

Here are the areas defined by the two paths:

pathintegral[x^2 + y^2, t,t,t,0,1] //N  (* path 1 *)
pathintegral[x^2 + y^2, t,t^6,t,0,1] //N   (* path 2 *)

Assignment 8. A much more interesting path

Still using  f(x,y)=(x^2)+(y^2), consider yet another path going from (0,0,0) to (1,1,2): 

Path 3:  x = sqrt(t);  y = sin(t 5π/2);  t goes from 0 to 1 

(a) Make a plot of the path represented by that parametrization
(b) Make a plot of the area represented by the the integral of f  along this path. To verify that this plot is correct, also Show it together with the above plot1. It should intersect the function in plot 1 everywhere along the path.
(c) Use my “pathintegral” function to calculate the area represented by your plot.

Contour integrals

Now back to complex numbers. A path integral involving the complex plane is typically called a “contour integral”. It often has a graphical representation, and like regular path integrals of functions of x and y, the value of the integral depends strongly on the path or “contour” that one takes. As the linked-to Wikipedia article describes, the typical way to compute such integrals is to use a parametrization exactly as we did for functions of x and y in the last section.

Since the function we studied in the last section, f=(x^2)+(y^2) closely corresponds to the complex number function f=|z((|)^2) (think about it), let’s look at what happens when we integrate |z((|)^2)  from the point   0 + 0i  to  1 + 1i.

Note: Unlike the plot command, the Integrate command can take complex values for the limits. Also, when you specify limits like this it seems that Mathematica assumes a straight line from A to B--like my “path 1” above.

Integrate[Abs[z]^2, {z,0,1+I}]

That answer is closely related to the path integral we obtained for path 1 above, as can be seen by taking the absolute value.

Abs[%] //N

They match exactly! 

Warning: Don’t take the area analogy too literally. Since the contour integral of the complex function yields a complex number, it doesn’t exactly correspond to an area... and even though in this case the area can be obtained by simply taking the magnitude of the complex number result, that’s not true in general. So, the correspondence between complex number integrals and areas under the curve of x-y path integrals is not 100%. In fact, as will be seen in the next section there are many cases where a complex number function will integrate to zero, even when  there is clearly some area present.

Residue theorem

A special type of contour integral is one where the start and end points are the same: a loop in the complex plane. An important theorem of complex analysis states that when we integrate a function along a closed right-handed loop, the result is equal to 2πi  times the sum of the Residues of any poles enclosed by the loop. Wow! 

That “residue theorem” is certainly not intuitive, but I’ll leave it to others to prove for you. This theorem turns out to be a powerful tool for solving really hard definite integrals. (It’s also the source of this really bad mathematical joke: “Q: What’s the value of a contour integral around Western Europe? A: Zero, because all the Poles are in Eastern Europe.”)

The purpose of this final assignment is to demonstrate the residue theorem for a particular function and help you visualize it, by comparing the residues of the function to the results of circular contour integrals around the poles.

Here’s the function: g(z) = 3/(z+(z^3)).

It has three poles in the complex plane: z = -i, z = 0, z = i .  I.e., if you plug in any of those three numbers, g(z) becomes infinite. This can be visualized with a 3D plot (top view).

g[z_] = 3/(z+z^3)
gplot=Plot3D[Abs[g[x+I y]],{x,-3,3},{y,-3,3},AxesLabel->{"x","y",None},Axes->{True,True,False},MeshStyle->Gray,PlotPoints->20,PlotRange->{0,10},LabelStyle->Medium,ViewPoint->{0,0,20},SphericalRegion->True]

Assignment 9. Residue theorem demonstration.

(a) Determine the residues for the function g at each of the three poles. Recall the residue is the coefficient of the (((z-(Subscript[z, 0])))^(-1)) term of the series expansion. You can do this via a Series expansion, via the SeriesCoefficient function, or via the Residue command. Or all three, if you wish. Record (on a piece of paper) the residue of each pole.

(b) Here’s some code which defines and displays a parametric path for the integral, then integrates it. You can specify the details of the circular path by adjusting z0 and r in the first line. (It should currently be set to integrate around a circle containing just the middle pole.) Perform three integrals, one around each pole, and verify that the result of each integral is 2πi  times the residue at that pole.

z0 = 0;  r= 1/2;  (* z0 is the expansion point, r is the integration-path radius *)
parametricpath[th_] := {r*Cos[th]+Re[z0],r*Sin[th]+Im[z0], Abs[g[r*(Cos[th]+I Sin[th])+z0]]}
Show[ gplot,
Plot3D[Abs[g[x+I y]],{x,-3,3},{y,-3,3},AxesLabel->{"x","y",None},Axes->{True,True,False},MeshStyle->Gray,PlotPoints->20,PlotRange->{0,10},LabelStyle->Medium,ViewPoint->{0,0,20},SphericalRegion->True],
ParametricPlot3D[parametricpath[th],{th,0,2Pi},PlotStyle->{Thickness[0.01],Black}]
]
zzz= z0 +  r*Exp[I θ];
integrand = g[zzz] D[zzz,θ];
Integrate[integrand ,{θ,0,2Pi}]

(c) Test out the residue theorem by defining your path to include no poles. The function should integrate to zero.

(d)  Finally, change (Subscript[z, 0]) back to 0 and make r large enough to include all three poles.  Verify that this larger path integral yields the 2πi  times the sum of all three residues.

(* part b, about the pole at z = 0. (Also done above.) Result is 2pi * residue at 0 *)

z0 = 0;  r= 1/2;  (* z0 is the expansion point, r is the integration-path radius *)
parametricpath[th_] := {r*Cos[th]+Re[z0],r*Sin[th]+Im[z0], Abs[g[r*(Cos[th]+I Sin[th])+z0]]}
Show[ gplot,
Plot3D[Abs[g[x+I y]],{x,-3,3},{y,-3,3},AxesLabel->{"x","y",None},Axes->{True,True,False},MeshStyle->Gray,PlotPoints->20,PlotRange->{0,10},LabelStyle->Medium,ViewPoint->{0,0,20},SphericalRegion->True],
ParametricPlot3D[parametricpath[th],{th,0,2Pi},PlotStyle->{Thickness[0.01],Black}]
]
zzz= z0 +  r*Exp[I θ];
integrand = g[zzz] D[zzz,θ];
Integrate[integrand ,{θ,0,2Pi}]

Wrap-up

The following functions were discussed in today's lab, or are very closely related:

Basic complex number functions: real part (Re), imaginary part (Im), absolute value (Abs), phase in radians (Arg), normalized value (Sign) and complex conjugate (Conjugate)
Random complex numbers: RandomComplex 
Exp/trig conversion: TrigToExp and  ExpToTrig 
Separate into real & imaginary parts: ComplexExpand 
Laurent series expansion: Series (same as Taylor series), SeriesCoefficient 
Residue of a complex function at a point: Residue 


Potentially useful summary/tutorial pages:
(I don’t think there were any mentioned today.)