1.09 Parametric Plotting
TUTORIALS

\(Off[General::spell]; 
Off[General::spell1]; 
Off[Plot::plnr]; 

    Off[ParametricPlot::ppcom]; 


    If[MemberQ[{"Global`", \ "Miscellaneous`SIUnits`", \ 
          "Graphics`Colors`"}, 
\ \ \ \ Context[Gray]], \ 
      Remove["*`Gray"]]; 
<< "Miscellaneous`Units`"; 

    Remove["*`Gray"]; 
<< "Graphics`Colors`"; 

    Graphics`Colors`GosiaGreen = RGBColor[0, \ 0.392187, \ 0]; 

    Clear[Derivative]; \)

ThreeAxes[u,v]

Axes3D::"usage" = 
      "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, \
y, and z running from -a/3 to a, and with axes labels b units beyond the tips \
of the axes.  Axes3D[a] is Axes3D[a, a/8]."; Axes3D[u_, v_] := 
      Graphics3D[{{Blue, Line[{{-(u/3), 0, 0}, {u, 0, 0}}]}, 
          Text["x", {u + v, 0, 0}], {Blue, 
            Line[{{0, -(u/3), 0}, {0, u, 0}}]}, 
          Text["y", {0, u + v, 0}], {Blue, 
            Line[{{0, 0, -(u/3)}, {0, 0, u}}]}, 
          Text["z", {0, 0, u + v}]}]; Axes3D[u_] := Axes3D[u, (u/8)];

CMView  =  {2.7,  1.6,  1.2};

T.1) Parametric plotting for projectile motion

A projectile is fired from a cannon at ground level.  The cannon is inclined at an angle a (here 0 < a < (π/2)) with the horizontal, with a given muzzle speed determined by the explosive charge.  If you neglect air resistance, it's not a big deal to come up with parametric formulas for the position {x[t], y[t]} of the projectile t seconds after firing. 
Here t measures time in seconds from the instant the cannon is fired.
x[t] measures the horizontal distance in feet of the projectile down range from the cannon at time t. 
y[t] measures the height of the projectile in feet at time t.
Gravity acts on the y[t] component of the position but does not act on the x[t] component. Thus
        x''[t]  =  0 and y''[t]  =  -32. 
The muzzle speed, s, splits into horizontal and vertical components as follows: 
        x'[0]  =  s Cos[a] 
        y'[0]  =  s Sin[a].
Saying that the cannon is fired at ground level is to say
       {x[0], y[0]}  ==  0.
See what {x[t], y[t]} look like:

