VC.03 The Gradient
 BASICS

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

B.1)  The gradient 
       gradf[x, y] = {D[f[x, y], x], D[f[x, y], y]} 
and the chain rule
       D[f[x[t], y[t]], t]=
	   gradf[x[t], y[t]].{x'[t],y'[t]}

Here is the chain rule for a function of one variable:

\(Clear[f, x, t]; 
D[f[x[t]], \ t]\)

This tells you that 
       df[x[t]]/dt=f'[x[t]] x'[t], 
just as you have known for some time.
Here is the chain rule for a function of two variables:

\(Clear[f, x, t]; 
D[f[x[t], y[t]], \ t]\)

The notation 
       (f^((1, 0)))[x, y] stands for D[f[x, y], x] 
and 
       (f^((0, 1)))[x, y] stands for D[f[x, y], y].
Some folks like to use the notation
       ∂f[x, y]/∂x=(f^((1, 0)))[x, y] = D[f[x, y], x]
       ∂f[x, y]/∂y=(f^((0, 1)))[x, y] = D[f[x, y], x]
So this output gives you three ways of writing the same thing:
       df[x[t], y[t]]/dt 
       =D[f[x[t], y[t]], x]x'[t]+ D[f[x[t], y[t]], y]y'[t]
       =∂f[x[t], y[t]]/∂xx'[t]+ ∂f[x[t], y[t]]/∂yy'[t]
       =(f^((1, 0)))[x[t], y[t]]x'[t]+ (f^((0, 1)))[x[t], y[t]]y'[t].
This is the chain rule for functions of two variables.

B.1.a.i) The gradient

How do you calculate the gradient of a function f[x, y]?

The gradient of f[x, y] is a 2D vector given by:

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

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

In other words, the gradient of a function f[x, y] is given by
       gradf[x, y] = {D[f[x, y], x], D[f[x, y], y]}
       ={∂f[x, y]/∂x,∂f[x, y]/∂x}
        = {(f^((1, 0)))[x, y], (f^((0, 1)))[x, y]}.
The gradient is important enough to carry the whole course from this point until its conclusion.

B.1.a.ii)

How do you calculate the gradient of a function f[x, y, z]?

The gradient of f[x, y, z] is a 3D vector given by:

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

    gradf[x_, y_, z_] = {D[f[x, y, z], \ x], D[f[x, y, z], \ y], 
        D[f[x, y, z], \ z]}\)

In other words, the gradient of a function f[x, y, z] is given by
       gradf[x, y, z] =  
       {D[f[x, y, z], x], D[f[x, y, z], y], D[f[x, y, z], z]}
       ={∂f[x, y, z]/∂x,∂f[x, y, z]/∂y,∂f[x, y, z]/∂z}
        = {(f^((1, 0, 0)))[x, y, z], (f^((0, 1, 0)))[x, y, z], 
          (f^((0, 0, 1)))[x, y, z]}.
 Lots of times you'll see the notation
       ∇f[x, y] = gradf[x, y], or
       ∇f[x, y, z] = gradf[x, y, z]

B.1.b)

Here is Mathematica's calculation of the derivative of f[x[t], y[t]] with respect to t:

\(Clear[f, x, y, t]; 
D[f[x[t], y[t]], \ t]\)

Most folks call this formula the "chain rule." 
You can think of this as a dot product 
       gradf[x[t], y[t]].{x'[t],y'[t]}:

Clear[gradf];
gradf[x_, y_] = {D[f[x, y],  x], D[f[x, y],  y]};
gradf[x[t], y[t]].{x'[t],y'[t]}

