{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ ], "source": [ "from itertools import product\n", "import numpy as np\n", "import plotly.graph_objects\n", "from ipywidgets import HBox, Layout, interact as _original_interact" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": { }, "output_type": "execute_result" } ], "source": [ "%%html\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "_original_interact = interact\n", "def interact(_function_to_wrap=None, _layout=\"horizontal\", **kwargs):\n", " \"\"\"interact, but with widgets laid out in a horizontal flexbox layout\n", "\n", " This function works exactly like 'interact' (from SageMath or ipywidgets), \n", " except that instead of putting all of the widgets into a vertical box \n", " (VBox), it uses a horizontal box (HBox) by default. The HBox uses a flexbox \n", " layout, so that if there are many widgets, they'll wrap onto a second row. \n", "\n", " Options:\n", " '_layout' - 'horizontal' by default. Anything else, and it will revert \n", " back to using the default layout of 'interact' (a VBox). \n", " \"\"\"\n", " def decorator(f):\n", " retval = _original_interact(f, **kwargs)\n", " if _layout == \"horizontal\":\n", " widgets = retval.widget.children[:-1]\n", " output = retval.widget.children[-1]\n", " hbox = HBox(widgets, layout=Layout(flex_flow=\"row wrap\"))\n", " retval.widget.children = (hbox, output)\n", " return retval\n", " if _function_to_wrap is None:\n", " # No function passed in, so this function must *return* a decorator\n", " return decorator\n", " # This function was called directly, *or* was used as a decorator directly\n", " return decorator(_function_to_wrap)\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "def find_zeros(field, *ranges, **options):\n", " \"\"\"Numerically approximate all zeros of a vector field within a box\n", "\n", " Each 'range' should have the form (x, xmin, xmax), where 'x' is one of the\n", " variables appearing in the vector field, and 'xmin' and 'xmax' are the\n", " bounds on that variables.\n", "\n", " Options:\n", " 'intervals' - How many subintervals to use in each dimension (default\n", " 20). If there are zeros that are very close together, you may need\n", " to increase this to find them. But be warned that the running time\n", " is roughly this to the n power, where n is the number of variables.\n", " 'tolerance' - How close to approximate roots, roughly (default 1e-4)\n", " 'maxiter' - Maximum number of iterations of Newton's Method to run for\n", " any one (potential) solution. This defaults to -2*log(tolerance).\n", " There usually isn't much harm in increasing this, unless there are\n", " many false hits.\n", " 'round' - Round the components of the solutions to this many decimal\n", " places (default None, meaning do not round them at all)\n", " \"\"\"\n", "\n", " # Initialization\n", " intervals = options.get(\"intervals\", 20)\n", " tolerance = options.get(\"tolerance\", 1e-5)\n", " maxiter = options.get(\"maxiter\", int(round(-2*log(tolerance))))\n", " roundto = options.get(\"round\", None)\n", " n = len(ranges)\n", " if len(field) != n:\n", " raise ValueError(\"Dimension of vector field is {}, but {} ranges \"\n", " \"given\".format(len(field), n))\n", "\n", " mins = [xmin - (xmax - xmin)/intervals*tolerance for x, xmin, xmax in ranges]\n", " maxes = [xmax + (xmax - xmin)/intervals*tolerance for x, xmin, xmax in ranges]\n", " deltas = [(xmax - xmin) / intervals for xmin, xmax in zip(mins, maxes)]\n", " powers = [1 << i for i in range(n)]\n", " J = jacobian(field, [x for x, xmin, xmax in ranges])\n", " def dist(v, w):\n", " return sqrt(sum(((a - b) / d)**2 for a, b, d in zip(v, w, deltas)))\n", "\n", " # Set up the array of positive/negative signs of the vector field\n", " signs = np.zeros((intervals + 1,) * n, dtype=int)\n", " sranges = [srange(m, m + d*(intervals + 0.5), d) for m, d in zip(mins, deltas)]\n", " for vertex, index in zip(product(*sranges), product(range(intervals + 1), repeat=n)):\n", " v = field(*vertex)\n", " signs[index] = sum(powers[i] for i in range(n) if v[i] > 0)\n", "\n", " # Now search through that array for potential solutions\n", " solutions = []\n", " mask = int((1 << n) - 1)\n", " for index in product(range(intervals), repeat=n):\n", " indexpowers = list(zip(index, powers))\n", " all0s = 0\n", " all1s = mask\n", " for k in range(mask + 1):\n", " newindex = [i + (1 if k & p else 0) for i, p in indexpowers]\n", " vertex = signs[tuple(newindex)]\n", " all0s |= vertex\n", " all1s &= vertex\n", " if all0s & ~all1s == mask:\n", " # Now do Newton's method!\n", " v = vector(m + d*(i + 0.5) for i, m, d in zip(index, mins, deltas))\n", " for i in range(maxiter):\n", " previous_v = v\n", " v = v - J(*v).solve_right(field(*v))\n", " if dist(v, previous_v) < tolerance:\n", " break\n", " else:\n", " warn(\"{} iterations reached without convergence for solution \"\n", " \"{}\".format(maxiter, v), RuntimeWarning)\n", " if solutions and min(dist(v, w) for w in solutions) < 2*tolerance:\n", " continue\n", " if not all(xmin <= x <= xmax for x, xmin, xmax in zip(v, mins, maxes)):\n", " continue\n", " solutions.append(vector(RDF, v))\n", "\n", " # A convenience: round to some number of decimal places, if requested\n", " if roundto is not None:\n", " solutions = [vector(round(x, roundto) for x in v) for v in solutions]\n", " return solutions\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# Our mixin class for Figure and FigureWidget. Should never be instantiated! \n", "class MyFigure(object):\n", " def add(self, *items, subplot=None):\n", " if subplot is not None:\n", " try:\n", " subplot[0]\n", " except:\n", " subplot = (1, subplot)\n", " row = int(subplot[0])\n", " col = int(subplot[1])\n", " retval = []\n", " text_indices = []\n", " text3d_indices = []\n", " for item in items:\n", " if isinstance(item, plotly.graph_objects.layout.Annotation):\n", " if subplot is None:\n", " self.add_annotation(item)\n", " else:\n", " self.add_annotation(item, row=row, col=col)\n", " text_indices.append(len(retval))\n", " retval.append(None)\n", " elif isinstance(item, plotly.graph_objects.layout.scene.Annotation):\n", " self.layout.scene.annotations += (item,)\n", " text3d_indices.append(len(retval))\n", " retval.append(None)\n", " else:\n", " if subplot is None:\n", " self.add_trace(item)\n", " else:\n", " self.add_trace(item, row=row, col=col)\n", " retval.append(self.data[-1])\n", " for i, pos in enumerate(text_indices, start=-len(text_indices)):\n", " retval[pos] = self.layout.annotations[i]\n", " for i, pos in enumerate(text3d_indices, start=-len(text3d_indices)):\n", " retval[pos] = self.layout.scene.annotations[i]\n", " if len(retval) == 1:\n", " return retval[0]\n", " return retval\n", "\n", " def __iadd__(self, item):\n", " if isinstance(item, plotly.graph_objects.layout.Annotation):\n", " self.add_annotation(item)\n", " elif isinstance(item, plotly.graph_objects.layout.scene.Annotation):\n", " self.layout.scene.annotations += (item,)\n", " else:\n", " self.add_trace(item)\n", " return self\n", "\n", " def axes_labels(self, *labels):\n", " if len(labels) == 2:\n", " self.layout.xaxis.title.text = labels[0]\n", " self.layout.yaxis.title.text = labels[1]\n", " elif len(labels) == 3:\n", " self.layout.scene.xaxis.title.text = labels[0]\n", " self.layout.scene.yaxis.title.text = labels[1]\n", " self.layout.scene.zaxis.title.text = labels[2]\n", " else:\n", " raise ValueError(\"You must specify labels for either 2 or 3 axes.\")\n", "\n", " def axes_ranges(self, *ranges, scale=None):\n", " if len(ranges) == 2:\n", " (xmin, xmax), (ymin, ymax) = ranges\n", " self.layout.xaxis.range = (xmin, xmax)\n", " self.layout.yaxis.range = (ymin, ymax)\n", " if scale is not None:\n", " x, y = scale\n", " self.layout.xaxis.constrain = \"domain\"\n", " self.layout.yaxis.constrain = \"domain\"\n", " self.layout.yaxis.scaleanchor = \"x\"\n", " self.layout.yaxis.scaleratio = y / x\n", " elif len(ranges) == 3:\n", " (xmin, xmax), (ymin, ymax), (zmin, zmax) = ranges\n", " self.layout.scene.xaxis.range = (xmin, xmax)\n", " self.layout.scene.yaxis.range = (ymin, ymax)\n", " self.layout.scene.zaxis.range = (zmin, zmax)\n", " if isinstance(scale, str):\n", " self.layout.scene.aspectmode = scale\n", " elif scale is not None:\n", " x, y, z = scale\n", " x *= xmax - xmin\n", " y *= ymax - ymin\n", " z *= zmax - zmin\n", " c = sorted((x, y, z))[1]\n", " self.layout.scene.aspectmode = \"manual\"\n", " self.layout.scene.aspectratio.update(x=x/c, y=y/c, z=z/c)\n", " else:\n", " raise ValueError(\"You must specify ranges for either 2 or 3 axes.\")\n", "\n", "\n", "class Figure(plotly.graph_objects.Figure, MyFigure):\n", " def __init__(self, *args, **kwargs):\n", " specs = kwargs.pop(\"subplots\", None)\n", " if specs is None:\n", " super().__init__(*args, **kwargs)\n", " else:\n", " try:\n", " specs[0][0]\n", " except:\n", " specs = [specs]\n", " rows = len(specs)\n", " cols = len(specs[0])\n", " fig = plotly.subplots.make_subplots(rows=rows, cols=cols, \n", " specs=specs, **kwargs)\n", " super().__init__(fig)\n", "\n", "\n", "class FigureWidget(plotly.graph_objects.FigureWidget, MyFigure):\n", " def __init__(self, *args, **kwargs):\n", " self._auto_items = []\n", " specs = kwargs.pop(\"subplots\", None)\n", " if specs is None:\n", " super().__init__(*args, **kwargs)\n", " else:\n", " try:\n", " specs[0][0]\n", " except:\n", " specs = [specs]\n", " rows = len(specs)\n", " cols = len(specs[0])\n", " fig = plotly.subplots.make_subplots(rows=rows, cols=cols, \n", " specs=specs, **kwargs)\n", " super().__init__(fig)\n", "\n", " def auto_update(self, *items):\n", " if self._auto_items:\n", " for item, newitem in zip(self._auto_items, items):\n", " if isinstance(newitem, tuple):\n", " newitem = newitem[0]\n", " item.update(newitem)\n", " else:\n", " for item in items:\n", " if isinstance(item, tuple):\n", " item, subplot = item\n", " self._auto_items.append(self.add(item, subplot=subplot))\n", " else:\n", " self._auto_items.append(self.add(item))\n", "\n", " def initialized(self):\n", " return len(self._auto_items) > 0\n", "\n", "\n", "def plotly_points3d(points, **options):\n", " options = options.copy()\n", " options.setdefault(\"marker_color\", options.pop(\"color\", \"black\"))\n", " options.setdefault(\"marker_size\", options.pop(\"size\", 2.5))\n", " options.setdefault(\"mode\", \"markers\")\n", " x, y, z = np.array(points, dtype=float).transpose()\n", " if options.pop(\"update\", False):\n", " return dict(x=x, y=y, z=z)\n", " return plotly.graph_objects.Scatter3d(x=x, y=y, z=z, **options)\n", "\n", "\n", "def plotly_function3d(f, x_range, y_range, **options):\n", " options = options.copy()\n", " plotpoints = options.pop(\"plotpoints\", 81)\n", " try:\n", " plotpointsx, plotpointsy = plotpoints\n", " except:\n", " plotpointsx = plotpointsy = plotpoints\n", " color = options.pop(\"color\", \"lightblue\")\n", " options.setdefault(\"colorscale\", (color, color))\n", " x, xmin, xmax = x_range\n", " y, ymin, ymax = y_range\n", " f = fast_float(f, x, y)\n", " x = np.linspace(float(xmin), float(xmax), plotpointsx)\n", " y = np.linspace(float(ymin), float(ymax), plotpointsy)\n", " x, y = np.meshgrid(x, y)\n", " xy = np.array([x.flatten(), y.flatten()]).transpose()\n", " z = np.array([f(x0, y0) for x0, y0 in xy]).reshape(plotpointsy, plotpointsx)\n", " if options.pop(\"update\", False):\n", " return dict(x=x, y=y, z=z)\n", " return plotly.graph_objects.Surface(x=x, y=y, z=z, **options)\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ ], "source": [ "x, y, z = var(\"x, y, z\")\n", "bump(x0, y0, h) = h^2/(h^2 + (x - x0)^2 + (y - y0)^2)\n", "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)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1.025749249792184, 4.924388410447029)\n", "(0.9850627229439552, 6.988915199423456)\n", "(2.5888506099044983, 1.966106259573833)\n", "(4.223844149137523, 2.0532217928719096)\n", "(4.278189612110433, 5.5133658104706305)\n", "(6.540423395003371, 2.883242346820021)\n" ] } ], "source": [ "crit_points = find_zeros(f.gradient(), (x, 0, 8), (y, 0, 8))\n", "for pt in crit_points:\n", " print(pt)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ ], "source": [ "hessian = jacobian(f.gradient(), (x, y))\n", "critpt_dict = {\"Maxima\": [], \"Minima\": [], \"Other\": []}\n", "for (x0, y0) in crit_points:\n", " e1, e2 = hessian(x0, y0).eigenvalues()\n", " if e1 < 0 and e2 < 0:\n", " critpt_dict[\"Maxima\"].append((x0, y0, f(x0, y0)))\n", " elif e1 > 0 and e2 > 0:\n", " critpt_dict[\"Minima\"].append((x0, y0, f(x0, y0)))\n", " else:\n", " critpt_dict[\"Other\"].append((x0, y0, f(x0, y0)))\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ ], "source": [ "darkorange = colors[\"darkorange\"].darker().html_color()\n", "def show_critpts_interactive(f, x_range, y_range, z_range, critpt_dict, **options):\n", " x, xmin, xmax = x_range\n", " y, ymin, ymax = y_range\n", " z, zmin, zmax = z_range\n", " axes_labels = options.get(\"axes_labels\", (str(x), str(y), str(z)))\n", " scale = options.get(\"scale\", None)\n", " colors = {\"Maxima\": \"darkgreen\", \"Minima\": \"darkred\", \"Other\": darkorange}\n", "\n", " figure = FigureWidget()\n", " figure.axes_labels(*axes_labels)\n", " figure.axes_ranges((xmin, xmax), (ymin, ymax), (zmin, zmax), scale=scale)\n", " figure.layout.margin = dict(t=10, b=10, r=10, l=10)\n", " figure.layout.showlegend = False\n", "\n", " figure.add(plotly_function3d(f, (x, xmin, xmax), (y, ymin, ymax)))\n", " points = {}\n", " planes = {}\n", " for kind, critpts in critpt_dict.items():\n", " points[kind] = figure.add(plotly_points3d(critpts, color=colors[kind]))\n", " planes[kind] = []\n", " for (x0, y0, z0) in critpts:\n", " newplane = plotly_function3d(z0, (x, x0-3, x0+3), (y, y0-3, y0+3), \n", " color=colors[kind], plotpoints=2, opacity=0.3)\n", " planes[kind].append(figure.add(newplane))\n", "\n", " choices = (\"None\",) + tuple(critpt_dict.keys()) + (\"All\",)\n", " @interact(which=selector(choices, label=\"Show: \"), \n", " show_tangent=checkbox(default=False, label=\"Show tangent plane(s)\"))\n", " def update(which, show_tangent):\n", " with figure.batch_update():\n", " for kind, critpts in points.items():\n", " critpts.visible = (which == kind or which == \"All\")\n", " for plane in planes[kind]:\n", " plane.visible = show_tangent and (which == kind or which == \"All\")\n", "\n", " return figure\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "# Learning goals:\n", "\n", "- Be able to define critical points for a function of two (or more) variables $f\\!: \\mathbb{R}^n \\to \\mathbb{R}$. \n", "- Be able to find all of the critical points of a function $f\\!: \\mathbb{R}^n \\to \\mathbb{R}$. \n", "- Be able to identify the ***three*** main types of critical points that can occur for a function $f\\!: \\mathbb{R}^n \\to \\mathbb{R}$. " ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "82f55230b8104e6da1345ea59541af2c", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Interactive function .update at 0x7f7ba52eb280> with 2 widgets\n", " wh…" ] }, "execution_count": 10, "metadata": { }, "output_type": "execute_result" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f5b931a213c5479f9a3bdee120c4638c", "version_major": 2, "version_minor": 0 }, "text/plain": [ "FigureWidget({\n", " 'data': [{'colorscale': [[0.0, 'lightblue'], [1.0, 'lightblue']],\n", " 'type': 'su…" ] }, "execution_count": 10, "metadata": { }, "output_type": "execute_result" } ], "source": [ "show_critpts_interactive(f, (x, 0, 8), (y, 0, 8), (z, 0, 6), critpt_dict, scale=(1,1,1))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "## Critical points of a function $f\\!: \\mathbb{R}^n \\to \\mathbb{R}$\n", "
\n", "\n", "**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. \n", "\n", "Note: We will only focus on the kind of critical points at which all of the partial derivatives are $0$. \n", "\n", "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.*** \n", "\n", "(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.) \n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "## Conclusions\n", "
\n", "\n", "For a function $f\\!:\\mathbb{R}^n \\to \\mathbb{R}$... \n", "\n", "- 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. \n", "\n", "- 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. \n", "\n", "- 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. \n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ ], "source": [ ] } ], "metadata": { "celltoolbar": "Slideshow", "hide_input": false, "kernelspec": { "display_name": "SageMath 9.3", "language": "sagemath", "metadata": { "cocalc": { "description": "Open-source mathematical software system", "priority": 10, "url": "https://www.sagemath.org/" } }, "name": "sage-9.3", "resource_dir": "/ext/jupyter/kernels/sage-9.3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }