1.09 Parametric Plotting
GIVE IT A TRY!

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

    Off[ParametricPlot::ppcom]; 
Off[ParametricPlot3D::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]; \)

Axes3D[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};

G.1)  Quick calculations

G.1.a)

Plot the curve {x[t], y[t]} with
     x[t]  =  ((3 t)/(1  +  (t^( 3)))) and  y[t]  =  ((3 (t^( 2)))/(1  +  (t^( 3))))
for 1  ≤  t  ≤  9.
If you regard y to be a function f[x] of x, then which x's give you f'[x]  =  0? 
Which x gives rise to the largest value of f[x]?

G.1.b)

When you go with
        x  =  x[t]  =  r t  -  r Sin[t]
        y  =  y[t]  =  r  -  r Cos[t], 
as in the Basics, you get a cycloid.
Calculating y in the form y  =  f[x] is almost out of the question.
Nevertheless, it's not a big deal to calculate f'[x] when you are given a specific value of x.
Here is a table in the form {x, f'[x]} for selected x's.

Clear[r, t, f]Table[{N[r  t - r  Sin[t]], f'[N[r  t - r  Sin[t]]]}, {t, 
        0, 2 π, (π/6)}]

Unfortunately f'[x] has not been calculated simply because no clean formula for f[x] is available.
Copy, paste, and edit the last instruction so that each second slot exhibits the actual value of f'[x] for the corresponding x in the first slot.

G.1.c)

A curve is given in parametric form by
       x  =  
        x[t]  =  (t^2)  +  5 t  -  (e^(- .3 (t^2)))
       y  =  y[t]  =   3 (t^2)  +  54 t  +  8.
Here's a plot of part of this curve:

Clear[x, y, t]\(x[t_]\  = \ t\^2 + 5\ t - E\^\(\(-0.3\)\ t\^2\); 

    y[t_]\  = \ 3\ t\^2 + 54\ t + 8; 


    ParametricPlot[{x[t], y[t]}, {t, \(-2\), 2}, 
\ \ \ \ 
      AspectRatio -> 1\/GoldenRatio, 
\ \ \ \ 
      PlotStyle -> {{Thickness[0.01], Blue}}, 
\ \ \ \ 
      AxesLabel -> {"x", "y"}]; \)

Plot more and more of the curve until you get a pretty good idea of how this baby looks.
Use information from x'[t] to give a good estimate of the smallest value of x on the curve.
Use information from y'[t] to give a good estimate of the smallest value of y on the curve.

G.2)  Parametric plotting of circles and ellipses in two dimensions*

A handy way to plot the circle
       (x^2)  +  (y^2)  =  1 
is to write
       x[t]  =  Cos[t]
and
       y[t]  =  Sin[t]
and then to plot the points
     {x[t],  y[t]}  =  {Cos[t], Sin[t]} 
as t advances from 0 to 2  π.

Clear[x, y, t]\({x[t_], y[t_]}\  = \ {Cos[t], \ Sin[t]}; 


    circle\  = \ 
\ \ \ \ \ 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2\ Pi}, 
\ \ \ \ \ \ \ \ 
        PlotStyle -> {{Blue, Thickness[0.015]}}, 
\ \ \ \ \ \ \ \ 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}]; \)

See what happens when you multiply x[t]  =  Cos[t] by 3:

Clear[x, y, t]\({x[t_], y[t_]}\  = \ {3\ Cos[t], \ Sin[t]}; 


    ParametricPlot[{x[t], y[t]}\ , {t, 0, 2\ Pi}, 
\ \ \ \ \ \ \ \ 
      PlotStyle -> {{Blue, Thickness[0.01]}}, 
\ \ \ \ \ \ \ \ 
      AspectRatio -> Automatic, AxesLabel -> {"x", "y"}]; \)

The plot looks quite a bit like the ellipse
       ((((x/3)))^2)  +  (y^2)  =  1.
Check this out:

Clear[x, y, t]Expand[((((x/3)))^2) + (y^2) /. {x -> 3 Cos[t], y -> Sin[t]}, 
      Trig -> True]

The plot IS the ellipse  
       ((((x/3)))^2)  +  (y^2)  =  1.
The upshot: 
       {x[t], y[t]}  =  {3 Cos[t],  Sin[t]} 
gives you a parametric formula for the ellipse  
       ((((x/3)))^2)  +  (y^2)  =  1.

Now see what happens when you also multiply y[t]  =  Sin[t] by 2:

Clear[x, y, t]\({x[t_], y[t_]}\  = \ {3\ Cos[t], \ 2\ Sin[t]}; 


    ParametricPlot[{x[t], y[t]}\ , {t, 0, 2\ Pi}, 
\ \ \ \ \ \ \ \ 
      PlotStyle -> {{Blue, Thickness[0.01]}}, 
\ \ \ \ \ \ \ \ 
      AspectRatio -> Automatic, AxesLabel -> {"x", "y"}]; \)

The plot looks quite a bit like the ellipse
       ((((x/3)))^2) + ((((y/2)))^2) = 1.
Check this out:

Clear[x, y, t]Expand[
      ((((x/3)))^2) + ((((y/2)))^2) /. {x -> 3 Cos[t], y -> 2 Sin[t]}, 
      Trig -> True]

The plot IS the ellipse  
      ((((x/3)))^2) + ((((y/2)))^2) = 1.
The upshot: 
      {x[t], y[t]}  =  {3 Cos[t],  2 Sin[t]} 
 gives you a parametric formula for the ellipse  
      ((((x/3)))^2) + ((((y/2)))^2) = 1.

G.2.a.i)

Give, as above, parametric formulas for the ellipse
      ((((x/4)))^2) + ((((y/3)))^2) = 1
and use your formulas to plot the curve.

G.2.a.ii)

Give, as above, parametric formulas for the ellipse
       ((((x/2)))^2) + ((((y/4)))^2) = 1
and use your formulas to plot the curve.

G.2.a.iii)

Given positive constants a and b, give as above parametric formulas for the ellipse
       ((((x/a)))^2) + ((((y/b)))^2) = 1.

G.2.a.iv)

Look at this graphic:

\(a\  = \ 5; 
b\  = \ 3; 

    Clear[acircleplotter, bcircleplotter, ellipseplotter, t]\)\(acircleplotter[t_]\  = \ {a\ Cos[t], a\ Sin[t]}; 

    ellipseplotter[t_]\  = \ {a\ Cos[t], b\ Sin[t]}; 

    bcircleplotter[t_]\  = \ {b\ Cos[t], b\ Sin[t]}; 


    CBS\  = \ 
\ \ \ 
      ParametricPlot[

        \ \ \ {acircleplotter[t], 
\ \ \ \ \ \ \ ellipseplotter[t], 

          \ \ \ \ \ \ \ \ \ \ \ \ bcircleplotter[t]}, 
\ \ \ {t, 0, 2\ Pi}, 
        
\ \ \ PlotStyle -> {{Blue, Thickness[0.02]}, 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Red, Thickness[0.02]}, 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Thickness[0.02]}}, 
\ \ \ \ 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}]; 
    \ \ \ \ \ \ \ \ \)

Here are the points plotted on each curve when you go with t  =  (π/6):

Clear[points, ray, t]\(points[t_]\  = \ 
      Graphics[{PointSize[0.05], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Point[{a\ Cos[t], a\ Sin[t]}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Point[{a\ Cos[t], b\ Sin[t]}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Point[{b\ Cos[t], b\ Sin[t]}]}]; 

    ray[t_]\  = \ 
      Graphics[
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        Line[{{0, 0}, {a\ Cos[t], a\ Sin[t]}}]]; 

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

    Show[CBS, points[π\/6], ray[π\/6]]; \)

Here is how the two circles and the ellipse plot out as t runs from 0 to (π/2):

Table[Show[CBS, points[t], ray[t]], {t, 0, (π/2), (π/12)}];

Your friend has no machine, but has a ruler, a compass, a T-square, and a triangle. Tell your friend how to plot a lot of points on the ellipse 
         ((((x/4)))^2) + ((((y/2)))^2) = 1
by drawing concentric circles and then moving horizontally or vertically from well-chosen points on the circles.

G.2.b)

A handy way to plot the circle of radius 1 centered at {2, 3}
       (((x - 2))^2) + (((y - 3))^2) = 1 
is to use:

Clear[x, y, t]\({x[t_], y[t_]}\  = \ {2, 3}\  + \ {Cos[t], Sin[t]}; 


    ParametricPlot[{x[t], y[t]}, {t, 0, 2\ Pi}, 
\ \ \ \ \ \ \ \ 
      PlotStyle -> {{Blue, Thickness[0.015]}}, 
\ \ \ \ \ \ \ \ 
      AspectRatio -> Automatic, AxesOrigin -> {2, 3}, 
\ \ \ \ \ \ \ \ 
      AxesLabel -> {"x", "y"}, 
\ \ \ \ \ \ \ \ 
      Epilog -> {Red, PointSize[0.06], Point[{2, 3}]}]; \)

Set up parametric formulas for the circle of radius r centered at {h, k}:
       (((x  -  h))^2)  +  (((y  -  k))^2)  =  (r^2).
Use your formulas to plot the circle  
       (((x  -  1))^2)  +  (((y  -  2.5))^2)  =  5.

G.2.c)

A handy way to plot the ellipse
       (((((x  -  2)/5)))^2) + (((((y  -  4)/2)))^2) = 1
centered on {2, 4} is to use:

Clear[x, y, t]\({x[t_], y[t_]}\  = \ {2, 4}\  + \ {5\ Cos[t], 2\ Sin[t]}; 


    ParametricPlot[{x[t], y[t]}, {t, 0, 2\ Pi}, 
\ \ \ \ \ \ \ \ 
      PlotStyle -> {{Blue, Thickness[0.015]}}, 
\ \ \ \ \ \ \ \ 
      AspectRatio -> Automatic, AxesOrigin -> {2, 4}, 
\ \ \ \ \ \ \ \ 
      AxesLabel -> {"x", "y"}, 
\ \ \ \ \ \ \ \ 
      Epilog -> {Red, PointSize[0.05], Point[{2, 4}]}]; \)

Use parametric formulas as above to plot the ellipse
        (((((x  +  1)/2)))^2) + (((((y  -  4)/5)))^2) = 1.
At what point is this ellipse centered?

G.3)  Elliptical orbits of planets and asteroids

G.3.a) Elliptical parameters

You can plot the circle (x^2)  +  (y^2)  =  (r^2) using the parametric equations
        x[t]  =  r Cos[t]
        y[t]  =  r Sin[t]
and running t from 0 to 2  π.
You can plot the ellipse 
       ((((x/a)))^2) + ((((y/b)))^2) = (r^2) 
using the parametric equations
        x[t]  =  a r Cos[t]
        y[t]  =  b r Sin[t]
and running t from 0 to 2  π.
Here you go with a true scale plot of the ellipse
       ((((x/9)))^2) + ((((y/3)))^2) = 1:

Clear[x, y, t]\(x[t_]\  = \ 9\ Cos[t]; 
y[t_]\  = \ 3\ Sin[t]; 


    ParametricPlot[{x[t], y[t]}, {t, 0, 2\ Pi}, 
\ \ \ \ \ 
      AspectRatio -> Automatic, 
\ \ \ \ \ 
      PlotStyle -> {Blue, Thickness[0.015]}, 
\ \ \ \ \ 
      AxesLabel -> {"x", "y"}]; \)

Give a true scale plot of the ellipse
       ((((x/4)))^2) + ((((y/9)))^2) = 1.

G.3.b) Eccentricity

Here is the ellipse ((((x/4)))^2) + ((((y/3)))^2) = 1 in true scale:

Clear[x, y, t]\(x[t_]\  = \ 4\ Cos[t]; 
y[t_]\  = \ 3\ Sin[t]; 


    ParametricPlot[{x[t], y[t]}, {t, 0, 2\ Pi}, 
\ \ \ \ \ 
      AspectRatio -> Automatic, 
\ \ \ \ \ 
      PlotStyle -> {Blue, Thickness[0.01]}, 
\ \ \ \ \ 
      AxesLabel -> {"x", "y"}]; \)

The long axis of this ellipse runs along the x - axis; its length is 8. The short axis of this ellipse runs along the y - axis; its length is 6.

To get what folks call the eccentricity of an ellipse, you take the ratio 
       ratio  =  
        ((length  of  short  axis)/(length  of  long  axis
            ))
and put 
       eccentricity  =  (Sqrt[1 - (ratio^2)])  . 
Here is the eccentricity of the ellipse above:

long = 8; short = 6; ratio = (short/long); eccentricity = N[(Sqrt[1 - (ratio^2)])]

G.3.b.i)

Plot the ellipse 
        ((((x/8)))^2) + ((((y/3)))^2) = 1 
in true scale and calculate its eccentricity.

G.3.b.ii)

Plot the ellipse 
       ((((x/2)))^2) + ((((y/5)))^2) = 1 
