## The Gradient in Matlab

In the activity Directional Derivatives in Matlab, we investigated the derivative in an arbitrary direction, called the *directional derivative*. Let's repeat some of that work here.

### Defining the Gradient

We begin by picking an arbitrary point `(a,b)` at which we wish to find the directional derivative. Then we pick a unit vector `vec u=langle u_1, u_2 rangle`, is whose direction we wish to find the derivative. To find the equation of the line in the `xy`-plane that passes through `(a,b)` in the direction of `vec u`, let `X(x,y)` be an arbitrary point on the line. Then the equation of the line in vector form is

`vec (PX)= t vec(u).`

This is equivalent to the vector equation

`langle x-a,y-b rangle=t langle u_1, u_2 rangle`,

which leads to the parametric equations of the line.

$\begin{eqnarray} x&=&a+u_1 t\\ y&=&b+u_2 t \end{eqnarray}$

In the activity Directional Derivatives in Matlab, we looked at points on the surface having `x`- and `y`-values restricted to the above line. There we had:

`z=f(x,y) \quad\text{where}\quad x=a+u_1t \quad\text{and}\quad y = b+u_2t`

Using the chain rule, we calculate `dz//dt`.

`frac(dz)(dt)=frac(partial f)(partial x) frac(dx)(dt)+frac(partial f)(partial y)frac(dy)(dt)`

Because `dx//dt=u_1` and `dy//dt=u_2`, this becomes

`frac(dz)(dt)=frac(partial f)(partial x) u_1+frac(partial f)(partial y) u_2.`

However, note that the right-hand side of this last result can be written as a dot product.

`frac(dz)(dt)=langle frac(partial f)(partial x),frac(partial f)(partial y) rangle cdot langle u_1,u_2 rangle.`

Time to define the gradient.

Given a function `f:R^2 \to R`, the gradient of `f` is defined to be

`nabla f=langle frac(partial f)(partial x), frac(partial f)(partial y) rangle`.

With this definition, the directional derivative can be written as follows:

`frac(dz)(dt)=langle frac(partial f)(partial x),frac(partial f)(partial y) rangle cdot langle u_1,u_2 rangle =nabla f cdot vec(u).`

Finally, if we adopt the notation `D_u f` to mean the derivative of `f` in the direction of the unit vector `vec u` (the directional derivative), then we have the following result.

Suppose that `f:R^2\to R` and `vec u` is a unit vector. Then the derivative of `f` in the direction of `vec(u)` is given by

`D_u f=nabla f cdot vec(u)`.

Thus, the slope of the tangent line to the surface at the point `(a,b,f(a,b))` is given by

`D_u f(a,b)=nabla f(a,b) cdot vec(u)`.

Let's make an example calculation.

Given `f(x,y)=9-x^2-y^2`, calculate the gradient, then evaluate the gradient at the point `(1,1)`.

**Solution: ** Because `f(x,y)=9-x^2-y^2`, we have `partial f//partial x = -2x` and `partial f//partial y=-2y`. Thus, the gradient, evaluated at `(x,y)`, is given by

`nabla f(x,y)=langle frac(partial f)(partial y), frac(partial f)(partial y) rangle=langle -2x,-2y rangle`.

Therefore the gradient, evaluated at the point `(1,1)`, is the vector

`nabla f(1,1)=langle -2,-2 rangle`.

Let's build on this last example and calculate a directional derivative.

Given the function `f(x,y)=9-x^2-y^2`, find the directional derivative of `f` at the point `(1,1)` in the direction of the unit vector `vec u=langle sqrt(2)//2,sqrt(2)//2 rangle`.

**Solution: **We've already calculated the gradient in Example 1. Therefore,

$\begin{eqnarray} D_f f(1,1) &=&\nabla f(1,1) \cdot \vec u\\ &=&\langle -2, -2 \rangle \cdot \langle \sqrt2/2, \sqrt2/2 \rangle\\ &=&-\sqrt2-\sqrt2\\ &=&-2\sqrt2. \end{eqnarray}$

Readers might recognize this as the slope of the tangent line from the activity Directional Derivatives in Matlab.

### The Gradient Field

In Example 1, for the function `f(x,y)=9-x^2-y^2`, the gradient vector at the point `(1,1)` was `nabla f(1,1)=langle -2,-2 rangle`. We've plotted this gradient vector at the point `(1,1)` in Figure 1. We've also calculated an additional gradient vector at the point `(5,-5)`, getting `nabla f(5,-5)=langle -10,10 rangle` and including it in Figure 1. Similar calculations show that the further we get from the origin, the longer the gradient vectors become.

We could continue by hand, calculating the gradient vector at each point of the grid in Figure 1 and plotting the result. The result would be a *field of vectors*, or equivalently, a *vector field*. However, a hand plot is tedious, particularly when Matlab does an excellent job of plotting gradient fields.

First, create a grid of `(x,y)` pairs at which we will plot gradient vectors.

[x,y]=meshgrid(-5:1:5);

Note that **-5:1:5** produces the vector **[-5,-4,-3,-2,-1,0,1,2,3,4,5]** and **meshgrid(-5:1:5)** uses this vector for both `x` and `y`. Next, calculate the gradient vectors at each point of the grid. Because we are using the function `f(x,y)=9-x^2-y^2`, we have `partial f//partial x = -2x` and `partial f//partial y = -2y`.

