<<VectorAnalysis`

<<Miscellaneous`RealOnly`

Calculus IV
Lab 2, Part 1

Non-Cartesian coordinate systems

Changing your own coordinate system

The basics

Sometimes, the coordinate system you need to use to work an integral is fairly obvious.  As an example, consider finding given integral over the region between the following sets of curves: 

(Underscript[∫∫, D])x ydA where D is the region between the graphs of:
	3x-y=1, y=3x-6,x y=1, and x y=2 in the first quadrant.

If you think about it, this is really the region:

	1≤3x-y≤6 and 1≤x y≤ 2

which looks like:

RegionPlot[1≤3x-y≤6&&1≤x y≤2,{x,-5,5},{y,-8,8},Axes->Automatic,PlotPoints->80,Frame->False]

Now, you could do this using 3 separate integrals.  However, by using a change of coordinates, with u=3x-y and v=x y, you could work it using a single integral.  So, let's set this up as a change of coordinates.  I'm going to define two transformation:

XYtoUV:(ℝ^2)->(ℝ^2) maps from D->(SuperStar[D]) (takes an input of {x,y} and returns a point in {u,v}) 
UVtoXY:(ℝ^2)->(ℝ^2) maps from (SuperStar[D])->D (takes an input of {u,v} and returns a point in {x,y}) 

(Remember that you take the Jacobian of UVtoXY; in the notation we use in class, XYtoUV would be (T^(-1)) and UVtoXY would be T.)

XYtoUV[{x_,y_}]:={3x-y,x y}

This says u=3x-y and v=x y.  Now, the function we want to integrate:

f[{x_,y_}]:=x y

Now, in order to find T  (or UVtoXY), we need to solve the system of equations for x and y (in terms of u and v).  We can do this with the solve command:

solveInverse=Solve[XYtoUV[{x,y}]=={u,v},{x,y}]

Notice that I gave the solution set a name (solveInverse), so I can refer back to it.  Also, notice that there are two sets of solutions.  Should this be a surprise?  (Hint:  Check out the graph earlier...)  If you look at these carefully, you will see that the second set is the one in the first quadrant (for u and v in the ranges we nned).

Further, notice that the solutions are not given in the form {x=blah,y=bleh}, they are given in the form {x->blah,y->bleh}.  These are called "replacement rules", which means that they are a set of instructions ("replace x with blah and replace y with bleh").  We can use the "replace all" operator /. to apply the rules.  So, for example, the following code:

(x^2)+2x-y/.x->3

replaces all the x's in the expression (x^2)+2x-y with x=3.  So, we could assign these solutions to set up our UVtoXY function as follows:

UVtoXY[{u_,v_}]={x,y}/.solveInverse[[2]]

This replaces {x,y} with the solutions we found above (called solveInverse).  Since we only wanted the solution in the first quadrant (the second one), we used solveInverse[[2]].

Now, we also need the Jacobian determinant.  In Mathematica, you can find the Jacobian of a function by:

D[function[var1,var2,va3,...],{{var1,var2,va3,...}}]

(Notice the double {]'s enclosing the variables.  So, for example:

D[{(x^2)(y^3),x+3Sin[y]},{{x,y}}]

Or, to make it look cooler:

D[{(x^2)(y^3),x+3Sin[y]},{{x,y}}]//MatrixForm

We then need to take the determinant; so, here we can do the following:

detJ=Det[D[UVtoXY[{u,v}],{{u,v}}]]//FullSimplify

(I use the //FullSimplify at the end of the line to make sure Mathematica has simplified it as much as possible before I go on.)  Notice, this will always be positive, so I don't need to take the absolute value.

Now, we need to take the composition of f and our transform T, i.e., f(T(u,v)):

fUV[{u_,v_}]=f[UVtoXY[{u,v}]]//FullSimplify

(Just a kindly word of warning here:  In this problem, f(T(u,v)) simplified out really nicely.  However, it is often the case that it will be much more complicated than f(x,y) was.  If it is too complicated, the transformed integral may be too hard for Mathematica to solve and you may have to go back and try something else...)

Finally, pull it all together:

(∫_1^2)(∫_1^6)fUV[{u,v}]detJ\[DifferentialD]u\[DifferentialD]v

You could certainly combine several of these steps together, but hopefully this gives you an idea how to work your way through one of these.  If you want to, you could do the integral without using the fUV and detJ expressions:

(∫_1^2)(∫_1^6)f[UVtoXY[{u,v}]]Det[D[UVtoXY[{u,v}],{{u,v}}]]\[DifferentialD]u\[DifferentialD]v

(The real problem with this is that you don't explicitly see what the determinant of the Jacobian looks like, so you don't really know whether you have the right sign for the Jacobian determinant.  If it is negative everywhere, then your answer will only be off by a sign, but if it is negative in some parts of your domain and positive in others, you have a much bigger problem.)

Variation on a theme

Sometimes, your region of integration is a variant on something you are already familiar with.  For example, if you wanted to integrate:

(Underscript[∫∫∫, D])((x )^2)(y^4)-(z^2)dV where D is the region between the graphs of:
	2(x^2)+4x+3(y^2)+(z^2)-3z=1 and 2(x^2)+4x+3(y^2)+(z^2)-3z=4
	
First, let's see what this actually looks like:

regionD=RegionPlot3D[1≤ 2(x^2)+4x+3(y^2)+(z^2)-3z≤4,{x,-5,5},{y,-5,5},{z,-5,5},PlotPoints->50]

Well, it looks like an elllipsoid (or the region between two ellipsoids, actually, if you were to cut it open).  In fact, if you stare at the eqution a bit, this should become obvious to you.  We can complete the squares to find out what is going on:

2(x^2)+4x+3(y^2)+(z^2)-3z
=2((x^2)+2x)+3(y^2)+((z^2)-3z)
=2(((x+1))^2)-2+3(y^2)+(((z-(3/2)))^2)-(9/4)
So we get:
1≤2(((x+1))^2)-2+3(y^2)+(((z-(3/2)))^2)-(9/4)≤4 
(21/4)≤2(((x+1))^2)+3(y^2)+(((z-(3/2)))^2)≤(33/4) 

(There is probably a way to make Mathematica complete the squares for you, but unless I'm missing something, it's actually simpler just to do it by hand...)  So, what we have is a set of nested ellipsoids, centered at (-1,0,(3/2)).  A slightly different way to think of this is as spheres shifted from the origin to (-1,0,(3/2)) and stretched by a factor of (Sqrt[1/2])in the x direction and (Sqrt[1/3])in the y direction (why?).  

So, one thing we could do here is to take the standard coordinate transform for spherical coordinates and "distort" it in exactly this way.  So, if our standard spherical transform is:

mySphere[{ρ_,θ_,ϕ_}]:={ρ Cos[θ]Sin[ϕ],ρ Sin[θ]Sin[ϕ],ρ Cos[ϕ]}

Our distorted ellipsoidal form would be:

myEllipsoid[{ρ_,θ_,ϕ_}]:={(ρ/(Sqrt[2])) Cos[θ]Sin[ϕ]-1,(ρ/(Sqrt[3])) Sin[θ]Sin[ϕ],ρ Cos[ϕ]+(3/2)}

Notice that, in this case:

(Sqrt[(21/4)])≤ρ≤(Sqrt[(33/4)]), 0≤θ≤2π, 0≤ϕ≤π

(Why did I take the square root of the numbers for ρ?)  Let's graph this for the outer ellipse to check that we haven't done anything wrong here:

testRegion=ParametricPlot3D[myEllipsoid[{(Sqrt[33/4]),θ,ϕ}],{θ,0,π},{ϕ,0,π}]

Since I have a clear domain for ρ, θ, and ϕ, I use ParametricPlot3D here instead of InequalityPlot3D.  (Notice that I only used half the domain for θ.  This is so when I draw the two graphs together, I can see if they really coincide.)

Show[regionD,testRegion]

From here on out, it's just the usual stuff:

jDet=Det[D[myEllipsoid[{ρ,θ,ϕ}],{{ρ,θ,ϕ}}]]//FullSimplify

jDet=-jDet

(This will now be positive everywhere in our domain...)

f[{x_,y_,z_}]:=((x )^2)(y^4)-(z^2)

(∫_0^π)(∫_0^(2π))(∫_(Sqrt[21/4])^(Sqrt[33/4]))f[myEllipsoid[{ρ,θ,ϕ}]]jDet\[DifferentialD]ρ\[DifferentialD]θ\[DifferentialD]ϕ

(If you were working this by hand, it might actually be simpler to set this up as several separate transformations, perhaps one to stretch and shift the coordinate system and another to convert to spherical coordinates.  That would make the computations, especially the Jacobian, much simpler.  However, since we are using Mathematica for all the heavy duty algebra lifting here, I like combining it all together into one change of coordinates; I just think it's a little cleaner...)

Problems

Use an appropriate change of coordinates to evaluate the following integrals (graph the region of integration as well):

(Underscript[∫∫, D])x (y^2)dA where D is the region between the graphs of:
	y=((2-2(x^2))^(1/3)), (x^2)=3-((y^3)/2), (y^3)=2x-1, and x=4+((y^3)/2) for x≥0.

(Underscript[∫∫∫, D])x((y^2)+(z^2))dV where D is the region between the graphs of:
	3≤(y^2)+2y+4(z^2)≤8, -x +2y+3z=15  and  x=(y^2)-3(z^2)
	
(Kindly hint:  The coordinate system here is given by the first compound inequality.)

Using Mathematica's built-in coordinate systems

Mathematica supports a very wide range of 3-dimensional coordinate systems.  (Double-click on VectorAnalysis and then hit function key F1 to go to the place in the Help Browser where these are discussed.)  This can be very useful, but it can also be very confusing.  How do you know what these different coordinate systems "look like" and when to use each one?

Well, this lab won't teach you that.  However, it will show you some techniques to investigate this question and how to work with these alternate coordinate systems in Mathematica.

First, we will work through an investigation of some properties of different coordinate systems using the two main non-Cartesian systems that we are most familiar with.  Then, you will apply these techniques to some less familiar ones.  (For all the work that follows, you must load the library module VectorAnalysis`.  In this notebook, it is done automatically when you execute the first command.)

