VC.04 2D Vector Fields and Their Trajectories
GIVE IT A TRY!

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

    \(Off[ParametricPlot3D::ppcom]; \)\)Off[ParametricPlot::ppcom]; Off[Plot::plnr]; SetOptions[Limit, Analytic -> True];

CMView = {2.7, 1.6, 1.2};

Vector Drawer

BeginPackage["Vector3D`", "Graphics`Colors`"];

\(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]."; 

    \(Perpend::"usage" = 
      "Perpend[a] returns a unit vector which is perpendicular to the vector a."; \)\)Arrow::"usage" = 
      "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. \
Arrow[a, Tail -> tail] gives a vector from tail to a."; Vector::"usage" = 
      "Vector[a, b] gives a 2D or 3D vector from a to b."; ArrowHead::"usage" = 
      "ArrowHead[a, vec] produces an arrowhead with its tip placed at the \
point a, pointing in the direction of the vector vec."; Tail::"usage" = 
      "Tail -> point puts the tail of the vector at point."; Norm::"usage" = 
      "Norm[vec] returns the Euclidean length of vec."; Normalize::"usage" = 
      "Normalize[vec] returns a unit vector in the direction of vec, \
provided vec is nonzero."; Aperture::"usage" = 
      "Aperture is the ratio of the radius of the base to the length of the \
head of a vector."; TipSize::"usage" = 
      "TipSize specifies an absolute size for the inner tip length of the \
head of a vector."; TipRatio::"usage" = 
      "TipRatio is the ratio of the inner length of the arrowhead to the \
length of the shaft of a vector."; HeadRatio::"usage" = 
      "HeadRatio is the ratio of the outer length of the arrowhead to the \
length of the shaft of a vector."; HeadSize::"usage" = 
      "HeadSize specifies an absolute size for the head of a vector."; 
    ScaleFactor::"usage" = 
      "ScaleFactor specifies the amount to scale a vector in length and may \
be either a number or a pure function."; ZeroVectorPointSize::"usage" = 
      "ZeroVectorPointSize is the size of the point used to represent a \
zero vector."; Options[Vector] = 
      Options[Arrow] = 
        Options[ArrowHead] = {HeadRatio -> 0.18, 
            TipRatio -> 0.14, Aperture -> 0.3, 
            ScaleFactor -> 1, HeadType -> Polygon, 
            ShaftQ -> True, EdgesQ -> True, 
            VectorColor -> RGBColor[0, 0, 1], ShaftWidth -> 0.005, 
            ZeroVectorPointSize -> 0.01}; Begin["`Private`"]; SetOptions[ParametricPlot, AspectRatio -> Automatic]; SetOptions[Plot, AspectRatio -> Automatic]; SetOptions[Graphics, AspectRatio -> Automatic]; 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)]; Norm[a_?VectorQ] := N[(Sqrt[a . a])]; Normalize[a_?VectorQ] := If[a . a == 0, a, (a/Norm)[a]]; Perpend[a_?VectorQ] := 
      Normalize[{a[[2]], 
            -a[[1]]}] /; 
        Length[a] == 2Perpend[a_?VectorQ] := 
      Normalize[
          If[a[[1]] == 0, {1, 0, 0}, {0, 
              a[[3]], 
              -a[[2]]}]] /; 
        Length[a] == 3Base3D = 
      Table[{Cos[((2 k π)/8.) ], Sin[((2 k π)/8.) ]}, {k, 0, 
          8}]; Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = 
            {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, 
                  TipRatio, Aperture, ScaleFactor, VectorColor, 
                  ZeroVectorPointSize} /. {opts} /. Options[Vector]; 
          If[from == to, 
            Return[Graphics[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]; tr = 0.75 hr]; 
          If[NumberQ[ts], tr = (ts/Norm)[diff]]; 
          Graphics[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              With[{trans = ap hr Norm[diff] Perpend[diff]}, {vc, 
                  Thickness[sw], 
                  ht[{tip, tip - hr diff + trans, tip - tr diff, 
                      tip - hr diff - trans, tip}]}]}]] /; 
        Length[from] == Length[to] == 2Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, eq, sw, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = 
            {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, 
                  ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts}
                   /. Options[Vector]; 
          If[from == to, 
            Return[Graphics3D[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]]; tr = (hr/2); 
          Graphics3D[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              If[eq, {}, EdgeForm[]], {SurfaceColor[vc], 
                With[{temp = Perpend[diff]}, 
                  (Polygon[Append[#1, tip]]&)/@
                    Partition[
                      (tip - hr diff + #1&)/@
                        ((ap hr Norm[diff] Base3D) . {temp, 
                              Normalize[temp\[Cross]diff]}), 2, 1]]}}]] /; 
        Length[from] == Length[to] == 3Arrow[a_, Tail -> b_, opts___] := Vector[b, a + b, opts]Arrow[a_, opts___] := Vector[Table[0, {Length[a]}], a, opts]\(ArrowHead[a_, b_, opts___] := 
      Vector[a - b, a, ShaftQ -> False, opts]; 

End[]; 

    \(EndPackage[]; \)\)

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

Spherical Setup

\(Clear[polar]; 

    polar[{x_, y_, z_}] := {\@\(x\^2 + y\^2 + z\^2\), ArcTan[y\/x], 
        ArcCos[z\/\@\(x\^2 + y\^2 + z\^2\)]}; 
Clear[X, length]; 

    length[X_] := \@\(X . X\); 
\(rr = 3; \)\)R = 5; scale = 0.8; red = RGBColor[1, 0, 0]; blue = RGBColor[0, 0, 1]; color1 = RGBColor[0.8, 0.2, 0.8]; color2 = RGBColor[0.2, 0.8, 0.2]; black = RGBColor[0, 0, 0]; \(Clear[axes]; 

    axes := Graphics3D[{Line[{X, {1.5\ rr, 0, 0}}], 
          Line[{X, {0, 1.5\ rr, 0}}], Line[{X, {0, 0, 2\ rr}}]}]; 

    Clear[xyplane]; 

    xyplane := 
      Graphics3D[
        Polygon[{{\(-2\)\ rr, \(-2\)\ rr, 0}, {\(-2\)\ rr, 2\ rr, 0}, {2\ rr, 
              2\ rr, 0}, {2\ rr, \(-2\)\ rr, 0}}]]; 
\(X = {0, 0, 0}; \)\)origin = Graphics3D[{black, Point[X]}]; Z = {2, 3, 5}; \(Clear[point]; 

    point := Graphics3D[{red, PointSize[0.03], Point[scale\ Z]}]; 

    \(Z0 = {0, 0, length[0.8\ Z]}; \)\)\(Clear[zpoint]; 
zpoint = Graphics3D[{black, Point[Z0]}]; 

    \(XY = scale\ Z\ {1, 1, 0}; \)\)\(Clear[xypoint]; 
xypoint := Graphics3D[{black, Point[XY]}]; 

    \(X0 = {length[scale\ Z\ {1, 1, 0}], 0, 0}; \)\)xpoint = Graphics3D[{black, Point[X0]}]; Y0 = {0, length[scale Z {1, 1, 0}], 0}; ypoint = Graphics3D[{black, Point[Y0]}]; \(Clear[ray]; 
ray := Graphics3D[{blue, Line[{X, scale\ Z}]}]; 

    Clear[ray1]; 
ray1 := Graphics3D[{color1, Line[{scale\ Z, Z}]}]; 

    Clear[xyray]; 
xyray := Graphics3D[{red, Line[{X, Z\ {1, 1, 0}}]}]; 

    \(rays := {ray, xyray}; \)\)xyaux = 
      (length[0.8 Z] (0.8 Z {1, 1, 0}))/
        length[0.8 Z {1, 1, 0}]; \(Clear[zarc]; 

    zarc = Graphics3D[{blue, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 0, 
                \(polar[0.8\ Z]\)[[3
                  ]],  .1}]]}]; 
Clear[zarc1]; 

    zarc1 = Graphics3D[{color1, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 
                \(polar[0.8\ Z]\)[[3]], 
                1,  .1}]]}]; 
