Mini-lecture 7-17 - The gradient vector field of a function of multiple variables.ipynb

Kernel: SageMath 9.3

from itertools import product
import numpy as np
import plotly.graph_objects
from ipywidgets import Button, Text, HTMLMath, HBox, VBox, Layout, interact as _original_interact

%%html
<!-- To disable the modebar IN ALL PLOT.LY PLOTS -->
<style>
.modebar {
    display: none !important;
}
</style>

_original_interact = interact
def interact(_function_to_wrap=None, _layout="horizontal", **kwargs):
    """interact, but with widgets laid out in a horizontal flexbox layout

    This function works exactly like 'interact' (from SageMath or ipywidgets), 
    except that instead of putting all of the widgets into a vertical box 
    (VBox), it uses a horizontal box (HBox) by default. The HBox uses a flexbox 
    layout, so that if there are many widgets, they'll wrap onto a second row. 
    
    Options:
        '_layout' - 'horizontal' by default. Anything else, and it will revert 
        back to using the default layout of 'interact' (a VBox). 
    """
    def decorator(f):
        retval = _original_interact(f, **kwargs)
        if _layout == "horizontal":
            widgets = retval.widget.children[:-1]
            output = retval.widget.children[-1]
            hbox = HBox(widgets, layout=Layout(flex_flow="row wrap"))
            retval.widget.children = (hbox, output)
        return retval
    if _function_to_wrap is None:
        # No function passed in, so this function must *return* a decorator
        return decorator
    # This function was called directly, *or* was used as a decorator directly
    return decorator(_function_to_wrap)

def find_zeros(field, *ranges, **options):
    """Numerically approximate all zeros of a vector field within a box

    Each 'range' should have the form (x, xmin, xmax), where 'x' is one of the
    variables appearing in the vector field, and 'xmin' and 'xmax' are the
    bounds on that variables.

    Options:
        'intervals' - How many subintervals to use in each dimension (default
            20). If there are zeros that are very close together, you may need
            to increase this to find them. But be warned that the running time
            is roughly this to the n power, where n is the number of variables.
        'tolerance' - How close to approximate roots, roughly (default 1e-4)
        'maxiter' - Maximum number of iterations of Newton's Method to run for
            any one (potential) solution. This defaults to -2*log(tolerance).
            There usually isn't much harm in increasing this, unless there are
            many false hits.
        'round' - Round the components of the solutions to this many decimal
            places (default None, meaning do not round them at all)
    """

    # Initialization
    intervals = options.get("intervals", 20)
    tolerance = options.get("tolerance", 1e-5)
    maxiter = options.get("maxiter", int(round(-2*log(tolerance))))
    roundto = options.get("round", None)
    n = len(ranges)
    if len(field) != n:
        raise ValueError("Dimension of vector field is {}, but {} ranges "
                         "given".format(len(field), n))

    mins  = [xmin - (xmax - xmin)/intervals*tolerance for x, xmin, xmax in ranges]
    maxes = [xmax + (xmax - xmin)/intervals*tolerance for x, xmin, xmax in ranges]
    deltas = [(xmax - xmin) / intervals for xmin, xmax in zip(mins, maxes)]
    powers = [1 << i for i in range(n)]
    J = jacobian(field, [x for x, xmin, xmax in ranges])
    def dist(v, w):
        return sqrt(sum(((a - b) / d)**2 for a, b, d in zip(v, w, deltas)))

    # Set up the array of positive/negative signs of the vector field
    signs = np.zeros((intervals + 1,) * n, dtype=int)
    sranges = [srange(m, m + d*(intervals + 0.5), d) for m, d in zip(mins, deltas)]
    for vertex, index in zip(product(*sranges), product(range(intervals + 1), repeat=n)):
        v = field(*vertex)
        signs[index] = sum(powers[i] for i in range(n) if v[i] > 0)

    # Now search through that array for potential solutions
    solutions = []
    mask = int((1 << n) - 1)
    for index in product(range(intervals), repeat=n):
        indexpowers = list(zip(index, powers))
        all0s = 0
        all1s = mask
        for k in range(mask + 1):
            newindex = [i + (1 if k & p else 0) for i, p in indexpowers]
            vertex = signs[tuple(newindex)]
            all0s |= vertex
            all1s &= vertex
        if all0s & ~all1s == mask:
            # Now do Newton's method!
            v = vector(m + d*(i + 0.5) for i, m, d in zip(index, mins, deltas))
            for i in range(maxiter):
                previous_v = v
                v = v - J(*v).solve_right(field(*v))
                if dist(v, previous_v) < tolerance:
                    break
            else:
                warn("{} iterations reached without convergence for solution "
                     "{}".format(maxiter, v), RuntimeWarning)
            if solutions and min(dist(v, w) for w in solutions) < 2*tolerance:
                continue
            if not all(xmin <= x <= xmax for x, xmin, xmax in zip(v, mins, maxes)):
                continue
            solutions.append(vector(RDF, v))

    # A convenience: round to some number of decimal places, if requested
    if roundto is not None:
        solutions = [vector(round(x, roundto) for x in v) for v in solutions]
    return solutions