Clear[x, y, t, a, s, Derivative]\(DSolve[{\(\(x'\)'\)[t]\  == \ 0, 
\ \ \ \ \ \ \ \ x[0]\  == \ 0, 

        \ \ \ \ \ \ \ \ \(x'\)[0]\  == \ s\ Cos[a]}, x[t], t]\)

\(DSolve[{\(\(y'\)'\)[t]\  == \ \(-32\), 
\ \ \ \ \ \ \ \ y[0]\  == \ 0, 
        
\ \ \ \ \ \ \ \ \(y'\)[0]\  == \ \ s\ Sin[a]}, y[t], t]\)

For a given muzzle speed s and cannon angle a:

x[t_, s_, a_] = s  t  Cos[a]

y[t_, s_, a_] = -16 (t^2) + s  t  Sin[a]

Now you can go to work.

T.1.a)

Plot the trajectory of the projectile given that the cannon is inclined at an angle (π/6) with the horizontal and with muzzle speed 200 ft/sec.

See when the projectile is back on the ground:

Solve[y[t, 200, (π/6)] == 0, t]

Here comes the plot:

Clear[trajectoryplotter]trajectoryplotter[t_] = {x[t, 200, (π/6)], y[t, 200, (π/6)]}

\(\(ParametricPlot[trajectoryplotter[t], {t, 0, 25\/4}, 
	
      AspectRatio -> 1\/GoldenRatio, 
	
      PlotStyle -> {{Red, Thickness[0.015]}}]; \)\)

There it is, a beautiful parabolic arch.

T.1.b)

Measure the horizontal range of the cannon as a function of muzzle speed and angle of inclination.
For a given muzzle speed, what angle of inclination maximizes the horizontal range of the cannon?

Go back to:

\(Clear[x, y, t, s, a, t]; 
x[t_, s_, a_] = s\ \ t\ \ Cos[a]\)

y[t_, s_, a_] = -16 (t^2) + s  t  Sin[a]

Now find when the projectile hits the ground:

Solve[y[t, s, a] == 0, t]

For a given muzzle speed s and angle a, the horizontal range of the cannon is:

x[(1/16) s Sin[a], s, a]

For a given muzzle speed, this is as big as it can be when the angle a is set so that 
        Cos[a] Sin[a] 
is as big as it can be.
Now examine the derivative of Cos[a] Sin[a] with respect to a:

Expand[(∂_a)(Cos[a] Sin[a]), Trig -> True]

So 
       D[Cos[a] Sin[a], a]  =  0 
for 
       Cos[2 a]  =  0;
which gives 
       2 a  =  (π/2);
this is the same as
        a  =  (π/4).
For a given muzzle speed, you get the greatest horizontal range by setting the firing angle a equal to (π/4). 
Not much of a surprise.

T.1.c)

How do you use the parametric formulas:

Clear[x, y, t, s, a, t]x[t_, s_, a_] = s  t  Cos[a]

y[t_, s_, a_] = -16 (t^2) + s  t  Sin[a]

to explain why, no matter what the muzzle speed and the angle are, the projectile moves on a parabola?

Look at the parametric formulas:

Clear[x, y, t, s, a, t]x[t_, s_, a_] = s  t  Cos[a]

y[t_, s_, a_] = -16 (t^2) + s  t  Sin[a]

If you put x  =  x[t, s, a] and y  =  y[t, s, a] and then you eliminate t by solving for t in terms of y and substituting the result into the formula for x, you get:

\(ExpandAll[
      Eliminate[{x\  == \ x[t, s, a], 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y\  == \ y[t, s, a]}, t]]\)

This tells you that when you fix a muzzle speed s and a firing angle a, then the projectile moves on the curve
       y  =  -((16 (x^2) Sec[a]^2)/(s^2))  +  x Tan[a].
For a given muzzle speed and angle, this curve is a parabola because no powers of x other than (x^2), (x^1), and (x^0) are present.

T.2) Parametric plotting for designing a cam

T.2.a)

A cam is to be built with a smooth shape and such that the following central displacements occur at the given angles with the horizontal:
        (30^∘) ->  2.03
        (85^∘) ->  1.75
        (175^∘) -> 2.27
        (270^∘) ->  1.58
        \(\(315\^∘\) -> 	1.38\)
        
To see what this means, feed in the given information:

Clear[displacement]\(\({displacement[1] = 2.03, 
	displacement[2] = 1.75, 
	
      displacement[3] = 2.07, 
	displacement[4] = 1.58, 
	
      displacement[5] = 1.38}; \)\)

Needs["Miscellaneous`Units`"]

Clear[angle]\({angle[1] = π\ \(Convert[30\ Degree, Radian]\)[\([1]\)], 
	
      angle[2] = π\ \(Convert[85\ Degree, Radian]\)[\([1]\)], 
	
      angle[3] = π\ \(Convert[175\ Degree, Radian]\)[\([1]\)], 
	
      angle[4] = π\ \(Convert[270\ Degree, Radian]\)[\([1]\)], 
	
      angle[5] = π\ \(Convert[315\ Degree, Radian]\)[\([1]\)]}\)

Using circular parmetrization, you get the following table of points the edge of the cam should pass through:

Clear[point, k]\(point[k_] := displacement[k]\ {Cos[angle[k]], Sin[angle[k]]}
\)points = Table[point[k], {k, 1, 5}]

And a plot:

\(\(pointplot\  = 

      ListPlot[points, PlotStyle -> {PointSize[0.04], Red}, 
\ \ \ \ 
        AspectRatio -> Automatic, 
\ \ \ \ AxesOrigin -> {0, 0}, 
\ \ \ \ 
        PlotRange -> {{\(-2.5\), 2.5}, {\(-2.5\), 2.5}}]; \)\)

Your problem is to find a parametric description of a cam whose edge passes through these points.

Answer:

Here's an idea first suggested by John Baptiste Fourier, a French mathematician from the 1800s who pioneered the notion that functions can be made from Sine and Cosine waves the same way that music is made from basic harmonics.
You want to find a function r[t] so that
         r[angle[k]]  =  displacement[k] 
for each k  =  1, 2, 3, 4, and 5.
Here's the data:

Clear[k]data = Table[{angle[k], displacement[k]}, {k, 1, 5}]

Because you have five data points, you have five degrees of freedom; so you use Fourier's idea, incorporating one constant function, two Sine waves, and two Cosine waves:

Clear[r, t]r[t_] = Fit[data, {1, Cos[t], Sin[t], Cos[2 t], Sin[2 t]}, t]

Here is the parametric formula for the curve you are after:

Clear[x, y]{x[t_], y[t_]} = r[t] {Cos[t], Sin[t]}

And a plot:

\(camcurve\  = \ 
\ \ \ 
      ParametricPlot[Evaluate[{x[t], y[t]}], 

        \ \ \ \ \ \ \ \ {t, 0, 2\ \ π}, PlotStyle -> Thickness[0.01], 

        \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> Identity]; 

    \ \ \ \ \ \ \ \ \ \ \ \ 

    Show[pointplot, camcurve, 
\ \ \ \ AspectRatio -> Automatic, 
\ \ \ \ 
      AxesOrigin -> {0, 0}, 
\ \ \ \ 
      PlotRange -> {{\(-2.5\), 2.5}, {\(-2.5\), 2.5}}, 
\ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \)

Damn good.
Fourier was no fool.

T.3) Parametric plotting of the predator-prey model

Remember the predator-prey model from the previous lesson.  In case you've forgotten about it, here's the scoop again:

Two species coexist in a closed environment. One species, the predator, feeds on the other, the prey. There is always plenty of food for the prey, but the predators eat nothing but the prey.
Put 
      pred[t]  =  population of predators at time t.
Put  
      prey[t]  =  population of prey at time t.
It's reasonable to assume that there are positive constants a and b such that:
     prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
because the abundance of food for the prey allows the birth rate of the prey to be proportional to their current number, and the death rate of prey is proportional to both the current number of prey and the current number of predators.

It also makes some sense to assume that there are positive constants c and d such that:
     pred'[t]  =  -c pred[t]  +  d pred[t] prey[t]
because it's reasonable to assume that the death rate of the predators is likely to be proportional to the current population of predators and that the birth rate of the predators is proportional to both the current number of the predators and the size of the food supply (the prey).

Here is what happens for a sample choice of a, b, c, d with the prey starting off with a population of 4 units and the predators starting out with a population of 1 unit.

Here is Mathematica's fake plot of the prey population as a function of time t for the first 50 time units:

\(endtime\  = \ 50; 

a\  = \ 0.7; 
b\  = \ 0.3; 
c\  = \ 0.3; 

    d\  = \ 0.1; 

Clear[pred, prey, t, Derivative]\)\(approxsolutions\  = \ 

      NDSolve[{\(prey'\)[t]\  == \ \ a\ prey[t]\  - \ b\ prey[t]\ pred[t], 

          \ \ \ \ \ \ \ \ \ 
          \(pred'\)[t]\  == \ \(-c\)\ pred[t]\  + \ d\ pred[t]\ prey[t], 

          \ \ \ \ \ \ \ \ \ prey[0]\  == \ 4, 
\ \ \ \ \ \ \ \ \ 
          pred[0]\  == \ 1}, 
\ \ \ {prey[t], pred[t]}, {t, 0, endtime}]; 

    \ \ \ 
Clear[fakepred, fakeprey]\)\(fakepred[t_]\  = \ pred[t] /. approxsolutions[\([1]\)]; 


    fakeprey[t_]\  = \ prey[t] /. approxsolutions[\([1]\)]; 


    preyplot\  = 

      Plot[fakeprey[t], {t, 0, endtime}, 
	\ \ \ 
        AspectRatio -> 1\/GoldenRatio, 
\ \ \ \ \ \ \ 
        PlotStyle -> {{Blue, Thickness[0.015]}}, 
\ \ \ \ \ \ \ 
        AxesLabel -> {"t", "prey"}]; \)

And the predator population as a function of time t:

\(\(predatorplot\  = 

      Plot[fakepred[t], {t, 0, endtime}, 
\ \ \ \ \ \ \ 
        AspectRatio -> 1\/GoldenRatio, 
\ \ \ \ \ \ \ 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
\ \ \ \ \ \ \ 
        AxesLabel -> {"t", "predators"}]; \)\)

Here they are together:

both = Show[predatorplot, preyplot];

Note the relationship between the predator population's crests and dips and prey population's crests and dips.  Just after the predators start growing, the prey is just about eaten up.

T.3.a)

Use parametric plotting to display the predator population as a function of the prey population and describe what you see.

Here it is:

\(\(
\(fullplot\  = 
\ \ \ 
      ParametricPlot[{fakeprey[t], fakepred[t]}, 

        \ \ \ \ \ \ \ \ \ {t, 0, endtime}, PlotStyle -> Thickness[0.01], 

        \ \ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ \ 
        AspectRatio -> 1\/GoldenRatio, 
\ \ \ \ \ \ \ \ \ 
        AxesLabel -> {"prey", "predator"}]; \)\)\)

A distorted ellipse that some folks call a closed curve.
This is a little hard to interpret. 
To get a grasp on it, see what happens as t goes from 0 to 5:

\(\(
\(five\  = \ 
      ParametricPlot[{fakeprey[t], fakepred[t]}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ {t, 0, 5}, PlotStyle -> Thickness[0.01], 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AspectRatio -> 1\/GoldenRatio, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        AxesLabel -> {"prey", "predator"}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, 9}, {0, 5}}]; 
    \)\)\)

Hmm. 
Now see what happens as t goes from 0 to 10:

\(\(ten\  = \ 
      ParametricPlot[{fakeprey[t], fakepred[t]}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ {t, 0, 10}, PlotStyle -> Thickness[0.01], 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AspectRatio -> 1\/GoldenRatio, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        AxesLabel -> {"prey", "predator"}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, 9}, {0, 5}}]; \)\)

Hmmm. 
Now see what happens as t goes from 0 to 15:

\(\(
\(fifteen\  = \ 
      ParametricPlot[{fakeprey[t], fakepred[t]}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ {t, 0, 15}, PlotStyle -> Thickness[0.01], 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AspectRatio -> 1\/GoldenRatio, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        AxesLabel -> {"prey", "predator"}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, 9}, {0, 5}}]; 
    \)\)\)