Clear[zarc2]; 

    zarc2 = Graphics3D[{color1, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 
                \(-0.8\), 0,  .1}]]}]; 
Clear[xyarc]; 

    xyarc = Graphics3D[{blue, 
          Line[Table[
              Cos[t]\ X0 + Sin[t]\ Y0, {t, 0, 
                \(polar[XY]\)[[2
                  ]],  .1}]]}]; 
Clear[xyarc1]; 

    xyarc1 = Graphics3D[{color1, 
          Line[Table[Cos[t]\ X0 + Sin[t]\ Y0, {t, \(- .5\), 0,  .1}]]}]; 

    Clear[xyarc2]; 

    xyarc2 = Graphics3D[{color1, 
          Line[Table[
              Cos[t]\ X0 + Sin[t]\ Y0, {t, 
                \(polar[XY]\)[[2]], 
                \(polar[XY]\)[[2]] + 
                  0.3,  .1}]]}]; 
Clear[stick]; 

    stick = Graphics3D[{red, Line[{scale\ Z, XY}]}]; 

    \(xaxislabel = Text["x", {0.363, 0.26}]; \)\)yaxislabel = Text["y", {0.73, 0.3}]; zaxislabel = Text["z", { .5, 0.67}]; stickplot = 
      Show[xyplane, axes, ray, point, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; rho = Text["r", {0.53, 0.45}]; rlabel = 
      Show[xyplane, axes, ray, point, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {rho, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; phi = Text["s", {0.55, 0.57}]; slabel = 
      Show[xyplane, zpoint, axes, ray, point, zarc, zarc1, zarc2, 
        Boxed -> False, ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {phi, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; theta = Text["t", { .5, 0.3}]; tlabel = 
      Show[xyplane, axes, xypoint, xpoint, ray, stick, point, xyarc, xyarc1, 
        xyarc2, Boxed -> False, ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {theta, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; all = 
      Show[xyplane, axes, ray, stick, point, zpoint, xpoint, xypoint, xyarc, 
        xyarc1, xyarc2, zarc, zarc1, zarc2, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {theta, phi, rho, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity];

Crisscross

\(Clear[x]; 

    crisscross = 
      Plot[{Cos[x], Sin[x]}, {x, π\/8, \(3\ π\)\/8}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.02]}}, 
        Axes -> None, AspectRatio -> 1\/3, 
        DisplayFunction -> Identity]; \)

G.1)  Looking for sinks (drains)*

G.1.a)

Here is a vector field and a plot:

\(Clear[Field, x, y]; 

    Field[x_, y_] = {\(-0.036\)\ x\ E\^\(x\^2\) + 0.123\ y, 
        0.123\ x - 0.087\ y\ E\^\(y\^2\)}; 

    fieldplot = 
      Table[Arrow[Field[x, y]\/\(1 + Norm[Field[x, y]]\), 
          Tail -> {x, y}], {x, \(-1.25\), 1.25, 2.5\/10}, {y, \(-1.25\), 
          1.25, 2.5\/10}]; 


    Show[fieldplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)

As you can see from the plot, if you drop a cork into this flow at any starting point within the plot other than {0, 0}, then it will float along to one of two special points and stop dead. 
Plot some trajectories in an effort to find the approximate locations of these two special points.

Copy, paste and edit.

G.1.b) Scaling a vector field

Here's an unscaled plot of a two dimensional vector field:

\(Clear[Field, x, y]; 

    Field[x_, y_] = {18.4 - 0.8\ x\^3 - 8.2\ y, \(-8.2\)\ x - 0.4\ y\^3}; 

    \(fieldplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}], {x, \(-5\), 5}, {y, 
          \(-5\), 5}]; 


    Show[fieldplot, Axes -> True, AxesLabel -> {"x", "y"}]; 
    \)\)

This plot sucks. You can't see what's happening.
Give a scaled plot of this vector field at the same points to show it off to good advantage.
Pretend that this field represents fluid flow and use your good plot to determine whether there are any drains within the plot.  If you can't spot them visually, run some trajectories and see where they lead you.

G.2)  Flow along and flow across*

Here's the vector field
     Field[x, y]={x+ (y/3),x- (y/2)} 
plotted at selected points on the circle of radius 2 centered at {1, 0}:

\(Clear[Field, m, n, x, y]; 
m[x_, y_] = x + y\/3; 

    \(n[x_, y_] = x - y\/2; \)\)Field[x_, y_] = {m[x, y], n[x, y]}; \(Clear[t]; 
{x[t_], y[t_]} = {1, 0} + 2\ {Cos[t], Sin[t]}; 

    \(curveplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2\ π}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; \)\)jump = (π/12); \(fieldoncurve = 
      Table[Arrow[Field[x[t], y[t]], Tail -> {x[t], y[t]}], {t, 0, 
          2\ π - jump, jump}]; 


    Show[curveplot, fieldoncurve, AxesOrigin -> {0, 0}, 
      AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

Here are the tangential components of the field vectors shown above:

Clear[tangent, tancomponent];
tangent[t_]={x'[t],y'[t]};
tancomponent[t_] = 
        Field[x[t], y[t]] . tangent[t]/tangent[t] . tangent[t] 
          tangent[t];

actualflowalong = 
        Table[Arrow[tancomponent[t], Tail -> {x[t], y[t]}], {t, 0, 
            2 π - jump, jump}];

flowalongplot = 
        Show[curveplot, actualflowalong, AxesOrigin -> {0, 0}, 
          AspectRatio -> Automatic, 
          AxesLabel -> {"x", "y"}, 
          DisplayFunction -> $DisplayFunction];

G.2.a.i)

You make the call:
Is the net flow of this vector field along this curve clockwise or counterclockwise?

G.2.a.ii)

Go with the same vector field and the same curve as above, but this time look at normal components of the vector field at points on the curve:

Clear[normal, normalcomponent];
normal[t_]={y'[t],-x'[t]};
normalcomponent[t_] = 
        Field[x[t], y[t]] . normal[t]/normal[t] . normal[t] normal[t];

actualflowacross = 
        Table[Arrow[normalcomponent[t], Tail -> {x[t], y[t]}], {t, 0, 
            2 π - jump, jump}];

\(flowacrossplot = 
        Show[curveplot, actualflowacross, 
		
          AspectRatio -> Automatic, 
          AxesLabel -> {"x", "y"}, 
          DisplayFunction -> $DisplayFunction]\);

You make the call:
Is the net flow of this vector field across this curve from inside to outside or from outside to inside?

G.2.b)

Here's a new vector field plotted at a selection of points on an ellipse:

\(Clear[Field, m, n, x, y, t]; 
m[x_, y_] = \(x - y\)\/3; 

    \(n[x_, y_] = \(x + y\)\/2; \)\)Field[x_, y_] = {m[x, y], n[x, y]}; {x[t_], y[t_]} = {2, 1} + {2 Cos[t], Sin[t]}; curveplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2 π}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; jump = (π/12); \(fieldoncurve = 
      Table[Arrow[Field[x[t], y[t]], Tail -> {x[t], y[t]}], {t, 0, 
          2\ π - jump, jump}]; 


    Show[curveplot, fieldoncurve, AxesOrigin -> {0, 0}, 
      AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

Plot, as in parts a.i) and a.ii), the tangential and normal components of the field vectors on this curve and then make the calls:
-> Is the net flow of this vector field along this curve clockwise or counterclockwise?
-> Is the net flow of this vector field across this curve from inside to outside or from outside to inside?

G.2.c)

Here's a new vector field plotted at a selection of points on the circle of radius 1 centered at {0, 0}:

\(Clear[Field, m, n, x, y, t]; 
m[x_, y_] = x\^3 - 3\ x\ y\^2; 

    \(n[x_, y_] = \(-3\)\ x\^2\ y + y\^3; \)\)Field[x_, y_] = {m[x, y], n[x, y]}; {x[t_], y[t_]} = {Cos[t], Sin[t]}; curveplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2 π}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; jump = (π/12); \(fieldoncurve = 
      Table[Arrow[Field[x[t], y[t]], Tail -> {x[t], y[t]}], {t, 0, 
          2\ π - jump, jump}]; 


    Show[curveplot, fieldoncurve, 
	Axes -> True, 
      AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

Plot, as in parts a.i) and a.ii), the tangential and normal components of these field vectors on this curve and then make the calls:
-> Is the net flow of this vector field along this curve clockwise, counterclockwise, or is the flow of this vector field along this curve 0 in the sense that the counterclockwise flow balances the clockwise flow?
-> Is the net flow of this vector field across this curve from inside to outside or from outside to inside, or is the flow of this vector field across this curve 0 in the sense that the inward flow balances the outward flow?

G.3)  Normals, tangents and dot plots*

G.3.a.i) Direction of the tangent vector {x'[t],y'[t]} at {x[t], y[t]}

When a closed curve is parameterized by {x[t], y[t]}, then as you advance along the curve in the direction of the parameterization, which way do the tangent vectors {x'[t],y'[t]} at {x[t], y[t]} point; in the direction you are going, or in the direction opposite to the direction you are going?

G.3.a.ii) Direction of the normal vector {y'[t],-x'[t]} at {x[t], y[t]} for a counterclockwise parameterization

When a closed curve is parameterized by {x[t], y[t]} in the counterclockwise way, then as you advance along the curve in the direction of the parameterization, which way do the normal vectors {y'[t],-x'[t]} at {x[t], y[t]} point, to the right to the outside of the curve, or to the left to the inside of the curve?

Take a look at this:

Clear[x, y, t, tangent, normal];
{x[t_], y[t_]} = {0.6 Cos[t], 1.7 Sin[t]};
tangent[t_]={x'[t],y'[t]};normal[t_]={y'[t],-x'[t]};curveplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2 π}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; jump = (π/3); tangentvectors = 
      Table[Arrow[tangent[t], Tail -> {x[t], y[t]}], {t, 0, 
          2 π - jump, jump}]; \(normalvectors = 
      Table[Arrow[normal[t], Tail -> {x[t], y[t]}], {t, 0, 
          2\ π - jump, jump}]; 


    Show[curveplot, tangentvectors, normalvectors, 
      AxesLabel -> {"x", "y"}, AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)

G.3.a.iii) Direction of the normal vector {y'[t],-x'[t]} at {x[t], y[t]}

When a closed curve is parameterized by {x[t], y[t]} in the clockwise way, then as you advance along the curve in the direction of the parameterization, which way do the normal vectors {y'[t],-x'[t]} at {x[t], y[t]} point, to the right to the inside of the curve, or to the left to the outside of the curve?

Take a look at this:

Clear[x, y, t, tangent, normal];
{x[t_], y[t_]} = {0.6 Sin[t], 1.7 Cos[t]};
tangent[t_]={x'[t],y'[t]};normal[t_]={y'[t],-x'[t]};curveplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2 π}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; jump = (π/3); tangentvectors = 
      Table[Arrow[tangent[t], Tail -> {x[t], y[t]}], {t, 0, 
          2 π - jump, jump}]; \(normalvectors = 
      Table[Arrow[normal[t], Tail -> {x[t], y[t]}], {t, 0, 
          2\ π - jump, jump}]; 


    Show[curveplot, tangentvectors, normalvectors, 
      AxesLabel -> {"x", "y"}, AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)

G.3.b.i)

Here's a vector field:

\(Clear[Field, m, n, x, y]; 
m[x_, y_] = 1; 
\(n[x_, y_] = Sin[x\ y]; 
    \)\)Field[x_, y_] = {m[x, y], n[x, y]}

Here's a curve:

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

    \(curveplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2\ π}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        AxesLabel -> {"x", "y"}]; \)\)

Look at this plot of some of the unit tangent vectors
       unittangent[t]=(({x'[t], y'[t]})/(Sqrt[((x'[t])^2) +  ((y'[t])^2)]))
pointing in the direction of the parameterization:

Clear[t, unittangent];
unittangent[t_]=(({x'[t],y'[t]})/(Sqrt[((x'[t])^2)+((y'[t])^2)]));
jump = (π/8); \(unittangents = 
      Table[Arrow[unittangent[t], Tail -> {x[t], y[t]}, 
          VectorColor -> Red], {t, 0, 2\ π - jump, jump}]; 

    \(Show[curveplot, unittangents]; \)\)

Now look at this plot of 
       Field[x[t], y[t]] . unittangent[t]:

\(\(dotplot = 
      Plot[Field[x[t], y[t]] . unittangent[t], {t, 0, 2\ π}, 
        PlotStyle -> {{Thickness[0.01], Blue}}, 
        AxesLabel -> {"t", \*""""}, AspectRatio -> 1\/2]; 
    \)\)

Interpret this plot.

G.3.b.ii)

Go with the same vector field and curve as in part i) and look at this plot of 
       Field[x[t], y[t]].{y'[t],-x'[t]}:

dotplot=Plot[Field[x[t], y[t]].{y'[t],-x'[t]},{t, 0, 2 π},PlotStyle -> {{Thickness[0.01], Blue}},\(AxesLabel -> {"t", \*""""}\),AspectRatio -> (1/2)];

Interpret this plot.  Does your answer change if you replace 
       {y'[t],-x'[t]}
by the normalized vector 
       (({y'[t], -x'[t]})/(Sqrt[((x'[t])^2) +  ((y'[t])^2)]))?

G.4)  The most important vector field of them all:
The gradient field*

One way to come up with a vector field is to take a function f[x, y] and put 
       Field[x, y] = gradf[x, y] = ∇f[x, y].
Try it out:

\(Clear[f, gradf, Field, x, y]; 

    f[x_, y_] = 1 - 1\/2\ \((x\^2 + 1.5\ y\^2)\); 

    gradf[x_, y_] = {D[f[x, y], \ x], D[f[x, y], \ y]}; 


    Field[x_, y_] = gradf[x, y]\)

This is the gradient field of 
       f[x, y] = 1 -  (1/2) ((x^2) +  1.5 (y^2)).
Note that {x, y} = {0, 0 maximizes f[x, y].
Here's a plot of the scaled gradient field of f[x, y] shown with the maximizer at {0, 0}:

maximizerplot = Graphics[{Red, PointSize[0.04], Point[{0, 0}]}]; 
    scalefactor = 0.2; \(gradfieldplot = 
      Table[Arrow[gradf[x, y], Tail -> {x, y}, 
          ScaleFactor -> scalefactor], {x, \(-2\), 2, 0.5}, {y, \(-2\), 
          2, 0.5}]; 


    Show[maximizerplot, gradfieldplot, Axes -> Automatic, 
      AxesLabel -> {"x", "y"}]; \)

G.4.a.i)

Why did this happen?
Where are the trajectories in this gradient field headed?
What do you learn about a function f[x, y] by looking at a plot of its gradient field?

G.4.a.ii)

This time go with the function
       f[x, y]=(2/(1  +  (((y  -  2))^2)  +  
                1.5 (((x  -  1))^2))).
Say why you are sure that {x, y} = {1, 2} maximizes f[x, y].  Then plot the (scaled) gradient field of f[x, y] near {1, 2}.
Describe what you see and explain why you see it.

G.4.a.iii)

You are given a certain function f[x, y] and the information that 
       {a, b, f[a, b]} 
sits on the top of a crest of the surface given by 
       z = f[x, y]. 
You center a small circle at {a, b} and plot the gradient field of f[x, y] on this circle.  
Is the net flow of the gradient field of f[x, y] across this circle from inside to outside or from outside to inside?

G.4.b.i) Minimizers

This time go with 
       f[x, y] = (e^((1/9) ((x^2) +  2  (y^2))))· .

\(Clear[f, gradf, Field, x, y]; 

    f[x_, y_] = E\^\(1\/9\ \((x\^2 + 2\ y\^2)\)\); 

    gradf[x_, y_] = {D[f[x, y], \ x], D[f[x, y], \ y]}; 


    Field[x_, y_] = gradf[x, y]\)

This is the gradient field of 
       f[x, y] = (e^((1/9) ((x^2) +  2  (y^2)))).
Note that {x, y} = {0, 0} minimizes f[x, y].
Here's a plot of the scaled gradient field of f[x, y] shown with the minimizer at {0, 0}:

minimizerplot = Graphics[{Red, PointSize[0.04], Point[{0, 0}]}]; 
    scalefactor = 0.8; \(gradfieldplot = 
      Table[Arrow[gradf[x, y], Tail -> {x, y}, 
          ScaleFactor -> scalefactor], {x, \(-1\), 1, 0.25}, {y, \(-1\), 
          1, 0.25}]; 


    Show[minimizerplot, gradfieldplot, Axes -> Automatic, 
      AxesLabel -> {"x", "y"}]; \)

Why did this happen?
What are the trajectories in this gradient field trying to get away from?

G.4.b.ii)

You are given a certain function f[x, y] and the information that 
       {a, b, f[a, b]} 
sits on the bottom of a dip of the surface given by 
       z = f[x, y]. 
You center a small circle at {a, b} and plot the gradient field of f[x, y] on this circle.  Is the net flow of the gradient field of f[x, y] across this circle from inside to outside or from outside to inside?

G.4.c) Crests and dips

Here's a plot of the surface
       z=f[x, y]=((Sin[2 x]  +  Sin[3 y])/(0.5  +  (x^2)  +  (y^2)))
for -1 ≤ x ≤ 1 and -1 ≤ y ≤ 1:

\(Clear[f, x, y]; 

    f[x_, y_] = \(Sin[2\ x] + Sin[3\ y]\)\/\(0.5 + x\^2 + y\^2\); 

    \(Plot3D[f[x, y], {x, \(-1\), 1}, {y, \(-1\), 1}, 
      AxesLabel -> {"x", "y", "z"}, 
      ViewPoint -> CMView]; \)\)

Here's a plot of the (scaled) gradient field of the same function for the same x's and the same y's:

\(Clear[gradf]; 
gradf[x_, y_] = {D[f[x, y], \ x], D[f[x, y], \ y]}; 

    scalefactor = 0.07; 

    gradfieldplot = 
      Table[Arrow[gradf[x, y], Tail -> {x, y}, \ 
          ScaleFactor -> scalefactor], {x, \(-1\), 1, 0.2}, {y, \(-1\), 
          1, 0.2}]; 


    Show[gradfieldplot, Axes -> Automatic, 
      AxesLabel -> {"x", "y"}]; \)

What information about the crests and dips on the surface does this plot of the gradient field of f[x, y] convey?
Use the information in this plot and the FindMinimum instruction to locate the position {a, b, f[a, b]} of the top of the highest crest, and the position {c, d, f[c, d]} of the bottom of the deepest dip on the part of the surface 
       z = f[x, y] 
plotted above.

G.4.d) The negative gradient field

Here's a contour plot of
       f[x, y] = 1.2  (e^((x^2))) +  2.9  (e^((y^2))) -  8.2  x (y^2).
for 
       0.8 ≤ x ≤ 1.3 and 0.7 ≤ y ≤ 1.2:

\(Clear[f, x, y]; 

    f[x_, y_] = 1.2\ E\^\(x\^2\) + 2.9\ E\^\(y\^2\) - 8.2\ x\ y\^2; 

    \(ContourPlot[f[x, y], {x, 0.8, 1.3}, {y, 0.7, 1.2}]; \)\)

There's a dip bottoming out at {a, b, f[a, b]} for some {a, b} near {1, 1}.
Look at a plot of the (scaled) negative gradient in the vicinity of {1, 1}.

\(Clear[negativegradient]; 

    negativegradf[x_, y_] = \(-{D[f[x, y], \ x], D[f[x, y], \ y]}\); 

    scalefactor = 0.01; 

    neggradfieldplot = 
      Table[Arrow[negativegradf[x, y], Tail -> {x, y}, 
          ScaleFactor -> scalefactor], {x, 0.9, 1.3, 0.4\/8}, {y, 0.8, 
          1.2, 0.4\/8}]; 


    Show[neggradfieldplot, Axes -> Automatic, 
      AxesLabel -> {"x", "y"}]; \)

Look at that flow.
Use what you see in the plot above, and then use the FindMinimum instruction to get a good estimate of the point {a, b, f[a, b]} at the bottom of the dip. 
Add a plot of {a, b} to the plot of the negative gradient field above. 
Describe what you see and say why you see it.
How do trajectories in this negative gradient field look when they pass through the part of the negative gradient field plotted above?

G.5)  Differential equations and their associated vector fields*

G.5.a)

Here's the vector field coming from the differential equation
       y'[x]=y[x]^2 -  x.

\(Clear[DEField, x, y]; 
DEField[x_, y_] = {1, y\^2 - x}; 

    scalefactor = 0.3; 


    DEFieldplot = 
      Table[Arrow[DEField[x, y], Tail -> {x, y}, 
          ScaleFactor -> scalefactor], {x, \(-2\), 6, 8\/10}, {y, 
          \(-2\), 2, 4\/10}]; 


    Show[DEFieldplot, Axes -> Automatic, 
      AxesLabel -> {"x", "y"}]; \)

Look at the slopes of those field vectors.
Add plots of the solutions of
       y'[x]=y[x]^2 -  x  with y[-1.5] = -2
and
       y'[x]=y[x]^2 -  x  with y[0] = -1
as trajectories of the vector field to the DEField plot above, and discuss the relation between the plots of the solutions and the DEField.
Finally, examine the flow of this vector field, and use it to describe the behavior of solutions of the differential equation
       y'[x]=y[x]^2 -  x
that pass through the part of the xy-plane plotted above.

DSolve cannot find a formula for the solution of either of these differential equations:

Clear[x, y, Derivative];
DSolve[{y'[x]==y[x]^2 - x,y[-1.5] == -2},y[x],x]

Clear[x, y, Derivative];
DSolve[{y'[x]==y[x]^2 - x,y[0] == -1},y[x],x]

This is not a deficiency of Mathematica.  No man, woman, or machine has ever been able to come up with a formula for the solution.  It's an unfortunate fact of science that most differential equations have solutions for which no clean formula is available.
The consequence:
With no formula available, the best way to get an idea of how the solutions of
       y'[x]=y[x]^2 -  x 
look is to look at the flow of the DEField coming from this differential equation.

G.5.b.i)

Here's the vector field coming from the differential equation
       y'[x]=y[x](1- y[x]/2):

\(Clear[DEField, x, y]; 
DEField[x_, y_] = {1, y\ \((1 - y\/2)\)}; 

    \(scalefactor = 0.3; \)\)\(DEFieldplot = 
      Table[Arrow[DEField[x, y], Tail -> {x, y}, 
          ScaleFactor -> scalefactor], {x, 0, 8, 8\/10}, {y, 0, 4, 
          4\/10}]; 

    cutoff = Graphics[{Red, Thickness[0.01], 
          Line[{{\(-0.5\), 2}, {8.5, 2}}]}]; 


    Show[cutoff, DEFieldplot, Axes -> True, 
      AxesLabel -> {"x", "y"}, 
      PlotLabel -> "Logistic field"]; \)

Use what you see to describe the fundamental difference between the behavior of a solution of            
       y'[x]=y[x](1- y[x]/2)
whose plot goes through a point {a, b} with b > 2
and the behavior of a solution of            
       y'[x]=y[x](1- y[x]/2)
whose plot goes through a point {a, b} with 0 < b < 2.

G.5.b.ii)

Here's more of the vector field coming from the logistic differential equation
       y'[x]=y[x](1- y[x]/2)

\(Clear[DEField, x, y]; 
DEField[x_, y_] = {1, y\ \((1 - y\/2)\)}; 

    \(scalefactor = 0.3; \)\)\(DEFieldplot = 
      Table[Arrow[DEField[x, y], Tail -> {x, y}, \ 
          ScaleFactor -> scalefactor], {x, 0, 10, 1}, {y, \(-2\), 4, 
          4\/10}]; 

    cutoff = Graphics[{Red, Thickness[0.01], 
          Line[{{\(-0.5\), 2}, {10.5, 2}}]}]; 


    Show[DEFieldplot, cutoff, Axes -> True, 
      AxesLabel -> {"x", "y"}]; \)

Now you see three distinct types of solutions of
       y'[x]=y[x](1- y[x]/2).
Describe the three distinct types of solutions in words.

G.5.c)

Plot the vector field coming from the differential equation 
       y'[x]=x y[x]
for a healthy selection of {x, y}'s with 
      -4 ≤ x ≤ 4 and -3 ≤ y ≤ 3.
Apply an appropriate scale factor so you can see what's going on.
Look at the resulting plot of the vector field, and describe how solutions with y[-4] > 0 are fundamentally different from solutions with y[-4] < 0.

G.6)  Trajectories: Can they cross?

G.6.a) Crossing of trajectories

Can two trajectories of a vector field ever cross over each other like this:

Show[crisscross, DisplayFunction -> $DisplayFunction];

Why or why not?

Think of it this way:
At the cross-over point, which way must the floating cork move?

G.6.b) Crossing of plots of solutions of differential equations

You have a differential equation 
       y'[x]=f[x, y[x]].
If you have a formula for f[x, y] and you start with 
       y[0] = 1, 
then you can use NDSolve to plot out a function y[x] with
       y[0] = 1 and y'[x]=f[x, y[x]].
On the other hand if you start with y[0] = 2, then you can use NDSolve to plot out a function y[x] with 
       y[0] = 2 and y'[x]=f[x, y[x]].
The question here is: Can the plots of these two solutions ever cross over each other like this:

Show[crisscross, DisplayFunction -> $DisplayFunction];

Why or why not?

G.7)  Drifting along with a tumbleweed

G.7.a)

Wind is swirling through the desert region
       -4 ≤ x ≤ 4, 0 ≤ y ≤ 8.
At a point {x, y} in this region, the velocity of the wind is
       wind[x, y]={0.18 +  0.6  Sin[y],((0.12 y)/(1  +  (x^2)))}.
Here is a picture of the wind as you look down from a helicopter at this desert region:

\(Clear[wind, x, y]; 

    wind[x_, y_] = {0.18 + 0.6\ Sin[y], \(0.12\ y\)\/\(1 + x\^2\)}; 

    \(windplot = 
      Show[Table[
          Arrow[wind[x, y], Tail -> {x, y}], {x, \(-4\), 4}, {y, 0, 8}], 
        Axes -> Automatic, AxesLabel -> {"x", "y"}]; \)\)

A tumbleweed drifting along was spotted at the point {-2, 2}.
Add a plot of the approximate path of this drifting tumbleweed to the plot above.

G.8)  Logistic harvesting revisited

Take the logistic equation
       y'[t]=a y[t] (1- y[t]/b) 
with 0 < a and 0 < b.
And start with 0 < y[0] < b. 
As t advances from 0, you are guaranteed that y[t] grows with some pep until y[t] gets near b.  
Once y[t] gets near b, then y[t] settles into global scale with
       (lim_(t  ->  ∞)) y[t] = b. 
Check it out with with a = 0.23 and b = 13:

a = 0.23; b = 13; \(Clear[DEField, t, y]; 
DEField[t_, y_] = {1, a\ y\ \((1 - y\/b)\)}; 

    \(scalefactor = 2.5; \)\)\(DEFieldplot = 
      Table[Arrow[DEField[t, y], Tail -> {t, y}, 
          ScaleFactor -> scalefactor], {t, 0, 52, 52\/12}, {y, 0, b, 
          b\/10}]; 
bline = Graphics[{Red, Line[{{0, b}, {52, b}}]}]; 


    Show[DEFieldplot, bline, Axes -> Automatic, 
      AxesLabel -> {"t", "y"}, AspectRatio -> 1\/2]; \)

You can see that the trajectories start low on the left and are sucked onto the line y = b plotted at the very top:

endtime = 52; starter = 2; Clear[t, y, Derivative, fakey];
solution=NDSolve[{y'[t]==a y[t] (1 - y[t]/b),y[0] == starter},y[t],{t, 0, endtime}];
fakey[t_] = 
        y[t] /. solution[[1]]; \(solutionplot = 
      Plot[{fakey[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; 


    Show[DEFieldplot, solutionplot, bline, Axes -> Automatic, 
      AxesLabel -> {"t", "y"}, AspectRatio -> 1\/2, 
      DisplayFunction -> $DisplayFunction]; \)

Try it again with y[0] = 5:

starter = 5; Clear[t, y, Derivative, fakey];
solution=NDSolve[{y'[t]==a y[t] (1 - y[t]/b),y[0] == starter},y[t],{t, 0, endtime}];
fakey[t_] = 
        y[t] /. solution[[1]]; \(solutionplot = 
      Plot[{fakey[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; 


    Show[DEFieldplot, solutionplot, bline, Axes -> Automatic, 
      AxesLabel -> {"t", "y"}, AspectRatio -> 1\/2, 
      DisplayFunction -> $DisplayFunction]; \)

You get a reasonable interpretation of the logistic differential equation     
       y'[t]=a y[t] (1- y[t]/b) with 0 < y[0] < b
by imagining that y[t] is number of catfish (in thousands) in a given lake on a catfish farm t weeks after the lake was stocked with y[0] catfish.  As time goes on, the catfish population increases until it reaches its steady-state population of b catfish.
But the catfish farmer doesn't grow catfish as pets; the farmer is in business to harvest catfish and to sell them so that hungry folks can fry them up and then wash them down with a couple of cold beers or iced teas.

G.8.a.i)

Measure time in weeks, assume that the farmer wants to harvest r fish per week, and explain why
       y'[t]=a y[t] (1- y[t]/b)- r
lays the base for a reasonable model.

G.8.a.ii)

Go with the model
       y'[t]=a y[t] (1- y[t]/b)- r,
where t measures the number of weeks that harvesting goes on, a = 0.23, b = 13 (thousand fish), and the weekly harvest is r = 0.5 (thousand fish).
Here is the associated DEField:

a = 0.23; b = 13; r = 0.5; \(Clear[DEField, t, y]; 

    DEField[t_, y_] = {1, a\ y\ \((1 - y\/b)\) - r}; 
\(scalefactor = 2.5; 
    \)\)\(DEFieldplot = 
      Table[Arrow[DEField[t, y], Tail -> {t, y}, 
          ScaleFactor -> scalefactor], {t, 0, 52, 52\/12}, {y, 0, b, 
          b\/10}]; 
bline = Graphics[{Red, Line[{{0, b}, {52, b}}]}]; 


    Show[DEFieldplot, bline, Axes -> Automatic, 
      AxesLabel -> {"t", "y"}, AspectRatio -> 1\/2]; \)

Examine the vector field and use what you see to estimate the smallest number of fish the lake must contain at the beginning (t = 0) to sustain the weekly harvest of r = 0.5 (thousand fish) for a full 52 consecutive weeks.

G.8.a.iii)

Here's what happens with a harvest of r = 1.5 (thousand fish) per week from the same lake:

r = 1.5; \(Clear[DEField, t, y]; 

    DEField[t_, y_] = {1, a\ y\ \((1 - y\/b)\) - r}; 
\(scalefactor = 2.5; 
    \)\)\(DEFieldplot = 
      Table[Arrow[DEField[t, y], Tail -> {t, y}, 
          ScaleFactor -> scalefactor], {t, 0, 52, 52\/12}, {y, 0, b, 
          b\/10}]; 
bline = Graphics[{Red, Line[{{0, b}, {52, b}}]}]; 


    Show[DEFieldplot, bline, Axes -> Automatic, 
      AxesLabel -> {"t", "y"}, AspectRatio -> 1\/2]; \)

Harvesting at the level r = 1.5 (thousand fish) per week is a very bad idea because all the trajectories that start on the far left are going negative before the end of the 52-week period under study.  This means that before the end of the 52-week period, all the fish will be gone and harvesting will have to stop no matter how many fish you start out with.
Your job is to use DEField plots to estimate the largest possible weekly harvest rate r (in thousands) and a number 
       y[0] = starter 
so that if the lake starts out containing starter fish, 
then the harvest can continue for the full 52-week period.
Once you have your estimates for r and y[0] = starter, show the resulting DEField together with the plot of the solution of 
       y'[t]=a y[t] (1- y[t]/b)- r 
with y[0] = starter.

If you can't come up with a grand theory, prospect a little bit.  What happens with r = 1.25? How about r = 0.95? 
Remember you are after an estimate, not an exact number.

G.9)  Water flow with spigots and drains

G.9.a)

Here is water flowing from left to right:

\(Clear[currentField, x, y]; 
currentstrength = 0.51; 

    currentField[x_, y_] = currentstrength\ {1, 0}; 


    currentfieldplot = 
      Table[Arrow[currentField[x, y], Tail -> {x, y}], {x, \(-4\), 4}, {
          y, \(-2\), 3}]; 

Show[currentfieldplot, Axes -> Automatic]; 
    \)

Here is water flowing into a drain at {0.5, 0.5}:

\(Clear[drainField, x, y]; 
{a, b} = {0.5, 0.5}; 

    \(drainstrength = 0.57; \)\)drainField[x_, y_] = 
      -((drainstrength {x - a, y - b})/((((x - a))^2) + 
              (((y - b))^2))); \(drainfieldplot = 
      Table[Arrow[drainField[x, y], Tail -> {x, y}], {x, \(-4\), 4}, {y, 
          \(-2\), 3}]; 

    drain = Graphics[{Blue, PointSize[0.04], Point[{a, b}]}]; 


    Show[drainfieldplot, drain, 
	AxesLabel -> {"x", "y"}, 
      Axes -> Automatic]; \)

Here is water flowing out of a spigot at {-1.5, 0.5}:

\(Clear[spigotField, x, y]; 
{aa, bb} = {\(-1.5\), 0.5}; 

    \(spigotstrength = 0.67; \)\)spigotField[x_, y_] = 
      ((spigotstrength {x - aa, y - bb})/((((x - aa))^2) + 
          (((y - bb))^2))); \(spigotfieldplot = 
      Table[Arrow[spigotField[x, y], Tail -> {x, y}], {x, \(-4\), 4}, {
          y, \(-2\), 3}]; 

    spigot = Graphics[{Red, PointSize[0.04], Point[{aa, bb}]}]; 


    Show[spigotfieldplot, spigot, 
	AxesLabel -> {"x", "y"}, 
      Axes -> Automatic]; \)

Here is water flowing under the influence of all three together:

\(Clear[combinationField, x, y]; 

    combinationField[x_, y_] = 
      currentField[x, y] + drainField[x, y] + spigotField[x, y]; 


    combinationfieldplot = 
      Table[Arrow[combinationField[x, y], Tail -> {x, y}], {x, \(-4\), 
          4}, {y, \(-2\), 3}]; 


    Show[combinationfieldplot, drain, spigot, 
	
      AxesLabel -> {"x", "y"}, Axes -> True]; \)

Experiment with different strengths and positions of the sources and sinks and report on your observations.
What happens you make the spigot strength equal to the drain strength, and you put the drain and the spigot close together?
How big does the current strength have to be in relation to the spigot strength and the drain strength to make the influence of the drain and the spigot negligible?

G.10)  2D Electrical fields

G.10.a)

Plot and describe the two dimensional electric field resulting from equal charges of the same sign at {1.5, 0} and {-1.5, 0}.
Include a plot of the unit vector field.

Make sure that your plot does not ask for the vectors at the points {1.5, 0} and {-1.5, 0}.

G.10.b) Dipoles

A dipole can be approximated by two charges of the same magnitude but opposite sign separated by a small distance.
Dipoles are especially important in atomic theory.
Richard Feynman explains it this way:

"Although an atom or molecule remains neutral in an external electronic field, there is a tiny separation of its positive and negative charges, and it becomes a microscopic dipole."

Plot and describe the electric field resulting from a positive charge at {0.01, 0} and the opposite charge of the same magnitude at {-0.01, 0}.
Include a plot of the unit vector field.

All hell breaks loose near {0, 0}, so make sure your plot does not ask for the vectors too close to {0, 0}.

G.11) Gradient fields for max-min, Hamiltonian fields for level curves, and implicitly defined functions

G.11.a.i) Gradient systems

Given a function f[x, y], you get the gradient of f[x, y] by going with
       gradf[x, y] = {(∂_x) f[x, y], (∂_y) f[x, y]}.
This is the same as saying 
       ∇f[x, y]  =  {∂f[x, y]/∂x, ∂f[x, y]/∂y}
        .

To get the gradient field you just go with
       Field[x, y] = gradf[x[t], y[t]], 
which is the same as putting 
       Field[x, y] = {∂f[x, y]/∂x, 
          ∂f[x, y]/∂y}.
 Here's the gradient field coming from 
       f[x, y] = x Cos[y] + y Sin[x]:

\(Clear[m, n, x, y, t, f, gradf, Field]; 

    f[x_, y_]\  = \ x\ Cos[y]\  + 0.3\ y\^2\ Sin[x]; 

    gradf[x_, y_]\  = {∂\_x f[x, y], ∂\_y f[x, y]}; 
    
{m[x_, y_], n[x_, y_]}\  = \ gradf[x, y]; 


    Field[x_, y_]\  = \ {m[x, y], n[x, y]}\)

Here's a trajectory in this gradient field of f[x, y] starting at a random point:

{a, b} = {Random[Real, {-1, 1}], Random[Real, {-1, 1}]};
starterpoint = {a, b};

Clear[x, y, t];
equationx=x'[t]==m[x[t], y[t]];
equationy=y'[t]==n[x[t], y[t]];starterx = x[0] == a; startery = y[0] == b; endtime = 7; \(approxsolutions = 
      NDSolve[{equationx, equationy, starterx, startery}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[trajectory]; 

    trajectory[t_] = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; 


    Show[starterplot, trajectoryplot, PlotRange -> All, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction, 
      AspectRatio -> 1\/GoldenRatio]; \)

Agree that {xtrajectory[t], ytrajectory[t]} specifies the point on this trajectory at time t. 
Here's what happens to f[xtrajectory[t], ytrajectory[t]] as t advances from 0:

Clear[xtrajectory, ytrajectory]; \({xtrajectory[t_], ytrajectory[t_]} = {x[t], y[t]} /. 
        approxsolutions[[1]]; 


    Plot[f[xtrajectory[t], ytrajectory[t]], {t, 0, endtime}, 
      PlotStyle -> {{Red, Thickness[0.015]}}, 
      AxesLabel -> {"t", "f[xtrajectory[t],ytrajectory[t]]"}]; 
    \)