# Our mixin class for Figure and FigureWidget. Should never be instantiated! 
class MyFigure(object):
    def add(self, *items, subplot=None):
        if subplot is not None:
            try:
                subplot[0]
            except:
                subplot = (1, subplot)
            row = int(subplot[0])
            col = int(subplot[1])
        retval = []
        text_indices = []
        text3d_indices = []
        for item in items:
            if isinstance(item, plotly.graph_objects.layout.Annotation):
                if subplot is None:
                    self.add_annotation(item)
                else:
                    self.add_annotation(item, row=row, col=col)
                text_indices.append(len(retval))
                retval.append(None)
            elif isinstance(item, plotly.graph_objects.layout.scene.Annotation):
                self.layout.scene.annotations += (item,)
                text3d_indices.append(len(retval))
                retval.append(None)
            else:
                if subplot is None:
                    self.add_trace(item)
                else:
                    self.add_trace(item, row=row, col=col)
                retval.append(self.data[-1])
        for i, pos in enumerate(text_indices, start=-len(text_indices)):
            retval[pos] = self.layout.annotations[i]
        for i, pos in enumerate(text3d_indices, start=-len(text3d_indices)):
            retval[pos] = self.layout.scene.annotations[i]
        if len(retval) == 1:
            return retval[0]
        return retval

    def __iadd__(self, item):
        if isinstance(item, plotly.graph_objects.layout.Annotation):
            self.add_annotation(item)
        elif isinstance(item, plotly.graph_objects.layout.scene.Annotation):
            self.layout.scene.annotations += (item,)
        else:
            self.add_trace(item)
        return self

    def axes_labels(self, *labels):
        if len(labels) == 2:
            self.layout.xaxis.title.text = labels[0]
            self.layout.yaxis.title.text = labels[1]
        elif len(labels) == 3:
            self.layout.scene.xaxis.title.text = labels[0]
            self.layout.scene.yaxis.title.text = labels[1]
            self.layout.scene.zaxis.title.text = labels[2]
        else:
            raise ValueError("You must specify labels for either 2 or 3 axes.")

    def axes_ranges(self, *ranges, scale=None):
        if len(ranges) == 2:
            (xmin, xmax), (ymin, ymax) = ranges
            self.layout.xaxis.range = (xmin, xmax)
            self.layout.yaxis.range = (ymin, ymax)
            if scale is not None:
                x, y = scale
                self.layout.xaxis.constrain = "domain"
                self.layout.yaxis.constrain = "domain"
                self.layout.yaxis.scaleanchor = "x"
                self.layout.yaxis.scaleratio = y / x
        elif len(ranges) == 3:
            (xmin, xmax), (ymin, ymax), (zmin, zmax) = ranges
            self.layout.scene.xaxis.range = (xmin, xmax)
            self.layout.scene.yaxis.range = (ymin, ymax)
            self.layout.scene.zaxis.range = (zmin, zmax)
            if isinstance(scale, str):
                self.layout.scene.aspectmode = scale
            elif scale is not None:
                x, y, z = scale
                x *= xmax - xmin
                y *= ymax - ymin
                z *= zmax - zmin
                c = sorted((x, y, z))[1]
                self.layout.scene.aspectmode = "manual"
                self.layout.scene.aspectratio.update(x=x/c, y=y/c, z=z/c)
        else:
            raise ValueError("You must specify ranges for either 2 or 3 axes.")


class Figure(plotly.graph_objects.Figure, MyFigure):
    def __init__(self, *args, **kwargs):
        specs = kwargs.pop("subplots", None)
        if specs is None:
            super().__init__(*args, **kwargs)
        else:
            try:
                specs[0][0]
            except:
                specs = [specs]
            rows = len(specs)
            cols = len(specs[0])
            fig = plotly.subplots.make_subplots(rows=rows, cols=cols, 
                                                specs=specs, **kwargs)
            super().__init__(fig)


