Calculus IV Lab 1: Functions, Graphs, and Integrals

Functions, operators, and graphs

Composition of functions

There are different ways to define a function of more than one variable in Mathematica.

f[x_,y_,z_,t_]:=((x y z Sin[2π t])/((x^2)+(y^2)+(z^2)))

g[x_,y_]:={y,(x^2)}

h[t_]:={Sin[t],Cos[t],t,2t}

Taking the composition of two functions is relatively easy (don't forget to evaluate the functions above first):

h[f[x,y,z,t]]

However, try:

g[g[x,y]]

This doesn't work.  It SHOULD work, but it doesn't. To get around this problem, we will often define functions of more than one variable in the following way:

g[{x_,y_}]:={y,(x^2)}

Now:

g[g[x,y]]

Will f[h[t]] work with the above definitions?  Should it?  If so, modify the definitions so it does.

Functions as arguments (Operators)

Notice that if you have a function of a single variable, Mathematica knows what you mean by:

h'[t]

Partial derivatives are also pretty straightforward:

(∂_x) g[x,y]

Sometimes, it is useful to define a function that operates on other functions.  So, for example, I could define a function to find the derivative matrix for functions from (ℝ^2)->(ℝ^2):

deriv2By2[expression_]:={(∂_x) expression,(∂_y) expression}//Transpose//MatrixForm

deriv2By2[{x Sin[y],x Cos[y]}]

deriv2By2[g[x,y]]

(Notice that I needed to transpose the matrix to get the rows/columns to work out like we want.  I used MatrixForm just because it looks cool.)  However, you have to be really careful with these.  For example, if I just make one minor little change when calling it:

deriv2By2[{u Cos[v],u Sin[v]}]

The problem is that, when we defined our function, we assumed that our variables would be called x and y, so any other choice of variables gets interpreted as a constant.  One way to fix this would be to explicitly specify your independent variables in the definition:

deriv2By2b[expression_,var1_,var2_]:={(∂_var1) expression,(∂_var2) expression}//Transpose//MatrixForm

deriv2By2b[g[u,v],u,v]

deriv2By2b[{u Cos[v],u Sin[v]},u,v]

This now works, though at the cost of a bit of redundancy in calling the function.  (There are fancier tricks we could use to generalize this, but this is a good starting point.)  

Actually, there is a built-in extension of the standard D (derivative) function in Mathematica that will take care of this particular problem much more easily:

D[g[{x,y}],{{x,y}}]//MatrixForm

D[{u Cos[v],u Sin[v]},{{u,v}}]//MatrixForm

Be sure you use the double curly brackets {{}} for the "with respect to" variables; if you only used a single pair, it would try to take a second derivative, once with respect to the first and the second time with respect to the other variable.  

If you need to access one of the components of a vector, you can use the notation:

g[x,y][[1]]

g[x,y][[2]]

Define your own functions to compute the gradient (call it myGrad), divergence (call it myDiv), and curl (myCurl) of an arbitrary function of 3 variables.  (You may not use any of the built-in functions that handle these.)  Demonstrate that each of these actually work for a few different functions, using different independent variable names.

ParametricPlot3D with spherical and cylindrical coordinates

Mathematica has special built-in functions for graphing cylindrical and polar coordinates (SphericalPlot3D).  Unfortunately, these aren't really as versatile as we might wish.  For example, they wouldn't be capable of plotting something as simple as ϕ=(π/6).  These commands are based on a much more flexible command, ParametricPlot3D.

ParametricPlot3D requires a function from either ℝ->(ℝ^3) or (ℝ^2)->(ℝ^3) as an input.  It will then graph the range of the function.  Thus:

curve[t_]:={t,Sin[t],Cos[2t]}

ParametricPlot3D[curve[t],{t,0,3π}]

If you like, you can use ParametricPlot3D to draw a surface:

surface[u_,v_]:={u Sin[v],u Cos[v],(u^2)}

ParametricPlot3D[surface[u,v],{u,-1,1},{v,0,2π}]

You can even use ParametricPlot3D to graph a function of 2 variables:

f[x_,y_]:=(x^2)+(y^2)

ParametricPlot3D[{x,y,f[x,y]},{x,-1,1},{y,-1,1}]

(If you had used the Plot3D function instead, it would have looked a little nicer...)

These two graphs look suspiciously alike.  Show that they are really the same surface.  Why do they look so different?

In order to use ParametricPlot3D to graph things in spherical and cylindrical coordinates, you must first develop a "change of variable" function. For example, a change of variable function for polar coordinates might look something like:

polarToRect[{r_,θ_}]:={r Cos[θ],r Sin[θ]}

rectToPolar[{x_,y_}]:={(Sqrt[(x^2)+(y^2)]),ArcTan[(y/x)]}

polarToRect[{2,(π/6)}]

rectToPolar[{(Sqrt[3]),1}]

rectToPolar[polarToRect[{2,(π/6)}]]

Warning:

rectToPolar[{-(Sqrt[3]),1}]

What is wrong with this and why?

Fortunately, we are mainly interested in polarToRect (actually the equivalents for cylindrical and spherical coordinates).

Create 2 new functions: cylinderToRect (accepts points in the form {r,θ,z}) and sphereToRect (accepts points in the form {ρ,θ,ϕ}) to make the appropriate conversions.  Notice that both of these functions are from (ℝ^3)->(ℝ^3).

You can use these new functions in conjunction with ParametricPlot3D to graph most equations in cylindrical and spherical coordinates. (You can even use it to graph cylindrical or spherical "parametric equations," something you don't see in many textbooks...)

For example, to graph something like ρ=Sin[θ]Cos[ϕ], you would do something like:

ρ=Sin[θ]Cos[ϕ]

ParametricPlot3D[sphereToRect[{ρ,θ,ϕ}],{θ,0,2π},{ϕ,0,2π}]

(Notice that you can only have 2 independent variables inside ParametricPlot3D, θ and ϕ in this case.)

Of course, the built-in SphericalPlot3D function could have done this (with less effort). However, you can also use your new functions to plot θ=(π/4):

ParametricPlot3D[sphereToRect[{ρ,(π/4),ϕ}],{ρ,0,5},{ϕ,0,π}]

Sometimes, it is important to impose limitations on your domain or range.  If they are simple limitiations, you can just take care of it in the ParametricPlot domain indications.  (In the previous example, we know that ρ≥0 and 0≤ϕ≤π, so we assigned those domains.)  However, sometimes finding the appropriate region to graph is sometimes more complicated.  Let's say you want to plot the graph of ϕ=ρ +sin(θ).  Well, it is certainly possible to compute a domain that will work for this, but it can get a bit tricky.  Mathematica has a built-in method to handle this, called RegionFunction.  You want to ensure that 0 ≤ ϕ ≤ π, so that means 0≤ρ +sin(θ)≤π.  You can graph this (using your sphereToRect function) as follows:

ParametricPlot3D[sphereToRect[{ρ,θ,ρ +Sin[θ]}],{ρ,0,5},{θ,0,2π}]

However, that doesn't satisfy the restrictions on ϕ in spherical coordinates.  To have Mathematica take that into account for you, try:

ParametricPlot3D[sphereToRect[{ρ,θ,ρ +Sin[θ]}],{ρ,0,5},{θ,0,2π},RegionFunction->Function[{θ,ρ},0≤ρ +Sin[θ]≤π]]

This looks weirder.  The reason is that I am using such a large domain for ρ (I get some disconnected sections).

ParametricPlot3D[sphereToRect[{ρ,θ,ρ +Sin[θ]}],{ρ,0,π},{θ,0,2π},RegionFunction->Function[{θ,ρ},0≤ρ +Sin[θ]≤π]]

Of course, it really depends on what you are doing whether this sort of restriction is necessary.

Use the functions you developed above to graph the following over domains appropriate to the coordinate system:

1. ϕ=ρ sin((θ/2))

2. r=sin(t), θ=(t^2), z=cos(t)  
This is a set of parametric equations in cylindrical coordinates; not something you see in textbooks very often, but not a problem in Mathematica.

3. (ρ^4)(sin^2)(θ)=(ϕ^3)

Stupid Integration Tricks

Viewing the region

There are several ways you can view regions over which you need to do double and triple integrals.  The RegionPlot and RegionPlot3D functions are theoretically the most useful.  So, for example, if you know a general description of the region:  "between the graph of (x^2)=(y^2)+2(z^2) and the YZ plane, and inside (y^2)+(z^2)=4", you could graph it as follows:

RegionPlot3D[(x^2)<(y^2)+2(z^2)&& (y^2)+(z^2)<4,{x,0,5},{y,-5,5},{z,-5,5},PlotPoints->100]

As you can see, even with using a relatively high PlotPoints value, you still get a pretty vague picture.  Also, it can often be helpful to see the region in the context of intersecting surfaces.  You could use ContourPlot3D to manage that:

ContourPlot3D[{(x^2)==(y^2)+2(z^2), (y^2)+(z^2)==4,x==0},{x,0,4},{y,-5,5},{z,-5,5}]

You can play around with some of the display options to make it a little clearer:

ContourPlot3D[{(x^2)==(y^2)+2(z^2), (y^2)+(z^2)==4,x==0},{x,0,4},{y,-5,5},{z,-5,5},ContourStyle->Opacity[0.6],Mesh->False,AxesLabel->Automatic]

Which of these versions is most useful really depends on the region involved.  

If you start with limits of integration in an integral, you can use these same methods:

(∫_0.5^1)(∫_(Sqrt[1-(x^2)])^(Sqrt[4-(x^2)]))(E^(y^3))\[DifferentialD]y\[DifferentialD]x

RegionPlot[0.5<x<1&&(Sqrt[1-(x^2)])<y<(Sqrt[4-(x^2)]),{x,0,3},{y,0,3},PlotPoints->100]

ContourPlot[{x==0.5,x==1,y==(Sqrt[1-(x^2)]),y==(Sqrt[4-(x^2)])},{x,0,3},{y,0,3}]

Graphing these regions can be helpful if you are trying to determine limits of integration (in the first case) or reverse the limits of integration (second case).  If you are working in a different coordinate system, you will need to combine a coordinate transform with these commands (as in the previous section of the lab).

Finding the intersection

Another thing Mathematica is good at is doing our algebra for us.  There are several Mathematica functions you can use to help out with this, but the ones probably most useful here are Solve, Reduce, and Eliminate.  So for example, if you want to integrate over the region between z=(x^2)+2(y^2) and (x^2)+(y^2)+(z^2)=1, you need to see what is going on:

ContourPlot3D[{z==(x^2)+2(y^2),(x^2)+(y^2)+(z^2)==1},{x,-1,1},{y,-1,1},{z,-1,1},ContourStyle->Opacity[0.6],Mesh->False,AxesLabel->Automatic]

So, it seems reasonable to let your inner limit be dz, and your outer limits of integration be dx dy or dy dx, so we need an equation of the trace of the curve of intersection projected into the XY plane.  Let's try:

trace1=Evaluate[Eliminate[{z==(x^2)+2(y^2),(x^2)+(y^2)+(z^2)==1},z]//FullSimplify]

Let's see what it look like:

ContourPlot[Evaluate[trace1],{x,-1,1},{y,-1,1},Axes->True]

(I'm not really sure why we need Evaluate here; there's probably some subtle Mathematica-type reason here - or a bug - but sometime you need to enclose the expression you are trying to graph inside Evaluate.)  So, let's say we want our outer limits of integration to be dy dx; we need to find a formula for y in terms of x:

Solve[trace1,y]//FullSimplify

To make it easier to work with these solutions, let's try this:

yLimits=y/.Solve[trace1,y]//FullSimplify

FOUR solutions?  Something odd is going on; let's graph each solution seperately:

Table[Plot[Evaluate[yLimits[[k]]],{x,-1,1}],{k,1,4}]

Notice, the first two don't give any graphs.

Explain why the first two solutions don't yield graphs.

Now, all we need are the limits for x.  Well, that's easy enough, just set one of the solutions we just found for y to 0:

xLimits=x/.Solve[yLimits[[4]]==0]

So, finally, our iterated integral might become:

(∫_(xLimits[[1]])^(xLimits[[2]]))(∫_(yLimits[[3]])^(yLimits[[4]]))(∫_((x^2)+2(y^2))^(Sqrt[1-(x^2)-(y^2)]))f[{x,y,z}]\[DifferentialD]z\[DifferentialD]y\[DifferentialD]x

You generally have to handle each problem seperately, but this gives you an idea of how these commands can be used.  Pay particular attention to how I combined looking at graphs with the algebra.  When I first got 4 solutions for y above, the graphs helped me figure out what was going on.

Find appropriate limits of integration for and compute the following integral.  You should do all computations within Mathematica and present appropriate graphs to support your choice.

(Underscript[∫∫∫, D])(x^2)y zdV,
where D is the region defined by (((y+1))^2)>(x^2)+(z^2) and (x^2)+(y^2)+(z^2)<1 and y + 2 z > -1.  (Hint:  You will need to chop this region up into several parts.  I did it with 3 seperate regions; you may need 3 or 4.)

Redefining the function

If you will recall, in class, we defined the double integral of a function f over a non-rectangular region D as the integral of the function (SuperStar[f]) over any rectangular region that contains D, where:

(SuperStar[f])(x,y)=(Piecewise[{{f(x,y)    , (x,y)∈D}, {0, Otherwise}}])

At the time, if you recall we discussed how this was more useful as a definition than as a useful method to actually do anything.  Well, now Mathematica can make me eat my words...  A much easier method to compute nasty double and triple integrals in Mathematica than finding iterated integrals and regions of integration is to define a new piecewise function that equals our integrand on the region of integration, but 0 everywhere else.  There are some limitations on this method, but it can save us quite a bit of work, especially if we are interested in a numerical approximation of an integral.  So, for example:

(Underscript[∫∫∫, D])(x^2)y(Sqrt[(z^2)+1])\[DifferentialD]V, where D is the region between the graph of(x^2)=(y^2)+2(z^2)and the YZ plane,and inside(y^2)+(z^2)=4

We could look back to the previous section where we graphed this and realize that this region can be described as:

0≤x≤ (Sqrt[(y^2)+2(z^2)]),-2≤y≤2,-(Sqrt[4-(y^2)])≤z≤(Sqrt[4-(y^2)])

So, we could integrate it as:

(∫_(-2)^2)(∫_(-(Sqrt[4-(y^2)]))^(Sqrt[4-(y^2)]))(∫_0^(Sqrt[(y^2)+2(z^2)]))x (y^2)((z^2)+1)\[DifferentialD]x\[DifferentialD]z\[DifferentialD]y

38(π/3)//N

If you had the urge, this could also be done in a pseudo-cylindrical coordinate system, where the cylinder circles the X axis.

Now, let's try this the old and useless/new and easy way.  First, let's define our integrand:

fRestriction[{x_,y_,z_}]:=Piecewise[{{(x (y^2)((z^2)+1)),(x^2)≤(y^2)+2(z^2)&&(y^2)+(z^2)≤4&&x≥0}},0]

Alternately, we could have used the fancy new way to write this as:

fRestriction[{x_,y_,z_}]:=(Piecewise[{{(x (y^2)((z^2)+1)), (x^2)≤(y^2)+2(z^2)&&(y^2)+(z^2)≤4&&x≥0}}]) 

using the Piecewise operator to define this function.  (To enter this, hit the <Escape> key on your keyboard (upper left corner), type "pw", then <Excape> again (i.e., <Esc>pw<Esc>).  To enter the condition, type <Control>, and <Control><Return> if you need more pieces.  Notice, this is not a standard curly bracket. If your input point doesn't satisfy the condition, the function will be assigned 0.)

So the above is equivalent to:

fRestriction[{x_,y_,z_}]:=(Piecewise[{{((x^2)(y^2)(Sqrt[(z^2)+1])), (x^2)≤(y^2)+2(z^2)&&(y^2)+(z^2)≤4&&x≥0}, {0, Otherwise}}])

Now, I need to know a box that is guaranteed big enough to include everything.  How about [-∞,∞]×[-∞,∞]×[-∞,∞]?  Now, we can integrate:

Integrate[fRestriction[{x,y,z}],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}]

This is, sadly, much slower than it was before, but does work eventually, if you go get some coffee...  (Reducing the box size doesn't really help.)  However, if you're looking for a quick, numerical result, how about trying this:

NIntegrate[fRestriction[{x,y,z}],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}]

So, the tradeoff here is much slower execution (for exact results, anyway) vs. much less work on our part.  Hmmm...  Now, before you get too excited, I really don't know the limitations on this method (i.e., will it work in as many situations as the first method and how does the execution speed depend on the complexity of the integral), so consider this a handy shortcut rather than a real replacement for the first method.

Rework this problem using the piecewise method and compare your answer to what you obtained above.  (Sob quietly into your coffee...)

(Underscript[∫∫∫, D])(x^2)y zdV,
where D is the region defined by (((y+1))^2)>(x^2)+z and (x^2)+(y^2)+(z^2)<1 and y + 2 z > -1.

Integrals and parameters

Sometimes, it is useful to look at how an integral's value depends on a parameter.  So, for example, we could examine how the arc-length of a helix depends on how elliptical it is (execute this command now, before you read further or things won’t work later):

helix[t_]:={Sin[t], k Cos[ t],t}

Manipulate[ParametricPlot3D[{Sin[t], k Cos[ t],t},{t,0,2π}],{k,0,10}]

Warning:  It would seem to make sense here to use helix[t] inside this ParametricPlot3D rather than retyping the function itself (easier to change, as well).  Unfortunately, since this is inside a Manipulate command, there is a bit of a problem with that.  (Try it and you'll see nothing gets graphed.)  I won't go into the details (it's a programming issue; if you're curious, ask me), but the fix is to replace the k with another constant (let's say a) and use Manipulate to change that.  So, this would look like:

Manipulate[ParametricPlot3D[Evaluate[helix[t]/.k->a],{t,0,2π}],{a,0,10}]

You can decide whether that is worth the effort or not...  Let's go ahead and compute the arc length:

arcHelix=(∫_0^(2π))(Sqrt[helix'[t].helix'[t]])\[DifferentialD]t

Since I suspect most of us aren't familiar with this special function, why don't we plot it:

Plot[arcHelix,{k,0,10},AxesLabel->{k,arc-length},AxesOrigin->{0,0}]

If we weren't as concerned about the exact formula for the integral, we could have just directly graphed the integral using NIntegrate:

Plot[NIntegrate[(Sqrt[helix'[t].helix'[t]]),{t,0,2π}],{k,0,10},AxesLabel->{k,arc-length},AxesOrigin->{0,0}]

Another thing we might be interested in examining is how the speed varies over the different curves (compare this graph to the graph of the curve for the same values of a, for example):

Manipulate[Plot[Evaluate[(Sqrt[helix'[t].helix'[t]])/.k->a],{t,0,2π},AxesOrigin->{0,0},PlotRange->{0,10}],{a,0,10}]

Examine how the given integral depends on the parameter k (-5≤k≤5).  Explain this behavior in terms of the regions of integration.  
	(Subscript[(Underscript[(∫_( )^( ))∫, D]),  ]^( ))x+(y^2)dA  
where D is the region between the graphs of (x^2)+2(y^2)==1 and (x^2)+(y^2)==1and above the graph of k x+y==1

Hints:  Don't waste your time trying to compute the integral exactly; skip straight to NIntegrate on this one.  Also, since the parameter is in the definition of your region of integration, you will need to pay close attention to that region.  I suggest using ContourPlot (or even Plot3D) within Manipulate to investigate that dependence.  You will also probably want to use the piecewise approach to integration mentioned earlier.  Graph the value of the integral vs. k.