Gradient

Version 2.2
9/19/91
by John B. Schneider

Introduction

In this NoteBook, the gradient will be examined.  All of the examples will be for two dimensional functions.  However, several of the Mathematica commands have been defined for three dimensions so that you can examine three dimensional functions.  The syntax for using all of these functions can be determined from the definitions given the section entitled Initialization Cells.

Cartesian Coordinates

Let us define a function of x and y.

fxy[x_,y_] := Exp[-(x^2 + y^2)]

This function has a maximum value of unity at the origin and decays exponentially away from the origin.  This is a scalar function Ð at every point (x,y), fxy(x,y) is a number.  The function fxy could represent such scalar functions as temperature or voltage on a plane.  We can plot this function over a section of the plane.

Plot3D[fxy[x,y],{x,-2,2},{y,-2,2},
   PlotPoints -> 21,
   PlotRange -> All,
   Lighting -> True,
   PlotLabel -> "Fig. 1"]

We can also draw the contour plot of this function.  The contour lines trace the points in the xy plane have a constant value (a constant temperature for instance).

ContourPlot[fxy[x,y],{x,-2,2},{y,-2,2},
   PlotPoints -> 21,
   PlotRange -> All,
   PlotLabel -> "Fig. 2"]

Even though we are considering a function in the Cartesian coordinate system, you should keep in mind that we are merely using that function to describe something that is physical.  Assume that the plot shown above is a representation of the temperature on a plane.  The temperature at a fixed physical point in that plane does not change with a change of coordinates system, or a shift or orgin, or a rotation of axes, etc.  The temperature is what it is.  How we choose to represent that temperature distribution is up to us, but the physical reality is not changed by that choice.

     The temperature at the point (x=1.5, y=1.0) is given by:

fxy[1.5,1.0]

In fact, given the mathematical description of the temperature, we can determine the temperature anywhere by plugging in the appropriate coordinates into the formula.  But, given the formula we can also determine other important aspect of this scalar field.  For example, we can determine how much the function changes with changes in x alone or y alone.
     
     We can take the derivative of fxy in the x direction with the assumption that y is fixed.

D[fxy[x,y], x]

This tells us how rapidly the function changes with changes in x.  A large (in magnitude) derivate means the function is changing rapidly and a small derivate means the function is relatively flat.  We can find the change in the y direction assuming that x is fixed.

D[fxy[x,y], y]

These two derivatives can be combined to indicate how the function changes with either a change in x or a change in y (or both) at a particular point.  To do this, we use a vector field.  The x component of the vector field indicates how the scalar field (the temperature) changes with changes in x, while the y component of the vector field indicates how the scalar field changes with changes in y.  This vector field is known as the gradient.
     
     For example, we can calculate the gradient of fxy at the point (1.5,1.0).

nGradient[fxy, (*x value*) 1.5, (*y value*) 1.0]

The gradient at the point (1.5,1.0) shows that increasing the value of x will decrease the temperature.  The decrease is indicated by the negative sign in x component of the gradient.  The slope at that point is negative, and hence an increase in x will bring a decrease in temperature.  Of course, you can see this in the plot of the temperature distribution shown in Fig. 1 above.  If x is positive and increasing, the temperature is decreasing.  (How about if x is negative?)  The y component of the gradient is also negative.  Again, this indicates that increasing the value of y at the point (1.5,1.0) will decrease the temperature.

     Let's assume that we are interested in moving from a point on the plane in such a way that we obtain the greatest increase in temperature for the smallest amount of movement.  If the x component of the gradient is negative, we know that the temperature will decrease if x is increased.  So, if the gradient is negative and x is decreased, the temperature will increase.  The same holds for y.  If the y component of the gradient is negative, temperature will decrease with an increase in y and the temperature will increase with a decrease in y.

    At the point (1.5,1.0), the gradient is (-0.116323, -0.0775484).  Since both the x and y components are negative, the temperature will increase the most by moving in the negative x direction and the negative y direction.  But, what is the correct ratio of  movement in the x direction to movement in the y direction for maximum change in temperature?   Let us say we will move in the direction ×.  We can find our rate of change in that direction by taking the dot product of × and the gradient, i.e. the vector (-0.116323, -0.0775484).  Since × has unit magnitude, the dot product is the magnitude of the gradient time cos(theta), where theta is the angle between the gradient and ×.  This is maximum when × points in the same direction as the gradient (when theta is zero).
     
     The gradient can be calculated symbolically.  This is merely the combination of the derivates with respect to x and y that were found earlier.