vx=-2*x; vy=-2*y;

We use the familiar **quiver** command, turning off all scaling.

quiver(x,y,vx,vy,0);

Adding the usual annotations produces the result in Figure 2.

Although the gradient field of Figure 2 is completely accurate, it is very difficult to interpret the result due to overcrowding. (Indeed, it appears that the gradient vectors are pointing outward, an optical illusion contradicted by the evidence in Figure 1.) A standard practice is to scale the vectors appropriately to aid in the interpretation. This is accomplished by letting Matlab's **quiver** command sketch the field with scaling turned on.

quiver(x,y,vx,vy);

The resulting *gradient field* is shown in Figure 3.

The vectors in Figure 3 have the inccorrect length. For example, the vector at the point `(1,1)` should have the length shown in Figure 1. Hoever, the vector is pointing in the correct direction and it has the *correct relative length*, as compared to the remaining vectors in the gradient field.

### Adding Contours

The gradient possesses several important properties, one of which is its relation to the set of level curves of the function. Let's add some contours to the gradient field in Figure 4.

First create a finer mesh for more accuracy.

[x,y]=meshgrid(linspace(-5,5,50));

Recalculate the function `f(x,y)=9-x^2-y^2` on this finer mesh.

z=9-x.^2-y.^2;

Hold the plot and add the level curves.

hold on contour(x,y,z)

Equalize the axes and add a grid.

grid on axis equal

The result is shown in Figure 4. Note that the gradient vectors all appear to be orthogonal to the level curves of the function.

### Direction of Gradient Vectors

One final question remains: If the gradient vectors are orthogonal to the level curves, there are two possible directions that they could point, inward or outward. Why do they point inward in our example? this is most easily vizualized by drawing a surface and adding both level curves and the gradient field to the image. We start a new image by closing all open figure windows.

close all

Create the mesh and recalculate the function values for`f(x,y)=9-x^2-y^2`.

[x,y]=meshgrid(linspace(-5,5,40)); z=9-x.^2-y.^2;

Now we draw the surface and contours with the **meshc** command. We also use Matlab's *handle graphics* to color the surface with magenta and set the mesh lines to black.

h1=meshc(x,y,z); set(h1,'FaceColor','m',... 'FaceAlpha',0.5,... 'EdgeColor','k')

Add a grid, turn on the box, and rotate into the standard orientation.

grid on box on view([130,30])

The result is shown in Figure 5.

Now, let's add the gradient field. Because our figure is three-dimensional, we will use the **quiver3** command to add the gradient field. Begin with a coarser grid where you want the gradient vectors to be plotted.

[x,y]=meshgrid(-5:5);

Calculate the positions where you want to plot the gradient vectors. In Figure 5, note that the bottom plane is set at `z=-50`. Thus, we want a matrix of `z`-values the same size as `x`- and `y`-matrices that is filled with the constant `-50`.

z=-50*ones(size(x));

Recall that `nabla f(x,y)=langle -2x, -2y rangle`. In addition, we want the `z`-component of the vectors to be zero, so we fill a matrix **vz** with zeros the same size as the matrix **vx**.

vx=-2*x; vy=-2*y; vz=zeros(size(vx));

Hold the plot and use the **quiver3** command to plot the gradient field in the desired position.

hold on quiver3(x,y,z,vx,vy,vz);

Add a grid and turn on the box, then rotate into the standard orientation.

grid on box on view([130,30])

Add the usual annotations.

xlabel('x-axis') ylabel('y-axis') title('Surface, Contours, and Gradient Field.')

The result is the surface, contours, and gradient field shown in Figure 6.

The question of direction now seems to be resolved. In class you will be given a more rigorous proof, but from the appearance in Figure 6, the gradient vectors seem to be pointing in the direction of greatest increase of the function. If you consider the contours to be a topographical map, then the gradient vectors are pointing in an "up-hill" direction.

### Matlab Files

Although the following file features advanced use of Matlab, we include it here for those interested in discovering how we generated the images for this activity. You can download the Matlab file at the following link. Download the file to a directory or folder on your system.

The file **gradient.m** is designed to be run in "cell mode." Open the file **gradient.m** in the Matlab editor, then enable cell mode from the **Cell Menu**. After that, use the entries on the **Cell Menu** or the icons on the toolbar to execute the code in the cells provided in the file. There are options for executing both single and multiple cells. After executing a cell, examine the contents of your folder and note that a PNG file was generated by executing the cell.

### Exercises

- Consider the function `f(x,y)=x^2+y^2-4`. Use Matlab to sketch a gradient field superimposed on contours as in Figure 6. Submit a printout of your result and the M-code that produced your result. Submit all mathematical computations required to produce your result.
Consider the function `f(x,y)=2-|x|-|y|`. Sketch a surface, level curve, and gradient plot such as that shown in Figure 6. Here's a little hint to help with the differentiation.

$\begin{eqnarray} D_x|x| &=&D_x\sqrt{x^2}\\ &=&D_x(x^2)^{1/2}\\ &=&\frac12(x^2)^{-1/2}(2x)\\ &=&\frac{x}{(x^2)^{1/2}}\\ &=&\frac{x}{\sqrt{x^2}}\\ &=&\frac{x}{|x|} \end{eqnarray}$

Submit a printout of your result and the M-code that produced your result. Submit all mathematical computations required to produce your result.