Holy smoke!
This is the same as:

Show[fullplot, PlotRange -> {{0, 9}, {0, 5}}];

And this plot shows what happens as t goes from 0 to 50.
Review the plots:

Show[GraphicsArray[{five, ten, fifteen}]];

Evidently as time advances the parametric points {prey[t], pred[t]} advance around the closed curve in a counterclockwise way.  As time advances on and on, the parametric points {prey[t], pred[t]} cycle around this curve over and over again. 
Take another look but this time setting the AxesOrigin -> {(c/d), (a/b)}:

Show[fullplot, PlotRange -> {{0, 9}, {0, 5}}, 
      AxesOrigin -> {(c/d), (a/b)}];

Again the closed curve tells you that predator and prey populations are periodic; they go through repeated cycles. And the time cycle for each is the same. The four sectors above defined by the axes indicate four phases depicting trends that reverse themselves as the curve crosses the axes at prey  =  (c/d) and predator  =  (a/b).

T.3.b.i)

Wait a minute!
 Take another look at the last plot.

\(\(Show[fullplot, PlotRange -> {{0, 9}, {0, 5}}, 
	
      AspectRatio -> 1\/GoldenRatio, AxesOrigin -> {c\/d, a\/b}]; \)\)

Those axes pierce the curve at the maximum and minimum values of prey and predator populations.
Is this just an accident?

Stupid question.

T.3.b.ii)

Explain why those axes pierce the plot at the maximum and minimum values of prey and predator populations.

Look at the original differential equations
        prey'[t]  =   a prey[t]  -  b prey[t] pred[t],
        pred'[t]  =  -c pred[t]  +  d pred[t] prey[t].
Take the first:
        prey'[t]  =   a prey[t]  -  b prey[t] pred[t].
At the times t at which the prey population is at its maximum or at its minimum, you gotta have 
        prey'[t]  =  0.
Consequently, at the times t at which the prey population is at its maximum or at its minimum,
        0  =  a prey[t]  -  b prey[t] pred[t],
            =  prey[t] (a   -  b pred[t]).
This tells you that at the times t at which the prey population is at its maximum or at its minimum,
        pred[t]  =  (a/b).
Similarly, use
        pred'[t]  =  -c pred[t]  +  d pred[t] prey[t] 
to see that for a maximum or minimum predator population, you have
        0  =  -c pred[t]  +  d pred[t] prey[t]
            =  pred[t] (-c   +  d  prey[t]),
so when the predator population is largest or smallest,
        prey[t]  =  (c/d).
This explains why when you use AxesOrigin -> {(c/d), (a/b)}, the  axes pierce the plot at the maximum and minimum values of prey and predator populations.