gradient[fxy,x,y]

The gradient of a function is a vector field.  Vector fields are difficult to draw.  At each point, a vector exists.  Obviously we cannot draw a vector at every point in the plane.  However, we can draw the vectors at points on a grid.  The following function plots the gradient vector at points on a grid.

showGradient[ fxy, (*x range*) {-2,2}, (*y range*) {-2,2},
   PlotPoints->10]

Feel free to change the PlotPoints, as well as the other parameters, in this function.  You should convince yourself that the gradient vectors that you see make sense.  
     
     It should be noted that the function showGradient normalizes the length of the vectors in the plot to a ªreasonableº value.  Therefore, you cannot compare the length of the vectors shown for one function with those of another.  However, the length of the vectors are consistent within any single plot Ð if one vector is twice as long at one point than the vector at another point, then the gradient is twice as large.

Cylindrical Coordinates

The temperature function given in the Cartesian Coordinate section has radial symmetry.  Therefore, it is a simple function in cylindrical coordinates.

frp[r_,p_] := Exp[-r^2]

We can plot this in cylindrical coordinates.

CylindricalPlot3D[frp[r,p],{r,0,2},{p,0,2Pi},
   BoxRatios -> {1, 1, 0.4},
   PlotRange-> {{-2,2},{-2,2},{0,1}},
   PlotPoints->15]

We can also calculate the gradient of this function symbolically.  Before clicking on the cell below, think about what the gradient should be.  What is the phi component of the gradient, i.e. what is the derivative in the phi direction?

gradientCylindrical[frp,r,p]

Because the function does not with changes in phi, the gradient in the phi direction is zero Ð the derivative in the phi direction is zero. 
     
     To display the gradient field in cylindrical coordinates, the gradient is evaluate at discrete points that are on concentric cirecles.  The gradient vector at each of these points is then drawn.

showGradientCylindrical[frp, (*r range*) {.25,2}, (*phi range*) {0,2Pi,Pi/6},
   PlotPoints->10]

You can see that the gradient at each point is directed towards the center.  As you know by looking at the function itself, this indicates that the way to increase the value of the function most rapidly is to move towards the origin.
     
     The functions showGradientCylindrical and showGradient do not normalize to vectors to the same value.  Therefore, you can not compare the magnitude of the vectors from the two different outputs.

Demonstration Functions

The Cartesian and cylindrical versions of two more functions are given in this section.  The first function, like the function given above, is only a function of the distance from the origin.  The second function varies with changes in phi.  Of course, you can change these functions in any manner you wish.  There is no text in the following cells, but all the commands that are pertinent to the gradient of these functions are given so that you just have to click on them.

Function 1

Cartesian

fxy1[x_,y_] := Sin[ Sqrt[x^2+y^2] ]

Plot3D[ fxy1[x,y], {x,-3,3},{y,-3,3},
   PlotPoints -> 21,
   Lighting -> True,
   PlotRange -> All]

ContourPlot[ fxy1[x,y], {x,-3,3},{y,-3,3},
   PlotPoints -> 21,
   PlotRange -> All]

gradient[fxy1,x,y]

showGradient[ fxy1, (*x range*) {-2,2}, (*y range*) {-2,2},
   PlotPoints -> 10]

Cylindrical

frp1[r_,p_] := Sin[r]

CylindricalPlot3D[ frp1[r,p], {r,0,3},{p,0,2Pi},
   BoxRatios -> {1, 1, 0.4},
   PlotRange -> All,
   PlotPoints -> 21]

gradientCylindrical[frp1,r,p]

showGradientCylindrical[frp1, (*r range*) {.25,3}, (*phi range*) {0,2Pi,Pi/6},
   PlotPoints -> 10]

Function 2

Cartesian

fxy2[x_,y_] := x^3 + y^2

Plot3D[ fxy2[x,y], {x,-2,2},{y,-2,2},
   PlotPoints -> 21,
   Lighting -> True,
   PlotRange -> All]