Does the chain rule formula
       D[f[x[t], y[t], z[t]], t] 
       =gradf[x[t], y[t], z[t]].{x'[t],y'[t],z'[t]}
work for functions of three variables?

Try it and see.
Here is gradf[x, y, z]:

\(Clear[f, x, y, z, t]; 

    gradf[x_, y_, z_] = {D[f[x, y, z], \ x], D[f[x, y, z], \ y], 
        D[f[x, y, z], \ z]}\)

Here is D[f[x[t], y[t], z[t]], t]:

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

Here is gradf[x[t], y[t], z[t]].{x'[t],y'[t],z'[t]}:

gradf[x[t], y[t], z[t]].{x'[t],y'[t],z'[t]}

Check whether  
       D[f[x[t], y[t], z[t]], t]=
	   gradf[x[t], y[t], z[t]].{x'[t],y'[t],z'[t]}:

D[f[x[t], y[t], z[t]],  t]==gradf[x[t], y[t], z[t]].{x'[t],y'[t],z'[t]}

You bet your sweet ear that the formula
       D[f[x[t], y[t], z[t]], t]=
	   gradf[x[t], y[t], z[t]].{x'[t],y'[t],z'[t]}
works for functions of three variables.
In fact, analogous versions of this formula work for functions of any number of variables.

B.1.c.i)

Here is Mathematica's calculation of 
       D[f[Sin[3  t], (e^(2  t))], t]
for a cleared function f[x, y]:

\(Clear[f, x, y, t]; 
x[t_] = Sin[3\ t]; 
\(y[t_] = E\^\(2\ t\); 

    D[f[x[t], y[t]], \ t]\)\)

Use the fact that
       D[f[x[t], y[t]], t]=gradf[x[t], y[t]].{x'[t],y'[t]} 
to explain the Mathematica output.

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

\(Clear[tan]; 
tan[t] = D[{x[t], y[t]}, \ t]\)

The chain rule says 
       D[f[x[t], y[t]], t] = gradf[x[t], y[t]] . tan[t]:

gradf[x[t], y[t]] . tan[t]

Check:

D[f[x[t], y[t]],  t]

Yes ma'am.

B.1.c.ii)

Here is Mathematica's calculation of 
       D[f[x[t], y[t]], t] 
for f[x, y] = Sin[(x^2) (y^3)] and cleared functions x[t] and y[t]:

\(Clear[f, x, y, t]; 
f[x_, y_] = Sin[x\^2\ y\^3]; 

    D[f[x[t], y[t]], \ t]\)

Use the fact that
       D[f[x[t], y[t], t]=gradf[x[t], y[t]].{x'[t],y'[t]} 
to explain the Mathematica output.

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

\(Clear[tan]; 
tan[t] = D[{x[t], y[t]}, \ t]\)

The chain rule says 
       D[f[x[t], y[t]], t] = gradf[x[t], y[t]] . tan[t]:

gradf[x[t], y[t]] . tan[t]

Check:

Expand[D[f[x[t], y[t]],  t]]

You can do these by hand with no sweat.

B.2)  Level curves, level surfaces and the gradient as normal vector

Anytime you are examining a function, you can learn a lot by plotting some of its level curves.

\(Clear[f, x, y]; 
f[x_, y_] = 5 - \((x\^2 + x\ y + y\^2)\)\)

Here is a plot of some level curves f[x, y] = c for various c's  as selected by Mathematica.  You are looking down at the surface z = f[x, y] from the positive z-direction.

a = -1; b = 1; levelcurves = 
      ContourPlot[Evaluate[f[x, y]], {x, a, b}, {y, a, b}, 
        ContourSmoothing -> Automatic, AxesLabel -> {"x", "y"}];

Here is the 3D plot of the surface z = f[x, y]:

surfaceplot = 
      ParametricPlot3D[{x, y, f[x, y]}, {x, a, b}, {y, a, b}, 
        ViewPoint -> CMView, AxesLabel -> {"x", "y", "f[x,y]"}];

The level curves are plots of the shapes of various trips on the surface that keep the height f[x, y] constant for the whole trip. 
Here is a plot of the shape of the trip that keeps 
       f[x, y] = 4.5:

c = 4.5; Show[levelcurves, Contours -> {c}, 
      ContourShading -> False];

Here is the actual plot of the same trip which lies on the intersection of the surface and the plane z = 4.5:

c = 4.5; \(plane = 
      ParametricPlot3D[{x, y, c}, {x, a, b}, {y, a, b}, 
        PlotPoints -> {2, 2}, DisplayFunction -> Identity]; 


    Show[surfaceplot, plane, ViewPoint -> CMView, 
      DisplayFunction -> $DisplayFunction]; \)

B.2.a) The gradients are perpendicular to level curves

Here's a new function 
       f[x, y] = 2  (x^2) +  x y +  (y^2) 
together with its gradient:

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

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

Here is a part of the level curve 
       f[x, y] = f[1.5, 2]
shown with gradf[1.5, 2] plotted with its tail at {1.5, 2}:

levelcurve = 
      ContourPlot[f[x, y], {x, 0, 3}, {y, 1, 4}, 
        Contours -> {f[1.5, 2]}, ContourShading -> False, 
        DisplayFunction -> Identity]; point = {1.5, 2}; \(gradient = 
      Arrow[gradf[1.5, 2], Tail -> point, VectorColor -> Red]; 


    Show[levelcurve, gradient, PlotRange -> All, 
      AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; \)

That gradient vector is perpendicular to the level curve.
Is this just an accident?

This is no accident.

Here's why:
If you parameterize the level curve 
       f[x, y] = c 
with a parameterization {x[t], y[t]} and then you put
       g[t] = f[x[t], y[t]], 
then you get g[t] = c no matter what t is.
Consequently 
       g'[t]=0 
no matter what t is.  But the chain rule says 
       g'[t]=gradf[x[t], y[t]].{x'[t],y'[t]}. 
So
       gradf[x[t], y[t]].{x'[t],y'[t]}=0 
no matter what t is.  This means that the gradient 
       gradf[x[t], y[t]] 
is perpendicular to the tangent vector 
       {x'[t],y'[t]} 
at {x[t], y[t]}.  
As a result, no matter what {x, y} you go to on a level curve 
       f[x, y] = c, 
the gradient gradf[x, y] is perpendicular to the level curve at {x, y}.
Read that again!
Check it out again, going with 
       f[x, y] = (x^2) +  4  (y^2):
Here is a true-scale plot of the level curve 
       f[x, y] = (x^2) +  4  (y^2) = 36 
and some of the gradients gradf[x, y] with tails at some points {x, y} on the level curve. 
The level curve 
       f[x, y] = (x^2) +  4  (y^2) = 36
is the ellipse             
          ((((x/6)))^2) +  ((((y/3)))^2) = 1.

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

    gradf[x_, y_] = {D[f[x, y], \ x], D[f[x, y], \ y]}; 
Clear[t]; 
    
{x[t_], y[t_]} = {6\ Cos[t], 3\ Sin[t]}; 

    \(levelcurve = 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2\ π}, 
        PlotStyle -> {{Blue, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; \)\)scalefactor = 0.4; \(gradients = 
      Table[Arrow[gradf[x[t], y[t]], Tail -> {x[t], y[t]}, 
          VectorColor -> Red, \ ScaleFactor -> scalefactor], {t, 0, 
          2\ π, \(2\ π\)\/12}]; 

    \(Show[levelcurve, gradients, AspectRatio -> Automatic, 
      AxesLabel -> {"x", "y"}, PlotRange -> All, 
      DisplayFunction -> $DisplayFunction]; \)\)

Prependicular as all get-out.
It will work any time and place.

B.2.b)

Does this work in three dimensions?

Check it out.
Here is part of the level surface
       f[x, y, z] = z -  x y = 1.

\(Clear[f, gradf, x, y, z]; 
f[x_, y_, z_] = z - x\ y; 

    \(gradf[x_, y_, z_] = {D[f[x, y, z], \ x], D[f[x, y, z], \ y], 
        D[f[x, y, z], \ z]}; 
Solve[f[x, y, z] == 1, z]\)\)

levelsurface = 
      ParametricPlot3D[{x, y, 1 + x y}, {x, -1, 1}, {y, -1, 1}, 
        ViewPoint -> CMView, BoxRatios -> Automatic, 
        AxesLabel -> {"x", "y", "z"}];

You can get one of the many curves on this surface by taking the surface plotting function {x, y, 1 +  x y} and setting 
       x = Cos[t]/2 and y = Sin[t]. 
Here come a lot of gradients with tails on this curve:

\(Clear[x, y, z, t]; 
    
{x[t_], y[t_], 
        z[t_]} = {x, y, 1 + x\ y} /. {x -> Cos[t]\/2, y -> Sin[t]}; 
    
\(gradients = 
      Table[Arrow[gradf[x[t], y[t], z[t]], Tail -> {x[t], y[t], z[t]}, 
          VectorColor -> Red], {t, 0, 2\ π, \(2\ π\)\/12}]; 

    \(Show[levelsurface, gradients]; \)\)\)

The gradients are all perpendicular to the surface–just as normal as Beaver, Wally, June, and Ward.

B.3)  The gradient points in the direction of greatest initial increase
The negative gradient points in the direction of greatest initial decrease

B.3.a)

Here is a function and its gradient:

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

    f[x_, y_] = \(-5\)\ x\^2 + \(x\^4\ y\)\/4 + 2\ y + 10; 

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

You are sitting at:

{a, b} = {1.5, 1.8}

At this point the function's value is:

f@@{a, b}

f[a, b]

Here is a plot of part of the level curve f[x, y] = f[a, b]:

levelcurve = 
      ContourPlot[f[x, y], {x, a - 2, a + 2}, {y, b - 2, b + 2}, 
        Contours -> {f[a, b]}, ContourShading -> False, 
        ContourSmoothing -> True, DisplayFunction -> Identity]; 
    \(labels = {Graphics[{PointSize[0.07], Point[{a, b}]}], 
        Graphics[Text["{a,b}", {a, b}, {\(-1.5\), 0}]]}; 


    Show[levelcurve, labels, AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)

You can leave the {a, b} in the direction of any point on the circle of radius 1 centered at {a, b}.

circle = Graphics[{Red, Thickness[0.01], Circle[{a, b}, 1]}]; Show[levelcurve, labels, circle, AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction];

Which direction should you go to get the greatest initial increase of the function?
Which direction should you go to get the greatest initial decrease of the function?

Two unit vectors come to mind.  They are unit tangent and unit normal vectors to the level curve at {a, b}.  One unit normal is the unit vector in the direction of the gradient.

unitgradnormal = gradf[a, b]/(Sqrt[gradf[a, b] . gradf[a, b]])

You can get a unit tangent vector by switching the components and hitting one of them with a minus one.

unittan = {unitgradnormal[[2]], 
        -unitgradnormal[[1]]}

setup = 
      Show[levelcurve, labels, circle, 
        Arrow[unitgradnormal, Tail -> {a, b}], 
        Arrow[unittan, Tail -> {a, b}], AspectRatio -> Automatic, 
        DisplayFunction -> $DisplayFunction];

You can leave the point in the direction of any vector with tail at {a, b} and tip on the circle of radius 1 centered at {a, b}.
Here are some possibilities:

Table[
      Show[setup, 
        Arrow[Cos[t] unitgradnormal + Sin[t] unittan, Tail -> {a, b}], 
        DisplayFunction -> $DisplayFunction, 
        AspectRatio -> Automatic], {t, 0, 2 π - (π/4), 
        (π/4)}];

In the movie, you saw the tips of the vectors
       Cos[t] unitnormal +  Sin[t] unittan
with tails at {a, b} sweep out the circle as t went from 0 to 2  π. The question is how to set t so that the vector
       Cos[t] unitnormal + Sin[t] unittan
with tail at {a, b} points in the direction of greatest initial change of f[x, y] as you leave {a, b}.
Think for a second and you will probably agree that going in the direction of unittan is a bad idea because it snuggles against the level curve and f[x, y] doesn't change at all on this level curve.  
In fact any direction that has a nonzero tangential component is a poor choice because the tangential component won't do much for changing the function.
To get the greatest initial change in f[x, y] you should leave {a, b} in the direction of 
       Cos[t] unitnormal +  Sin[t] unittan 
with t set so that the unittan term is zeroed out.  This means you set t = 0 or π and this leaves you the choice of directions:
       unitnormal 
or in the direction of 
       -unitnormal.
Because the unitnormal points in the same direction as the gradient, to get the greatest initial change in f[x, y] you should leave {a, b} in the direction of either  
       gradf[a, b]  or -gradf[x, y].
Let's see which one:

\(Clear[s]; 
    
{f@@\(({a, b} + s\ gradf[a, b])\), f@@{a, b}, 
        f@@\(({a, b} + s\ \((\(-gradf[a, b]\))\))\)} /. s -> 0.1\)

{f@@({a, b} + s gradf[a, b]), f@@{a, b}, 
        f@@({a, b} + s (-gradf[a, b]))} /. s -> 0.01

{f@@({a, b} + s gradf[a, b]), f@@{a, b}, 
        f@@({a, b} + s (-gradf[a, b]))} /. s -> 0.001

{f@@({a, b} + s gradf[a, b]), f@@{a, b}, 
        f@@({a, b} + s (-gradf[a, b]))} /. s -> 0.0001

The gradient points in the direction of greatest initial increase. 
The negative gradient points in the direction of the greatest initial decrease.
This will happen for any function at any point at which the gradient is not {0, 0}.

B.4)  Using linearizations to help to explain the chain rule 
       D[f[x[t], y[t]], t]=
	  gradf[x[t], y[t]].{x'[t],y'[t]}

Here is a simple curve in two dimensions:

\(Clear[f, x]; 
f[x_] = Sin[x]; 

    \(fplot = 
      Plot[f[x], {x, \(-3\), 3}, PlotStyle -> {{Blue, Thickness[0.01]}}, 
        AxesLabel -> {"x", "y"}]; \)\)

Here's the linearization of f[x] at {a, f[a]}:

Clear[linearf, a];
linearf[x_, a_]=f[a]+f'[a] (x - a)

And a plot in the case a = 0.5:

\(Clear[a, linplot]; 

    linplot[a_] := 
      Plot[linearf[x, a], {x, a - 0.5, a + 0.5}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; 

    pointplot[a_] = Graphics[{PointSize[0.03], Point[{a, f[a]}]}]; 


    Show[fplot, linplot[0.5], pointplot[0.5], 
      DisplayFunction -> $DisplayFunction]; \)

Here's a little movie:

Table[
      Show[fplot, linplot[a], pointplot[a], 
        DisplayFunction -> $DisplayFunction], {a, -2.5, 2.5}];

That's right; linearf[x, a] is nothing but the tangent line at x = a.
To linearize a function f[x, y], you can use the gradient.
The gradient of a function f[x, y] is given by:

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

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

Just as the linearization of a function f[x] at a point x = a is given by
       linearf[x, a]=f[a]+ f'[a](x -  a),
the linearization of a function f[x, y] at a 2D point {a, b} is given by
       linearf[x, y, {a, b}] = 
        f[a, b] +  gradf[a, b] . {x -  a, y -  b}.
To help you stick your mental tongs into the idea behind using the gradient for linearization, look at the following plot of a function f[x, y] and a point {a, b}:

\(Clear[f, x, y]; 
f[x_, y_] = 3 - x\^2 - 0.5\ y\^2; 

    \(actualsurfaceplot = 
      ParametricPlot3D[{x, y, f[x, y]}, {x, \(-1\), 1}, {y, \(-1\), 1}, 
        Axes -> Automatic, 
        AxesLabel -> {"x", "y", "f[x,y]"}, 
        PlotRange -> All, ViewPoint -> CMView, 
        DisplayFunction -> $DisplayFunction]; \)\)

Now throw in the plot of 
       linearf[x, y, {a, b}] 
at the indicated point {a, b, f[a, b]}:

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

    linearf[x_, y_, {a_, b_}] := f[a, b] + gradf[a, b] . {x - a, y - b}; 

    \(linearizedplot[a_, b_] := 
      ParametricPlot3D[{x, y, linearf[x, y, {a, b}]}, {x, a - 0.5, a + 0.5}, {
          y, b - 0.5, b + 0.5}, PlotPoints -> {2, 2}, 
        DisplayFunction -> Identity]; \)\)\(point[a_, b_] := 
      Graphics3D[{Red, PointSize[0.02], Point[{a, b, f[a, b]}]}]; 

    
{a, b} = {\(-0.1\), 0.4}; 

    Show[actualsurfaceplot, linearizedplot[a, b], point[a, b], 
      Axes -> Automatic, AxesLabel -> {"x", "y", "z"}, 
      PlotRange -> All, ViewPoint -> CMView, 
      DisplayFunction -> $DisplayFunction]; \)

The plot of linearf[x, y, {a, b}] is that little bit of plane tangent to the surface at the plotted point {a, b, f[a, b]}.
See what happens for other some other points:

{a, b} = {0.5, 0}; Show[actualsurfaceplot, linearizedplot[a, b], point[a, b], 
      Axes -> Automatic, AxesLabel -> {"x", "y", "z"}, 
      PlotRange -> All, ViewPoint -> CMView, 
      DisplayFunction -> $DisplayFunction];

{a, b} = {-0.25, 0.2}; Show[actualsurfaceplot, linearizedplot[a, b], point[a, b], 
      Axes -> Automatic, AxesLabel -> {"x", "y", "f[x,y]"}, 
      PlotRange -> All, ViewPoint -> CMView, 
      DisplayFunction -> $DisplayFunction];

{a, b} = {0.5, -0.4}; Show[actualsurfaceplot, linearizedplot[a, b], point[a, b], 
      Axes -> Automatic, AxesLabel -> {"x", "y", "z"}, 
      PlotRange -> All, ViewPoint -> CMView, 
      DisplayFunction -> $DisplayFunction];

{a, b} = {0, 0}; Show[actualsurfaceplot, linearizedplot[a, b], point[a, b], 
      Axes -> Automatic, AxesLabel -> {"x", "y", "z"}, 
      PlotRange -> All, ViewPoint -> CMView, 
      DisplayFunction -> $DisplayFunction];

Play with some other choices of a and b until you are satisfied.

B.4.a) Why the linearized plot is tangent to the actual plot

In all the experiments, the plot of linearf[x, y, {a, b}] is tangent to the surface at {a, b, f[a, b]}.
Explain why this will always happen, no matter what function f[x, y] and point {a, b, f[a, b]} you go with.

One possible explanation is that at the point {a, b, f[a, b]} the surface 
       z = f[x, y] 
and the plane
       z = 
        linearf[x, y, {a, b}] = 
          f[a, b] +  gradf[a, b] . {x - a, y - b}.  
have have parallel normal vectors.
To investigate this, run with cleared functions to see what these normal vectors are:

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

    linearf[x_, y_, {a_, b_}] = f[a, b] + gradf[a, b] . {x - a, y - b}\)

An xyz-equation of the plot of 
       z = linearf[x, y, {a, b}] 
is
       (f^((1, 0)))[a, b] (x -  a) +  
          (f^((0, 1)))[a, b] (y -  b) -  (z -  f[a, b]) = 0
        .
This is the equation of a plane with normal vector
       {(f^((1, 0)))[a, b], (f^((0, 1)))[a, b], -1}.
Now check out the normal to the actual surface at a point {a, b, f[a, b]}:

\(Clear[surfacepoint, xcurve, ycurve, normal]; 

    surfacepoint[x_, y_] = {x, y, f[x, y]}; 

    \(xcurve[x_] = surfacepoint[x, b]; \)\)ycurve[y_] = surfacepoint[a, y]; \(tanxcurve[x_] = D[xcurve[x], \ x]; 
tanycurve[x_] = D[ycurve[y], \ y]; 
    
normal = tanxcurve[a]\[Cross]tanycurve[b]\)

