Mini-lecture 7-16 - Critical points of functions of multiple variables.ipynb

Kernel: SageMath 9.3

from itertools import product
import numpy as np
import plotly.graph_objects
from ipywidgets import HBox, 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


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_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)

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": [], "Other": []}
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["Other"].append((x0, y0, f(x0, y0)))

darkorange = colors["darkorange"].darker().html_color()
def show_critpts_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", "Other": 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

    figure.add(plotly_function3d(f, (x, xmin, xmax), (y, ymin, ymax)))
    points = {}
    planes = {}
    for kind, critpts in critpt_dict.items():
        points[kind] = figure.add(plotly_points3d(critpts, color=colors[kind]))
        planes[kind] = []
        for (x0, y0, z0) in critpts:
            newplane = plotly_function3d(z0, (x, x0-3, x0+3), (y, y0-3, y0+3), 
                    color=colors[kind], plotpoints=2, opacity=0.3)
            planes[kind].append(figure.add(newplane))

    choices = ("None",) + tuple(critpt_dict.keys()) + ("All",)
    @interact(which=selector(choices, label="Show: "), 
              show_tangent=checkbox(default=False, label="Show tangent plane(s)"))
    def update(which, show_tangent):
        with figure.batch_update():
            for kind, critpts in points.items():
                critpts.visible = (which == kind or which == "All")
                for plane in planes[kind]:
                    plane.visible = show_tangent and (which == kind or which == "All")

    return figure

Learning goals:

• Be able to define critical points for a function of two (or more) variables $f\!: \mathbb{R}^n \to \mathbb{R}$.

• Be able to find all of the critical points of a function $f\!: \mathbb{R}^n \to \mathbb{R}$.

• Be able to identify the three main types of critical points that can occur for a function $f\!: \mathbb{R}^n \to \mathbb{R}$.

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

Critical points of a function $f\!: \mathbb{R}^n \to \mathbb{R}$

Definition: Let $f\!: \mathbb{R}^n \to \mathbb{R}$. A point $(x^*, y^*, \dotsc)$ is a critical point of $f$ if every partial derivative of $f$ (i.e. $\frac{\partial f}{\partial x}$, $\frac{\partial f}{\partial y}$, etc...) is either $0$ or undefined at that point.

Note: We will only focus on the kind of critical points at which all of the partial derivatives are $0$.

Just as in the single-variable case, there is a theorem that says that any local maximum or minimum of $f$ must occur at a critical point.

(If $f$ is defined on some closed domain, then local maxima/minima can also occur at the boundary points of the domain. This can be much more complicated in the multivariable setting, since the domain isn't just an interval. So we will not consider this case either.)

Conclusions

For a function $f\!:\mathbb{R}^n \to \mathbb{R}$...

• A critical point of $f$ is a point in the domain of $f$ at which every partial derivative of $f$ is either $0$ or undefined. (We will only focus on the kind where the partial derivatives are $0$.) Each local maximum and minimum of $f$ must occur at a critical point.

• Therefore, to find the critical points of $f$, take all of its partial derivatives $\frac{\partial f}{\partial x}$, $\frac{\partial f}{\partial y}$, etc, set them all to $0$, and solve simultaneously. Note that this is very similar to finding the equilibrium points of a system of differential equations.

• A critical point of $f$ can be a local maximum, or a local minimum, or a saddle point. So, just as in the single-variable case, we need a way to classify critical points.