\( \DeclareMathOperator{\abs}{abs} \newcommand{\ensuremath}[1]{\mbox{$#1$}} \)

Vector Fields, Gradients, and Potentials

 1 Set-Up

Once again we will be working with vectors, and so we must load the vect library.  We will also build our usual helper functions mag() and hat() as well as redefine the cross product to  the infix operator ##:

(%i5) load("vect")$
mag(u):=sqrt(u.u)$
hat(u):=u/mag(u)$
infix("##")$
u##v:=express(u~v)$

Also note that kill(all) is now off-limits if we do not want to undo everything that we've done above.

 2 Vector Fields as Expressions and Functions

Vector fields are simply functions from either R^2 or R^3 into R^2 or R^3.  That is for a given point (x,y) [or (x,y,z)],
a vector field generates a new point in R^2 or R^3.  We can encode such an object either as an expression:

--> E:[x^2, y^2, z^2];
\[\tag{E}[{{x}^{2}},{{y}^{2}},{{z}^{2}}]\]

And then extract values using subst():

--> subst([x=1,y=2,z=3],E);
\[\tag{%o7} [1,4,9]\]

--- this works but is inefficient --- or as a function:

--> F(x,y,z):= [x^2,y^2,z^2];
F(1,2,3);
\[\tag{%o8} \operatorname{F}\left( x,y,z\right) :=[{{x}^{2}},{{y}^{2}},{{z}^{2}}]\] \[\tag{%o9} [1,4,9]\]

In general, we will prefer the latter approach as it is much more streamlined and flexible.

--> kill(E,F);
\[\tag{%o10} \mathit{done}\]

 3 Gradients

The vect library provides a grad() function that produces 3D gradients of scalar functions of either 2 or 3 variables,
but it's a little goofy.  Let's take a look at it:

--> f(x,y,z):=x^2+y^2+z^2;
grad(f(x,y,z));
\[\tag{%o11} \operatorname{f}\left( x,y,z\right) :={{x}^{2}}+{{y}^{2}}+{{z}^{2}}\] \[\tag{%o12} \operatorname{grad}\left( {{z}^{2}}+{{y}^{2}}+{{x}^{2}}\right) \]

Ok, so we defined f(x,y,z) and it looks like grad() didn't throw an error, but it really doesn't look like what we were expecting.  We've seen that the vect library often doesn't produce exactly what we want (recall ~ and the express() in the cross product?)  So, let's try to express() the gradient:

--> express(grad(f(x,y,z)));
\[\tag{%o13} [\frac{d}{d x} \left( {{z}^{2}}+{{y}^{2}}+{{x}^{2}}\right) ,\frac{d}{d y} \left( {{z}^{2}}+{{y}^{2}}+{{x}^{2}}\right) ,\frac{d}{d z} \left( {{z}^{2}}+{{y}^{2}}+{{x}^{2}}\right) ]\]

This is promising!  So if we express the gradient, we get a symbolic expression that is correct, but not evaluated for the derivatives.  We can get this to evaluate using the ev(,diff) command that we've seen before:

--> ev(express(grad(f(x,y,z))),diff);
\[\tag{%o14} [2 x,2 y,2 z]\]

Again the reason behind why this is so obtuse has to do with the different symbolic levels that expressions can occupy in Maxima.  We have functions, expressions, operations, and compound versions of all three of these things.  Since the programmer of the vect library doesn't know what every user expects to get when they compute grad(), they opt for the safest -- highest level -- object that they can: an unexpressed symbolic expression involving the derivative.

We as users have to pull that expression down to the level of a vector field.  This will be the same issue that we will face when we consider the curl() and the div(), and so to help us, we introduce:

 3.1 The "@" Helper Function

We will write our own helper function called @ that will automatically do the ev(express(),diff) part for us, so that we can streamline our work.  So instead of writing ev(express(grad(f(x,y,z))),diff) to get the gradient, we will instead just write:

@grad(f(x,y,z))

You have to admit, that's much better!  The @ symbol will be something called a "prefix" operator since it will prefix whatever
we want to operate on.  Let's define it here:

(%i7) prefix("@");
@expression:=ev(express(expression),diff);
\[\tag{%o6} @ \] \[\tag{%o7} \operatorname{@ }\left( \mathit{expression}\right) :=\operatorname{ev}\left( \operatorname{express}\left( \mathit{expression}\right) ,\mathit{\neq \{Lisp function\}}\right) \]

Now let's see how we use this:

(%i9) @grad(x^2+y^2);
@grad(x^2+y^2+z^2);
\[\tag{%o8} [2 x,2 y,0]\] \[\tag{%o9} [2 x,2 y,2 z]\]

@ works just fine with gradients of functions too:

(%i13) f(x,y):=x*sin(y);
@grad(f(x,y));
g(x,y,z):=x*exp(y*z);
@grad(g(x,y,z));
\[\tag{%o10} \operatorname{f}\left( x,y\right) :=x \sin{(y)}\] \[\tag{%o11} [\sin{(y)},x \cos{(y)},0]\] \[\tag{%o12} \operatorname{g}\left( x,y,z\right) :=x \operatorname{exp}\left( y z\right) \] \[\tag{%o13} [{{\% e}^{y z}},x z\, {{\% e}^{y z}},x y\, {{\% e}^{y z}}]\]

Note: we could have named @ anything as long as it didn't conflict with operators (like +,-,*,^, etc) that already exist in Maxima.  Further, we will continue to use @ when we apply the curl() and the div() operators later on.

 3.2 Catching the Output of the Gradient

Let's look at f(x,y) from the previous section:

(%i14) f(x,y):=x*sin(y);
\[\tag{%o14} \operatorname{f}\left( x,y\right) :=x \sin{(y)}\]

Let's say that we want to take the gradient and then store it in the vector field gf(x,y).  We do it in the usual way:

(%i15) define(gf(x,y),@grad(f(x,y)));
\[\tag{%o15} \operatorname{gf}\left( x,y\right) :=[\sin{(y)},x \cos{(y)},0]\]

Now we can compute maximum rates of increase -- for example:

(%i16) mag(gf(x,y));
\[\tag{%o16} \sqrt{{{\sin{(y)}}^{2}}+{{x}^{2}}\, {{\cos{(y)}}^{2}}}\]

For a gradient involving three variables:

(%i18) g(x,y,z):=x*exp(y*z);
define(gg(x,y,z),@grad(g(x,y,z)));
\[\tag{%o17} \operatorname{g}\left( x,y,z\right) :=x \operatorname{exp}\left( y z\right) \] \[\tag{%o18} \operatorname{gg}\left( x,y,z\right) :=[{{\% e}^{y z}},x z\, {{\% e}^{y z}},x y\, {{\% e}^{y z}}]\]

We also note that we can "chain" these operations so that they aren't so cumbersome to type and keep track of all the parentheses.  Let's redo the above example by way of chaining:

(%i21) g(x,y,z):=x*exp(y*z); /*Define our function*/
@grad(g(x,y,z))$     /*Take the gradient of the function*/
define(gg(x,y,z),%); /*Define a new function that uses the output of the last calculation (grad)*/
\[\tag{%o19} \operatorname{g}\left( x,y,z\right) :=x \operatorname{exp}\left( y z\right) \] \[\tag{%o21} \operatorname{gg}\left( x,y,z\right) :=[{{\% e}^{y z}},x z\, {{\% e}^{y z}},x y\, {{\% e}^{y z}}]\]

 3.3 Notes About grad()

 3.3.1 Only x,y,z Allowed (sort of)

The way that the vect constructs the gradient requires that we use functions (or expressions) with x, y, and z as variables.  We note that grad() works just fine on expressions:

(%i22) @grad(x^3-2*y^2+exp(z));
\[\tag{%o22} [3 {{x}^{2}},-4 y,{{\% e}^{z}}]\]

But the operator chokes when the variables are not x,y, or z:

(%i23) @grad(s+t+u);
\[\tag{%o23} [0,0,0]\]

If we look at the symbolic underpinnings of grad():

(%i24) express(grad(s+t+u));
\[\tag{%o24} [\frac{d}{d x} \left( u+t+s\right) ,\frac{d}{d y} \left( u+t+s\right) ,\frac{d}{d z} \left( u+t+s\right) ]\]

We can easily see why it fails!  There is a cheap way around this limitation however.  Imagine that we want to compute the gradient of s^2+t^2+u^2.  Then we have:

(%i28) express(grad(s^2+t^2+u^2));
subst([s=x,t=y,u=z],%);
ev(%,diff);
subst([x=s,y=t,z=u],%);
\[\tag{%o25} [\frac{d}{d x} \left( {{u}^{2}}+{{t}^{2}}+{{s}^{2}}\right) ,\frac{d}{d y} \left( {{u}^{2}}+{{t}^{2}}+{{s}^{2}}\right) ,\frac{d}{d z} \left( {{u}^{2}}+{{t}^{2}}+{{s}^{2}}\right) ]\] \[\tag{%o26} [\frac{d}{d x} \left( {{z}^{2}}+{{y}^{2}}+{{x}^{2}}\right) ,\frac{d}{d y} \left( {{z}^{2}}+{{y}^{2}}+{{x}^{2}}\right) ,\frac{d}{d z} \left( {{z}^{2}}+{{y}^{2}}+{{x}^{2}}\right) ]\] \[\tag{%o27} [2 x,2 y,2 z]\] \[\tag{%o28} [2 s,2 t,2 u]\]

So we simply swap out the variables for x,y,z, take the gradient, and then reswap back in.

 3.3.2 Always 3D Output

Note that if we have a function of two variables like:

(%i30) f(x,y):=x*y;
@grad(f(x,y));
\[\tag{%o29} \operatorname{f}\left( x,y\right) :=x y\] \[\tag{%o30} [y,x,0]\]

The output is still a 3D vector.  If we insist on a 2D vector as output, we can do the following:

(%i33) f(x,y):=x*y;
@grad(f(x,y))$
[%[1],%[2]];
\[\tag{%o31} \operatorname{f}\left( x,y\right) :=x y\] \[\tag{%o33} [y,x]\]

The above chain works by simply dropping the 3rd component from the gradient.

 3.3.3 grad() from Scratch to Fix These Issues

You could always compute a gradient from scratch using partial derivatives like so:

(%i35) f(x,y):=x*y^2;
define(
   gradf(x,y),
   [diff(f(x,y),x),diff(f(x,y),y)]
   );
\[\tag{%o34} \operatorname{f}\left( x,y\right) :=x\, {{y}^{2}}\] \[\tag{%o35} \operatorname{gradf}\left( x,y\right) :=[{{y}^{2}},2 x y]\]
(%i37) gradf(x,y);gradf(1,1);
\[\tag{%o36} [{{y}^{2}},2 x y]\] \[\tag{%o37} [1,2]\]

As you can see, the above gets around the limitation of grad() producing only 3D vectors.  Also:

(%i39) G:s^3+sqrt(s-t^2);
define(
   gradG(s,t),
   [diff(G,s),diff(G,t)]
);
\[\tag{G}\sqrt{s-{{t}^{2}}}+{{s}^{3}}\] \[\tag{%o39} \operatorname{gradG}\left( s,t\right) :=[\frac{1}{2 \sqrt{s-{{t}^{2}}}}+3 {{s}^{2}},-\frac{t}{\sqrt{s-{{t}^{2}}}}]\]

So doing things from scratch removes the requirement of x,y,z variables as well!

 3.4 Finding Critical Points and the Discriminant

Now that we know how to compute gradients, let's try to find some critical points.   Suppose we want to find the critical points of the paraboloid: z=x^2+y^2 --- which we know, of course, to be a minimum at (0,0).

(%i40) z(x,y):=x^2+y^2;
\[\tag{%o40} \operatorname{z}\left( x,y\right) :={{x}^{2}}+{{y}^{2}}\]

This function only has two variables, so our gradient is only the first two components of grad(z(x,y)):

(%i42) @grad(z(x,y));
define(gz(x,y),[%[1],%[2]]);
\[\tag{%o41} [2 x,2 y,0]\] \[\tag{%o42} \operatorname{gz}\left( x,y\right) :=[2 x,2 y]\]
(%i43) solve(gz(x,y))[1];
\[\tag{%o43} [y=0,x=0]\]

Recall that solve(): (1) defaults to setting everything to 0 unless told otherwise and (2) always gives a list of solutions, which in this case is a list of points.  We only want the first point, so we ask for solve()[1].

Let's compute the discriminant and test our point.  This is going to require some partial derivatives:

(%i44) define(
   D(x,y),
   diff(z(x,y),x,2)*diff(z(x,y),y,2)-diff(diff(z(x,y),x),y)^2
   );
\[\tag{%o44} \operatorname{D}\left( x,y\right) :=4\]

We see that D(0,0)>0 is clearly true, so that the concavity is consistent around the critical point of (0,0).  Let's get that concavity:

(%i45) subst([x=0,y=0],diff(z(x,y),x,2));
\[\tag{%o45} 2\]

Note that the concavity is up at (0,0) and we see that the point we found was a minimum.

 3.5 A More Complicated Extrema Problem

Now consider z(x,y)=x^2 + y^2 - 1/x -1/y

(%i46) z(x,y):=x^2+y^2+1/x+1/y;
\[\tag{%o46} \operatorname{z}\left( x,y\right) :={{x}^{2}}+{{y}^{2}}+\frac{1}{x}+\frac{1}{y}\]
(%i48) @grad(z(x,y))$
define(gz(x,y),[%[1],%[2]]);
\[\tag{%o48} \operatorname{gz}\left( x,y\right) :=[2 x-\frac{1}{{{x}^{2}}},2 y-\frac{1}{{{y}^{2}}}]\]
(%i49) S:solve(gz(x,y));
\[\tag{S}[[y=0.79370049732915,x=-0.6873648184993 \% i-0.39685026299204],[y=0.79370049732915,x=0.6873648184993 \% i-0.39685026299204],[y=0.79370049732915,x=0.79370049732915],[y=0.6873648184993 \% i-0.39685026299204,x=-0.6873648184993 \% i-0.39685026299204],[y=0.6873648184993 \% i-0.39685026299204,x=0.6873648184993 \% i-0.39685026299204],[y=0.6873648184993 \% i-0.39685026299204,x=0.79370049732915],[y=-0.6873648184993 \% i-0.39685026299204,x=-0.6873648184993 \% i-0.39685026299204],[y=-0.6873648184993 \% i-0.39685026299204,x=0.6873648184993 \% i-0.39685026299204],[y=-0.6873648184993 \% i-0.39685026299204,x=0.79370049732915]]\]

Yikes!  Let's try to_poly_solve() to see if that's any better.  To note, to_poly_solve() requires we list the variables as well as the equations:

(%i50) S2:to_poly_solve(gz(x,y),[x,y]);
\[\mbox{}\\\mbox{to\_ poly\_ solve: to\_ poly\_ solver.mac is obsolete; I'm loading to\_ poly\_ solve.mac instead.}\mbox{}\\\mbox{Loading maxima-grobner \$ Revision: 1.6 \$ \$ Date: 2009-06-02 07:49:49 \$ }\] \[\tag{S2}\]

To grab the real solutions out of the %union() we use part():

(%i52) part(S2,1,1);part(S2,1,2);
\[\tag{%o51} x=\frac{1}{{{2}^{\frac{1}{3}}}}\] \[\tag{%o52} y=\frac{1}{{{2}^{\frac{1}{3}}}}\]

The part(S2,1,1) means first element of S2, first part, while part(S2,1,2) means first element of S2, second part.  It should be noted that:

(%i54) float(part(S2,1,1));float(part(S2,1,2));
\[\tag{%o53} x=0.79370052598409\] \[\tag{%o54} y=0.79370052598409\]

Were exactly the values that were given in solve(), they were just numerical approximations.  to_poly_solve() was able to extract exact values.  Now let's test the discriminant:

(%i55) define(
   D(x,y),
   diff(z(x,y),x,2)*diff(z(x,y),y,2)-diff(diff(z(x,y),x),y)^2
);
\[\tag{%o55} \operatorname{D}\left( x,y\right) :=\left( \frac{2}{{{x}^{3}}}+2\right) \, \left( \frac{2}{{{y}^{3}}}+2\right) \]
(%i56) D(1/2^(1/3),1/2^(1/3));
\[\tag{%o56} 36\]

We see that this is positive, so that the critical point has consistent concavity.  Now to determine the behavior:

(%i57) subst([part(S2,1,1),part(S2,1,2)],diff(z(x,y),x,2));
\[\tag{%o57} 6\]

We see that the second partial is positive and so the graph is concave up at the critical point, making it a minimum.

 4 Vector Potentials

Often it is of importance to us to find a potential function for a given vector field. Thankfully, Maxima's vect library
has a command called potential() that does exactly this, and it works very straightforwardly.  Let's illustrate with
a few examples:

--> potential([x,y],[x,y]);
\[\tag{%o111} \frac{{{y}^{2}}+{{x}^{2}}}{2}\]
--> potential([x,exp(y)],[x,y]);
\[\tag{%o119} \frac{2 {{\% e}^{y}}+{{x}^{2}}-2}{2}\]

Recall that the potential function is unique only up to a scalar constant, and the above has an additional -1 that could be dropped from the potential.

--> potential([x,y,0],[x,y,z]);
\[\tag{%o120} \frac{{{y}^{2}}+{{x}^{2}}}{2}\]
--> potential([x,y,z],[x,y,z]);
\[\tag{%o121} \frac{{{z}^{2}}+{{y}^{2}}+{{x}^{2}}}{2}\]
--> potential([x^2,y^2,z^2],[x,y,z]);
\[\tag{%o130} \frac{{{z}^{3}}+{{y}^{3}}+{{x}^{3}}}{3}\]

CAUTION: only use potential with x,y, and z!  Using other variables will reset the variables used in grad() leading to difficult
to diagnose errors in your calculations!

So we see that potential() takes in a vector-field along with a list of variables for that field and then it proceeds to uncover the potential for that field.  Let's put this through its paces a little bit by taking a gradient of a pretty nast function, and trying to
reverse it with potential():

(%i59) f(x,y,z):=x*exp(y*z);
define(F(x,y,z), @grad(f(x,y,z)));
\[\tag{%o58} \operatorname{f}\left( x,y,z\right) :=x \operatorname{exp}\left( y z\right) \] \[\tag{%o59} \operatorname{F}\left( x,y,z\right) :=[{{\% e}^{y z}},x z\, {{\% e}^{y z}},x y\, {{\% e}^{y z}}]\]
(%i60) potential(F(x,y,z),[x,y,z]);
\[\tag{%o60} x\, {{\% e}^{y z}}\]

Nice!  Let's ratchet up the difficulty one more notch:

(%i63) f(x,y):=(x^2+y*cos(x))*sin(x*sqrt(y));
@grad(f(x,y))$
define(gf(x,y),[%[1],%[2]]);
\[\tag{%o61} \operatorname{f}\left( x,y\right) :=\left( {{x}^{2}}+y \cos{(x)}\right) \sin{\left( x\, \sqrt{y}\right) }\] \[\tag{%o63} \operatorname{gf}\left( x,y\right) :=[\sin{\left( x\, \sqrt{y}\right) } \left( 2 x-\sin{(x)} y\right) +\cos{\left( x\, \sqrt{y}\right) } \sqrt{y}\, \left( \cos{(x)} y+{{x}^{2}}\right) ,\frac{x \cos{\left( x\, \sqrt{y}\right) } \left( \cos{(x)} y+{{x}^{2}}\right) }{2 \sqrt{y}}+\cos{(x)} \sin{\left( x\, \sqrt{y}\right) }]\]
(%i64) P:potential(gf(x,y),[x,y]);
\[\mbox{}\\unable to prove that the following difference between the input and the gradient of the returned result is zero [-\frac{\left( \cos{\left( x\, \sqrt{y}+x\right) }+\cos{\left( x\, \sqrt{y}-x\right) }-2 \cos{(x)} \cos{\left( x\, \sqrt{y}\right) }\right) \, {{y}^{\frac{3}{2}}}+\left( \cos{\left( x\, \sqrt{y}+x\right) }-\cos{\left( x\, \sqrt{y}-x\right) }+2 \sin{(x)} \sin{\left( x\, \sqrt{y}\right) }\right) y}{2},-\frac{\left( x \cos{\left( x\, \sqrt{y}+x\right) }+x \cos{\left( x\, \sqrt{y}-x\right) }-2 x \cos{(x)} \cos{\left( x\, \sqrt{y}\right) }\right) \, \sqrt{y}+2 \sin{\left( x\, \sqrt{y}+x\right) }+2 \sin{\left( x\, \sqrt{y}-x\right) }-4 \cos{(x)} \sin{\left( x\, \sqrt{y}\right) }}{4}] \] \[\tag{P}\frac{\left( \sin{\left( x\, \sqrt{y}+x\right) }+\sin{\left( x\, \sqrt{y}-x\right) }\right) y+2 {{x}^{2}} \sin{\left( x\, \sqrt{y}\right) }}{2}\]

Maxima confides in us that it's having a little difficulty confirming that the above is genuinely correct, but if we trigexpand() the last result:

(%i65) trigexpand(P);
\[\tag{%o65} \frac{2 \cos{(x)} \sin{\left( x\, \sqrt{y}\right) } y+2 {{x}^{2}} \sin{\left( x\, \sqrt{y}\right) }}{2}\]

and then ratsimp() the above:

(%i67) ratsimp(%o65);
\[\tag{%o67} \cos{(x)} \sin{\left( x\, \sqrt{y}\right) } y+{{x}^{2}} \sin{\left( x\, \sqrt{y}\right) }\]

We see that Maxima has nailed it!


Created with wxMaxima.