class FigureWidget(plotly.graph_objects.FigureWidget, MyFigure):
    def __init__(self, *args, **kwargs):
        self._auto_items = []
        specs = kwargs.pop("subplots", None)
        if specs is None:
            super().__init__(*args, **kwargs)
        else:
            try:
                specs[0][0]
            except:
                specs = [specs]
            rows = len(specs)
            cols = len(specs[0])
            fig = plotly.subplots.make_subplots(rows=rows, cols=cols, 
                                                specs=specs, **kwargs)
            super().__init__(fig)

    def auto_update(self, *items):
        if self._auto_items:
            for item, newitem in zip(self._auto_items, items):
                if isinstance(newitem, tuple):
                    newitem = newitem[0]
                item.update(newitem)
        else:
            for item in items:
                if isinstance(item, tuple):
                    item, subplot = item
                    self._auto_items.append(self.add(item, subplot=subplot))
                else:
                    self._auto_items.append(self.add(item))

    def initialized(self):
        return len(self._auto_items) > 0


# Below are all the actual plotting methods. First, the 2D graphics:
def plotly_text(text, location, **options):
    options = options.copy()
    x, y = np.array(location, dtype=float)
    if options.pop("update", False):
        return dict(text=text, x=x, y=y)
    options.setdefault("font_color", options.pop("color", "black"))
    size = options.pop("size", None)
    if size is not None:
        options.setdefault("font_size", size)
    arrow = options.pop("arrow", None)
    if arrow:
        options.setdefault("ax", float(arrow[0]))
        options.setdefault("ay", float(arrow[1]))
        options.setdefault("showarrow", True)
    else:
        options.setdefault("showarrow", False)
    if options.pop("paper", False):
        options.update(xref="paper", yref="paper")
        options.setdefault("xanchor", "left")
        options.setdefault("yanchor", "bottom")
    return plotly.graph_objects.layout.Annotation(
            text=text, x=x, y=y, **options)


def plotly_points(points, **options):
    options = options.copy()
    x, y = np.array(points, dtype=float).transpose()
    if options.pop("update", False):
        return dict(x=x, y=y)
    options.setdefault("marker_color", options.pop("color", "black"))
    options.setdefault("marker_size", options.pop("size", 8))
    options.setdefault("mode", "markers")
    return plotly.graph_objects.Scatter(x=x, y=y, **options)


def plotly_lines(points, **options):
    options = options.copy()
    x, y = np.array(points, dtype=float).transpose()
    if options.pop("update", False):
        return dict(x=x, y=y)
    color = options.pop("color", "blue")
    options.setdefault("line_color", color)
    options.setdefault("mode", "lines")
    return plotly.graph_objects.Scatter(x=x, y=y, **options)


def plotly_vector(vec, start=(0, 0), axes_scale=(1, 1), **options):
    options = options.copy()
    options["mode"] = "lines"
    options["fill"] = "toself"
    color = options.pop("color", "black")
    options.setdefault("line_color", color)
    options.setdefault("fillcolor", color)
    options.setdefault("line_width", 1)

    tipfactor = 0.2  # Arrowhead length, as a fraction of the total arrow length
    tipslope1 = 2.5  # For an arrow pointing up, the slope of the leading edge of the left side of the arrowhead
    tipslope2 = 1.0  # For an arrow pointing up, the slope of the trailing edge of the left side of the arrowhead

    vec = np.array(vec, dtype=float)
    start = np.array(start, dtype=float)
    end = start + vec
    NaN = np.array((np.nan, np.nan), dtype=float)
    axes_scale = axes_scale[1] / axes_scale[0]
    tip = tipfactor * vec
    tip_perp = np.array(( -axes_scale/tipslope1 * tip[1], 
                         1/axes_scale/tipslope1 * tip[0]), dtype=float)
    tipleft = end - tip + tip_perp
    tipmiddle = end - float(1 - tipslope2/tipslope1) * tip
    tipright = end - tip - tip_perp
    xy = np.array((start, end, NaN, end, tipleft, tipmiddle, tipright, end))
    x, y = xy.transpose()
    if options.pop("update", False):
        return dict(x=x, y=y)
    return plotly.graph_objects.Scatter(x=x, y=y, **options)