in true scale and calculate its eccentricity.

Remember:
       ratio  =  
        (length of short axis)/(length of long axis)
and 
       eccentricity  =  (Sqrt[1 - (ratio^2)]) .

G.3.b.iii)

Given positive numbers a and b, how do you know that the eccentricity of the ellipse 
       ((((x/a)))^2) + ((((y/b)))^2) = 1 
is a number no smaller than 0 but less than 1?
What do you call an ellipse whose eccentricity is 0?

G.3.b.iv)

Plot in true scale an ellipse whose long axis is 10 units long and whose eccentricity is 0.9.

Put a  =  5 and b  =   (short/2).
      
       ratio  =  
        ((length  of  short  axis)/(length of  long  axis
            )) 
                   =  ((2 b)/(2 a)) = (b/a)
      
       eccentricity  =  (Sqrt[1 - (ratio^2)]) 
So 
       0.9  =  (Sqrt[1 - ((((b/a)))^2)]))  =  
        (Sqrt[1 - ((((b/5)))^2)]).

G.3.b.v)

Plot in true scale an ellipse whose long axis is 10 units and whose eccentricity is 0.1.

G.3.c) Focuses

If a  ≥  b, then the focuses of the ellipse
       ((((x/a)))^2) + ((((y/b)))^2) = 1 
are located at  the points 
       {-eccentricity a, 0} and {eccentricity a, 0}.
But if a  <  b, then the focuses of the ellipse
       ((((x/a)))^2) + ((((y/b)))^2) = 1 
are located at  the points 
       {0, -eccentricity b} and {0, eccentricity b}. 
Here is a true scale plot of the ellipse 
       ((((x/6)))^2) + ((((y/3)))^2) = 1 
shown along with its two focuses:

Clear[x, y, a, b, t]\(x[t_]\  = \ a\ \ Cos[t]; 
y[t_]\  = \ b\ \ Sin[t]; 
a\  = \ 6; 

    b\  = \ 3; 
eccentricity\  = \ Sqrt[1\  - \((b\/a)\)\^2]; 

    ellipse\  = 

      ParametricPlot[{x[t], y[t]}, 
\ \ \ \ \ {t, 0, 2\ Pi}, 
        PlotStyle -> {Blue, Thickness[0.01]}, 
\ \ \ \ \ 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 

        \ \ \ \ \ 
        Epilog -> 
          \ {{Red, PointSize[0.04], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
              Point[{\(-eccentricity\)\ a, 0}]}, 

            \ \ \ \ \ \ \ \ \ \ \ \ {Red, PointSize[0.04], 

              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[{eccentricity\ a, 0}]}}]; 
    \)

Get it?

Plot the ellipse 
       ((((x/2)))^2) + ((((y/9)))^2) = 1  
and show it along with its two focuses.

If your focuses did not plot out inside the ellipse, then you screwed up.

G.3.d.i) Planets, asteroids and comets

Johannes Kepler (1571-1630) earned a permanent place in history by perfecting Copernicus's idea that the planets revolve about the sun.  Greatly influenced by Copernicus and getting his start by working with other giants like Brahe and Galileo, he determined that planets, asteroids, and comets move in elliptical orbits around the sun with the sun at one of the focuses of the elliptical orbit.

Kepler's first law tells you that to give high-quality plots of the orbits of planets, asteroids, and comets, you are well advised to plot with the sun at a focus.
Here is the usual plot of 
          ((((x/7)))^2) + ((((y/3)))^2) = 1 
with the center of the ellipse at the origin:

Clear[x, y, t]\(x[t_]\  = \ a\ \ Cos[t]; 
y[t_]\  = \ b\ \ Sin[t]; 
a\  = \ 7; 

    b\  = \ 3; 
eccentricity\  = \ Sqrt[1\  - \ \((b\/a)\)\^2]; 

    ellipse\  = 

      ParametricPlot[{x[t], y[t]}, 
\ \ \ \ \ {t, 0, 2\ Pi}, 
        PlotStyle -> {Blue, Thickness[0.01]}, 
\ \ \ \ \ \ 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 

        \ \ \ \ \ \ 
        Epilog -> 
          \ {{Red, PointSize[0.04], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
              Point[{\(-eccentricity\)\ a, 0}]}, 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Red, PointSize[0.04], 

              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[{eccentricity\ a, 0}]}}]; 
    \)

Here is how you plot the same curve shifted to the left so that the right focus is at the origin:

Clear[x, y, t, a, eccentricity]\(x[t_]\  = \ a\ Cos[t]\  - \ a\ eccentricity; 
y[t_]\  = \ b\ Sin[t]; 

    
a\  = \ 7; 
b\  = \ 3; 

    eccentricity\  = \ Sqrt[1\  - \ \((b\/a)\)\^2]; 


    ellipse\  = 

      ParametricPlot[{x[t], y[t]}, 
\ \ \ \ \ {t, 0, 2\ Pi}, 
        PlotStyle -> {Blue, Thickness[0.01]}, 
\ \ \ \ \ 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 

        \ \ \ \ \ 
        Epilog -> 

          \ \ \ \ \ \ \ \ \ {{Red, PointSize[0.04], 

              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
              Point[{\(-2\)\ eccentricity\ a\ , 0}]}, 

            \ \ \ \ \ \ \ \ \ \ {Red, PointSize[0.04], 

              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[{0, 0}]}}]; \)

Get it?
Plot the ellipse 
          ((((x/2)))^2) + ((((y/6)))^2) = 1. 
Then plot the same curve adjusted so that the top focus is at the origin. Show the focuses in both plots.

G.3.d.ii)

Plot an ellipse whose long axis is 20 units long and whose eccentricity is 0.65 with its center at the origin.  Then plot the same curve adjusted so that the right focus is at the origin. Show the focuses in both plots.

G.3.d.iii)

Kepler's first law says that planets, asteroids, and comets move in elliptical orbits around the sun with the sun at one of the focuses of the elliptical orbit.
One big mystery is that all these orbits lie roughly in the same plane.
This mystery is intriguing, but the fact that all these orbits lie in one plane makes it easy to give reasonable plots of the orbits of the planets, asteroids, and comets.
Here are the data for all the planets in the form
{planet name, length of long axis of elliptical orbit, eccentricity}

{Mercury, 0.774, 0.206}          {Venus, 1.446, 0.007}{Earth, 2.00, 0.017}               
       {Mars, 3.04, 0.093}{Jupiter, 10.4, 0.048}           
       {Saturn, 19.08, 0.056}{Uranus, 38.36, 0.05}             
       {Neptune, 60.14, 0.01}{Pluto, 78.88, 0.25}