This is just the negative of the normal vector  
       {(f^((1, 0)))[a, b], (f^((0, 1)))[a, b], -1}
to the plot of linearf[x, y, {a, b}] found above.
The upshot: 
The plot of linearf[x, y], {a, b}] has no choice but to be tangent to the surface at the point {a, b, f[a, b]}.

B.4.b.i)

You are given a curve {x[t], y[t]} in two dimensions and a function f[x, y]. You can put 
       z[t] = f[x[t], y[t]]
You can also linearize f[x, y] at a point 
     {x[(Subscript[t, 0])], y[(Subscript[t, 0])], f[x[(Subscript[t, 0])], y[(Subscript[t, 0])]]}
and look at the new function
       znew[t] = linearf[x[t], y[t], {x[(Subscript[t, 0])], y[(Subscript[t, 0])]}].
Here's what happens in the specific case with 
       f[x, y] = Sin[x y] 
and 
       {x[t], y[t]} = {t, (t^2)} 
for some sample selections of (Subscript[t, 0]) :

\(Clear[f, x, y, z, t]; 
f[x_, y_] = Sin[x\ y]; 

    \({x[t_], y[t_]} = {t, t\^2}; \)\)z[t_] = f[x[t], y[t]]; actualcurve = 
      Plot[z[t], {t, -0.5, 1.5}, 
        PlotStyle -> {{Blue, Thickness[0.01]}}, 
        AxesLabel -> {"t", "z[t]"}];