def plotly_vector_field(f, x_range, y_range, **options):
    options = options.copy()
    axes_scale = options.pop("axes_scale", None)
    axes_aspect = options.pop("axes_aspect", (4, 3))
    scalefactor = options.pop("scalefactor", 1)
    region = options.pop("region", None)
    plotpoints = options.pop("plotpoints", 21)
    try:
        plotpointsx, plotpointsy = plotpoints
    except:
        plotpointsx = plotpointsy = plotpoints
    options["mode"] = "lines"
    options["fill"] = "toself"
    color = options.pop("color", "limegreen")
    options.setdefault("line_color", color)
    options.setdefault("fillcolor", color)
    options.setdefault("line_width", 1)

    # Arrow size/shape parameters
    zero_threshold = 1e-12
    minlength = 0.06 # Tip size won't be smaller than for an arrow of this length, as a fraction of the graph diagonal
    maxlength = 0.16 # Tip size won't be larger than for an arrow of this length, as a fraction of the graph diagonal
    tipfactor = 0.2  # Arrowhead length, as a fraction of the total arrow length
    tipslope1 = 2.5  # For an arrow pointing up, the slope of the leading edge of the left side of the arrowhead
    tipslope2 = 1.0  # For an arrow pointing up, the slope of the trailing edge of the left side of the arrowhead

    # Set the axes ranges, and the axes_scale and axes_aspect
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    if axes_scale is None:
        axes_scale = axes_aspect[1] / axes_aspect[0] * (xmax - xmin) / (ymax - ymin)
    else:
        axes_scale = axes_scale[1] / axes_scale[0]
    axes_aspect = np.array((xmax - xmin, axes_scale * (ymax - ymin)), dtype=float)
    diagonal = np.linalg.norm(axes_aspect)

    # Choose the start and step size, for both x and y, for the grid points
    def grid_params(min1, max1, min2, max2, ratio):
        step1 = (max1 - min1) / (plotpoints + 0.75)
        start1 = min1 + 0.375 * step1
        step2 = step1 / ratio
        remainder = divmod(float((max2 - min2) / step2), int(1))[1] / 2
        if remainder < 0.125:
            remainder += 0.5
        start2 = min2 + remainder * step2
        return start1, step1, start2, step2
    try:
        plotpoints[1]
    except: # A single value was given for plotpoints. Set grid automatically. 
        if axes_aspect[0] > axes_aspect[1]:
            x0, xstep, y0, ystep = grid_params(xmin, xmax, ymin, ymax, axes_scale)
        else:
            y0, ystep, x0, xstep = grid_params(ymin, ymax, xmin, xmax, 1/axes_scale)
    else: # A 2-tuple was given for plotpoints. Use exactly that many. 
        xstep = (xmax - xmin) / (plotpoints[0] + 0.75)
        x0 = xmin + 0.375 * xstep
        ystep = (ymax - ymin) / (plotpoints[1] + 0.75)
        y0 = ymin + 0.375 * ystep

    # The actual computations
    f1, f2 = [fast_float(f_, x, y) for f_ in f]
    xvals = np.arange(float(x0), float(xmax), float(xstep))
    yvals = np.arange(float(y0), float(ymax), float(ystep))
    xy = np.array(np.meshgrid(xvals, yvals)).reshape(2, -1).transpose()
    if region:
        region = fast_float(region, x, y)
        xy = xy[[region(x0, y0) > 0 for x0, y0 in xy]]
    vec = np.array([(f1(x0, y0), f2(x0, y0)) for x0, y0 in xy]).transpose()
    keep = np.linalg.norm(vec, axis=0) > zero_threshold
    xy = xy[keep].transpose()
    vec = vec[:,keep]
    grid_scale = np.array((xstep, ystep), dtype=float).reshape(2, 1)
    longest_vector = np.linalg.norm(vec / grid_scale, axis=0, ord=np.inf).max()
    vec *= scalefactor / longest_vector
    paper_scale = np.array((1, axes_scale), dtype=float).reshape(2, 1)
    length = np.linalg.norm(vec * paper_scale, axis=0) / diagonal
    tip = (tipfactor * np.clip(length, minlength, maxlength) / length) * vec
    tip_perp = np.array((float( -axes_scale/tipslope1) * tip[1], 
                         float(1/axes_scale/tipslope1) * tip[0]))
    NaN = np.full_like(xy, np.nan)
    end = xy + vec
    tipleft = end - tip + tip_perp
    tipmiddle = end - float(1 - tipslope2/tipslope1) * tip
    tipright = end - tip - tip_perp
    xy = np.array((xy, end, NaN, end, tipleft, tipmiddle, tipright, end, NaN))
    x, y = np.moveaxis(xy, 0, 2).reshape(2, -1)
    if options.pop("update", False):
        return dict(x=x, y=y)
    return plotly.graph_objects.Scatter(x=x, y=y, **options)