All but Mercury and Pluto have small eccentricities; so their orbits are nearly circular. No fun plotting these. 
Plot the orbit of Mercury with the sun at one focus.

G.3.e.iv)

To get a better kick out of this, give a plot that shows the orbits of Neptune and Pluto.
Is Pluto always the outmost planet?

The eccentricity of Neptune's orbit is so small that you can safely plot it as a circle of radius 30.07 centered at the Sun.

G.3.e.v)

To get an even better kick out of this, give a plot that shows the orbits of Earth and a typical member of the Apollo asteroid group. 
In the form of the data above the Apollo data are {Apollo, 3.0, 0.56}.

The eccentricity of Earth's orbit is so small that you can safely plot it as a circle of radius 1.00 centered at the Sun.

G.4)  Parametric plotting of circles, tubes and horns in three dimensions*

Here are the three coordinate axes in three dimensions:

\(\(h\  = \ 5; 
spacer\  = \ h\/10; 

    threedims\  = \ 
\ \ \ \ \ \ 
      Graphics3D[{

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{\(-h\), 0, 0}, {h, 0, 0}}]}, 
          
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["x", {h\  + \ spacer, 0, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, \(-h\), 0}, {0, h, 0}}]}, 
          
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["y", {0, \ h\  + \ spacer, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, 0, \(-h\)}, {0, 0, h}}]}, 
          
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["z", {0, 0, h\  + \ spacer}]}]; \ 

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

    CMView\  = \ {2.7, \ 1.6, \ 1.2}; 

    Show[threedims, PlotRange -> All, 
ViewPoint -> CMView, Boxed -> False]; 
    \ \ \ \ \ \ \ \ \)\)

The xy-plane is the plane that you get by laying down a rigid flat sheet on the girders defined by the x-axis and the y-axis.
Here is a piece of it:

\(\(xyplane\  = 

      Graphics3D[
\ \ \ 
        Polygon[{{0, 0, 0}, {h, 0, 0}, {h, h, 0}, {0, h, 0}}]]; 


    Show[threedims, xyplane, PlotRange -> All, 
\ \ \ \ \ \ 
      ViewPoint -> CMView, 
\ \ \ \ \ \ Boxed -> False]; \ \ \ \ \ \ \ \ 
    \)\)

The xz-plane is the plane that you get by laying down a rigid flat sheet on the girders defined by the x-axis and the z-axis.
Here is a piece of it:

\(\(xzplane\  = \ 
      Graphics3D[
	Polygon[{{0, 0, 0}, {0, 0, h}, {h, 0, h}, {h, 0, 0}}]]; 
    

Show[threedims, xzplane, PlotRange -> All, 
\ \ ViewPoint -> CMView, 
      Boxed -> False]; \ \ \ \ \ \ \ \ \)\)

Show a piece of the yz-plane together with the coordinate axes.
Continue to use the option ViewPoint -> CMView.

G.4.a)

Here is the circle in the xy-plane of radius 1.5 centered at {2.5, 3, 0}:

Clear[t]\(circle\  = \ 
      ParametricPlot3D[
\ \ \ \ \ \ 
        Evaluate[{2.5, 3, 0}\  + \ 1.5\ {Cos[t], Sin[t], 0}], 

        \ \ \ \ \ {t, 0, 2\ Pi}, DisplayFunction -> Identity]; 

h\  = \ 5; 
    
xyplane\  = 

      Graphics3D[
\ \ \ \ \ 
        Polygon[{{0, 0, 0}, {h, 0, 0}, {h, h, 0}, {0, h, 0}}]]; 


    Show[threedims, xyplane, circle, PlotRange -> All, 
\ \ \ \ \ \ 
      ViewPoint -> CMView, Boxed -> False, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Here is the circle parallel to the xz-plane of radius 1.5 centered at {2.5, 3, 2}:

Clear[t]\(circle\  = \ 
      ParametricPlot3D[
\ \ \ \ \ \ 
        Evaluate[{2.5, 3, 2}\  + \ 1.5\ {Cos[t], 0, Sin[t]}], 

        \ \ \ \ \ \ \ {t, 0, 2\ Pi}, DisplayFunction -> Identity]; 

    h\  = \ 5; 

    xzplane\  = \ 
      Graphics3D[
	Polygon[{{0, 0, 0}, {0, 0, h}, {h, 0, h}, {h, 0, 0}}]]; 
    

Show[threedims, xzplane, circle, PlotRange -> All, 
\ \ \ \ \ \ 
      ViewPoint -> CMView, Boxed -> False, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Plot the circle parallel to the yz-plane of radius 1.5 centered at {2, 3, 2}.

G.4.b)

Here are several circles of radius 3 parallel to the xz-plane centered on the y-axis:

Clear[t]\(circle1\  = \ 
      ParametricPlot3D[
Evaluate[{0, 1, 0}\  + \ 3\ {Cos[t], 0, Sin[t]}], 

        \ \ \ \ \ \ \ \ \ \ \ {t, 0, 2\ Pi}, DisplayFunction -> Identity]; 

    
circle2\  = \ 
      ParametricPlot3D[
Evaluate[{0, 2, 0}\  + \ 3\ {Cos[t], 0, Sin[t]}], 

        \ \ \ \ \ \ \ \ \ \ \ {t, 0, 2\ Pi}, DisplayFunction -> Identity]; 

    
circle3\  = \ 
      ParametricPlot3D[
Evaluate[{0, 3, 0}\  + \ 3\ {Cos[t], 0, Sin[t]}], 

        \ \ \ \ \ \ \ \ \ \ \ {t, 0, 2\ Pi}, DisplayFunction -> Identity]; 

    
Show[threedims, circle1, circle2, circle3, 
\ \ \ \ \ \ 
      ViewPoint -> CMView, Boxed -> False, 
\ \ \ \ \ \ PlotRange -> All, 

      \ \ \ \ \ \ \ \ DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ 
    \)

It's no fun to settle for a few crummy circles when you can get a whole tube composed of circles of radius 3 all parallel to the xz-plane and all centered on the y-axis:

Clear[t, y]\(tube\  = \ 
      ParametricPlot3D[
Evaluate[{0, y, 0}\  + \ 3\ {Cos[t], 0, Sin[t]}], 

        \ \ \ \ \ \ \ \ \ \ \ {t, 0, 2\ Pi}, {y, 0, 3}, 
\ \ \ \ \ \ \ \ \ 
        DisplayFunction -> Identity]; 

    Show[threedims, tube, PlotRange -> All, 
\ \ \ \ \ \ 
      ViewPoint -> CMView, Boxed -> False, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Plot a tube composed of circles of radius 3.5 all parallel to the xy-plane and all centered on the z-axis.

G.4.c)

Here is a modest curve in three dimensions:

Clear[t, curveplotter]\(curveplotter[t_]\  = \ {Cos[t], 3\ t\ Sin[t], 3\ t}; 


    curve\  = \ 
      ParametricPlot3D[Evaluate[curveplotter[t]], {t, 0, 2}, \ 
        DisplayFunction -> Identity]; 


    Show[threedims, curve, ViewPoint -> CMView, 
\ \ \ \ \ \ \ 
      Boxed -> False, 
\ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Here is what you get when you make a tube consisting of all circles parallel to the xy - plane of radius 2 centered at all points on this curve:

Clear[s, t]\(tube\  = \ 
      ParametricPlot3D[
\ \ \ \ \ 
        Evaluate[curveplotter[t]\  + \ 2\ {Cos[s], Sin[s], 0}], 

        \ \ \ \ \ \ \ {t, 0, 2}, {s, 0, 2\ Pi}, \ 
\ \ \ \ \ \ \ 
        DisplayFunction -> Identity]; 


    Show[threedims, tube, ViewPoint -> CMView, 
\ \ \ \ \ \ \ \ 
      Boxed -> False, 
\ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Here is another curve:

Clear[t, curveplotter]\(curveplotter[t_]\  = \ {t\ Sin[t], 2\ t, Cos[2\ t]}; 


    curve\  = \ 
      ParametricPlot3D[Evaluate[curveplotter[t]], 

        \ \ \ \ \ \ \ \ \ \ {t, 0, 4}, \ DisplayFunction -> Identity]; 


    Show[threedims, curve, ViewPoint -> CMView, 
\ \ \ \ \ \ \ 
      Boxed -> False, 
\ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Plot the tube composed of all circles of radius 1.5 parallel to the xz - plane and centered on this curve.

G.4.d)

Look at this plot:

Clear[s, t, curveplotter, rad]\(curveplotter[t_]\  = \ {t\ Sin[t], 2\ t, Cos[2\ t]}; 


    rad[t_]\  = \ t\^2\/6; 


    horn\  = \ 
      ParametricPlot3D[
\ \ \ \ \ 
        Evaluate[
          curveplotter[t]\  + \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            rad[t]\ {Cos[s], 0, Sin[s]}], 
{t, 0, 4}, {s, 0, 2\ Pi}, \ 
        DisplayFunction -> Identity]; 


    Show[threedims, horn, ViewPoint -> CMView, 
\ \ \ \ \ Boxed -> False, 

      \ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Describe the circles that this horn is composed of.

G.4.e)

Do something that's either artistic or crazy.

G.5)  Surfaces you can make by rotating curves

Here is a curve in two dimensions:

Clear[f, x]f[x_] = ((Cos[π x] + 2)/(Sqrt[x + 1])); Plot[f[x], {x, 0, 3}, PlotStyle -> Blue, 
      AxesLabel -> {"x", ""}];

You can look at the same curve in three dimensions as a plot in the xz - plane:

Tip:

We'll use a shorthand way of getting the threedims axes:

?Axes3D

Clear[x, curveplotter]\(curveplotter[x_]\  = \ {x, 0, f[x]}; 


    curve\  = \ 
      ParametricPlot3D[curveplotter[x], {x, 0, 3}, \ 
        DisplayFunction -> Identity]; 
threedims\  = \ Axes3D[4]; 

    curveplot\  = 
\ \ 
      Show[threedims, curve, PlotRange -> All, 
\ \ \ \ \ \ \ 
        ViewPoint -> CMView, Boxed -> False, 
\ \ \ \ \ \ \ 
        DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

When you rotate this curve about the x - axis, you get the surface composed of all circles parallel to the yz - plane centered at {x, 0, 0} with radius f[x]:

Clear[s, x, surfaceplotter]\(surfaceplotter[s_, x_]\  = \ {x, 0, 0}\  + \ f[x] {0, Cos[s], Sin[s]}; 
    

surface\  = \ 
      ParametricPlot3D[
\ \ \ \ \ \ \ Evaluate[surfaceplotter[s, x]], 

        \ \ \ \ \ \ \ {x, 0, 3}, {s, 0, 2\ Pi}, \ 
\ \ \ \ \ \ \ 
        DisplayFunction -> Identity]; 


    Show[threedims, surface, ViewPoint -> CMView, 
\ \ \ \ \ \ \ \ 
      Boxed -> False, 
\ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Take another look at the curve:

Show[curveplot];

When you rotate the curve about the z - axis, you get all circles parallel to the xy - plane of radius x and centered at the point {0, 0, f[x]} on the z - axis:

Clear[s, x, surfaceplotter]\(surfaceplotter[s_, x_]\  = \ 

      \ \ \ \ \ \ {0, 0, f[x]}\  + \ x\ {Cos[s], \ Sin[s], 0}; \ 

    surface\  = \ 
      ParametricPlot3D[
\ \ \ \ \ \ \ 
        Evaluate[surfaceplotter[s, x]], {x, 0, 3}, {s, 0, 2\ Pi}, \ 

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


    Show[threedims, surface, ViewPoint -> CMView, 
\ \ \ \ \ \ \ \ 
      Boxed -> False, 
\ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Beautiful tool, this Mathematica is.

G.5.a)

Here's another curve in two dimensions:

Clear[f, y]\(f[y_] = E\^\(\(-0.5\)\ y\)\ \((Sin[π\ y] + 1)\); 


    Plot[f[y], {y, 0, 3.5}, PlotStyle -> Blue, 
      AxesLabel -> {"y", ""}]; \)

You can look at the same curve in three dimensions as a plot in the yz - plane:

Clear[y, curveplotter]\(curveplotter[y_]\  = \ {0, y, f[y]}; 


    curve\  = \ 
      ParametricPlot3D[curveplotter[y], {y, 0, 3.5}, \ 
        DisplayFunction -> Identity]; 
threedims\  = \ Axes3D[4]; 


    curveplot\  = 
\ \ 
      Show[threedims, curve, PlotRange -> All, 
\ \ \ \ \ \ \ 
        ViewPoint -> CMView, Boxed -> False, 
\ \ \ \ \ \ \ 
        DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \)

Plot the surface you get when you rotate this curve about the y - axis.
Follow it up with a plot of what you get when you rotate this curve about the z - axis.
Clearly identify each plot.

G.6)  Projectile marksmanship

G.6.a) Off the wall

On flat terrain, a projectile with an initial velocity of 190 ft/sec is fired at a 20 - foot high vertical wall 1000 feet from the tip of the cannon inclined at an angle a with the horizontal.

G.6.a.i)

If a = (π/6),  does the projectile hit the wall?

The projectile falls short of the wall; you give the calculation.
For some help, see the Tutorials.

G.6.a.ii)

For approximately what angles of inclination will the projectile hit the wall or soar over the wall?

Plot the horizontal range as a function of a.
For more on the horizontal range, see the Tutorials.

G.6.a.iii)

For approximately what angles of inclination will the projectile hit the wall?

G.6.b.i) Trick shots

You have a cannon with a muzzle speed of 190 ft/sec.
At time t  =  0, you fire the cannon with an angle of inclination 
      (π/4)  +  z.
Then you wait d seconds and fire again with an angle of inclination
      (π/4)  -  z.
Here is a plot which reveals what the second projectile is doing while the first is in flight in the case z  =  (π/8) and d  =  1.

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

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


    z\  = \ π\/8; 
d\  = \ 1; 


    ParametricPlot[

      \ \ \ \ {{x[t, 190, π\/4 + z], y[t, 190, π\/4 + z]}, {
          x[t - d, 190, π\/4 - z], y[t - d, 190, π\/4 - z]}}, 

      \ \ \ \ \ \ \ \ \ \ \ {t, 0, 190\/16\ Sin[π\/4 + z]}, 

      \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {\(-100\), 500}, 
\ \ \ \ \ \ 
      PlotStyle -> {{Thickness[0.015], Blue}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Thickness[0.010], Red}}]; \)

That its plot is below the ground on the right indicates that the second shot hit the ground before the first shot hit the ground.
That its plot is below the ground on the left indicates that the second shot was fired after the first shot was fired. 
Your mission, if you accept it, is to find a new delay time d so that both projectiles hit the ground at the same instant.

G.6.b.ii)

You have a cannon with a muzzle speed of s ft/sec.
At time t  =  0, you fire a projectile with an angle of inclination 
      (π/4)  +  z.
Then you wait d[z] seconds and fire another projectile with an angle of inclination
      (π/4)  -  z.
This time your mission is:
Determine a formula for the delay time
        d  =  d[s, z] (for 0  <  z  <  (π/4))  
as a function of z and s so that both projectiles hit the ground at the same spot and at the same instant.
Show off your formula with a couple of sharp plots.

You want:

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


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


    equation1 = y[t, s, π\/4 + z] == 0\)

\(Clear[d]; 
equation2 = y[t - d, s, π\/4 - z] == 0\)

Solve these two equations for t and d.

Solve[{equation1, equation2}, {t, d}]

Now it's up to you to decide which of these to go with.

G.7)  More cams

G.7.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^∘)    ->   3.82
        (85^∘)    ->   5.00
        (175^∘)  ->   1.87
        (270^∘)  ->   5.00
        (315^∘)  ->   4.36
        
To see what this means, feed in the given information:

Clear[displacement]\({displacement[1]\  = \ 3.82, 
\ displacement[2]\  = \ 5.00, 
\ 
      displacement[3]\  = \ 1.87, 
\ displacement[4]\  = \ 5.00, 
\ 
      displacement[5]\  = \ 4.36}\)

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 parametrization, 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}]; \)\)

Your problem is to find a parametric description of a cam whose edge passes through these points.
Show off your work by plotting your cam together with the given points.
Are you totally satisfied with the results you got?

G.7.b)

