1021-Gradients.sagews

Gradients

Based on material from Ben Woodruff

It is simple to calculate the gradient of a function.

f(x,y)=9-x^2-y^2
xbounds=(x,-3,3)
ybounds=(y,-3,3)
plot3d(f(x,y), xbounds,ybounds)

h=f.gradient()
h(x,y)

(-2*x, -2*y)

We can plot the gradient field easily as well.

f(x,y)=9-x^2-y^2
xbounds=(x,-3,3)
ybounds=(y,-3,3)

h=f.gradient()
plot_vector_field(h(x,y), xbounds, ybounds,aspect_ratio=1)

We can also plot the level curves (i.e., a contour plot) on top of the gradient field, so that we can see that the gradient is always normal (i.e., perpendicular) to the level curve.

f(x,y)=9-x^2-y^2
xbounds=(x,-3,3)
ybounds=(y,-3,3)

h=f.gradient()
contour_plot(f(x,y), xbounds, ybounds,fill=False,aspect_ratio=1)+plot_vector_field(h(x,y), xbounds, ybounds)

This all holds for functions from $\RR^3 \to \RR$ as well.  However, in this case, the function has "level surfaces", and the gradient is a 3-dimensional vector.

f(x,y,z)=x^2-y^2+z
xbounds=(x,-3,3)
ybounds=(y,-3,3)
zbounds=(z,-3,3)

h=f.gradient()
h(x,y,z)

(2*x, -2*y, 1)

Evaluate the cell below to get the plot_vector_field3d command, which allows us to plot this 3d gradient field.

%auto
from matplotlib.cm import get_cmap
from matplotlib.colors import LinearSegmentedColormap
from sage.plot.misc import setup_for_eval_on_grid

def plot_vector_field3d(functions, xrange, yrange, zrange, 
                        plot_points=5, colors='jet', center_arrows=False,**kwds):
    r"""
    Plot a 3d vector field

    INPUT:

    - ``functions`` - a list of three functions, representing the x-,
      y-, and z-coordinates of a vector

    - ``xrange``, ``yrange``, and ``zrange`` - three tuples of the
      form (var, start, stop), giving the variables and ranges for each axis

    - ``plot_points`` (default 5) - either a number or list of three
      numbers, specifying how many points to plot for each axis

    - ``colors`` (default 'jet') - a color, list of colors (which are
      interpolated between), or matplotlib colormap name, giving the coloring
      of the arrows.  If a list of colors or a colormap is given,
      coloring is done as a function of length of the vector

    - ``center_arrows`` (default False) - If True, draw the arrows
      centered on the points; otherwise, draw the arrows with the tail
      at the point

    - any other keywords are passed on to the plot command for each arrow

    EXAMPLES::
    
        sage: x,y,z=var('x y z')
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi))
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),colors=['red','green','blue'])
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),colors='red')
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),plot_points=4)
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),plot_points=[3,5,7])
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),center_arrows=True)
    """
    (ff,gg,hh), ranges = setup_for_eval_on_grid(functions, [xrange, yrange, zrange], plot_points)
    xpoints, ypoints, zpoints = [srange(*r, include_endpoint=True) for r in ranges]
    points = [vector((i,j,k)) for i in xpoints for j in ypoints for k in zpoints]
    vectors = [vector((ff(*point), gg(*point), hh(*point))) for point in points]
    
    try:
        cm=get_cmap(colors)
        assert(cm is not None)
    except (TypeError, AssertionError):
        if isinstance(colors, (list, tuple)):
            cm=LinearSegmentedColormap.from_list('mymap',colors)
        else:
            cm = lambda x: colors
                
    max_len = max(v.norm() for v in vectors)
    scaled_vectors = [v/max_len for v in vectors]
    
    if center_arrows:
        return sum([plot(v,color=cm(v.norm()),**kwds).translate(p-v/2) for v,p in zip(scaled_vectors, points)])
    else:
        return sum([plot(v,color=cm(v.norm()),**kwds).translate(p) for v,p in zip(scaled_vectors, points)])

You'll notice that the gradient vectors in the plot below are all perpendicular to the level surfaces.

implicit_plot3d(f(x,y,z), xbounds, ybounds, zbounds, contour=[-4,0,4],opacity=0.8)+plot_vector_field3d(h(x,y,z), xbounds, ybounds, zbounds)

Directional Derivatives

A directional derivative is just the gradient dotted with a unit vector in the direction. Remember to evaluate at the point.

f(x,y)=9-x^2-y^2
u=vector([sqrt(3)/2, -1/2])
pt=(1/2, sqrt(3)/2)
directional_derivative = f.gradient()*(u/norm(u))
directional_derivative(x,y)

-sqrt(3)*x + y

To evaluate it at a point, just do the following.  Literally, this is taking the first thing in point (the $x$-coordinate) and making that the first argument of the function, and so on.

directional_derivative(*pt)

0

In the interact below, the directional derivative is the slope of the green tangent line.  The cross-section in the direction is the red curve, and the gradient vector is the orange vector.

var('x,y,t,z')
f(x,y)=sin(x)*cos(y)

pif = float(pi)

line_thickness=3
surface_color='blue'
plane_color='purple'
line_color='red'
tangent_color='green'
gradient_color='orange'

@interact
def myfun(location=input_grid(1, 2, default=[0,0], label = "Location (x,y)", width=2), angle=slider(0, 2*pif, label = "Angle"), 
show_surface=("Show surface", True)):
    location3d = vector(location[0]+[0])
    location = location3d[0:2]
    direction3d = vector(RDF, [cos(angle), sin(angle), 0])
    direction=direction3d[0:2]
    cos_angle = math.cos(angle)
    sin_angle = math.sin(angle)
    df = f.gradient()
    direction_vector=line3d([location3d, location3d+direction3d], arrow_head=True, rgbcolor=line_color, thickness=line_thickness)
    curve_point = (location+t*direction).list()
    curve = parametric_plot(curve_point+[f(*curve_point)], (t,-3,3),color=line_color,thickness=line_thickness)
    plane = parametric_plot((cos_angle*x+location[0],sin_angle*x+location[1],t), (x, -3,3), (t,-3,3),opacity=0.8, color=plane_color)
    pt = point3d(location3d.list(),color='green', size=10)

    tangent_line = parametric_plot((location[0]+t*cos_angle, location[1]+t*sin_angle, f(*location)+t*df(*location)*(direction)), (t, -3,3), thickness=line_thickness, color=tangent_color)
    picture3d = direction_vector+curve+plane+pt+tangent_line

    picture2d = contour_plot(f(x,y), (x,-3,3),(y,-3,3), plot_points=100)
    picture2d += arrow(location.list(), (location+direction).list()) 
    picture2d += point(location.list(),rgbcolor='green',pointsize=40)
    if show_surface:
        picture3d += plot3d(f(x,y), (x,-3,3),(y,-3,3),opacity=0.7)
        
    dff = df(location[0], location[1])
    dff3d = vector(RDF,dff.list()+[0])
    picture3d += line3d([location3d, location3d+dff3d], arrow_head=True, rgbcolor=gradient_color, thickness=line_thickness)
    picture2d += arrow(location.list(), (location+dff).list(), rgbcolor=gradient_color, width=line_thickness)
    show(picture3d,aspect=[1,1,1], axes=True)
    show(picture2d, aspect_ratio=1)

 

The 2nd Derivative test

The second derivative test looks at the eigenvalues of the Hessian matrix evaluated at critical points.

