Application Center - Maplesoft

App Preview:

Intro to complex numbers and functions with Maple for physics students

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Chapter4.mws

Chapter 4: Complex numbers and functions

You first learned about complex numbers when an algebra teacher tried to convince you that sqrt(-1) made sense. It doesn't just make sense; it is a powerful way to model the real world we live in, in spite of i being a so-called "imaginary number". Maple knows all about the complex plane, starting with a special variable to represent sqrt(-1) : Maple calls it I ; yes, capital I . So in Maple the complex number 5 + 3i would be written this way

> z1:=5 + 3*I;

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;

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 J[0](z) and Gamma(z) . 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 z = x+I*y .

Complex arithmetic

Debugging

Let's review what you can do with complex numbers. First, there is simple addition and subtraction. And since the rule is that real parts add and imaginary parts add, if you think about the real numbers being the x-axis and the imaginary numbers being the y-axis, together forming the complex xy-plane, then it is not too hard to see that complex addition and subtraction are just addition and subtraction of vectors in this 2-dimensional plane. And Maple knows how to do it:

> z1:=5+3*I;

z1 := 5+3*I

> z2:=1-6*I;

z2 := 1-6*I

> z3:=z2+z1;

z3 := 6-3*I

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]);

[Maple Plot]

Ok: addition and subtraction are easy. Multiplication and division, however, are not anything like the vector products you know about from your physics classes. The complex multiplication rule is simple enough: just treat the real and imaginary parts like variables in algebra and multiply everything out, remembering that I^2 = -1 . Here, for example is a multiplication:

(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;

23-27*I

Division is a little trickier because I appears in the denominator. We take care of this problem with the idea of the complex conjugate. The complex conjugate of z1 is almost the same as z1, except that the sign of the imaginary part of z1 has its sign changed, like this:

> conjugate(6+4*I);

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

52

> assume(x,real);assume(y,real);

> conjugate(x+y*I)*(x+y*I);

(x-I*y)*(x+I*y)

> expand(%);

x^2+y^2

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 6^2+4^2 = 52 . 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 (3+4*I)/(2-5*I) , then do it with Maple and show that the answers agree.

---------------------------------------------------------------------------------------------

Ok, this is simple enough to do, but what does it mean? Is there a picture in the complex plane that allows us to visualize what multiplication and division mean? Well, there is; but to understand it we need to have another way of thinking about complex numbers. When we write 5+4I we are using Cartesian coordinates in the complex plane. But it is useful to use polar coordinates in the complex plane as well. In polar coordinates the number 5+4I would be represented by a radius of magnitude sqrt(5^2+4^2) = sqrt(41) , and an angle theta measured counterclockwise from the real axis, with theta = arctan(4/5) , like this

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

[Maple Plot]

>

Polar coordinates are a wonderful alternate way to visualize complex numbers. The connection between the complex number represented by x + iy and the polar form ( r, theta ) is just a simple exercise in trigonometry: x = r*cos*theta and y = r*sin(theta) . The Swiss mathematician Leonhard Euler discovered that the exponential function applied to a purely imaginary number gives the remarkable result e^(i*theta) = cos(theta)+i*sin(theta) . This means that the connection between the Cartesian and polar forms of complex numbers can be written in this compact form: x+i*y = r*e^(i*theta) . And this form, combined with the properties of the exponential function under multiplication and division, allows us to interpret what complex multiplication and division mean. Let's look at multiplication first. Writing z[1] = r[1]*e^(i*theta[1]) and z[2] = r[2]*e^(i*theta[2]) allows us to see that their product can be written as z[1]*z[2] = r[1]*r[2]*e^(i*(theta[1]+theta[2])) . So this is what multiplication means: the product has a magnitude which is the product of the magnitudes of z[1] and z[2] and an angle which is the sum of the angles of z[1] and z[2] , as shown in the picture below. The two complex numbers z[1] and z[2] are in green and their product is in red.

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

[Maple Plot]

>

Division works similarly, but with the two magnitudes being divided and the two angles subtracting: z[1]/z[2] = r[1]*e^(i*(theta[1]-theta[2]))/r[2] .

Well, that's about it for the simple arithmetic. But before proceeding onto more interesting stuff I should mention that Maple has some commands that come in handy when working with complex numbers. Re(z) picks off the real part of z, Im(z) gives you the imaginary part, abs(z) gives you the polar magnitude, and arctan(Im(z),Re(z)) gives you the angle, like this

> z:=5-3*I;

z := 5-3*I

> Re(z);Im(z);abs(z);arctan(Im(z),Re(z));

5

-3

sqrt(34)

-arctan(3/5)

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

polar(sqrt(34),-arctan(3/5))

And finally, remember that Maple will take the complex conjugate for you.

> conjugate(3+4*I);

3-4*I

>

Problem 4.2

Do the operations indicated below on the complex numbers z[1] and z[2] 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;

z1 := 3-5*I

z2 := -2+3*I

(a) z[1]+z[2] (b) z[1]-z[2] (c) z[1]*z[2] (d) z[1]/z[2]

>

----------------------------------------------------------------------------------------------

Problem 4.3

Show graphically, by making a Maple arrow plot like the one shown above, that rotating a complex number through angle theta in the complex plane is easily accomplished by multiplying by the factor e^(i*theta) . Do this by defining a complex number z1 and plotting it as an arrow, then multiplying z1 by exp(I*theta) for some angle theta 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)] .