Going up.  Use your knowledge of the fact that when you plot gradf[x, y] with tail at {x, y}, the direction of gradf[x, y] is the direction of greatest initial increase in f[x, y] to explain why this curve had no choice but to go up.
Explain this bold statement: 
No matter what function f[x, y] you go with, if 
       {xtrajectory[t], ytrajectory[t]} 
parameterizes a trajectory in the gradient field of f[x, y], then you are guaranteed that f[xtrajectory[t], ytrajectory[t]] goes up as t goes up.

G.11.a.ii)

Here's a new vector field:

\(Clear[m, n, x, y, t, Field]; 
m[x_, y_]\  = \ \ 0.2\ x\  + \ 1.5\ y; 

    n[x_, y_]\  = \ \(-0.5\)\ x\ \  - \ 0.2\ y; 


    Field[x_, y_]\  = \ {m[x, y], n[x, y]}; \)

And one of its trajectories:

{a, b} = {-4, -3}; starterpoint = {a, b};

Clear[x, y, t];
equationx=x'[t]==m[x[t], y[t]];
equationy=y'[t]==n[x[t], y[t]];starterx = x[0] == a; startery = y[0] == b; endtime = 8; \(ndssol = 
      NDSolve[{equationx, equationy, starterx, startery}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[trajectory]; 

    trajectory[t_] = {
        x[t] /. ndssol[[1]], 
        y[t] /. ndssol[[1]]}; 

    trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; 


    Show[starterplot, trajectoryplot, PlotRange -> All, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

An elliptical closed curve.
How does the shape of this trajectory signal (with no extra calculation) that the given system is not a gradient system?

Remember from part i) above:
No matter what function f[x, y] you go with, if 
       {xtrajectory[t], ytrajectory[t]} 
parameterizes a trajectory in the gradient field of f[x, y], then you are guaranteed that 
       f[xtrajectory[t], ytrajectory[t]]
goes up as t goes up.

G.11.b.i)

Here's the gradient of 
       f[x, y] = (x^3) +  (E^(3 y)) - 3 x (E^y):

\(Clear[f, gradf, x, y]; 
f[x_, y_] = x\^3 + E\^\(3\ y\) - 3\ x\ E\^y; 

    
gradf[x_, y_] = {∂\_x f[x, y], ∂\_y f[x, y]}\)

Here's the corresponding gradient field:

\(\(Clear[m, n, x, y, t, Field]; 
    
{m[x_, y_], n[x_, y_]}\  = \ gradf[x, y]; 


    Field[x_, y_]\  = \ gradf[x, y]\ \ 
\)\)

You can see that {1, 0} is a point at which gradf[x, y] = {0, 0}:

\({xcritical, ycritical}\  = \ {1, 0}; 
gradf[xcritical, ycritical]\)

Look at the level curves of f[x, y] in the vicinity of {xcritical, ycritical}:

h = 0.2; levelcurves = 
      ContourPlot[
        Evaluate[f[x, y]], {x, xcritical - h, xcritical + h}, {y, 
          ycritical - h, ycritical + h}, ContourSmoothing -> Automatic, 
        AxesLabel -> {"x", "y"}];

Throw in the flow of this gradient field:

\(flowplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Red, 
          ScaleFactor -> 0.1, HeadSize -> 0.02], {x, xcritical - h, 
          xcritical + h, h\/8}, {y, ycritical - h, ycritical + h, h\/8}]; 

    equilibplot = 
      Graphics[{PointSize[0.04], Gold, Point[{xcritical, ycritical}]}]; 


    Show[levelcurves, flowplot, equilibplot, Axes -> True, 
      AxesOrigin -> {xcritical, ycritical}, 
      AxesLabel -> {"x", "y"}, 
      PlotRange 
        -> {{xcritical - h, xcritical + h}, {ycritical - h, 
            ycritical + h}}, 
      PlotLabel -> "Actual flow with level curves"]; \)

What does all this graphical information give you about what happens to f[x, y] at the plotted critical point {xcritical, ycritical}?

G.11.b.ii)

Go with this function:

\(Clear[f, gradf, x, y]; 

    f[x_, y_] = \(x + 2  y\)\/\(1 + x\^2 + y\^2\)\)

Say how you know in advance that this function has a global maximizer and a global minimizer.

G.11.b.iii)

Stay with the same function and calculate its gradient:

\(Clear[f, gradf, x, y]; 
f[x_, y_] = \(x + 2  y\)\/\(1 + x\^2 + y\^2\); 
    

gradf[x_, y_] = {∂\_x f[x, y], ∂\_y f[x, y]}\)

Look at this:

Chop[N[Solve[gradf[x,y] == {0,0}]]]

This gives you two candidates for global maximizers and minimizers. Study in the manner of part i) each of these points in effort to determine which is the global maximizer and which is the global minimizer.
Can you come up with an alternate way of making this determination?

G.11.c) Fishing for possible maximizers

Just as you are working on gradients and max-min, an old retired prof (sort of a geezer) sits down at a computer next to you and says, "I don't use Solve commands or anything like that when I am looking for maximizers in a certain region. Instead I just go fishing in the gradient flow. For instance if I am trying to find local max's near {0, 0} of
             f[x, y] = 
        (((E^(-0.3 (((x - 0.4))^2))) Sin[x y  +  0.4])/(1  +  
            (y^2))), 
I throw four corks in like this:

\(Clear[xstarter, ystarter, k]; 
h = 2; 
    
{xstarter[1], ystarter[1]} = {h, 0}; 
    
{xstarter[2], ystarter[2]} = {0, h}; {xstarter[3], ystarter[3]} = {
        \(-h\), 0}; {xstarter[4], ystarter[4]} = {0, \(-h\)}; 
    starterplots = 
      Table[Graphics[{Red, PointSize[0.04], 
            Point[{xstarter[k], ystarter[k]}]}], {k, 1, 4}]; 
    starterlabels = 
      Table[Graphics[{Red, 
            Text["start", {xstarter[k], ystarter[k]}, {0, 1}]}], {k, 1, 
          4}]; 

Show[starterplots, starterlabels, PlotRange -> All, 
      Axes -> True, AxesLabel -> {"x", "y"}]; \)

"Then I enter the function and plot the trajectories in the gradient field that start at the plotted points and then I see where they go:"

Clear[f, gradf, m, n, x, y, t];
f[x_, y_] = 
        (((E^(-0.3 (((x - 0.4))^2))) Sin[x y + 0.4])/(1 + 
            (y^2)));gradf[x_, y_] = {(∂_x) f[x, y], (∂_y) f[x, y]};
{m[x_, y_], n[x_, y_]} = gradf[x, y];
equationx=x'[t]==m[x[t], y[t]];
equationy=y'[t]==n[x[t], y[t]];
endtime = 25;

ndssol1 = 
        NDSolve[{equationx, equationy, x[0] == h, y[0] == 0}, {x[t], y[t]}, {
            t, 0, endtime}];
ndssol2 = 
        NDSolve[{equationx, equationy, x[0] == 0, y[0] == h}, {x[t], y[t]}, {
            t, 0, endtime}];
ndssol3 = 
        NDSolve[{equationx, equationy, x[0] == -h, y[0] == 0}, {x[t], 
            y[t]}, {t, 0, endtime}];
ndssol4 = 
        NDSolve[{equationx, equationy, x[0] == 0, y[0] == -h}, {x[t], 
            y[t]}, {t, 0, endtime}];

Clear[trajectory1, trajectory2, trajectory3, trajectory4];
trajectory1[t_] = {x[t], y[t]} /. 
          ndssol1[[1]];
trajectory2[t_] = {x[t], y[t]} /. 
          ndssol2[[1]];
trajectory3[t_] = {x[t], y[t]} /. 
          ndssol3[[1]];
trajectory4[t_] = {x[t], y[t]} /. 
          ndssol4[[1]];