def illustrate_test(f, fvals=(-5,5), tvals=(-3,3)):
    """
    Return a table of critical points, eigenvalues, and eigenvectors, plus a 3D graph illustrating the curvature at each critical point.
    
    INPUTS:
      f -- a function.  Usually this should be a callable function, defined like f(x,y)=x^2+y^2, for example.
      fvals (default (-5,5)) -- minimum and maximum values for the function plot
      tvals (default (-3,3)) -- minimum and maximum values for the parameter in the curves at each critical point

    EXAMPLES:
        sage: var('x,y')
        sage: f(x,y)=-x^2+y^2
        sage: table, graph=illustrate_test(f)
        sage: print table
        sage: show(graph)
        
        sage: var('x,y')
        sage: f(x,y)=x*y*exp((-x^2-y^2)/3)
        sage: table, graph=illustrate_test(f)
        sage: print table
        sage: show(graph)
    """
    line_thickness=3
    surface_color='blue'
    line_color = ['red','blue']
    x,y,t=var('x,y,t')
    tmin, tmax = tvals
    fmin, fmax = fvals
    f_parametric = vector([x,y,f(x,y)])
    picture = plot3d(f,(x,fmin,fmax),(y,fmin,fmax), opacity=0.3, surface_color=surface_color)
    hessian=f(x,y).hessian()
    solns=solve([eqn(x,y)==0 for eqn in f.gradient()],[x,y], solution_dict=True)

    table=[['Color','Critical Point', 'Eigenvalue', 'Eigenvector']]
    pictures=[]
    i=0
    for solution in [s for s in solns if s[x] in RR and s[y] in RR]:
        soln = dict([(key,RDF(val)) for key,val in solution.items()])
        critical_point = f_parametric(soln)
        hessian_at_point = hessian(soln).change_ring(RDF)
        evals,evecs=hessian_at_point.eigenvectors_right()
        picture += point3d(critical_point.list(), color='green', size=10)
        for eval,evec in zip(evals,evecs.columns()):
            direction_plot = evec.plot(width=line_thickness,color=line_color[i%2]).plot3d().translate(critical_point)    
            curve_formula = f_parametric(x=soln[x]+t*evec[0],y=soln[y]+t*evec[1])
            curve = parametric_plot(curve_formula, (t,tmin,tmax),color=line_color[i%2],thickness=line_thickness)
            picture+=curve+direction_plot
            pictures += [curve]
            table.append([line_color[i%2],(critical_point), eval, (evec)])
            i+=1

    return table, picture

f(x,y)=-x^2+y^2
table, graph=illustrate_test(f)
html.table(table,header=True)
show(graph)

Color	Critical Point	Eigenvalue	Eigenvector
red	\left(0.0,0.0,0\right)	-2.0	\left(1.0,0.0\right)
blue	\left(0.0,0.0,0\right)	2.0	\left(0.0,1.0\right)

f(x,y)=x*y*exp((-x^2-y^2)/3)
table, graph=illustrate_test(f)
html.table(table,header=True)
show(graph)

Color	Critical Point	Eigenvalue	Eigenvector
red	\left(0.0,0.0,0\right)	1.0	\left(0.707106781187,0.707106781187\right)
blue	\left(0.0,0.0,0\right)	-1.0	\left(-0.707106781187,0.707106781187\right)
red	\left(-1.22474487139,-1.22474487139,0.551819161757\right)	-0.735758882343	\left(1.0,0.0\right)
blue	\left(-1.22474487139,-1.22474487139,0.551819161757\right)	-0.735758882343	\left(0.0,1.0\right)
red	\left(1.22474487139,-1.22474487139,-0.551819161757\right)	0.735758882343	\left(1.0,0.0\right)
blue	\left(1.22474487139,-1.22474487139,-0.551819161757\right)	0.735758882343	\left(0.0,1.0\right)
red	\left(-1.22474487139,1.22474487139,-0.551819161757\right)	0.735758882343	\left(1.0,0.0\right)
blue	\left(-1.22474487139,1.22474487139,-0.551819161757\right)	0.735758882343	\left(0.0,1.0\right)
red	\left(1.22474487139,1.22474487139,0.551819161757\right)	-0.735758882343	\left(1.0,0.0\right)
blue	\left(1.22474487139,1.22474487139,0.551819161757\right)	-0.735758882343	\left(0.0,1.0\right)