A new 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\^▫		\  -> 			3.82\)
        \(85\^▫		\  -> 		\ \ 5.00\)
        (125^▫)    ->       2.05
        \(175\^▫	\  -> 		\ \ 1.87\)
        (215^▫)    ->       6.02
        (265^▫)    ->       1.86
        \(335\^▫	\  -> 		\ \ 4.36\)
        
To see what this means, feed in the given information:

Clear[displacement]\({displacement[1]\  = \ 3.82, 
\ displacement[2]\  = \ 5.00, 
\ 
      displacement[3]\  = \ 2.05, 
\ displacement[4]\  = \ 1.87, 
\ 
      displacement[5]\  = \ 6.02, 
\ displacement[6]\  = \ 0.86, 
\ 
      displacement[7]\  = \ 0.36}\)

Clear[angle]\({angle[1]\  = \ π\ \ \(Convert[30\ Degree, Radian]\)[\([1]\)], 
\ 
      angle[2]\  = \ \ π\ \ \(Convert[85\ Degree, Radian]\)[\([1]\)], 
\ 
      angle[3]\  = \ π\ \ \(Convert[125\ Degree, Radian]\)[\([1]\)], 
\ 
      angle[4]\  = \ π\ \ \(Convert[175\ Degree, Radian]\)[\([1]\)], 
\ 
      angle[5]\  = \ π\ \ \(Convert[215\ Degree, Radian]\)[\([1]\)], 
\ 
      angle[6]\  = \ π\ \ \(Convert[265\ Degree, Radian]\)[\([1]\)], 
\ 
      angle[7]\  = \ π\ \ \(Convert[335\ Degree, Radian]\)[\([1]\)]}\)

Using circular parametrization, 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, 7}]

And a plot:

\(\(pointplot\  = 

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

Your problem is to find a parametric description of a cam whose edge passes through these points.
Show off your work by plotting your cam together with the given points.
Are you totally satisfied with the results you got this time?
Does the positioning of the data have something to do with the outcome?

Tip:

You have seven points so you can fit with 
      1, Cos[t], Sin[t], Cos[2 t], Sin[2 t], Cos[3 t], and Sin[3 t].

G.8)  Parametric plotting of a predator-prey model in which 
         the prey don't reproduce and the predators don't die

The basic predator-prey model is:
     prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
     pred'[t]  =  -c pred[t]  +  d pred[t] prey[t].
In the case in which the predators don't die and the prey don't reproduce, you have a  =  c  =  0.
Here are some plots related to this set-up:

\(endtime\  = \ 30; 

a\  = \ 0; 
b\  = \ 0.2; 
c\  = \ 0; 

    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}, 
\ \ \ \ \ \ \ 
        PlotStyle -> {{Blue, Thickness[0.015]}}, 
\ \ \ \ \ \ \ 
        PlotRange -> All, 
\ \ \ \ \ \ \ 
        AxesLabel -> {"t", "prey"}]; \)