trajectoryplots = 
        ParametricPlot[{trajectory1[t], trajectory2[t], trajectory3[t], 
            trajectory4[t]}, {t, 0, endtime}, 
          PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
          DisplayFunction -> Identity];
endplots = {
          Graphics[{Blue, PointSize[0.04], Point[trajectory1[endtime]]}], 
          Graphics[{Blue, PointSize[0.04], Point[trajectory2[endtime]]}], 
          Graphics[{Blue, PointSize[0.04], Point[trajectory3[endtime]]}], 
          Graphics[{Blue, PointSize[0.04], Point[trajectory4[endtime]]}]};

Show[starterplots, starterlabels, trajectoryplots, endplots, 
        PlotRange -> All, Axes -> True, 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction];

The old guy goes on to say, "Those blue dots at trajectory1[endtime], trajectory2[endtime], trajectory3[endtime] and trajectory4[endtime] are the points at which the trajectories stalled. So, goldurn it, I'm pretty well convinced that these are worth investigating as maximizers.
Is this guy still with it? 
Does his idea have any merit?. 
If so, explain why.

G.11.d.i) Hamiltonian systems

Given a function f[x, y], you get the gradient system coming from f[x, y] by going with
       m[x, y] = (∂_x) f[x, y] (which is the same as ∂f[x, y]/∂x) 
       n[x, y] = (∂_y) f[x, y] (which is the same as ∂f[x, y]/∂y)

and putting
       {x'[t],y'[t]}={m[x[t], y[t]], n[x[t], y[t]]}.

Given a function f[x, y], you get the Hamiltonian system coming from f[x, y] by going with
       m[x, y] = -(∂_y) f[x, y] (which is the same as -∂f[x, y]/∂y) 
       n[x, y] = (∂_x) f[x, y] (which is the same as ∂f[x, y]/∂x) 
and putting
       {x'[t],y'[t]}={m[x[t], y[t]], n[x[t], y[t]]}.

Here's the Hamiltonian system coming from
       f[x, y] = (y^2) - y Sin[(x/3)]^2:

Clear[m, n, x, y, t, f];
f[x_, y_]  = (y^2) - y Sin[(x/3)]^2;

{m[x_, y_], n[x_, y_]}={-(∂_y) f[x, y],(∂_x) f[x, y]};
hamiltoniansystem  =  
        ({x'[t], y'[t]}  ==  {m[x[t], y[t]], n[x[t], y[t]]});

ColumnForm[Thread[hamiltoniansystem]]

Here's a random trajectory in this Hamiltonian flow:

{a, b} = {Random[Real, {-2, 2}], Random[Real, {-2, 2}]}; starterpoint = {a, b};

Clear[x, y, t];
equationx=x'[t]==m[x[t], y[t]];
equationy=y'[t]==n[x[t], y[t]];starterx = x[0] == a; startery = y[0] == b; \(endtime = 7; 


    ndssol = NDSolve[{equationx, equationy, starterx, startery}, {x[t], 
          y[t]}, {t, 0, endtime}]; 

Clear[trajectory]; 

    trajectory[t_] = {
        x[t] /. ndssol[[1]], 
        y[t] /. ndssol[[1]]}; 

    \(trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; \)\)\(starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; 


    Show[starterplot, trajectoryplot, PlotRange -> All, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction, 
      AspectRatio -> 1\/GoldenRatio]; \)

Agree that {xtrajectory[t], ytrajectory[t]} specifies the point on this trajectory at time t. 
Here's what happens to f[xtrajectory[t], ytrajectory[t]]:

Clear[xtrajectory, ytrajectory]; \({xtrajectory[t_], ytrajectory[t_]} = {x[t], y[t]} /. 
        ndssol[[1]]; 


    Plot[f[xtrajectory[t], ytrajectory[t]], {t, 0, endtime}, 
      PlotStyle -> {{Red, Thickness[0.015]}}, 
      AxesLabel -> {"t", "f[xtrajectory[t],ytrajectory[t]]"}, 
      PlotRange -> {f[a, b] - 2, f[a, b] + 2}]; \)

Got any idea why that happened?

G.11.d.ii)

Here's the Hamiltonian system coming from 
       f[x, y] = (x^2) - 5 Sin[x] y + 4 (y^2):

Clear[f, gradf, m, n, x, y];
f[x_, y_]  = (x^2) - 5 Sin[x] y + 4 (y^2);
{m[x_, y_], n[x_, y_]} = {-(∂_y) f[x, y],(∂_x) f[x, y]};
hamiltoniansystem  =  
        ({x'[t], y'[t]}  ==  {m[x[t], y[t]], n[x[t], y[t]]});

ColumnForm[Thread[hamiltoniansystem]]

Here's a trajectory for this system shown with a contour plot of some of the level curves of f[x, y]:

{a, b} = {Random[Real, {-5, 5}], Random[Real, {-3, 3}]}; starterpoint = {a, b};

Clear[x, y, t];
equationx=x'[t]==m[x[t], y[t]];
equationy=y'[t]==n[x[t], y[t]];starterx = x[0] == a; startery = y[0] == b; \(endtime = 5; 


    ndssol = NDSolve[{equationx, equationy, starterx, startery}, {x[t], 
          y[t]}, {t, 0, endtime}]; 

Clear[trajectory]; 

    trajectory[t_] = {
        x[t] /. ndssol[[1]], 
        y[t] /. ndssol[[1]]}; 

    \(trajectoryplot = 
      ParametricPlot[trajectory[t], {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; \)\)starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}]; \(levelcurves = 
      ContourPlot[Evaluate[f[x, y]], {x, \(-7\), 7}, {y, \(-4\), 4}, 
        ContourSmoothing -> Automatic, 
        DisplayFunction -> Identity]; 


    Show[levelcurves, starterplot, trajectoryplot, PlotRange -> All, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

Give your opinion about why each plot turns out the way it does.

G.11.d.iii) Defining y in terms of x implicitly via f[x,y] = f[a,b]

Go with:

\(Clear[f, x, y]; 
f[x_, y_] = \(-2.3\)\ Sin[y] + 0.9\ Cos[3.0\ x]\)

When you set 
       f[x, y] = f[0, 0]
and vary x, then y also has to vary to maintain the equality
       f[x, y] = f[0, 0].
So, when you set 
       f[x, y] = f[0, 0],
you are implicitly defining y in terms of of x. 
Mathematica can give you the formula for y in terms of x:

sols = Solve[f[x, y] == f[0, 0], y]; y[x_] = N[y /. sols[[1]], 6]

And you can plot at least part of y[x] as a function of x:

ParametricPlot[{x, y[x]}, {x, 0, 5}, 
      PlotStyle -> {{Magenta, Thickness[0.01]}}, 
      AxesLabel -> {"x", "y[x]"}, AspectRatio -> (1/2)];

At every point {x, y} on this curve, you are guaranteed that 
       f[x, y] = f[0, 0]:

x = Random[Real, {0, 5}]; f[x, y[x]] == f[0, 0]

Plot[f[x, y[x]], {x, 0, 6}, 
      PlotStyle -> {{Blue, Thickness[0.01]}}, 
      AxesLabel -> {"x", "f[x,y[x]]"}, AspectRatio -> (1/2), 
      PlotRange -> {f[0, 0] - 1, f[0, 0] + 1}];

Now go with:

\(Clear[f, x, y]; 

    f[x_, y_] = \(-2.3\)\ Sin[y] + E\^Cos[y]\ Cos[3.0\ x]\)

When you set 
       f[x, y] = f[0, 0], 
you are implicitly defining y in terms of x. But when you go to Mathematica for a formula for y  in terms of x, you get:

sols = Solve[f[x,y] == f[0,0],y]

No dice. No formula is available.
In spite of this apparent setback, you can use a trajectory plot in an appropriate Hamiltonian system to come up with a nice plot of y in terms of x through the implicit definition 
       f[x, y] = f[0, 0].
Do it.