Divergence

Version 0.4
7/8/92
John B. Schneider

Introduction

This Notebook provides an introduction to the divergence vector operator.  The following commands are defined within this notebook:

   divergence
   cylindricalDivergence
   sphericalDivergence
   nDivergence
   nCylindricalDivergence
   nSphericalDivergence
 plotScalarField
 plotVectorField
 plotVectorFieldAll

Divergence

Let us define a vector function of two variables.  This function might represent any number of things, e.g. velocity and direction of particles in a fluid, the magnitude and direction of the electric field in a planar region of space, or the flow of heat in a sheet of metal.

vF1[x_,y_] := {Sin[x],Cos[y]}

We can examine this function in various ways.  We can plot the x component of the field, the y component, or we can even draw vectors representing the field at points on a uniform grid.  Here is a command to do all that.

plotVectorFieldAll[vF1,{-4,4,.5},{-4,4,.5}]

Notice that there are points in this field where the vectors seem to be converging (such as the point x ≈ 3.1, y ≈ 1.5)  and  other points where the field appears to be diverging (such as the point x ≈ 0, y ≈ ±1.5).  If this field represents fliud flow, then the areas where vectors are converging can be considered ™sinks∫ and the areas where the lines are diverging can be thought of as ™sources.∫  In other words, particles flow from the sources to the sinks.  But how do we measure the strength of a source or a sink?  That's where the divergence is useful.

The divergence of a vector function f(x,y) is written
   ∇• f(x,y) = (◇x  ∂/∂x + ◇y  ∂/∂y)∑(◇x  fx + ◇y  fy)= ∂fx / ∂x + ∂fy / ∂y.

We can calculate the divergence of vF1 symbolically.

divF1 = divergence[vF1,x,y]

Note that this is a scalar function of x and y.  Now let's plot the divergence over the same region we examined above.

plotScalarField[divF1,{-4,4},{-4,4}]

The more positive the divergence, the greater the source strength.  The more negative, the stronger the sink.  When the divergence is zero, the point is neither a source nor a sink.  We can say that there is no change in vector magnitude about a point with zero divergence.

Another example

This function is interesting in that its divergence may not be what you think it is.

vF2[x_,y_] := {x Cos[y], -Sin[y]}

We can look at the x component, the y component, and the ™vector map.∫

plotVectorFieldAll[vF2,{-4,4,.5},{-4,4,.5}]

Where are the sources and sinks?  Looking at the plots above, it is rather hard to tell.  There is obviously a lot of ™circulation∫ in this field, but there aren't any obvious sources or sinks.  Let's compute the divergence.

divergence[vF2,x,y]

Ah ha!  The divergence is zero everywhere.  The question is why?  For this function, 
   ∂fx / ∂x  = ∂(x Cos[y]) / ∂x = Cos[y]
   ∂fy / ∂y   = ∂(±Sin[y]) / ∂y  = ±Cos[y]
This means that at every point, the amount of change in the function for an infinitesimal change in the x direction is exactly opposite to the amount of change that occurs with an equal change in the y direction.  Imagine that we are measuring the vector fields entering and leaving a small box.  We would find that if the vectors decreased as we moved from the left side of the box to the right side (so that it is behaving as something like a sink), the vectors would increase as we moved from the bottom to the top (so that here it is behaving something like a source).  In the limit as this box gets small, the increase and decrease in the field exactly cancel and we can say that the amount of ™flux∫ entering the box is the amount leaving.  Hence it has zero divergence.

Other Examples

Here are some other illustrative examples with which you can explore divergence.  They are provided without comment but with most of the commands you need to start your explorations.  You should spend some time thinking about the results for each example.  Of course, you should feel free to create your own functions.

Function vF3

vF3[x_,y_] := {Exp[-(x/4)^2], 0.5 - (y/4)^2}

plotVectorFieldAll[vF3,{-4,4,.5},{-4,4,.5}]

divF3 = divergence[vF3,x,y]

plotScalarField[divF3,{-4,4},{-4,4}]

Function vF4

vF4[x_,y_] := {2^Sqrt[(x/2)^2+(y/4)^2]-1, Sin[(x y)/3]}

plotVectorFieldAll[vF4,{-4,4,.5},{-4,4,.5}]

divF4 = divergence[vF4,x,y]

plotScalarField[divF4,{-4,4},{-4,4}]

Function vF5

vF5[x_,y_] := {y^3, x y^2}

plotVectorFieldAll[vF5,{-4,4,.5},{-4,4,.5}]

divF5 = divergence[vF5,x,y]

plotScalarField[divF5,{-4,4},{-4,4}]

Initialization Cells

Turn the spell checker off.

Off[General::spell1]

divergence[vF,x,y,(z)] calculates ∇∑(vF(x,y,(z))) symbolically.  Define this in cartesian, cylindrical, and spherical coordinates.

Cartesian coordinates.

divergence[vFunction_,x_,y_] :=
   D[vFunction[x,y],x][[1]] + D[vFunction[x,y],y][[2]];
divergence[vFunction_,x_,y_,z_] :=
   D[vFunction[x,y,z],x][[1]] + 
   D[vFunction[x,y,z],y][[2]] + 
   D[vFunction[x,y,z],z][[3]]

Cylindrical coordinates.

divergenceCylindrical[vFunction_,r_,p_] :=
   1/r D[r vFunction[r,y],r][[1]] + 1/r D[vFunction[r,p],p][[2]];