ContourPlot[ fxy2[x,y], {x,-2,2},{y,-2,2},
   PlotPoints -> 21,
   PlotRange -> All]

gradient[ fxy2, x,y]

showGradient[ fxy2, (*x range*) {-2,2}, (*y range*) {-2,2},
   PlotPoints -> 10]

Cylindrical

frp2[r_,p_] := (r Cos[p])^3 + (r Sin[p])^2

CylindricalPlot3D[ frp2[r,p], {r,0,2},{p,0,2Pi},
   BoxRatios -> {1, 1, 0.4},
   PlotRange -> All,
   PlotPoints -> 15]

gradientCylindrical[frp2,r,p]

showGradientCylindrical[frp2, (*r range*) {.25,2}, (*phi range*) {0,2Pi,Pi/10},
   PlotPoints -> 10]

Initialization Cells

The definitions of the functions used in this NoteBook are given here.  From these definitions, it should not be difficult to determine the calling syntax for the commands not discussed in the text above.  Some of these commands use some sophisticated tricks in Mathematica.  Therefore, you may not be interested in exploring this section in any detail.

Off[General::spell1]

Gradient Functions in Cartesian Coordinates

gradient::usage = "gradient[f,x,y(,z)] calculates the symbolic gradient
 of the function f(x,y) or f(x,y,z).  Cartesian coordinates are assumed.";

gradient[function_,x_,y_] := 
   {D[function[x,y],x], D[function[x,y],y]};
gradient[function_,x_,y_,z_] := 
   {D[function[x,y,z],x], D[function[x,y,z],y], D[function[x,y,z],z]}

nGradient::usage = "nGradient[f,x,y] calculates the value grad f(x,y) where
x and y are numbers.  The gradient is first calculated symbolically with dummy
variables, and then those variables are assigned the values of x and y.  
nGradient[f,x,y,z] does the same thing for three variables.  Cartesian coordinates
are assumed.";

nGradient[function_,x_,y_] := 
   gradient[function,nGrad`x,nGrad`y] /. {nGrad`x -> x, nGrad`y -> y};
nGradient[function_,x_,y_,z_] := 
   gradient[function,nGrad`x,nGrad`y,nGrad`z] /. {nGrad`x -> x, nGrad`y -> y, nGrad`z ->z}

Gradient Functions in Cylindrical Coordinates

gradientCylindrical::usage = "gradientCylindrical[f,r,p(,z)] calculates the
symbolic gradient of the function f(r,p) or f(r,p,z).  Cylindrical coordinates
are assumed.";

gradientCylindrical[function_,r_,p_] := 
   {D[function[r,p],r], 1/r D[function[r,p],p]};
gradientCylindrical[function_,r_,p_,z_] := 
   {D[function[r,p,z],r], 1/r D[function[r,p,z],p], D[function[r,p,z],z]}

nGradientCylindrical::usage = "nGradientCylindrical[f,r,p(,z)] calculates the
numerical gradient of the function f at the point (r,p) or (r,p,z).  Cylindrical
coordinates are assumed.";

nGradientCylindrical[function_,r_,p_] := 
   gradientCylindrical[function,nGrad`r,nGrad`p] /. {nGrad`r -> r, nGrad`r -> p};
nGradientCylindrical[function_,r_,p_,z_] := 
   gradientCylindrical[function,nGrad`r,nGrad`p,nGrad`z] /.
   {nGrad`r -> r, nGrad`r -> p, nGrad`z ->z}

Gradient Functions in Fpherical Coordinates

gradientSpherical::usage = "gradientSpherical[f,r,t,p] calculates the 
symbolic gradient of the function f(r,t,p).  Spherical coordinates are
assumed.";

gradientSpherical[function_,r_,t_,p_] := 
   {D[function[r,t,p],r],
    1/r D[function[r,t,p],t],
    1/(r Sin[t]) D[function[r,t,p],p]}

nGradientSpherical::usage = "nGradientSpherical[f,r,t,p] calculates the 
numerical gradient of the function f at the point (r,t,p).  Spherical
coordinates are assumed.";

nGradientSpherical[function_,r_,t_,p_] := 
   gradientSPherical[function,nGrad`r,nGrad`t,nGrad`p] /.
   {nGrad`r -> r, nGrad`t -> t, nGrad`p -> p}

Functions to Display Vector Fields.

Cartesian coordinates

showGradient::usage = "showGradient[f,{x0,x1(,dx)},{y0,y1(,dy)}] calculates
the gradient of the function f over a grid of points.  The vector at each
of these grid points is shown.  The x coordinate starts at x0 and incrementally
increases by dx until reaching x1.  The y coordinate starts at y0 and
incrementally increases by dy until reaching y1.  The vector length is
normalized so that the maximum length is 1.2 times the grid spacing.  Arrows
are placed at the tips of the vectors.  Two dimensional Cartesian coordinates
are assumed.";

showGradient[function_,{x0_,x1_,dx_},{y0_,y1_,dy_}] :=
Block[{ix,jy,iDum,grad,max,l1,l2,lines,gridSpace},
   (* Calculate the gradient. *)
   grad[x_,y_] = gradient[function,x,y];
   (* Find the gradient over a grid of points. *)
   l1 = Table[N[grad[x,y]],{x,x0,x1,dx},{y,y0,y1,dy}];
   {ix,jy,iDum} = Dimensions[l1];
   (* Remove the indeterminate entries from the grid. *)
   l2 = 
   If[NumberQ[ Flatten[#][[1]] ] && NumberQ[ Flatten[#][[2]] ],
      #,{0.,0.}]& /@ Flatten[l1,1];
   (* Normalize by the max vector length. *)
   max = Max[Map[vectorLength,l2]];
   gridSpace = vectorLength[{dx,dy}];
   l2 = 1.5 gridSpace l2 / max;
   (* Repartition the list into the proper form. *)
   l2 = Partition[l2,jy];
   (* Define the lines for plotting. *)
   Off[ArcTan::ind]; 
   lines = 
   Table[
      vectorLines[{x,y},l2[[Round[(x-x0)/dx]+1,Round[(y-y0)/dy]+1]]],
      {x,x0,x1,dx},{y,y0,y1,dy}];
   On[ArcTan::ind]; 
   (* Find PlotRange coefficients *)
   (* Show the beast. *)
   Show[Graphics[{lines,axes[x0,x1,y0,y1]}],
      PlotRange -> {{x0,x1},{y0,y1}},
      Frame->True,
      AspectRatio->(y1-y0)/(x1-x0)]
   ];
   
showGradient[function_,{x0_,x1_,dx_},{y0_,y1_},opts___] :=
Block[{dy,plotPoints},
   plotPoints = PlotPoints /. Flatten[{opts,Options[Plot3D]}];
   dy = (y1-y0)/(plotPoints-1);
   showGradient[function,{x0,x1,dx},{y0,y1,dy}]
   ];

showGradient[function_,{x0_,x1_},{y0_,y1_,dy_},opts___] :=
Block[{dx,plotPoints},
   plotPoints = PlotPoints /. Flatten[{opts,Options[Plot3D]}];
   dx = (x1-x0)/(plotPoints-1);
   showGradient[function,{x0,x1,dx},{y0,y1,dy}]
   ];

showGradient[function_,{x0_,x1_},{y0_,y1_},opts___] :=
Block[{dx,dy,plotPoints},
   plotPoints = PlotPoints /. Flatten[{opts,Options[Plot3D]}];
   dx = (x1-x0)/(plotPoints-1);
   dy = (y1-y0)/(plotPoints-1);
   showGradient[function,{x0,x1,dx},{y0,y1,dy}]
   ];

Cylindrical coordinates

showGradientCylindrical::usage = 
"showGradientCylindrical[f,{r0,r1(,dr)},{p0,p1(,dp)}] displays the gradient
vector of f at points on concentric circles.  The vector length is normalized
so that the maximum length is 1.5 times the minimum of the radius spacing and
the phi spacing at the maximum radius.  Arrows are placed at the tips of the
vectors.  Cylindrical coordinates are assumed.";

showGradientCylindrical[function_,{r0_,r1_,dr_},{p0_,p1_,dp_}] :=
Block[{i,j,ir,jy,iDum,grad,max,gridSpace,l1,l2,lines,r,p,runit,punit,point},
   (* Calculate the gradient. *)
   grad[r_,p_] = gradientCylindrical[function,r,p];
   (* Find the gradient over a grid of points. *)
   l1 = Table[N[grad[r,p]],{r,r0,r1,dr},{p,p0,p1,dp}];
   {ir,jp,iDum} = Dimensions[l1];
   (* Remove the indeterminate entries from the grid. *)
   l2 = 
   If[NumberQ[ Flatten[#][[1]] ] && NumberQ[ Flatten[#][[2]] ],
      #,{0.,0.}]& /@ Flatten[l1,1];
   (* Normalize by the max vector length. *)
   max = Max[Map[vectorLength,l2]];
   gridSpace = Min[N[dr],N[r1 dp]];
   l2 = 1.5 gridSpace l2 / max;
   (* Repartition the list into the proper form. *)
   l2 = Partition[l2,jp];
   (* Define the lines for plotting. *)
   Off[ArcTan::ind]; 
   xMin = Infinity;
   xMax = -Infinity;
   yMin = Infinity;
   yMax = -Infinity;
   lines = 
   Table[
      point = N[{r Cos[p],r Sin[p]}];
      If [point[[1]] < xMin, xMin = point[[1]]];
      If [point[[1]] > xMax, xMax = point[[1]]];
      If [point[[2]] < yMin, yMin = point[[2]]];
      If [point[[2]] > yMax, yMax = point[[2]]];
      runit = point/r;
      punit = {-Sin[p],Cos[p]};
      i = Round[(r-r0)/dr]+1;
      j = Round[(p-p0)/dp]+1;
      vectorLines[point,
         Flatten[{l2[[i,j,1]] runit + l2[[i,j,2]] punit}]],
      {r,r0,r1,dr},{p,p0,p1,dp}];
   On[ArcTan::ind]; 
   (* Show the beast. *)
      Show[Graphics[{lines,axes[xMin,xMax,yMin,yMax]}],
      PlotRange -> {{xMin,xMax},{yMin,yMax}},
      Frame->True,
      AspectRatio->(yMax-yMin)/(xMax-xMin)
      ]
   ];

showGradientCylindrical[function_,{r0_,r1_},{p0_,p1_,dp_},opts___] :=
Block[{dr,plotPoints},
   plotPoints = PlotPoints /. {opts} /. Options[Plot3D];
   dr = (r1-r0)/(plotPoints-1);
   showGradientCylindrical[function,{r0,r1,dr},{p0,p1,dp}]
   ];

showGradientCylindrical[function_,{r0_,r1_,dr_},{p0_,p1_},opts___] :=
Block[{dp,plotPoints},
   plotPoints = PlotPoints /. {opts} /. Options[Plot3D];
   dp = (p1-p0)/(plotPoints-1);
   showGradientCylindrical[function,{r0,r1,dr},{p0,p1,dp}]
   ];

showGradientCylindrical[function_,{r0_,r1_},{p0_,p1_},opts___] :=
Block[{dr,dp,plotPoints},
   plotPoints = PlotPoints /. {opts} /. Options[Plot3D];
   dr = (r1-r0)/(plotPoints-1);
   dp = (p1-p0)/(plotPoints-1);
   showGradientCylindrical[function,{r0,r1,dr},{p0,p1,dp}]
   ];

Functions used in both coordinate systems

vectorLines::usage = "vectorLines[p1,p2] generates lines that are used
to draw a vector.  p2 (a list of two numbers) gives the horizontal and
vertical dimensions of the vector and p1 defines its origin.  An arrow is
drawn at the tip of the vector.";

vectorLines[p1_,p2_] := 
Block[{angle,arrow},
   arrow = vectorLength[p2];
   arrow = Min[.08,.3 arrow];
   angle = N[ArcTan[-p2[[1]],-p2[[2]]]];
   {Line[{p1,p2+p1}],
    Line[{p1+p2,p1+p2+{arrow Cos[angle+Pi/6],arrow Sin[angle+Pi/6]}}],
    Line[{p1+p2,p1+p2+{arrow Cos[angle-Pi/6],arrow Sin[angle-Pi/6]}}]
    }
   ];

vectorLength::usage = "vectorLength[v] returns the length of the vector v.";

vectorLength[{x_,y_}] := Sqrt[x^2 + y^2];

axes::usage = "axes[x0,x1,y0,y1]  generates the graphics commands needed
to draw xy axes.";

axes[x0_,x1_,y0_,y1_] := 
Block[{tick},
   tick = .02 Min[{x1-x0,y1-y0}];
   {(* x axis *) Line[{{0,y0},{0,y1}}],
   (* y axis *) Line[{{x0,0},{x1,0}}],
   (* x axis tick marks *) Table[Line[{{x,-tick},{x,tick}}],{x,x0,x1}],
   (* y axis tick marks *) Table[Line[{{-tick,y},{tick,y}}],{y,y0,y1}]}
   ];

Cylindrical Plotting Commands

The following commands are needed to draw a three-dimensional plot of a function that uses cylindrical coordinates.  CylindricalPlot3D is a modified version of a routine published by Roman Maeder.

Needs["Graphics`ParametricPlot3D`"]

BeginPackage["CylindricalPlot3D`","Graphics`ParametricPlot3D`"]

CylindricalPlot3D::usage = "CylindricalPlot[f,{r,r0,r1(,dr)},{p,p0,p1(,dp)}]
generates a three-dimensional plot of f as a function of r and p.";

Begin["`Private`"]

Unprotect[CylindricalPlot3D];

CylindricalPlot3D[z_, rlist:{r_,__}, plist:{p_,__}, opts___] :=
Block[{newopts,r4list,p4list,plotPoints,pointer,j},
   plotPoints = PlotPoints /. {opts} /. Options[Plot3D];
   pointer = Position[{opts},PlotPoints->__][[1,1]];
   newopts = 
      Flatten[
         {Table[{opts}[[j]],{j,1,pointer-1}],
          Table[{opts}[[j]],{j,pointer+1,Length[{opts}]}]}
      ];
   r4list = Flatten[{rlist,(rlist[[3]]-rlist[[2]])/plotPoints}];
   p4list = Flatten[{plist,(plist[[3]]-plist[[2]])/plotPoints}];
   ParametricPlot3D[{r Cos[p], r Sin[p], z}, r4list, p4list, newopts]
   ] /; Length[rlist] == 3 && Length[plist] == 3;

CylindricalPlot3D[z_, rlist:{r_,__}, plist:{p_,__}, opts___] :=
Block[{newopts,r4list,p4list,plotPoints,pointer,j},
   plotPoints = PlotPoints /. {opts} /. Options[Plot3D];
   pointer = Position[{opts},PlotPoints->__][[1,1]];
   newopts = 
      Flatten[
         {Table[{opts}[[j]],{j,1,pointer-1}],
          Table[{opts}[[j]],{j,pointer+1,Length[{opts}]}]}
      ];
   r4list = Flatten[{rlist,(rlist[[3]]-rlist[[2]])/plotPoints}];
   ParametricPlot3D[{r Cos[p], r Sin[p], z}, r4list, plist, newopts]
   ] /; Length[rlist] == 3 && Length[plist] == 4;

CylindricalPlot3D[z_, rlist:{r_,__}, plist:{p_,__}, opts___] :=
Block[{newopts,r4list,p4list,plotPoints,pointer,j},
   plotPoints = PlotPoints /. {opts} /. Options[Plot3D];
   pointer = Position[{opts},PlotPoints->__][[1,1]];
   newopts = 
      Flatten[
         {Table[{opts}[[j]],{j,1,pointer-1}],
          Table[{opts}[[j]],{j,pointer+1,Length[{opts}]}]}
      ];
   p4list = Flatten[{plist,(plist[[3]]-plist[[2]])/plotPoints}];
   ParametricPlot3D[{r Cos[p], r Sin[p], z}, rlist, p4list, newopts]
   ] /; Length[rlist] == 4 && Length[plist] == 3;

CylindricalPlot3D[z_, rlist:{r_,__}, plist:{p_,__}, opts___] :=
   ParametricPlot3D[{r Cos[p], r Sin[p], z}, rlist, plist, opts] /;
   Length[rlist] == 4 && Length[plist] == 4;
   
Protect[CylindricalPlot3D];

End[]
EndPackage[];

On[General::spell1]