Go to top of section

Elementary complex functions

Debugging

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 e^(i*theta) = cos(theta)+i*sin(theta) and the properties of the exponential function allow us to write for e^z (with z = x+i*y ), e^(x+i*y) = e^x*(cos(y)+i*sin(y)) . 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");

[Maple Plot]

> 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");

[Maple Plot]

Notice that at negative x it is very small and for positive x it is very big--this is the behavior of e^x . 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

Use Maple to verify the following functional identities. Do this problem by going to ?evalc , reading it over and looking at the examples, and then using this command to check these identities.

(a) cos(x) = (e^(i*x)+e^(-i*x))/2 (b) sin(x) = (e^(i*x)-e^(-i*x))/(2*i) (c) cos(i*x) = cosh(x) (d) sin(i*x) = i*sinh(x)

>

-----------------------------------------------------------------------------------------------------

Problem 4.5

Make 3-D surface plots over the complex plane of the real parts of cos(z) and sin(z) , and of the imaginary part of tan(z) , where z = x+i*y . 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.

>

-----------------------------------------------------------------------------------------------------

Now let's move on to power-law functions, i.e., we want to study the behavior of z^p in the complex plane. It is easier to see what this function does if we use the polar representation: z^p = r^p*e^(i*p*theta) . This behavior is not too hard to visualize for integer powers p . The magnitude is just raised to the p power and the imaginary exponential just wiggles faster as we go around in polar angle.

Problem 4.6

Make a 3-D plot of the real parts of the functions z^2 and z^5 using the range -3..3 in both x and y for z^2 and a smaller range for z^5 .

----------------------------------------------------------------------------------------------------

For non-integer powers things get a lot more complicated. To see why, let's consider the square root function, p = 1/2 . 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]);

[Maple Plot]

> 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]);

[Maple Plot]

> 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]);

[Maple Plot]

Surely ordinary functions can't behave this way!? Well, they can in the complex plane. Here's the problem. Maple defines its polar angle theta to be in the range -Pi .. Pi . Let's see what the square root function does to these two different values of z: z[1] = e^(i*pi) and z[2] = e^(-i*pi) . The first thing to notice is that z[1] and z[2] are exactly the same number, namely -1. The second thing to notice is that the square root function turns the first one into i and the second one into -i ! So the square root of -1 has different answers depending on what angle we use to get to -1 in the complex plane. This is the origin of the cliff face in the imaginary part plot, and the line in the complex plane where the cliff occurs is called a branch cut. Branch cuts are a pain, and you will study them in more detail when you take a course in complex analysis. All fractional powers have this problem, as does the plain old ordinary log function, ln(z), as you will soon discover.

>

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 ln(x+I*y) 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 f(z) = z^2/(1+z^2) (remember that z = x+I*y ) . 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.

Go to top of section

Special complex functions

Debugging

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

Use Maple to show that the Bessel functions J[0](z) and I[0](z) are related in the complex plane by asking Maple to evaluate J[0](i*x) and I[0](i*x) . Also do it for a few more orders of Bessel functions, say n=1,2,3. Then try the same thing for the Bessel function Y[0] . In this second case Maple will not tell you anything interesting, even though there is something interesting to be told. It would be nice if it would tell you that Y[0](i*x) = i*I[0](x)-2*K[0](x)/Pi , but I don't know how to get it to do it. Just as a check, at least, try verifying this identity numerically to see if Maple knows what it's doing. Use z = I*x for a few real values of x.