Here's what happens when you take (Subscript[t, 0]) = 0.5 and plot the new curve
       znew[t] = linearf[x[t], y[t], {x[(Subscript[t, 0])], y[(Subscript[t, 0])]}]
and show the two plots together:

Clear[a, b, gradf, linearf, point, znew]; \(gradf[x_, y_] = {D[f[x, y], \ x], D[f[x, y], \ y]}; 

    linearf[x_, y_, {a_, b_}] := f[a, b] + gradf[a, b] . {x - a, y - b}; 

    \(znew[t_] = linearf[x[t], y[t], {a, b}]; \)\)\(point[t_] = {
        Graphics[
          Text["{\!\(t\_0\),z[\!\(t\_0\)]}", {t, z[t]}, {0, \(-3\)}]], 
        Graphics[{Red, PointSize[0.03], Point[{t, znew[t]}]}]}; 
t\_0 = 0.5; 
    
{a, b} = {x[t\_0], y[t\_0]}; 

    newcurve = 
      Plot[znew[t], {t, t\_0 - 0.5, t\_0 + 0.5}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; 

    Show[actualcurve, newcurve, point[t\_0], PlotRange -> All, 
      AxesLabel -> {"t", ""}, 
      DisplayFunction -> $DisplayFunction]; \)

Try it again for different selections of (Subscript[t, 0]):

\(t\_0 = 0.2; 
{a, b} = {x[t\_0], y[t\_0]}; 

    newcurve = 
      Plot[znew[t], {t, t\_0 - 0.5, t\_0 + 0.5}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; 

    Show[actualcurve, newcurve, point[t\_0], PlotRange -> All, 
      AxesLabel -> {"t", \*""""}, 
      DisplayFunction -> $DisplayFunction]; \)

Another selection of (Subscript[t, 0]):

\(t\_0 = 1.3; 
{a, b} = {x[t\_0], y[t\_0]}; 

    newcurve = 
      Plot[znew[t], {t, t\_0 - 0.5, t\_0 + 0.5}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; 

    Show[actualcurve, newcurve, point[t\_0], PlotRange -> All, 
      AxesLabel -> {"t", \*""""}, 
      DisplayFunction -> $DisplayFunction]; \)

Another selection of (Subscript[t, 0]):

\(t\_0 = \(-0.2\); 
{a, b} = {x[t\_0], y[t\_0]}; 

    newcurve = 
      Plot[znew[t], {t, t\_0 - 0.5, t\_0 + 0.5}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; 

    Show[actualcurve, newcurve, point[t\_0], PlotRange -> All, 
      AxesLabel -> {"t", \*""""}, 
      DisplayFunction -> $DisplayFunction]; \)

Play with some more selections of (Subscript[t, 0]).
Describe what you see and explain why you are seeing it.

Do it again:

\(t\_0 = 0.8; 
{a, b} = {x[t\_0], y[t\_0]}; 

    newcurve = 
      Plot[znew[t], {t, t\_0 - 0.5, t\_0 + 0.5}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; 

    Show[actualcurve, newcurve, point[t\_0], PlotRange -> All, 
      AxesLabel -> {"t", \*""""}, 
      DisplayFunction -> $DisplayFunction]; \)

The two curves are tangent to each other at their point of contact at {(Subscript[t, 0]), z[(Subscript[t, 0])]}.  When you think about it, you'll have to agree this outcome is, as Lewis Carroll put it, "large as life and twice as natural." 
Reason:
One curve is 
       z[t] = f[x[t], y[t]]; 
the other is
       znew[t] = linearf[x[t], y[t], {x[(Subscript[t, 0])], y[(Subscript[t, 0])]}].
They are tangent at {(Subscript[t, 0]), z[(Subscript[t, 0])]}, as above, because the plane
       z = linearf[x, y, {x[(Subscript[t, 0])], y[(Subscript[t, 0])]}] 
is tangent to the surface 
       z = f[x, y] 
at the point 
        {x[(Subscript[t, 0])], y[(Subscript[t, 0])], f[x[(Subscript[t, 0])], y[(Subscript[t, 0])]]} = {x[(Subscript[t, 0])], y[(Subscript[t, 0])], 
          z[(Subscript[t, 0])]}.
And that's all there is to it!
This will work for any choice of f[x, y], x[t], and y[t] you want.

B.4.b.i) Why the chain rule works

Given 
       z[t] = f[x[t], y[t]], 
Mathematica calculates z'[(Subscript[t, 0])] as follows:

\(Clear[x, y, z, t, f]; 
z[t_] = f[x[t], y[t]]; 

    D[z[t], \ t] /. t -> to\)

Most folks call this formula the "chain rule." 
Explain where the chain rule comes from.

As you saw above, the two curves 
       z[t] = f[x[t], y[t]] 
and         
       znew[t] = linearf[x[t], y[t], {x[(Subscript[t, 0])], y[(Subscript[t, 0])]}] 
are tangent at {(Subscript[t, 0]), z[(Subscript[t, 0])]}.
This tells you that z'[(Subscript[t, 0])]=znew'[(Subscript[t, 0])].
Look at znew[t] for cleared functions f[x, y], x[t], and y[t]:

Clear[a, b, x, y, z, gradf, linearf, znew]; \(gradf[x_, y_] = {D[f[x, y], \ x], D[f[x, y], \ y]}; 

    linearf[x_, y_, {a_, b_}] := f[a, b] + gradf[a, b] . {x - a, y - b}; 

    znew[t_] = 
      linearf[x[t], y[t], {a, b}] /. {a -> x[to], b -> y[to]}\)

You don't need Mathematica to see that znew'[t] is:

D[znew[t],  t]

But
       z'[(Subscript[t, 0])]=znew'[(Subscript[t, 0])]; 
so z'[(Subscript[t, 0])] is given by:

D[znew[t],  t] /. t -> to

This explains Mathematica's calculation:

\(z[t_] = f[x[t], y[t]]; 
D[z[t], \ t] /. t -> to\)

D[f[x[t], y[t]],  t] /. t -> to

This explains where the chain rule comes from.

B.4.b.ii) Why the chain rule works for functions of three variables

Given 
       w[t] = f[x[t], y[t], z[t]],
Mathematica calculates w'[(Subscript[t, 0])] as follows:

\(Clear[x, y, z, w, t, f, to]; 
w[t_] = f[x[t], y[t], z[t]]; 

    D[w[t], \ t] /. t -> to\)

This is the chain rule for functions of three variables. 
Explain where it comes from.

You just linearize
       f[x, y, z] 
at a point {a, b, c} with
            \(linearf[x, y, z, {a, b, c}] = 
	\ \ \ \ 
        f[a, b, c] + \ gradf[a, b, c] . {x - \ a, y - \ b, z - \ c}\) 
and proceed as above.
This time you will not be able to do the plots, but the calculations will hold up.