def plotly_points3d(points, **options):
    options = options.copy()
    options.setdefault("marker_color", options.pop("color", "black"))
    options.setdefault("marker_size", options.pop("size", 2.5))
    options.setdefault("mode", "markers")
    x, y, z = np.array(points, dtype=float).transpose()
    if options.pop("update", False):
        return dict(x=x, y=y, z=z)
    return plotly.graph_objects.Scatter3d(x=x, y=y, z=z, **options)


def plotly_lines3d(points, **options):
    options = options.copy()
    color = options.pop("color", "blue")
    options.setdefault("line_color", color)
    options.setdefault("mode", "lines")
    x, y, z = np.array(points, dtype=float).transpose()
    if options.pop("update", False):
        return dict(x=x, y=y, z=z)
    return plotly.graph_objects.Scatter3d(x=x, y=y, z=z, **options)


def plotly_function3d(f, x_range, y_range, **options):
    options = options.copy()
    plotpoints = options.pop("plotpoints", 81)
    try:
        plotpointsx, plotpointsy = plotpoints
    except:
        plotpointsx = plotpointsy = plotpoints
    color = options.pop("color", "lightblue")
    options.setdefault("colorscale", (color, color))
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    f = fast_float(f, x, y)
    x = np.linspace(float(xmin), float(xmax), plotpointsx)
    y = np.linspace(float(ymin), float(ymax), plotpointsy)
    x, y = np.meshgrid(x, y)
    xy = np.array([x.flatten(), y.flatten()]).transpose()
    z = np.array([f(x0, y0) for x0, y0 in xy]).reshape(plotpointsy, plotpointsx)
    if options.pop("update", False):
        return dict(x=x, y=y, z=z)
    return plotly.graph_objects.Surface(x=x, y=y, z=z, **options)


def plotly_vector_field_on_surface3d(f, field, u_range, v_range, **options):
    options = options.copy()
    plotpoints = options.pop("plotpoints", 21)
    try:
        plotpointsu, plotpointsv = plotpoints
    except:
        plotpointsu = plotpointsv = plotpoints
    color = options.pop("color", "limegreen")
    options.setdefault("colorscale", (color, color))
    u, umin, umax = u_range
    v, vmin, vmax = v_range
    f1, f2, f3 = [fast_float(f_, u, v) for f_ in f]
    if len(field) == 2:
        field = f.derivative() * field
    fx, fy, fz = [fast_float(f_, u, v) for f_ in field]
    u = np.linspace(float(umin), float(umax), plotpointsu)
    v = np.linspace(float(vmin), float(vmax), plotpointsv)
    uv = np.array(np.meshgrid(u, v)).reshape(2, -1).transpose()
    points = np.array([(f1(u0, v0), f2(u0, v0), f3(u0, v0)) for u0, v0 in uv])
    vectors = np.array([(fx(u0, v0), fy(u0, v0), fz(u0, v0)) for u0, v0 in uv])
    x, y, z = points.transpose()
    u, v, w = vectors.transpose()
    if options.pop("update", False):
        return dict(x=x, y=y, z=z, u=u, v=v, w=w)
    return plotly.graph_objects.Cone(x=x, y=y, z=z, u=u, v=v, w=w, **options)

x, y, z = var("x, y, z")
bump(x0, y0, h) = h^2/(h^2 + (x - x0)^2 + (y - y0)^2)
f(x, y) = 1 + 3*bump(2.5, 2, 1.5) + 5*bump(6.5, 3, 2) - 2.5*bump(4.5, 5, 2) + 2*bump(1, 7, 1)

crit_points = find_zeros(f.gradient(), (x, 0, 8), (y, 0, 8))
for pt in crit_points:
    print(pt)

(1.025749249792184, 4.924388410447029)
(0.9850627229439552, 6.988915199423456)
(2.5888506099044983, 1.966106259573833)
(4.223844149137523, 2.0532217928719096)
(4.278189612110433, 5.5133658104706305)
(6.540423395003371, 2.883242346820021)

hessian = jacobian(f.gradient(), (x, y))
critpt_dict = {"Maxima": [], "Minima": [], "Saddle points": []}
for (x0, y0) in crit_points:
    e1, e2 = hessian(x0, y0).eigenvalues()
    if e1 < 0 and e2 < 0:
        critpt_dict["Maxima"].append((x0, y0, f(x0, y0)))
    elif e1 > 0 and e2 > 0:
        critpt_dict["Minima"].append((x0, y0, f(x0, y0)))
    else:
        critpt_dict["Saddle points"].append((x0, y0, f(x0, y0)))