divergenceCylindrical[vFunction_,r_,p_,z_] :=
   1/r D[r vFunction[r,p,z],r][[1]] + 
   1/r D[vFunction[r,p,z],p][[2]] + 
   D[vFunction[r,p,z],z][[3]]

Spherical coordinates.

divergenceSpherical[vFunction_,r_,t_,p_] :=
   1/r^2 D[r^2 vFunction[r,t,p],r][[1]] + 
   1/(r Sin[t]) D[Sin[t] vFunction[r,t,p],t][[2]] + 
   1/(r Sin[t]) D[vFunction[r,t,p],p][[3]]

nDivergence[vF,x,y,(z)] calculates the value of ∇∑(vF(x,y,(z))) where x, y, and z are given numbers.  The divergence is first calucalated symbolically with dummy variables and then the variables are assigned the values of x, y, and z.

Cartesian coordinates.

nDivergence[vFunction_,x_,y_] := 
   divergence[vFunction,nDiv`x,nDiv`y] /. {nDiv`x -> x, nDiv`y -> y};
nDivergence[vFunction_,x_,y_,z_] := 
   divergence[vFunction,nDiv`x,nDiv`y,nGrad`z] /.
   {nDiv`x -> x, nDiv`y -> y, nDiv`z -> z}

Cylindrical coordinates.

nDivergenceCylindrical[vFunction_,r_,p_] :=
   divergenceCylindrical[vFunction,nDiv`r,nDiv`p] /. {nDiv`r -> r, nDiv`p -> p};
nDivergenceCylindrical[vFunction_,r_,p_,z_] :=
   divergenceCylindrical[vFunction,nDiv`r,nDiv`p,nDiv`z] /. 
   {nDiv`r -> r, nDiv`p -> p, nDiv`z -> z};

Spherical coordinates.

nDivergenceSpherical[vFunction_,r_,t_,p_] :=
   divergenceSpherical[vFunction,nDiv`r,nDiv`t,nDiv`p] /.
   {nDiv`r -> r, nDiv`t -> t, nDiv`p -> p}

Command to plot a vector field.  The function is evaluated over a grid and the value at each grid point is displayed as an arrow where the direction is that of the vector and the length is proportional to the vector magnitude.

plotVectorField[function_,{x0_,x1_,dx_},{y0_,y1_,dy_},opts___] :=
Block[{ix,jy,iDum,grad,max,l1,l2,lines,gridSpace},
   (*
    * Turn off the warning messages that are par for the course.
    *)
   Off[Power::infy,ArcTan::indet,Infinity::indet];
   (*
    * Evaluate the function over a grid of points. 
    *)
   l1 = Table[N[function[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 = N[1.5 gridSpace l2 / max];
   (* 
    * Repartition the list into the proper form. 
    *)
   l2 = Partition[l2,jy];
   (*
    * Define the lines for plotting. 
    *)
   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[Power::infy,ArcTan::indet,Infinity::indet];
   (* Show the beast. *)
   Show[Graphics[{lines,axes[x0,x1,y0,y1]}],
      PlotRange -> {{x0,x1},{y0,y1}},
	  Axes -> True,
	  AxesLabel -> {"X","Y"},
      Frame->True,
      AspectRatio->(y1-y0)/(x1-x0),opts]
   ];

plotVectorFieldAll[function_,{x0_,x1_,dx_},{y0_,y1_,dy_}] :=
Block[{x,y,xComponent,yComponent,xyComponent},
   xComponent = Plot3D[function[x,y][[1]],{x,x0,x1},{y,y0,y1},
                       PlotLabel -> "X Component",
                       AxesLabel -> {"X","Y",""}, 
                       DisplayFunction->Identity];
   yComponent = Plot3D[function[x,y][[2]],{x,x0,x1},{y,y0,y1},
                       PlotLabel -> "Y Component",
                       AxesLabel -> {"X","Y",""},
                       DisplayFunction->Identity];
   xyVectors  = plotVectorField[function,{x0,x1,dx},{y0,y1,dy},
                       DisplayFunction->Identity];
   Show[GraphicsArray[{{xComponent,yComponent},{xyVectors}}],
        DisplayFunction->$DisplayFunction]
   ]

plotScalarField[function_,{x0_,x1_},{y0_,y1_},opts___] :=
Block[{densPlot,surfPlot,contPlot},
   densPlot = DensityPlot[function,{x,x0,x1},{y,y0,y1},
                 Axes->True,
                 AxesLabel -> {"X","Y"},
                 DisplayFunction->Identity,opts];
   surfPlot = Plot3D[function,{x,x0,x1},{y,y0,y1},
                 Axes->True,
                 AxesLabel -> {"X","Y",""},
                 DisplayFunction->Identity];
   contPlot = ContourPlot[function,{x,x0,x1},{y,y0,y1},
                 Axes->True,
                 AxesLabel -> {"X","Y"},
                 ContourShading->False,
                 DisplayFunction->Identity,opts];
   Show[GraphicsArray[{{contPlot,densPlot},{surfPlot}}],
      DisplayFunction->$DisplayFunction]
   ]

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_List,p2_List] := {} /; N[p2[[1]]] == 0. && N[p2[[2]]] == 0.;
vectorLines[p1_List,p2_List] := 
Block[{angle,arrow},
   arrow = vectorLength[p2];
   arrow = Min[.08,.3 arrow];
   angle = N[ArcTan[-p2[[1]],-p2[[2]]]];
   {Line[N[{p1,p2+p1}]],
    Line[N[{p1+p2,p1+p2+{arrow Cos[angle+Pi/6],arrow Sin[angle+Pi/6]}}]],
    Line[N[{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}]}
   ];