Chapter 4: Complex numbers and functions
> z1:=5 + 3*I;
>
You need to be careful with this syntax because Maple will allow you to use i, like this:
> z2:=5+3*i;
>
This looks fine, but when you do arithmetic with "i" it won't work right because i is just another variable that doesn't have a value yet. Only I is the real thing. (Actually, you can redefine I to be i if you want by using the interface command, but I suggest you don't.)
In this chapter we will review the various things you can do in the complex plane, including arithmetic, calculus, and even special functions like
and
. Note: in this chapter whenever you encounter a variable with a name like z, or z1, or z2, etc., it is a complex variable of the form
.
> z1:=5+3*I;
> z2:=1-6*I;
> z3:=z2+z1;
Just for fun, here's how this addition problem looks using arrows in the complex plane. (Note: I am showing you some of Maple's tools for drawing pictures in the commands below, so there is more useful information here than just the discussion of complex arithmetic). In the pictures z1 is green, z2 is blue, and z3 is red.
>
with(plottools):
l1 := arrow([0,0], [5,3], .2, .5, .1, color=green):
l2 := arrow([0,0], [1,-6], .2, .5, .1, color=blue):
l3 := arrow([0,0], [6,-3],.2,.5,.1,color=red):
>
plots[display](l1,l2,l3, axes=normal,view=[-10..10,-10..10]);
(5+3I)*(1-6I) = (5 + 3I -30I + 18) = 23-27I . You just use the distributive law and work it out. And Maple knows how to do it:
> z1*z2;
> conjugate(6+4*I);
The conjugate is useful because of what the product of a complex number and its conjugate turn out to be:
> conjugate(6+4*I)*(6+4*I);
>
Their product is just a real number, equal to the square of the magnitude of the complex number, treated as a vector in the usual Pythagorean way: square the components and add. The answer is
. This is useful for doing division because if you multiply the top and bottom of a complex fraction by the conjugate of the bottom, the bottom becomes real and you don't have to worry about I in the denominator anymore. Just multiply the numerator out using the distributive law and divide both real and imaginary parts by the real number in the bottom.
Problem 4.1
Work out by hand the complex division problem
, then do it with Maple and show that the answers agree.
---------------------------------------------------------------------------------------------
>
with(plottools):
l1 := arrow([0,0], [5,4], .2, .5, .1, color=green):
l2 := arc([0,0],4,0..arctan(.8),color=red):
>
plots[display](l1,l2, axes=normal,view=[-5..5,-5..5],scaling=constrained);
>
>
with(plottools):
l1 := arrow([0,0], [1,2], .1, .3, .1, color=green):
l1a := arc([0,0],1.5,0..arctan(2),color=green):
>
l2 := arrow([0,0], [1,.8], .1, .3, .1, color=green):
l2a := arc([0,0],.75,0..arctan(.8),color=green):
>
l3 := arrow([0,0], [-.6,2.8], .1, .3, .1, color=red):
l3a := arc([0,0],2.5,0..arctan(2.8,-.6),color=red):
>
plots[display](l1,l2,l3,l1a,l2a,l3a, axes=normal,view=[-3..3,-3..3],scaling=constrained);
>
Division works similarly, but with the two magnitudes being divided and the two angles subtracting:
.
> z:=5-3*I;
> Re(z);Im(z);abs(z);arctan(Im(z),Re(z));
Maple also has a command that converts z in Cartesian form into polar form, but I find it hard to use its results in subsequent calculations, so I usually just use the arctan form given above. But here it is anyway.
> readlib(polar):
> polar(z);
And finally, remember that Maple will take the complex conjugate for you.
> conjugate(3+4*I);
Problem 4.2
Do the operations indicated below on the complex numbers
and
and express your answers in both Cartesian (x+iy) form and polar form (
readlib(polar):
and the
polar
command make this problem easy).
Here are the definitions of z1 and z2 that you are supposed to use in this problem:
> z1:=3-5*I;z2:=-2+3*I;
>
(a)
(b)
(c)
(d)
----------------------------------------------------------------------------------------------
Problem 4.3
Show graphically, by making a Maple arrow plot like the one shown above, that rotating a complex number through angle
in the complex plane is easily accomplished by multiplying by the factor
. Do this by defining a complex number z1 and plotting it as an arrow, then multiplying z1 by
for some angle
that you choose, and then by plotting this new complex number as an arrow. As you do this you need to be careful to only use floating point numbers in the arguments of the
arrow
command. Maple will probably give you answers back that have square roots of integers in them and if you just copy these numbers into the arrow command no arrow will appear. Our old friend
evalf
will help you here, as will this syntax for putting the real and imaginary parts of a complex number into the arrow command:
[Re(z),Im(z)]
.
All of the functions you have ever heard of make sense for complex numbers as well as for real numbers. And in fact, the complex plane can give a more unified picture of these functions than is possible when we are confined only to the real axis.
We will begin with the exponential function. Euler's formula
and the properties of the exponential function allow us to write
. You may not be able to visualize what this function looks like just by staring at this formula, so to see the behavior of this function over the entire complex plane, let's make a 3-D surface plot of its real part and its imaginary part using Maple's
plot3d
command (see
3d plot).
> restart;
> plot3d(Re(exp(x+I*y)),x=-6..6,y=-6..6,shading=z,axes=framed,labels=["x","y","Re(f)"],orientation=[-60,60],title="Real Part");
> plot3d(Im(exp(x+I*y)),x=-6..6,y=-6..6,shading=z,axes=framed,labels=["x","y","Im(f)"],orientation=[-60,60],title="Imaginary Part");
>
Notice that at negative x it is very small and for positive x it is very big--this is the behavior of
. In y, the functions just oscillate like cosines and sines, which is, of course, exactly what they are.
This mixing together of the exponential function and sines and cosines gives rise to several interesting formulas which you are to verify in the next problem.
Problem 4.4
(a)
(b)
(c)
(d)
-----------------------------------------------------------------------------------------------------
Problem 4.5
Make 3-D surface plots over the complex plane of the real parts of
and
, and of the imaginary part of
, where
. Use x in the range -9..9 and y in the range -5..5 for sin and cos; use -5..5 and -3..3 for tan. Stare at them until they make sense. When you do the tangent function you will have to deal with its singularities. Do this by using the view option in the plot, in the form
view=-5..5
. This vertical size control keeps Maple from trying to plot infinity. It may also help to have a finer grid, which you can get with
grid=[30,30].
To make sense of it you will also have to remember what the tanh function looks like.
-----------------------------------------------------------------------------------------------------
Problem 4.6
Make a 3-D plot of the real parts of the functions
and
using the range -3..3 in both x and y for
and a smaller range for
.
----------------------------------------------------------------------------------------------------
For non-integer powers things get a lot more complicated. To see why, let's consider the square root function,
. Below I have coded three 3-D plots of this function. The first is of its magnitude, and it just looks like the normal square root function rotated about the z-axis, so this looks like we expected. But there is something funny going on along the negative x-axis in the real part plot, and it gets worse in the imaginary part plot.
> plot3d(abs((x+I*y)^(1/2)),x=-1..1,y=-1..1,shading=z,axes=framed,labels=["x","y","|f|"],orientation=[-60,45],title="Magnitude",grid=[40,40]);
> plot3d(Re((x+I*y)^(1/2)),x=-1..1,y=-1..1,shading=z,axes=framed,labels=["x","y","Re(f)"],orientation=[-60,35],title="Real Part",grid=[40,40]);
> plot3d(Im((x+I*y)^(1/2)),x=-1..1,y=-1..1,shading=z,axes=framed,labels=["x","y","Im(f)"],orientation=[-60,35],title="Imaginary Part",grid=[40,40]);
>
Problem 4.7
Do the same thing to the function ln(z) that we just did to square root. Locate the branch cut and explain why it's there by looking at the same problem discussed in the paragraph above. Note: be sure to use the function
in your plot commands.
------------------------------------------------------------------------------------------------------
Problem 4.8
Finally, do a set of 3-D surface plots of this rather ordinary looking and perfectly well-behaved-along-the-real-axis function
(remember that
) . First make a plot of f(x) along the real axis so you know what it looks like there. Use the x-range -4..4. Then make surface plots of the real part, the imaginary part, and the magnitude (use
abs
), again using -4..4 in both dimensions. Also use the option view=-3..3. Hopefully this will give you a sense of the surprise and beauty of the complex plane.
It will probably come as no surprise that the more advanced functions of mathematical physics do the same kinds of things in the complex plane that the elementary ones do. Rather than talking about it I'll let you discover a few of the interesting things that happen through these exercises.
Problem 4.9
-----------------------------------------------------------------------------------------------------
Problem 4.10
Do a 3-D plotting exploration of the behavior of the elliptic integral K(z) in the complex plane. This function is obtained with the Maple command EllipticK(z). Look at its real part and imaginary part with both x and y in the range -2..2.
> restart;
> s1:=Sum(cos(omega*t+2*Pi*'n'/13),'n'=1..12);
> s1:=value(s1);
> simplify(s1);
>
Well, that's an answer, but it's not very pretty, and I can't get Maple to simplify it. Now let's try doing it using the complex plane. Instead of adding 12 cosines together we'll add 12 complex exponentials together. Since the real part of each one is a cosine we can do the complex sum, then take the real part at the end and it will be as if we had added 12 cosines together. The only difference is that the answer will be a lot prettier. Ok, let's see how this goes. Rather than giving the problem to Maple right away, let's do a little pre-analysis first. Here's the sum whose real part we will take at the end
.
We will take out the common factor of
( notice: time factored out! this is the whole reason for using the complex exponential)
and then rewrite the remaining complex exponential this way
. So now our sum looks like this
where
. Now we remember from algebra class
> s:=Sum(z^n,n=1..N);
> s:=simplify(value(s));
> N:=12;s;
so that in our case we have
. But now look --
, but since
can be rewritten as
! Now Maple knows what to do
> simplify(z*(1-1/z)/(1-z));
>
Now we tack on the
we factored out front, take the real part, and discover that the big mess of cosine functions in the original Maple sum is just:
!
The lesson here is hopefully clear; if you are doing an interference problem where you have to add up a whole bunch of sines or cosines, do it with complex exponentials instead.
Problem 4.11
> restart:with(plots):
> l1:= listplot([[0,0],[2,0]]): # distance to screen
> t1:= textplot([1,-.1,"L"]): # L label
> l2 := listplot([[-.08,-.2],[-.1,-.2],[-.1,.2],[-.08,.2]]): # d bracket
> t2:= textplot([-.15,0,"d"]): # d label
> p1:= pointplot([0,-.2],symbol=circle): # lower source
> p2:= pointplot([0,.2],symbol=circle): # upper source
> l3:= listplot([[2,-1.25],[2,1.25]]): # screen
> t3:= textplot([2.08,.5,"y"]): # y label
> l4:= listplot([[0,.2],[2,.6]]): # upper ray
> t4:= textplot([1.5,.6,"r1"]): # r1 label
> l5:= listplot([[0,-.2],[2,.6]]): # lower ray
> t5:= textplot([1.5,.3,"r2"]): # r2 label
> plots[display]([l1,t1,l2,t2,p1,p2,l3,t3,l4,t4,l5,t5],axes=none,thickness=3,view=[-.5..2.1,-1.25..1.25],font=[HELVETICA,14]);
Warning, the name changecoords has been redefined
>
To study the interference of these two sources with as much accuracy as possible we will let the field from each one fall off like
from the source. The two distances are given by
and
so when we add the two fields together to get the total field at the screen we get the kind of mess Maple was invented to handle:
Define functions for r1(y) and r2(y)
> R1:=y-> sqrt((y-d/2)^2+L^2);
> R2:=y-> sqrt((y+d/2)^2+L^2);
The straightforward thing to do now would just be to add the wavefunctions from the two sources, take the real part, and square it to get the intensity, like this:
> Intensity:= y-> Re( exp(I*k*R1(y)-I*omega*t)/R1(y) + exp(I*k*R2(y)-I*omega*t)/R2(y))^2;
>
Unfortunately, it's not that simple. The problem is with the time. If we just choose some convenient time, like t=0, then we will be comparing instantaneous intensities at different places and at each place the intensity varies between zero and some maximum. We don't want to see fake zeros in our intensity pattern on the screen just because we look at the wrong time. To get rid of this problem we want to find the time-averaged intensity .
The time average of an oscillating quantity is defined this way: <f(t)> =
, where T is the period of the oscillation. This looks a little complicated, but is really just a dressed-up version of what you normally think of as the average: (i) add up all the data points (integration is addition of a whole bunch of f*dt's) and (ii) divide by the number of of points (divide by the integration interval). The connection works like this
but
, the approximate number of data points, so we have
, the usual average.
Since we usually work with sines and cosines it is useful to know what the time averages of sines, cosines, and their products are:
> Int(sin(t),t=0..2*Pi)/(2*Pi);value(%);
>
> Int(cos(t),t=0..2*Pi)/(2*Pi);value(%);
>
> Int(cos(t)*sin(t),t=0..2*Pi)/(2*Pi);value(%);
>
> Int(cos(t)^2,t=0..2*Pi)/(2*Pi);value(%);
>
> Int(sin(t)^2,t=0..2*Pi)/(2*Pi);value(%);
>
You need to memorize these results: the time averages of
,
,
and
are all zero. The time averages of
and
are both 1/2.
They even make sense: sin, cos, and sin*cos all spend equal amounts of time being positive and negative, so they average to zero. Their squares oscillate symmetrically between 0 and 1, so their averages are 1/2.
OK, with this straight we can tackle the intensity again. Here's what Maple gives us for the intensity in expression form:
> assume(k,real,r1,real,r2,real,omega,real,t,real);
> Intense:= Re( exp(I*k*r1-I*omega*t)/r1 + exp(I*k*r2-I*omega*t)/r2)^2;
To make the calculation simpler, let's redefine the arguments of these cosine functions to be x1 and x2, respectively:
> Intense:=(cos(x1)/r1+cos(x2)/r2)^2;
Now expand it like this
> Intense:=expand(Intense);
and then time average each term.
>
The two squared terms are easy: both x1 and x2 have
in them so they just time average to 1/2 (the phase shift due to
only shifts the function; it still spends equal amounts of time between 0 and 1.) The mixed term we have to handle with a trig identity:
.
The x1+x2 term goes like
and time averages to zero, while the x1-x2 term has no time dependence at all, so its time average is just itself. So we get
> Intense:= 1/(2*r1^2)+cos(x1-x2)/(r1*r2)+1/(2*r2^2);
With time averaging complete we can redefine our intensity function of position y
> Intensity:=y->1/(2*R1(y)^2)+cos(k*R1(y)-k*R2(y))/(R1(y)*R2(y))+1/(2*R2(y)^2);
Now define the details of physical situation with d=.01mm, L=1m, and red light of wavelength 600 nm
> d:=1e-5;L:=1;k:=2*Pi/600e-9;
Finally, we can plot the intensity as a function of y to see what the pattern should look like on the screen. We will just look in a window 1 m wide near the center of the pattern and see a beautiful interference pattern.
> plot(Intensity(y),y=-0.5..0.5);
>
---------------------------------------------------------------------------------------------------
Problem 4.12
Serway says that the first minimum away from y=0 is at the place where
, where
is the angle indicated below (
). Verify that the first interference minimum in the plot above is at the correct place by measuring y-values on the plot in 4.11 and converting y to angle
. Use the same values of the variables that were used in Problem/Example 4.11:
> restart:with(plots):
Warning, the name changecoords has been redefined
> l1:= listplot([[0,0],[2,0]]): # distance to screen
> t1:= textplot([1,-.1,"L"]): # L label
> l3:= listplot([[2,-1.25],[2,1.25]]): # screen
> t3:= textplot([2.08,.5,"y"]): # y label
> l4:= listplot([[0,0],[2,.6]]): # hypotenuse
> t4:= textplot([.6,.1,"q"]): # theta label
> p1:=plots[display]([l1,t1,l3,t3,l4],axes=none,thickness=3,view=[-.5..2.1,-1.25..1.25],font=[HELVETICA,14]):
> p2:=plots[display]([t4],axes=none,thickness=3,view=[-.5..2.1,-1.25..1.25],font=[SYMBOL,20]):
> display(p1,p2);
>
------------------------------------------------------------------------------------------------
Problem 4.13
Something that Serway doesn't do is show what happens when the screen gets closer and closer to the source. This is complicated because the two fields no longer have about the same magnitude and the distance functions
and
vary more rapidly with y. Rerun the calculation in Problem 4.11 above twice more, once with L=10d and once with L=2d and comment on what these near-field patterns look like. In particular, does the formula given in Problem 4.12 still work? (Be sure to reduce the window in y over which you plot to an appropriate size; -20d..20d is interesting.) There is a lesson here: things are really complicated when you are close to radiation sources.
Electrostatics with line charges
The complex plane is a wonderful place to visualize electric fields and equipotential surfaces. Since the complex plane is only two-dimensional, the only kinds of fields you can visualize there are the fields of infinitely long line charges. It turns out (you will just have to trust me until you take a course in complex analysis) that the complex function
has the property that its real part is the electrostatic potential of a line charge with charge/length
located in the xy-plane at complex position
, and that the lines in the complex plane where the imaginary part is constant are the electric field lines! This function combined with Maple's plotting ability makes it possible to quickly visualize equipotential surfaces and field lines for a variety of line charge arrangements. And in the process we will review how Maple does contour plots (
contour plot). Here is a
contourplot
that shows equipotential surfaces in blue and electric field lines in red for a single line charge. (Since we are just visualizing things I have left off all of the physical constants out front.)
> restart:with(plots):
> p1:=contourplot(Re(-ln(x+I*y)),x=-5..5,y=-5..5,grid=[30,30],contours=30,coloring=[cyan,blue]):
> p2:=contourplot(Im(-ln(x+I*y)),x=-5..5,y=-5..5,grid=[30,30],contours=30,coloring=[red,red]):
> plots[display]([p1,p2],scaling=constrained);
Warning, the name changecoords has been redefined
>
This looks quite nice except for the mess in the electric field along the negative x-axis. This is the infamous branch cut of the log function which plots as a cliff. If you can just ignore this lousy feature, the plots look nice. (Another way to visualize the electric field without this annoying cliff is to use Maple's fieldplot command, discussed in Chapter 2. A quick way to find it in this book is to click Edit on the toolbar, click Find, then search on fieldplot.)
Problem 4.14
Make contour plots just like the one above, but for (a) two line charges of equal density, one at
and the other at
, (b) two line charges with equal density magnitudes but opposite sign located as in part (a), and (c) put a line charge of strength 2 at
and one of strength 1 at
. In this part you will need a finer grid and more contours to see what's going on. See if you can identify a place where the electric field is zero.
---------------------------------------------------------------------------------------------------
>
> restart:with(plots):
> p1:=contourplot(Re(I*(x+I*y)^2),x=0..5,y=0..5,grid=[30,30],contours=40,coloring=[cyan,blue]):
> p2:=contourplot(Im(I*(x+I*y)^2),x=0..5,y=0..5,grid=[30,30],contours=40,coloring=[red,red]):
> plots[display]([p1,p2],scaling=constrained);
Warning, the name changecoords has been redefined
>
Again, the blue contours are equipotentials and the red lines are field lines. You are looking at the way electric fields behave in a grounded right-angle corner. The field is zero in the corner and the potential is proportional to
, so that if either x or y is zero, then the voltage is zero. But x=0 or y=0 defines the corner itself, so it is clearly grounded. And how does the electric field behave as we approach the corner? Remember from Physics 122 that
and
, so
is proportional to y and
is proportional to x, so as we approach the corner where x=0 and y=0, both fields go linearly to zero.
We can also find a function that represents the way fields behave in the neighborhood of a convex corner instead of a concave one by using the function
.
> restart:with(plots):
> p1:=contourplot(Re(exp(-I*Pi/6)*(x+I*y)^(2/3)),x=-5..5,y=-5..5,grid=[30,30],contours=40,coloring=[blue,blue]):
> p2:=contourplot(Im(exp(-I*Pi/6)*(x+I*y)^(2/3)),x=-5..5,y=-5..5,grid=[30,30],contours=40,coloring=[red,red]):
> plots[display]([p1,p2],scaling=constrained);
Warning, the name changecoords has been redefined
>
To see how this is the field of a corner, ignore the third quadrant and just focus on quadrants 1, 2, and 4. The blue lines are equipotentials and they get pulled around the corner, and the red field lines radiate out from the corner. In cylindrical coordinates this potential is proportional to
so the electric field, which goes like the derivative of the potential, is proportional to
.
You may be asking yourself where these weird formulas are coming from. They are called conformal mappings and you will study them in a course on complex analysis. Something fun to look forward to.
Problem 4.15
Now I'll let you do a very sharp needle-like corner. The function to use is
. Your job is to make a picture of potential and E-lines. This time the branch cut will be your friend, because it will be right where the needle is. By thinking about how this potential varies with
in cylindrical coordinates and realizing that the electric field goes like the r-derivative, find out how the electric field blows up as you approach the tip.