darkorange = colors["darkorange"].darker().html_color()
def gradient3d_interactive(f, x_range, y_range, z_range, critpt_dict, **options):
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    z, zmin, zmax = z_range
    axes_labels = options.get("axes_labels", (str(x), str(y), str(z)))
    scale = options.get("scale", None)
    colors = {"Maxima": "darkgreen", "Minima": "darkred", "Saddle points": darkorange}

    figure = FigureWidget()
    figure.axes_labels(*axes_labels)
    figure.axes_ranges((xmin, xmax), (ymin, ymax), (zmin, zmax), scale=scale)
    figure.layout.margin = dict(t=10, b=10, r=10, l=10)
    figure.layout.showlegend = False

    surface = figure.add(plotly_function3d(f, x_range, y_range))
    points = {}
    xypoints = {}
    xylines = {}
    for kind, critpts in critpt_dict.items():
        points[kind] = figure.add(plotly_points3d(critpts, color=colors[kind]))
        xypts = [(x0, y0, 0) for (x0, y0, z0) in critpts]
        lines = [(a, b, (np.nan, np.nan, np.nan)) for a, b in zip(critpts, xypts)]
        lines = np.array(lines).reshape(-1, 3)
        xypoints[kind] = figure.add(plotly_points3d(xypts, color=colors[kind]))
        xylines[kind] = figure.add(plotly_lines3d(lines, color=colors[kind], line_dash="dash"))
    xyplane(x, y) = (x, y, 0)
    field = figure.add(plotly_vector_field_on_surface3d(xyplane, f.gradient(), 
            x_range, y_range, color="limegreen"))

    choices = ("None",) + tuple(critpt_dict.keys()) + ("All",)
    @interact(which=selector(choices, label="Show: "), 
              show_xypoints=checkbox(default=False, label="Show critical points"), 
              show_gradient=checkbox(default=False, label="Show gradient vector field"))
    def update(which, show_xypoints, show_gradient):
        with figure.batch_update():
            surface.opacity = 0.6 if (show_xypoints or show_gradient) else 1
            field.visible = show_gradient
            for kind in critpt_dict:
                visible = (which == kind or which == "All")
                points[kind].visible = visible
                xypoints[kind].visible = show_xypoints and visible
                xylines[kind].visible = show_xypoints and visible

    return figure

def gradient2d_interactive(f, x_range, y_range, critpt_dict, **options):
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    axes_labels = options.get("axes_labels", (f"${x}$", f"${y}$"))
    scale = options.get("scale", None)
    zoom_factor = options.get("zoom_factor", 10)
    zoom_xrange = (xmax - xmin) / zoom_factor
    zoom_yrange = (ymax - ymin) / zoom_factor
    colors = {"Maxima": "green", "Minima": "red", "Saddle points": "darkorange"}

    figure = FigureWidget()
    figure.axes_labels(*axes_labels)
    figure.axes_ranges((xmin, xmax), (ymin, ymax), scale=scale)
    figure.layout.margin = dict(t=30, b=40, r=10, l=10)
    figure.layout.dragmode = "select" # For this figure, we want select, not pan

    field = figure.add(plotly_vector_field(f.gradient(), x_range, y_range, 
            color="limegreen", plotpoints=17, axes_scale=scale, showlegend=False))
    points = []
    for kind, critpts in critpt_dict.items():
        points.append(figure.add(plotly_points([(x0, y0) for (x0, y0, z0) in critpts], 
                name=kind, color=colors[kind])))

    def zoomin(widget):
        selected = []
        for critpts in points:
            if critpts.selectedpoints:
                for i in critpts.selectedpoints:
                    selected.append((critpts.x[i], critpts.y[i]))
        if len(selected) != 1:
            return
        x0, y0 = selected[0]
        new_xrange = (x, x0 - zoom_xrange/2, x0 + zoom_xrange/2)
        new_yrange = (y, y0 - zoom_yrange/2, y0 + zoom_yrange/2)
        button.description = "Zoom out"
        button.on_click(zoomin, remove=True)
        button.on_click(zoomout)
        with figure.batch_update():
            figure.axes_ranges(new_xrange[1:], new_yrange[1:], scale=scale)
            field.update(plotly_vector_field(f.gradient(), new_xrange, new_yrange, 
                    plotpoints=11, axes_scale=scale, update=True))

    def zoomout(widget):
        button.description = "Zoom in"
        button.on_click(zoomout, remove=True)
        button.on_click(zoomin)
        with figure.batch_update():
            figure.axes_ranges((xmin, xmax), (ymin, ymax), scale=scale)
            field.update(plotly_vector_field(f.gradient(), x_range, y_range, 
                    plotpoints=17, axes_scale=scale, update=True))

    button = Button(description="Zoom in")
    button.on_click(zoomin)

    return HBox((button, figure))