Since we are in 3 dimensions here, we need 3 coordinates to specify a point in space.  To find out what names Mathematica gives to these coordinates by defaulty, use:

Coordinates[Cylindrical]

Notice, Mathematica has built-in variable names for these, but they aren't very convenient to use.  I can tell Mathematica to use my coordinate variables directly, however.  So, what are the transformation equations to convert from Cylindrical to Cartesian coordinates?

CoordinatesToCartesian[{r,θ,z},Cylindrical]

(Notice that I used my own variable names here.  This makes the formulae look nicer. It can get you in trouble on occasion, however.  More on that later.)

So, what about Spherical coordinates?

Coordinates[Spherical]

We can get some further info on the possible ranges of these coordinates:

CoordinateRanges[Spherical]

Hmmm, something isn'g quite right here.  Notice that Mathematica reverses the roles of Theta and Phi here (and uses a slightly different domain for what we know as Theta).  We can fix that by specifying what variables Mathematica should use here as follows:

CoordinatesToCartesian[{ρ,ϕ,θ},Spherical]

So, an important warning here:  For some reason, Mathematica reverses the standard roles of θ and ϕ in spherical coordinates.  In Mathematica, θ stands for the angle from the Z axis, while ϕ is the angle of the projection into the XY plane from the X axis.  This isn't all that hard to deal with (just always remember to enter points in {ρ,ϕ,θ} format, rather than the more usual {ρ,θ,ϕ} form, using our standard θ and ϕ), but you do have to be careful to remember it.  Also, notice that they have chosen the range for their ϕ (our θ) to be [-π,π] rather than [0,2π].