\(\(predatorplot\  = 

      Plot[fakepred[t], {t, 0, endtime}, 
\ \ \ \ \ \ \ 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
\ \ \ \ \ \ \ 
        PlotRange -> All, 
\ \ \ \ \ \ \ 
        AxesLabel -> {"t", "predators"}]; \)\)

Here they are together:

both = Show[predatorplot, preyplot];

G.8.a.i)

Interpret the plots.

G.8.a.ii)

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

G.8.a.iii)   Just for your information

The model is:
        prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
        pred'[t]  =  -c pred[t]  +  d pred[t] prey[t]
with a  =  c  =  0. 
You can use the chain rule to explain the linear relationship of the number of prey as a function of the number of predators.
Here's how it goes:
Put 
        prey  =  f[pred]. 
This gives
        prey[t]  =  f[pred[t]].
The chain rule says
        prey'[t]  =  f'[pred[t]] pred'[t].
Consequently,
        f'[pred[t]] = ((prey'[t])/(pred'[t])).
But 
        prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
        pred'[t]  =  -c pred[t]  +  d pred[t] prey[t]
        with a  =  c  =  0. 
So 
        prey'[t]  =  - b prey[t] pred[t]
        pred'[t]  =   d pred[t] prey[t].
Consequently, 
        f'[pred[t]]  =  
        prey'[t]/pred'[t]  =  -(b/d).
This works for any t so
        f'[pred]  =  -(b/d).
This tells you and everyone else that 
        prey  =  f[pred] 
is a line function of pred and because prey[0]  =  4 and pred[0]  =  1, 
     prey  =  f[pred] 
is the line function with growth rate -(b/d) that goes through {1, 4}.
Try it out:

Clear[f, pred]f[pred_] = 4 - ((b (pred - 1))/d)

Solve[f[pred] == 0]

\(\(\ \ \ 
    \(Plot[f[pred], {pred, 1, 3}, PlotStyle -> Thickness[0.01], 

      \ \ \ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ \ \ 
      AxesLabel -> {"predator", "prey"}]; \)\)\)

This is the same as:

\(\(ParametricPlot[{fakepred[t], fakeprey[t]}, 
\ \ \ {t, 0, endtime}, 
      PlotRange -> All, 
		\ \ PlotStyle -> Thickness[0.01], 
		\ \ \ 
      AxesLabel -> {"predator", "prey"}]; \)\)

You will get 
      prey  =  prey[0]  -  (b/d) (pred  -  pred[0])
for any predator-prey model with a  =  c  =  0.

G.9)  Politics and the environment

G.9.a.i)

Study the Tutorial on the predator-prey model based on the differential equations
     prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
     pred'[t]  =  -c pred[t]  +  d pred[t] prey[t].
Use the same a, b, c, and d as in the Tutorial but use starting data with 
          prey[0]   = (c/d) + h,
          pred[0]  =  (a/b) + h
 for a relatively small value of h:

\(h\  = \ 0.6; 
endtime\  = \ 20; 

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]\  == \ c\/d\  + \ h, 
\ \ \ \ \ \ \ \ \ 
          pred[0]\  == \ a\/b\  + \ h}, 

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

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


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


    cycleplot\  = 
\ \ \ 
      ParametricPlot[{fakeprey[t], fakepred[t]}, 

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

        \ \ \ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ \ \ 
        AxesLabel -> {"prey", "predator"}, 
\ \ \ \ \ \ \ \ \ \ 
        AxesOrigin -> {c\/d, a\/b}]; \)

Why was it a good idea to set AxesOrigin -> {(c/d), (a/b)}?
Reducing the size of h has what effect on the fluctuations of the two populations? What happens when you make h  =  0.01?
How about h  =  0?
Try to explain, in your own terms, why this turned out the way it did.

As you plot, pay special attention to the numbers along the axes.

G.9.a.ii)

Most folks call the populations 
      prey  = (c/d) and predators  =  (a/b) 
for the predator-prey model
      prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
      pred'[t]  =  -c pred[t]  +  d pred[t] prey[t] 
that you worked with above by the name "equilibrium populations."
Imagine that the prey are insects and that the predators are birds or fish. The government goes into action, spreading an insecticide like DDT all over the place.  The insecticide kills both predator and prey in proportion to their current population, so the new model becomes
     prey'[t]  =  
        a prey[t]  -  b prey[t] pred[t]  -  u prey[t]                pred'[t]  =  
        -c pred[t]  +  d pred[t] prey[t]  -  v pred[t]
where a, b, c, d, u, and v are all positive constants.
This is the same as
    prey'[t]  =  
        (a  -  u) prey[t]  -  b prey[t] pred[t] 
    pred'[t]  =  
        -(c  +  v) pred[t]  +  d pred[t] prey[t] 
What are the new equilibrium populations?
How do these new equilibrium populations tell you that the net effect of the government action was to raise the equilbrium population of the prey (insects) but to lower the equilibrium population of the predators?
Is that what the government had in mind?

You don't need the machine for this one.

G.9.b)

In the Tutorials, you saw a rationale for looking at the differential equations
     prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
     pred'[t]  =  -c pred[t]  +  d pred[t] prey[t]
for modeling the predator-prey relationship. These differential equations are not carved in marble. In fact, current literature pays attention to trying to modify these equations to make them more realistic.

For example, you can decide to replace the second equation by
     pred'[t]  =  
        -((c pred[t])/prey)[t] + d pred[t] prey[t].
Rationale:
It's reasonable to assume that the birth rate of the predators is proportional both to the current number of the predators and the size of the food supply, and the death rate of the predators is likely to be directly proportional to the current population of predators and inversely proportional to the number of prey.
Here's a try:

\(endtime\  = \ 30; 
a\  = \ 0.5; 
b\  = \ 0.3; 
c\  = \ 0.6; 

    d\  = \ 0.2; 

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

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

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

          \ \ \ \ \ \ \ \ \ prey[0]\  == \ 5, 
\ \ \ \ \ \ \ \ \ 
          pred[0]\  == \ 3}, 
\ \ \ {prey[t], pred[t]}, {t, 0, endtime}]; 

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


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


    ParametricPlot[{fakeprey[t], fakepred[t]}, 

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

      \ \ \ \ \ \ \ \ \ \ AspectRatio -> 1\/GoldenRatio, 

      \ \ \ \ \ \ \ \ \ \ PlotRange -> All, 
\ \ \ \ \ \ \ \ \ \ 
      AxesLabel -> {"prey", "predator"}, 
\ \ \ \ \ \ \ \ \ \ 
      AxesOrigin -> {c\/d, a\/b}]; \)