def vectorfield_interactive(f, x_range, y_range, **options):
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    axes_labels = options.get("axes_labels", (f"${x}$", f"${y}$"))
    axes_scale = options.get("axes_scale", None)
    scalefactor = options.get("scalefactor", None)
    funcname = options.get("label", "f")

    figure = FigureWidget()
    figure.axes_labels(*axes_labels)
    figure.axes_ranges((xmin, xmax), (ymin, ymax), scale=axes_scale)
    field = figure.add(plotly_vector_field(f, x_range, y_range, name="Vector field", 
            plotpoints=11, color="limegreen", axes_scale=axes_scale, visible="legendonly"))
    point = figure.add(plotly_points([(xmin, ymin)], showlegend=False))

    def update_point(change):
        value = change["new"].strip()
        if value[0] != "(" or value[-1] != ")":
            return
        try:
            x0, y0 = [float(num) for num in value[1:-1].split(",")]
        except:
            return
        with figure.batch_update():
            point.x = [x0]
            point.y = [y0]

    def add_vector(widget):
        x0, y0 = point.x[0], point.y[0]
        vec = f(x0, y0) * scalefactor
        with figure.batch_update():
            figure.add(plotly_vector(vec, start=(x0, y0), color="purple", line_width=2, 
                    axes_scale=axes_scale, showlegend=False))

    point_input = Text(description="Point:")
    point_input.observe(update_point, names="value")
    point_input.value = f"({round(xmin)},{round(ymin)})"
    button = Button(description="Add vector")
    button.on_click(add_vector)
    label = fr"{latex(f(x,y)[0])} \\ {latex(f(x,y)[1])}"
    label = fr"$\qquad {funcname}({x}, {y}) = \begin{{bmatrix}} {label} \end{{bmatrix}}$"
    label = HTMLMath(label)

    return VBox((HBox((point_input, button, label)), figure))

Learning goals:

• Be able to compute the gradient vector field of a function $f\!: \mathbb{R}^n \to \mathbb{R}$.

• Be able to explain how the direction of the gradient vector field of $f$ relates to the “slope” of the function $f$.

• Be able to use the gradient vector field of $f$ to determine whether a critical point of $f$ is a local maximum, local minimum, or saddle point.

gradient3d_interactive(f, (x, 0, 8), (y, 0, 8), (z, 0, 6), critpt_dict, scale=(1,1,1))

The gradient of a function $f\!: \mathbb{R}^n \to \mathbb{R}$

Consider a function $f$ of two (or more) variables, say $f(x, y)$.

We know $f$ will have two partial derivatives: $$\frac{\partial f}{\partial x} \qquad \text{and} \qquad \frac{\partial f}{\partial y}$$

And generally, each of these partial derivatives will also be a function from $\mathbb{R}^2$ to $\mathbb{R}$. So if we put these two functions together into a vector:

$$\begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{bmatrix}$$

the result will be a function from from $\mathbb{R}^2$ to $\mathbb{R}^2$... a vector field!

This vector field is called the gradient of $f$.

g(x, y) = 5*x^2 + 20*y^2 - 1/4*x^3 - 2*y^3 - 1/2*x^2*y^2
H = g.derivative(2)
for x0, y0 in solve(list(g.derivative()(x, y)), (x, y)):
    x0, y0 = x0.rhs(), y0.rhs()
    if not (x0 in RR and y0 in RR):
        continue
    print(f"({x0}, {y0}): Eigenvalues {H(x0,y0).eigenvalues()}")

(0, 0): Eigenvalues [10, 40]
(40/3, 0): Eigenvalues [-1240/9, -10]
(0, 20/3): Eigenvalues [-310/9, -40]
(1/2*sqrt(133) + 3/2, -1/4*sqrt(133) + 3/4): Eigenvalues [9/16*sqrt(133) - 1/16*sqrt(-810*sqrt(133) + 276670) - 45/16, 9/16*sqrt(133) + 1/16*sqrt(-810*sqrt(133) + 276670) - 45/16]
(-1/2*sqrt(133) + 3/2, 1/4*sqrt(133) + 3/4): Eigenvalues [-9/16*sqrt(133) - 1/16*sqrt(810*sqrt(133) + 276670) - 45/16, -9/16*sqrt(133) + 1/16*sqrt(810*sqrt(133) + 276670) - 45/16]
(-8, -4): Eigenvalues [-sqrt(4177) + 15, sqrt(4177) + 15]
(5, 5/2): Eigenvalues [-35, 65/4]