Visualization

So, how can you visualize what these alternate coordinate systems "look like"?  This can be a bit tricky, so I will show you two ways to do this.  

Consider a coordinate system {u,v,w}.  If you fix w and let u and v vary, you will get a surface (why?).  If you do this for different values of w, you will get a bunch of surfaces that are "parallel" in some sense (literally, in Cartesian coordinates).  These can be thought of as "level surfaces" for the coordinate system:  everything on one of these surfaces has the same w coordinate.  You can do the same thing for u and v by letting each of them stay constant.  Putting all these together give a pretty good picture of how the coordinate system behaves.  Let's take a look at this for Cartesian and Spherical coordinates:

First, you use the Table function to generate your "level surfaces" for each combination of variables.  Notice that the Table function works like:

Table[ Stuff to do for different values of k, {k, kStart, kEnd, kStep}]

where:
kStart = The beginning value to plug in for k
kEnd = The final value to plug in for k
kStep = The amount to increment k by each time

So, when I use the command:
Table[ParametricPlot3D[ stuff...], {k,-4,4,2}]

I am telling it to make a series of ParametricPlot3D's, letting k = -4, -2, 0, 2, 4.

zSlices=Table[ParametricPlot3D[{u,v,k},{u,-5,5},{v,-5,5}],{k,-4,4,2}]

ySlices=Table[ParametricPlot3D[{u,k,v},{u,-5,5},{v,-5,5}],{k,-4,4,2}]

xSlices=Table[ParametricPlot3D[{k,u,v},{u,-5,5},{v,-5,5}],{k,-4,4,2}]