-----------------------------------------------------------------------------------------------------

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.

Go to top of section

Wave interference

Debugging

There are many ways that the complex plane helps us to do physics; perhaps the most important is in the study of phase shifts and interference patterns. This is probably a good time to give you a simple definition of the technical term phase. Phase is anything that appears as the argument of a sine or cosine function. For instance, if an oscillation in time is described by the function cos(omega*t) , then we say that the phase of the oscillation is changing linearly in time and is simply given by omega*t . If the oscillation is shifted in time so that it is given instead by cos(omega*t+phi) , then the new phase of the oscillation is omega*t+phi , and we say that relative to the previous oscillation, this one has a phase shift phi . Phase shifts can be a real pain to handle because they are stuck inside cosines and sines and you have to know lots of trig identities to get simple answers out. But if we remember that cos(omega*t) = Re(e^(i*omega*t)) (that funny script R is Maple's pretty representation of Re , the real part) then lots of formulas can be simplified. Here's an example. Suppose we have 12 different oscillations going like cos(omega*t+phi[n]) , all interfering with each other by adding together and suppose also that the nth wave has a phase shift given by phi[n] = 2*n*Pi/13 . Using the power of Maple we can add these oscillations together to get

> restart;

> s1:=Sum(cos(omega*t+2*Pi*'n'/13),'n'=1..12);

s1 := Sum(cos(omega*t+2/13*Pi*n),n = 1 .. 12)

> s1:=value(s1);

s1 := cos(omega*t+2/13*Pi)+cos(omega*t+4/13*Pi)+cos...
s1 := cos(omega*t+2/13*Pi)+cos(omega*t+4/13*Pi)+cos...
s1 := cos(omega*t+2/13*Pi)+cos(omega*t+4/13*Pi)+cos...
s1 := cos(omega*t+2/13*Pi)+cos(omega*t+4/13*Pi)+cos...

> simplify(s1);

cos(omega*t+2/13*Pi)+cos(omega*t+4/13*Pi)+cos(omega...
cos(omega*t+2/13*Pi)+cos(omega*t+4/13*Pi)+cos(omega...
cos(omega*t+2/13*Pi)+cos(omega*t+4/13*Pi)+cos(omega...
cos(omega*t+2/13*Pi)+cos(omega*t+4/13*Pi)+cos(omega...

Well, that's an answer, but it's not very pretty, and it's tough to 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

sum(e^(i*omega*t+i*phi[n]),n = 1 .. 12) .

We will take out the common factor of e^(i*omega*t) ( notice: time factored out! this is the whole reason for using the complex exponential) and then rewrite the remaining complex exponential this way

e^(i*2*Pi*n/13) = (e^(i*2*Pi/13))^n . So now our sum looks like this

e^(i*omega*t)*sum(z^n,n = 1 .. 12) where z = e^(i*2*Pi/13) . Now Maple can finish it off.

> z:=exp(I*2*Pi/13);

z := exp(2/13*I*Pi)

> s:=Sum(z^n,n=1..12);

s := Sum(exp(2/13*I*Pi)^n,n = 1 .. 12)

> s:=simplify(value(s));

s := -1

Finally, we tack on the e^(i*omega*t) we factored out front at the beginning, take the real part, and discover that the big mess of cosine functions in the original Maple sum is just:

sum(cos(omega*t+2*pi*n/13),n = 1 .. 12) = -cos(omeg...

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. Also, don't be discouraged if Maple doesn't give you a simple answer right away. Use what you know about complex numbers to give Maple some extra help.

>

Problem 4.11

(Note: this isn't really a problem. Just follow along, execute all of the Maple commands, and pay close attention; some real problems will follow.) Let's do a simple interference problem using complex analysis. The two-slit diffraction pattern is always taught in Physics 221, but only in the far field approximation. We'll use Maple to do it right. We'll have two sources of light a distance d apart and perpendicular to the line between them at distance L is a screen. We want to find the intensity of a wave on the screen as a function of distance y on the screen, measured upward from the midline from between the sources to the screen, as shown in the plot below. (Note that this shows you how to draw figures in Maple. Use ?plottools to find more drawing stuff.)

> 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

[Maple Plot]

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 e^(i*(k*r-omega*t))/r from the source. The two distances are given by r[1] = sqrt((y-d/2)^2+L^2) and r[2] = sqrt((y+d/2)^2+L^2) 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);

R1 := proc (y) options operator, arrow; sqrt((y-1/2...

R2 := proc (y) options operator, arrow; sqrt((y+1/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;

Intensity := proc (y) options operator, arrow; Re(e...

>

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 a special 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)> = int(f(t),t = 0 .. T)/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

Int(f(t),t = 0 .. T)/T = Sum(f(t[i])*dt)/T but T/dt = N , the approximate number of data points, so we have Int(f(t),t = 0 .. T)/T = Sum(f(t[i]))/N , 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(%);

1/2*Int(sin(t),t = 0 .. 2*Pi)/Pi

0

>

> Int(cos(t),t=0..2*Pi)/(2*Pi);value(%);

1/2*Int(cos(t),t = 0 .. 2*Pi)/Pi

0

>

> Int(cos(t)*sin(t),t=0..2*Pi)/(2*Pi);

1/2*Int(cos(t)*sin(t),t = 0 .. 2*Pi)/Pi

> value(%);

0

> Int(cos(t)^2,t=0..2*Pi)/(2*Pi);

1/2*Int(cos(t)^2,t = 0 .. 2*Pi)/Pi

> value(%);

1/2

> Int(sin(t)^2,t=0..2*Pi)/(2*Pi);value(%);

1/2*Int(sin(t)^2,t = 0 .. 2*Pi)/Pi

1/2

> value(%);

1/2

You need to memorize these results: the time averages of sin(t) , cos(t) , and sin(t)*cos(t) are all zero. The time averages of cos(t)^2 and sin(t)^2 are both 1/2. These answers 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, having learned this 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;

Intense := (1/r1*cos(k*r1-omega*t)+1/r2*cos(k*r2-om...

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;

Intense := (cos(x1)/r1+cos(x2)/r2)^2

Now expand it like this

> Intense:=expand(Intense);

and then time average each term.

Intense := cos(x1)^2/r1^2+2*cos(x1)/r1*cos(x2)/r2+c...

The two squared terms are easy: both x1 and x2 have omega*t in them so they just time average to 1/2 (the phase shift due to k*r1 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:

> combine(cos(x1)*cos(x2));

1/2*cos(x1-x2)+1/2*cos(x1+x2)

(Note: to get trig identities try expand, simplify, and combine .) The x1+x2 term goes like 2*omega*t 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);

Intense := 1/2*1/(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);

Intensity := proc (y) options operator, arrow; 1/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;

d := .1e-4

L := 1

k := 3333333.333*Pi

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

[Maple Plot]

>

---------------------------------------------------------------------------------------------------

Problem 4.12

Serway says that the first minimum away from y=0 is at the place where k*d*sin(theta) = Pi , where theta is the angle indicated below ( y = L*tan(theta) ). 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 theta . 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);

[Maple Plot]

>

------------------------------------------------------------------------------------------------

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 r[1] and r[2] 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.

Go to top of section

Electrostatics with line charges

Debugging

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 F(z) = -lambda*ln(z-z[0])/(2*Pi*epsilon[0]) has the property that its real part is the electrostatic potential of a line charge with charge/length lambda located in the xy-plane at complex position z[0] , 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

[Maple Plot]

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 z = i and the other at z = -i , (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 z = -i and one of strength 1 at z = i . 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.

---------------------------------------------------------------------------------------------------

It is also possible to use complex functions to study the behavior of potential and voltage near conducting corners. Consider the following contour plot of the function i*z^2 .

>

> 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

[Maple Plot]

Recalling that the potential is given by the real part and the field lines by the imaginary part we have

i*z^2 = -x*y+i*(x^2-y^2)

and in the plot the blue contours are equipotentials and the red lines are field lines. With V = -x*y both the x-axis and the y-axis are grounded, so you are looking at the way electric fields behave in a grounded right-angle corner. You can see the red field lines far apart from each other in the corner, indicating that the electric field is zero there. And how does the electric field behave as we approach the corner? Remember from electricity and magnetism that E[x] = -diff(V,x) and E[y] = -diff(V,y) , so E[x] is proportional to y and E[y] 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 e^(-i*Pi/6)*z^(2/3) .

> 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

[Maple Plot]

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 tightly around the corner, and the red field lines radiate out from the corner. In cylindrical coordinates this potential is proportional to r^(2/3) so the electric field, which goes like the derivative of the potential, is proportional to 1/(r^(1/3)) , which goes to infinity at the corner.

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. It's 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 sqrt(z) . 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 r 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.

Go to top of section

Go to top of chapter