An example of the gradient of a function

Define $f\!: \mathbb{R}^2 \to \mathbb{R}$ by $$f(x, y) = 5x^2 + 20y^2 - \tfrac{1}{4} x^3 - 2y^3 - \tfrac{1}{2} x^2 y^2$$

The partial derivative of $f$ with respect to $x$ is$\quad \frac{\partial f}{\partial x} = 10x - \tfrac{3}{4} x^2 - x y^2$

And the partial derivative of $f$ with respect to $y$ is $\quad \frac{\partial f}{\partial y} = 40y - 6y^2 - x^2 y$

So the gradient of $f$ is a function, which we'll write as $\mathrm{grad}f$, defined by

$$\mathrm{grad}f (\begin{bmatrix} x \\ y \end{bmatrix}) = \begin{bmatrix} 10x - \tfrac{3}{4} x^2 - x y^2 \\ 40y - 6y^2 - x^2 y \end{bmatrix}$$

options = dict(scalefactor=0.005, axes_scale=(1,1), label=r"\textrm{grad}f")
vectorfield_interactive(g.derivative(), (x, 3.6, 6.4), (y, 1.8, 3.2), **options)

Notation and a few details about the gradient

Notation: The gradient of $f$ is written as $\ \mathrm{grad}f \$ or $\ \vec{\nabla}f$.

So for our previous example, we could write

$$\mathrm{grad}f (x, y) = \begin{bmatrix} 10x - \tfrac{3}{4} x^2 - x y^2 \\ 40y - 6y^2 - x^2 y \end{bmatrix} \qquad \text{or} \qquad \vec{\nabla}f (x, y) = \begin{bmatrix} 10x - \tfrac{3}{4} x^2 - x y^2 \\ 40y - 6y^2 - x^2 y \end{bmatrix}$$

Dimensions: If $f\!: \mathbb{R}^n \to \mathbb{R}$, then $\mathrm{grad} f\!: \mathbb{R}^n \to \mathbb{R}^n$ is defined by

$$\mathrm{grad}f = \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \vdots \end{bmatrix}$$

The gradient, and critical points of $f$

Recall that to find the critical points of $f$, we set both partial derivatives to $0$, and solve simultaneously.

Note that this is the same way we find equilibrium points of a system of differential equations. So now we can say that finding the critical points of $f$ is the same as finding the equilibrium points of the $\,\mathrm{grad}f \,$ vector field!

How to visualize the gradient vector field

gradient3d_interactive(f, (x, 0, 8), (y, 0, 8), (z, 0, 6), critpt_dict, scale=(1,1,1))

How to visualize the gradient vector field

gradient2d_interactive(f, (x, 0, 8), (y, 0, 8), critpt_dict, scale=(1,1))

The direction of the gradient vector field

At any point in the domain of $f$, there is a gradient vector. It points in the direction in which the graph of $f$ has its steepest slope.

In other words, in short, it points in the “uphill” direction.

This will allow us to classify critical points of $f$:

• A local maximum of $f$ will be a stable equilibrium point (sink) of the $\,\mathrm{grad}f \,$ vector field.
• A local minimum of $f$ will be an unstable equilibrium point (source) of the $\,\mathrm{grad}f \,$ vector field.
• A saddle point of $f$ will be a saddle point of the $\,\mathrm{grad}f \,$ vector field.

Conclusions:

• For a function $f\!: \mathbb{R}^n \to \mathbb{R}$, its gradient is the vector field $\mathrm{grad}f\!: \mathbb{R}^n \to \mathbb{R}^n$ defined by $$\mathrm{grad}f (x, y, \dotsc) = \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \vdots \end{bmatrix}$$

• Therefore the critical points of $f$ are the same as the equilibrium points of the gradient vector field of $f$.

• At any point in the domain of $f$, the direction of $\mathrm{grad}f$ at that point is the direction in which $f$ increases the fastest. (greatest rate of change, highest slope)

• Therefore we can use the gradient of $f$ to classify the critical points of $f$:

  • A local maximum of $f$ will be a stable equilibrium point (sink) of the gradient vector field of $f$.

  • A local minimum of $f$ will be an unstable equilibrium point (source) of the gradient vector field of $f$.

  • A saddle point of $f$ will be a saddle point of the gradient vector field of $f$.