var('x,y')
f(x,y)=x^3 - y^3 - 2*x*y + 6
table, graph=illustrate_test(f, fvals=(-1,1), tvals=(-1,1))
html.table(table,header=True)
show(graph)

Color	Critical Point	Eigenvalue	Eigenvector
red	\left(-0.666666666667,0.666666666667,6.2962962963\right)	-2.0	\left(0.707106781187,-0.707106781187\right)
blue	\left(-0.666666666667,0.666666666667,6.2962962963\right)	-6.0	\left(0.707106781187,0.707106781187\right)
red	\left(0.0,0.0,6\right)	2.0	\left(0.707106781187,-0.707106781187\right)
blue	\left(0.0,0.0,6\right)	-2.0	\left(0.707106781187,0.707106781187\right)

Lagrange Multipliers

We'll work out an example problem.

Find the maximum and minimum temperatures on the circle $x^2+y^2=1$ if the temperature of a point $(x,y)$ is $T(x,y)=4+x^2-y$.

var('x,y,z,t')
T(x,y)=4+x^2-y
g(x,y)=x^2+y^2-1

p=parametric_plot([cos(t), sin(t), T(cos(t), sin(t))], (t, 0, 2*pi),thickness=8,color='black',aspect_ratio=1)
p+=plot3d(T(x,y), (x,-2,2), (y,-2,2),opacity=0.6)
p+=implicit_plot3d(g(x,y)==0,(x,-2,2),(y,-2,10),(z,0,10),aspect_ratio=1,color='red',opacity=0.6) # The y-value is set to work around a bug in Sage 4.2
show(p)

Another way to imagine it is just as the constraint function (the circle), drawn on the temperature function.

p=parametric_plot([cos(t), sin(t), T(cos(t), sin(t))], (t, 0, 2*pi),thickness=5)
p+=plot3d(T(x,y), (x,-2,2), (y,-2,2),opacity=0.4)
p+=points([[sqrt(3/4),-1/2,21/4],[-sqrt(3/4),-1/2,21/4]],size=10,color='green')
p+=point([0,1,3],size=10,color='red')
p+=circle((0,0),1,color='black').plot3d()
show(p,aspect_ratio=1)

p=contour_plot(T(x,y), (x,-2,2),(y,-2,2),contours=40,fill=False)
p+=implicit_plot(g(x,y)==0, (x,-2,2), (y,-2,2),cmap=['blue'])
p+=points([[sqrt(3/4),-1/2],[-sqrt(3/4),-1/2]],pointsize=50,color='green')
p+=point([0,1],pointsize=50,color='red')
show(p,aspect_ratio=1)

Calculate the Lagrangian $L(x,y,\lambda)=T(x,y)-\lambda*g(x,y)$.  We'll use "lambda_" for $\lambda$ because we can't use the word "lambda" in python/Sage (it means something different than what we want).

var('lambda_')
L(x,y,lambda_)=T(x,y)-lambda_*g(x,y)

gradL=L.gradient()(x,y,lambda_)
gradL

(-2*lambda_*x + 2*x, -2*lambda_*y - 1, -x^2 - y^2 + 1)

We get four solutions to the three equations given by $\nabla L=\vec 0$.

solve([gradL[0]==0,gradL[1]==0,gradL[2]==0], [x,y,lambda_])

[[x == 0, y == -1, lambda_ == (1/2)], [x == 0, y == 1, lambda_ == (-1/2)], [x == 1/2*sqrt(3), y == (-1/2), lambda_ == 1], [x == -1/2*sqrt(3), y == (-1/2), lambda_ == 1]]

Now let's plug each of these $(x,y)$ values in to $T(x,y)$ to find the global maxima and minima.

T(0,1)

3

T(0,-1)

5

T(sqrt(3/4), -1/2)

21/4

T(-sqrt(3/4), -1/2)

21/4

The maximum is therefore at both $(\sqrt{3/4},-1/2)$ and $(-\sqrt{3/4},-1/2)$ (max value $21/4$) and the minimum is at $(0,1)$ (min value 3).