Notice that in each of these, u and v are used as variables for graphing the surface, while k always stands for the variable that will get held constant for each surface.  Now, let's put all this together in a useful picture:

Show[xSlices,ySlices,zSlices,PlotRange->All]

Not surprisingly, we see that Cartesian coordinates can be thought of as a bunch of intersecting planes.  (I'm not really sure why you need the PlotRange->All here, but you can always give this a try if part of your graph is missing...)  You locate points at the intersection of the three perpindicular planes.  This is boring, but makes for a good start.

A couple of things you should notice about how I set this up:  

I didn't let k run over as large a domain as I did u and v.  This was so that when I put everything together at the end (using the Show command), we could actually see the intersecting planes.  (If you let k run from -5 to 5, like u and v, all you see is a cube.  Perfectly correct, but not very informative.)  

You will probably have to experiment with the domain for u and v, as well as the range for k (and the number of slices) in different coordinate systems to come up with an informative final graph.

Let's do the same thing for Spherical coordinates.  Notices that it is very similar, though we have to convert to Cartesian coordinates before we can plot anything (ParametricPlot3D plots only in Cartesian).  Also, notice the use of Evaluate before the change of coordinates; this is just to stop Mathematica complaining (it compiles the function before plotting).

ρSlices=Table[ParametricPlot3D[Evaluate[CoordinatesToCartesian[{k,ϕ,θ},Spherical]],{ϕ,0,π},{θ,-((2π)/3),π}],{k,1,5,2}]

Notice that I didn't complete the spheres (I restricted θ).  I did this so that when we put the graphs together, you can get a "cut-away" view inside (see below).

θSlices=Table[ParametricPlot3D[Evaluate[CoordinatesToCartesian[{ρ,ϕ,k},Spherical]],{ϕ,0,π},{ρ,0,6}],{k,-((3π)/6),((5π)/6),π/3}]

ϕSlices=Table[ParametricPlot3D[Evaluate[CoordinatesToCartesian[{ρ,k,θ},Spherical]],{θ,-((2π)/3),π},{ρ,0,6}],{k,0,π,π/6}]

Show[ρSlices,θSlices,ϕSlices,PlotRange->All]

Again, you locate points at the intersection of the three types of surfaces (planes, spheres, and cones).

You can make this a little more interactive by using the Manipulate command.  Let's set it up so we have sliders for each of the 3 coordinates ρ, θ, and ϕ:

Manipulate[Column[
{"{ρ,θ,ϕ}=",
{rho,theta,phi},
Show[
ParametricPlot3D[Evaluate[CoordinatesToCartesian[{ρ,phi,θ},Spherical]],{θ,-π,π},{ρ,0,6}],
ParametricPlot3D[Evaluate[CoordinatesToCartesian[{rho,ϕ,θ},Spherical]],{θ,-π,π},{ϕ,0,π}],ParametricPlot3D[Evaluate[CoordinatesToCartesian[{ρ,ϕ,theta},Spherical]],{ϕ,0,π},{ρ,0,6}],PlotRange->All]
}],
{rho,0,6},{theta,-π,π},{phi,0,π}]

By moving rho, theta, and phi around, you can get a good idea of the effect each of these coordinates have (your point is at the intersection of the 3 surfaces, of course).  This is kind of overkill for spherical coordinates, but can be very helpful for some of the other more bizarre coordinate systems.  (Notice that you will need to change the domains of all the coordinates to fit your other coordinate system.  You may need to play around with that a bit to get the best picture/performance.)

Problem:

Pick two other non-standard coordinate systems that Mathematica supports (listed in the Help Browser page I mention above; you may not use Cartesian, Cylindrical, or Spherical) and generate visualizations of those systems like I did above.  Also, have Mathematica give you the formula to translate from these other coordinate systmes into Cartesian and find the ranges Mathematica expects for each coordinate.  Choose your domains, etc., so that your final picture is easy to interpret.  Try to explain how points are located in space (i.e., at the intersection of a cylinder around the Z axis, a plane parallel to the XY plane, and a plane containing the Z axis would be an explanation for Cylindrical).  Set up a Manipulate expression for it and see if that agrees (or clarifies) you explanation.

Warning:  Some of the coordinate systems in the list seem to have more than 3 coordinates.  This isn't really true; the first 3 variables listed are always the coordinates.  The "extra" variables are actually constant parameters; you get a slightly different coordinate system for each choice of those parameters.  You can use the default values for each one for this exercise.