Come up with the equillibrium populations in this model and reset the axes origin accordingly.

Tip:

Look at
       0  =  a prey[t]  -  b prey[t] pred[t]
           =  prey[t] (a   -  b pred[t]).  
And look at      
       0  = -(c/prey)[t] + d prey[t]
           =  ((-c  +  d prey[t]^2)/prey)[t].
Don't be afraid to take a square root.

G.10)  Epidemics

The March 4, 1978 issue of the British Medical Journal reported all the details about a flu epidemic that swept through a boarding school in early 1978. Scientists were delighted that the actual data fell in line with those predicted by the SIR model.
The model for this epidemic is:
        Sus'[t]  =  -a Sus[t] Inf[t],
        Inf'[t]  =  a Sus[t] Inf[t]  -  b Inf[t]
        Recov'[t]  =  b Inf[t] 
with
        Sus[0]  =  762,  Inf[0]  =  1, Recov[0]  =  0, 
        a  =  0.00218 and b  = 0.44036.

Here is a three-dimensional replay of this epidemic:

\(endtime\  = \ 30; 

a\  = \ 0.00218; 
b\  = \ 0.44036; 


    Clear[Sus, Inf, Recov, t, Derivative]\)\(approxsolutions\  = \ 

      NDSolve[{\(Sus'\)[t]\  == \ \(-\ a\)\ Sus[t]\ Inf[t], 

          \ \ \ \ \ \ \ \ \ 
          \(Inf'\)[t]\  == \ a\ Sus[t]\ Inf[t]\  - \ b\ Inf[t], 

          \ \ \ \ \ \ \ \ \ \(Recov'\)[t]\  == \ b\ Inf[t], \ 

          \ \ \ \ \ \ \ \ \ Sus[0]\  == \ 763, 
\ \ \ \ \ \ \ \ \ 
          Inf[0]\  == \ 1, 
\ \ \ \ \ \ \ \ \ Recov[0]\  == \ 0}, 

        \ \ \ {Sus[t], Inf[t], Recov[t]}, {t, 0, endtime}]; 


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


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


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


    epidemicplot\  = \ 
      ParametricPlot3D[{Sus[t], Inf[t], Recov[t]}, 
\ \ \ \ {t, 0, endtime}, 
        DisplayFunction -> Identity]; 
threedims\  = \ Axes3D[500, 20]; 

    Show[threedims, epidemicplot, ViewPoint -> CMView, 
\ \ 
      PlotRange -> All, Boxed -> False, 
\ \ 
      DisplayFunction -> $DisplayFunction]; \)

You can trace the progress of the epidemic starting at the lower left with lots of susceptibles and finishing at the top with lots of recovereds. In the midst of the full-blown epidemic, there are lots of infecteds.
Here's how to use a two-dimensional plot to see how the number of infecteds varies as a function of the susceptibles.

\(\(ParametricPlot[{Sus[t], Inf[t]}, {t, 0, endtime}, 
\ \ \ \ \ \ 
      PlotStyle -> {{Blue, Thickness[0.015]}}, 
\ \ \ \ \ \ 
      AxesLabel -> {"Sus", "Inf"}]; \)\)

G.10.a)

The model says 
      Inf'[t]  =  a Sus[t] Inf[t]  -  b Inf[t]
                      =  (a Sus[t]  -  b) Inf[t] . 
Why does this tell you that the number of infecteds is the greatest when there are (b/a) susceptibles?
Check this out by adding the vertical line through {(b/a), 0} and {(b/a), 300} to the plot immediately above.

G.10.b)

Plot the susceptibles as a function of the infecteds.

G.10.c.i)

Plot the susceptibles as a function of the recovereds.

G.10.c.ii)   Just for your information

The curve you got  in part c.i) above should look like a dead ringer for exponential decay. This can be explained as follows:
Put 
       Sus  =  f[recov]. 
This gives
       Sus[t]  =  f[recov[t]].
The chain rule says
       Sus'[t]  =  f'[recov[t]] recov'[t].
Consequently,
       f'[recov[t]] = ((Sus'[t])/(recov'[t])).
But in the model, 
       Sus'[t]  =  - a Sus[t] Inf[t],
       Recov'[t]  =  b Inf[t].
So  
      f'[recov[t]]=-((a Sus[t] Inf[t])/(b Inf[t]))
                              =  - ((a Sus[t])/b)
                              =  - ((a f[recov[t]])/b),
because 
       f[recov[t]]  =  Sus[t]
       f'[recov]  =  -((a f[recov])/b).
This explains the exponential decay you saw above because the solution of this little differential equation is given by:

Clear[y, f, recov, a, b]\(DSolve[{\(y'\)[recov]\  == \ \(-\(\(a\ y[recov]\)\/b\)\), 

        \ \ \ \ \ \ \ \ y[0]\  == \ f[0]}, y[recov], recov]\)

So
        Sus  =  
        f[recov]  =  f[0] (e^(-[a/b] recov))  . 
Beause Sus  =  762 when recov  =  0,this gives
        Sus  =  762 (e^(-[a/b] recov))  . 
And this confirms the exponential decay.

G.11)  Collision?

G.11.a)

A Spartan missile and a Trojan missle are both flying at the same constant altitude.
At time t, the Spartan missile is at the point:

Clear[spartan, t]spartan[t_] = {16.1 - 7 t + (t^2), 13 t - 2 (t^2)}

At the same time t, the Trojan missile is at the point:

Clear[trojan]trojan[t_] = {26 - 13 t + 2 (t^2), 23 - 5 t + (t^2)}

Here is a plot of the paths of the two missiles:

\(spartanpath\  = 
\ \ \ \ \ 
      ParametricPlot[spartan[t], {t, 0, 6}, 
\ \ \ \ \ 
        PlotStyle -> {{Blue, Thickness[0.02]}}, 
\ \ \ \ \ 
        DisplayFunction -> Identity]; 


    trojanpath\  = 
\ \ 
      ParametricPlot[trojan[t], {t, 0, 6}, 
\ \ \ \ \ 
        PlotStyle -> {{Red, Thickness[0.02]}}, 
\ \ \ \ \ 
        DisplayFunction -> Identity]; 


    Show[spartanpath, trojanpath, 
\ \ \ \ \ AspectRatio -> 1\/GoldenRatio, 
      
\ \ \ \ \ DisplayFunction -> $DisplayFunction]; \)

Their paths cross, but do they crash?

In order to crash, both missiles must be at the same point at the same time.