\(\(Show[fullplot, 
	PlotRange -> {{0, 9}, {0, 5}}, 
      AxesOrigin -> {c\/d, a\/b}]; \)\)

Math happens again.

T.4) Quick calculations

T.4.a)

What is the instantaneous growth rate of the parametric curve
     x  =  x[t]  =  (e^t)  +  t
     y  =  y[t]  =  t (3  -  t) Log[t]
at the point at which x  =  5?

Answer:

Take a look:

Clear[x, y, t]\(x[t_]\  = \ E\^t + t; 
y[t_]\  = \ t\ \((3\  - \ t)\)\ Log[t]; 


    curveplot\  = 
\ \ \ \ \ 
      ParametricPlot[{x[t], y[t]}, {t, 0.01, 5}, 
\ \ \ \ \ 
        PlotStyle -> {Thickness[0.01], Red}, 
\ \ \ \ \ 
        AspectRatio -> Automatic, 
\ \ \ \ \ 
        AxesLabel -> {"x", "y"}]; \)

Now find the t for which x[t]  =  5:

Plot[{x[t], 5}, {t, 0, 4}];

We start the search near t  =  1.2:

FindRoot[x[t] == 5, {t, 1.2}]

The instantaneous growth rate at x  =  5 is:

y'[t]/(x')[t] /. t -> 1.30656

T.4.b)

What is the highest point on the parametric curve
     x  =  x[t]  =  (e^t)  +  t
     y  =  y[t]  =  t (3  -  t) Log[t]?

Answer:

Evidently the curve stays negative for larger t's because 
     y[t]  =  t (3  -  t) Log[t] 
is negative for t  ≥  3.
We search for a point t with 0  <  t  ≤  3 at which (dy/dx) =  0.
At {x[t],  y[t]}, we know (dy/dx)is given by ((y'[t])/(x'[t])).
So we have to find the t's with 0  <  t  ≤  3 at which y'[t]  =  0.

Plot[y'[t], {t, 0.01, 3}];

FindRoot[y'[t] == 0, {t, 0.3}]

FindRoot[y'[t] == 0, {t, 2.1}]

Look at:

y[0.32105]

y[2.10316]

The high point:

high = {x[2.10316], y[2.10316]}

is the highest point on the curve.
Check:

\(pointplot = Graphics[RGBColor[1, 0, 0], PointSize[0.02], Point[high]]; 
    
\(Show[curveplot, pointplot]; \)\)

Nailed it.

T.4.c)

If you have x[t] and y[t]as follows,
     x'[t]  =  Cos[t] y[t]  with  x[0]  =  5
     y'[t]  =  4 Sin[t] y[t] with  y[0]  =  9
then what is (dy/dx)when t  =  π?

Answer:

When t  =  π, then
    (dy/dx)=((y'[π])/(x'[π])) 
      
          = 
        ((4 y[π] Sin[π])/((y[π] Cos[π])))
       
          = ((4 Sin[π])/Cos)[π]
      
          = -((4 0)/1) = 0