Divergence

Multivariable Calculus, Math 211, Carleton College, Fall 2023, Joshua R. Davis, jdavis@carleton.edu

Basic divergence examples

vecF = {x, y};
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

vecF = {-x, -y};
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

vecF = {1, 0};
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

vecF = {-y / (x^2 + y^2), x / (x^2 + y^2)};
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

vecF = {y, 0};
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

Harder divergence examples

The vectors in this vector field are all parallel to each other. Does that mean that the divergence is 0?

vecF = {x^2, 0};
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

The vectors in this vector field all point away from the origin. Does that mean that the divergence is positive?

vecF = {x, y} / Sqrt[x^2 + y^2];
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

The vectors in this vector field all point away from the origin. Does that mean that the divergence is positive?

vecF = {x, y} / (x^4 + y^4);
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

The vectors in this vector field all point toward the origin. Does that mean that the divergence is negative?

vecF = -{x, y} / (x^2 + y^2);
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}]
divVecF = Simplify[Div[vecF, {x, y}]]
ContourPlot[divVecF, {x, -2, 2}, {y, -2, 2}]

Divergence of gradient

One way to make a vector field is to take the gradient of a scalar field, right? So let’s do some examples like that.

f = Exp[-2.0 (x - 1)^2 - (y - 1)^2] + 1.5 Exp[-0.4 (x + 1)^2 - (y + 1)^2];
ContourPlot[f, {x, -2, 2}, {y, -2, 2}]
gradF = Grad[f, {x, y}]
VectorPlot[gradF, {x, -2, 2}, {y, -2, 2}]
divGradF = Div[gradF, {x, y}]
ContourPlot[divGradF, {x, -2, 2}, {y, -2, 2}]

And here’s Newtonian gravity with simplified constants. The level sets of f are spheres centered on the origin.

f = (x^2 + y^2 + z^2)^(-1 / 2);
gradF = Grad[f, {x, y, z}]
VectorPlot3D[gradF, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]
divGradF = Simplify[Div[gradF, {x, y, z}]]

Geology again

Here’s a transpressional shear zone, like the deforming rock along the San Andreas fault. I’m being terse because we saw this example earlier in the course.

l11 = 1;
l12 = 4;
l21 = 0;
l22 = -l11;
l = {{l11, l12}, {l21, l22}};
eig = Eigensystem[l];
If[eig[[1, 1]] < eig[[1, 2]],
	u = eig[[2, 1]],
	u = eig[[2, 2]]];
u = l.{u[[1]] / u[[2]], 1};
vecF = If[y > 1, u,
	If[y < -1, -u,
		l.{x, y}]];
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}, VectorPoints->20]

To compute the divergence, Mathematica insists that we break the vector field into its three cases. All have divergence 0 everywhere. This means that, at each point, the mass flux into the point balances the mass flux out of the point. It’s a form of conservation of mass.

vecFNorth = Simplify[vecF, {1 < y}]
Div[vecFNorth, {x, y}]
vecFMiddle = Simplify[vecF, {-1 < y < 1}]
Div[vecFMiddle, {x, y}]
vecFSouth = Simplify[vecF, {y < -1}]
Div[vecFSouth, {x, y}]

Here’s a second example. Is the divergence still 0? If not, then what could that mean, physically?

l11 = 1;
l12 = 4;
l21 = 0;
l22 = -1.5;
l = {{l11, l12}, {l21, l22}};
eig = Eigensystem[l];
If[eig[[1, 1]] < eig[[1, 2]],
	u = eig[[2, 1]],
	u = eig[[2, 2]]];
u = l.{u[[1]] / u[[2]], 1};
vecF = If[y > 1, u,
	If[y < -1, -u,
		l.{x, y}]];
VectorPlot[vecF, {x, -2, 2}, {y, -2, 2}, VectorPoints->20]
vecFTop = Simplify[vecF, {1 < y}]
Div[vecFTop, {x, y}]
vecFMiddle = Simplify[vecF, {-1 < y < 1}]
Div[vecFMiddle, {x, y}]
vecFBottom = Simplify[vecF, {y < -1}]
Div[vecFBottom, {x, y}]

Lotka-Volterra predator-prey model again

You can compute the divergence of any smooth vector field. But sometimes it’s hard to attach practical meaning to the answer. Talk to me if you want to know more.

a = 1.5;
b = 1.0;
c = 3.0;
d = 1.0;
vecF = {a x - b x y, d x y - c y}
VectorPlot[vecF, {x, 0, 10}, {y,0, 10}, VectorPoints->20]
Div[vecF, {x, y}]