{ "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 Button, Text, HTMLMath, HBox, VBox, 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", "# Below are all the actual plotting methods. First, the 2D graphics:\n", "def plotly_text(text, location, **options):\n", " options = options.copy()\n", " x, y = np.array(location, dtype=float)\n", " if options.pop(\"update\", False):\n", " return dict(text=text, x=x, y=y)\n", " options.setdefault(\"font_color\", options.pop(\"color\", \"black\"))\n", " size = options.pop(\"size\", None)\n", " if size is not None:\n", " options.setdefault(\"font_size\", size)\n", " arrow = options.pop(\"arrow\", None)\n", " if arrow:\n", " options.setdefault(\"ax\", float(arrow[0]))\n", " options.setdefault(\"ay\", float(arrow[1]))\n", " options.setdefault(\"showarrow\", True)\n", " else:\n", " options.setdefault(\"showarrow\", False)\n", " if options.pop(\"paper\", False):\n", " options.update(xref=\"paper\", yref=\"paper\")\n", " options.setdefault(\"xanchor\", \"left\")\n", " options.setdefault(\"yanchor\", \"bottom\")\n", " return plotly.graph_objects.layout.Annotation(\n", " text=text, x=x, y=y, **options)\n", "\n", "\n", "def plotly_points(points, **options):\n", " options = options.copy()\n", " x, y = np.array(points, dtype=float).transpose()\n", " if options.pop(\"update\", False):\n", " return dict(x=x, y=y)\n", " options.setdefault(\"marker_color\", options.pop(\"color\", \"black\"))\n", " options.setdefault(\"marker_size\", options.pop(\"size\", 8))\n", " options.setdefault(\"mode\", \"markers\")\n", " return plotly.graph_objects.Scatter(x=x, y=y, **options)\n", "\n", "\n", "def plotly_lines(points, **options):\n", " options = options.copy()\n", " x, y = np.array(points, dtype=float).transpose()\n", " if options.pop(\"update\", False):\n", " return dict(x=x, y=y)\n", " color = options.pop(\"color\", \"blue\")\n", " options.setdefault(\"line_color\", color)\n", " options.setdefault(\"mode\", \"lines\")\n", " return plotly.graph_objects.Scatter(x=x, y=y, **options)\n", "\n", "\n", "def plotly_vector(vec, start=(0, 0), axes_scale=(1, 1), **options):\n", " options = options.copy()\n", " options[\"mode\"] = \"lines\"\n", " options[\"fill\"] = \"toself\"\n", " color = options.pop(\"color\", \"black\")\n", " options.setdefault(\"line_color\", color)\n", " options.setdefault(\"fillcolor\", color)\n", " options.setdefault(\"line_width\", 1)\n", "\n", " tipfactor = 0.2 # Arrowhead length, as a fraction of the total arrow length\n", " tipslope1 = 2.5 # For an arrow pointing up, the slope of the leading edge of the left side of the arrowhead\n", " tipslope2 = 1.0 # For an arrow pointing up, the slope of the trailing edge of the left side of the arrowhead\n", "\n", " vec = np.array(vec, dtype=float)\n", " start = np.array(start, dtype=float)\n", " end = start + vec\n", " NaN = np.array((np.nan, np.nan), dtype=float)\n", " axes_scale = axes_scale[1] / axes_scale[0]\n", " tip = tipfactor * vec\n", " tip_perp = np.array(( -axes_scale/tipslope1 * tip[1], \n", " 1/axes_scale/tipslope1 * tip[0]), dtype=float)\n", " tipleft = end - tip + tip_perp\n", " tipmiddle = end - float(1 - tipslope2/tipslope1) * tip\n", " tipright = end - tip - tip_perp\n", " xy = np.array((start, end, NaN, end, tipleft, tipmiddle, tipright, end))\n", " x, y = xy.transpose()\n", " if options.pop(\"update\", False):\n", " return dict(x=x, y=y)\n", " return plotly.graph_objects.Scatter(x=x, y=y, **options)\n", "\n", "\n", "def plotly_vector_field(f, x_range, y_range, **options):\n", " options = options.copy()\n", " axes_scale = options.pop(\"axes_scale\", None)\n", " axes_aspect = options.pop(\"axes_aspect\", (4, 3))\n", " scalefactor = options.pop(\"scalefactor\", 1)\n", " region = options.pop(\"region\", None)\n", " plotpoints = options.pop(\"plotpoints\", 21)\n", " try:\n", " plotpointsx, plotpointsy = plotpoints\n", " except:\n", " plotpointsx = plotpointsy = plotpoints\n", " options[\"mode\"] = \"lines\"\n", " options[\"fill\"] = \"toself\"\n", " color = options.pop(\"color\", \"limegreen\")\n", " options.setdefault(\"line_color\", color)\n", " options.setdefault(\"fillcolor\", color)\n", " options.setdefault(\"line_width\", 1)\n", "\n", " # Arrow size/shape parameters\n", " zero_threshold = 1e-12\n", " minlength = 0.06 # Tip size won't be smaller than for an arrow of this length, as a fraction of the graph diagonal\n", " maxlength = 0.16 # Tip size won't be larger than for an arrow of this length, as a fraction of the graph diagonal\n", " tipfactor = 0.2 # Arrowhead length, as a fraction of the total arrow length\n", " tipslope1 = 2.5 # For an arrow pointing up, the slope of the leading edge of the left side of the arrowhead\n", " tipslope2 = 1.0 # For an arrow pointing up, the slope of the trailing edge of the left side of the arrowhead\n", "\n", " # Set the axes ranges, and the axes_scale and axes_aspect\n", " x, xmin, xmax = x_range\n", " y, ymin, ymax = y_range\n", " if axes_scale is None:\n", " axes_scale = axes_aspect[1] / axes_aspect[0] * (xmax - xmin) / (ymax - ymin)\n", " else:\n", " axes_scale = axes_scale[1] / axes_scale[0]\n", " axes_aspect = np.array((xmax - xmin, axes_scale * (ymax - ymin)), dtype=float)\n", " diagonal = np.linalg.norm(axes_aspect)\n", "\n", " # Choose the start and step size, for both x and y, for the grid points\n", " def grid_params(min1, max1, min2, max2, ratio):\n", " step1 = (max1 - min1) / (plotpoints + 0.75)\n", " start1 = min1 + 0.375 * step1\n", " step2 = step1 / ratio\n", " remainder = divmod(float((max2 - min2) / step2), int(1))[1] / 2\n", " if remainder < 0.125:\n", " remainder += 0.5\n", " start2 = min2 + remainder * step2\n", " return start1, step1, start2, step2\n", " try:\n", " plotpoints[1]\n", " except: # A single value was given for plotpoints. Set grid automatically. \n", " if axes_aspect[0] > axes_aspect[1]:\n", " x0, xstep, y0, ystep = grid_params(xmin, xmax, ymin, ymax, axes_scale)\n", " else:\n", " y0, ystep, x0, xstep = grid_params(ymin, ymax, xmin, xmax, 1/axes_scale)\n", " else: # A 2-tuple was given for plotpoints. Use exactly that many. \n", " xstep = (xmax - xmin) / (plotpoints[0] + 0.75)\n", " x0 = xmin + 0.375 * xstep\n", " ystep = (ymax - ymin) / (plotpoints[1] + 0.75)\n", " y0 = ymin + 0.375 * ystep\n", "\n", " # The actual computations\n", " f1, f2 = [fast_float(f_, x, y) for f_ in f]\n", " xvals = np.arange(float(x0), float(xmax), float(xstep))\n", " yvals = np.arange(float(y0), float(ymax), float(ystep))\n", " xy = np.array(np.meshgrid(xvals, yvals)).reshape(2, -1).transpose()\n", " if region:\n", " region = fast_float(region, x, y)\n", " xy = xy[[region(x0, y0) > 0 for x0, y0 in xy]]\n", " vec = np.array([(f1(x0, y0), f2(x0, y0)) for x0, y0 in xy]).transpose()\n", " keep = np.linalg.norm(vec, axis=0) > zero_threshold\n", " xy = xy[keep].transpose()\n", " vec = vec[:,keep]\n", " grid_scale = np.array((xstep, ystep), dtype=float).reshape(2, 1)\n", " longest_vector = np.linalg.norm(vec / grid_scale, axis=0, ord=np.inf).max()\n", " vec *= scalefactor / longest_vector\n", " paper_scale = np.array((1, axes_scale), dtype=float).reshape(2, 1)\n", " length = np.linalg.norm(vec * paper_scale, axis=0) / diagonal\n", " tip = (tipfactor * np.clip(length, minlength, maxlength) / length) * vec\n", " tip_perp = np.array((float( -axes_scale/tipslope1) * tip[1], \n", " float(1/axes_scale/tipslope1) * tip[0]))\n", " NaN = np.full_like(xy, np.nan)\n", " end = xy + vec\n", " tipleft = end - tip + tip_perp\n", " tipmiddle = end - float(1 - tipslope2/tipslope1) * tip\n", " tipright = end - tip - tip_perp\n", " xy = np.array((xy, end, NaN, end, tipleft, tipmiddle, tipright, end, NaN))\n", " x, y = np.moveaxis(xy, 0, 2).reshape(2, -1)\n", " if options.pop(\"update\", False):\n", " return dict(x=x, y=y)\n", " return plotly.graph_objects.Scatter(x=x, y=y, **options)\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_lines3d(points, **options):\n", " options = options.copy()\n", " color = options.pop(\"color\", \"blue\")\n", " options.setdefault(\"line_color\", color)\n", " options.setdefault(\"mode\", \"lines\")\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", "\n", "def plotly_vector_field_on_surface3d(f, field, u_range, v_range, **options):\n", " options = options.copy()\n", " plotpoints = options.pop(\"plotpoints\", 21)\n", " try:\n", " plotpointsu, plotpointsv = plotpoints\n", " except:\n", " plotpointsu = plotpointsv = plotpoints\n", " color = options.pop(\"color\", \"limegreen\")\n", " options.setdefault(\"colorscale\", (color, color))\n", " u, umin, umax = u_range\n", " v, vmin, vmax = v_range\n", " f1, f2, f3 = [fast_float(f_, u, v) for f_ in f]\n", " if len(field) == 2:\n", " field = f.derivative() * field\n", " fx, fy, fz = [fast_float(f_, u, v) for f_ in field]\n", " u = np.linspace(float(umin), float(umax), plotpointsu)\n", " v = np.linspace(float(vmin), float(vmax), plotpointsv)\n", " uv = np.array(np.meshgrid(u, v)).reshape(2, -1).transpose()\n", " points = np.array([(f1(u0, v0), f2(u0, v0), f3(u0, v0)) for u0, v0 in uv])\n", " vectors = np.array([(fx(u0, v0), fy(u0, v0), fz(u0, v0)) for u0, v0 in uv])\n", " x, y, z = points.transpose()\n", " u, v, w = vectors.transpose()\n", " if options.pop(\"update\", False):\n", " return dict(x=x, y=y, z=z, u=u, v=v, w=w)\n", " return plotly.graph_objects.Cone(x=x, y=y, z=z, u=u, v=v, w=w, **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\": [], \"Saddle points\": []}\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[\"Saddle points\"].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 gradient3d_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\", \"Saddle points\": 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", " surface = figure.add(plotly_function3d(f, x_range, y_range))\n", " points = {}\n", " xypoints = {}\n", " xylines = {}\n", " for kind, critpts in critpt_dict.items():\n", " points[kind] = figure.add(plotly_points3d(critpts, color=colors[kind]))\n", " xypts = [(x0, y0, 0) for (x0, y0, z0) in critpts]\n", " lines = [(a, b, (np.nan, np.nan, np.nan)) for a, b in zip(critpts, xypts)]\n", " lines = np.array(lines).reshape(-1, 3)\n", " xypoints[kind] = figure.add(plotly_points3d(xypts, color=colors[kind]))\n", " xylines[kind] = figure.add(plotly_lines3d(lines, color=colors[kind], line_dash=\"dash\"))\n", " xyplane(x, y) = (x, y, 0)\n", " field = figure.add(plotly_vector_field_on_surface3d(xyplane, f.gradient(), \n", " x_range, y_range, color=\"limegreen\"))\n", "\n", " choices = (\"None\",) + tuple(critpt_dict.keys()) + (\"All\",)\n", " @interact(which=selector(choices, label=\"Show: \"), \n", " show_xypoints=checkbox(default=False, label=\"Show critical points\"), \n", " show_gradient=checkbox(default=False, label=\"Show gradient vector field\"))\n", " def update(which, show_xypoints, show_gradient):\n", " with figure.batch_update():\n", " surface.opacity = 0.6 if (show_xypoints or show_gradient) else 1\n", " field.visible = show_gradient\n", " for kind in critpt_dict:\n", " visible = (which == kind or which == \"All\")\n", " points[kind].visible = visible\n", " xypoints[kind].visible = show_xypoints and visible\n", " xylines[kind].visible = show_xypoints and visible\n", "\n", " return figure\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ ], "source": [ "def gradient2d_interactive(f, x_range, y_range, critpt_dict, **options):\n", " x, xmin, xmax = x_range\n", " y, ymin, ymax = y_range\n", " axes_labels = options.get(\"axes_labels\", (f\"${x}$\", f\"${y}$\"))\n", " scale = options.get(\"scale\", None)\n", " zoom_factor = options.get(\"zoom_factor\", 10)\n", " zoom_xrange = (xmax - xmin) / zoom_factor\n", " zoom_yrange = (ymax - ymin) / zoom_factor\n", " colors = {\"Maxima\": \"green\", \"Minima\": \"red\", \"Saddle points\": \"darkorange\"}\n", "\n", " figure = FigureWidget()\n", " figure.axes_labels(*axes_labels)\n", " figure.axes_ranges((xmin, xmax), (ymin, ymax), scale=scale)\n", " figure.layout.margin = dict(t=30, b=40, r=10, l=10)\n", " figure.layout.dragmode = \"select\" # For this figure, we want select, not pan\n", "\n", " field = figure.add(plotly_vector_field(f.gradient(), x_range, y_range, \n", " color=\"limegreen\", plotpoints=17, axes_scale=scale, showlegend=False))\n", " points = []\n", " for kind, critpts in critpt_dict.items():\n", " points.append(figure.add(plotly_points([(x0, y0) for (x0, y0, z0) in critpts], \n", " name=kind, color=colors[kind])))\n", "\n", " def zoomin(widget):\n", " selected = []\n", " for critpts in points:\n", " if critpts.selectedpoints:\n", " for i in critpts.selectedpoints:\n", " selected.append((critpts.x[i], critpts.y[i]))\n", " if len(selected) != 1:\n", " return\n", " x0, y0 = selected[0]\n", " new_xrange = (x, x0 - zoom_xrange/2, x0 + zoom_xrange/2)\n", " new_yrange = (y, y0 - zoom_yrange/2, y0 + zoom_yrange/2)\n", " button.description = \"Zoom out\"\n", " button.on_click(zoomin, remove=True)\n", " button.on_click(zoomout)\n", " with figure.batch_update():\n", " figure.axes_ranges(new_xrange[1:], new_yrange[1:], scale=scale)\n", " field.update(plotly_vector_field(f.gradient(), new_xrange, new_yrange, \n", " plotpoints=11, axes_scale=scale, update=True))\n", "\n", " def zoomout(widget):\n", " button.description = \"Zoom in\"\n", " button.on_click(zoomout, remove=True)\n", " button.on_click(zoomin)\n", " with figure.batch_update():\n", " figure.axes_ranges((xmin, xmax), (ymin, ymax), scale=scale)\n", " field.update(plotly_vector_field(f.gradient(), x_range, y_range, \n", " plotpoints=17, axes_scale=scale, update=True))\n", "\n", " button = Button(description=\"Zoom in\")\n", " button.on_click(zoomin)\n", "\n", " return HBox((button, figure))\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ ], "source": [ "def vectorfield_interactive(f, x_range, y_range, **options):\n", " x, xmin, xmax = x_range\n", " y, ymin, ymax = y_range\n", " axes_labels = options.get(\"axes_labels\", (f\"${x}$\", f\"${y}$\"))\n", " axes_scale = options.get(\"axes_scale\", None)\n", " scalefactor = options.get(\"scalefactor\", None)\n", " funcname = options.get(\"label\", \"f\")\n", "\n", " figure = FigureWidget()\n", " figure.axes_labels(*axes_labels)\n", " figure.axes_ranges((xmin, xmax), (ymin, ymax), scale=axes_scale)\n", " field = figure.add(plotly_vector_field(f, x_range, y_range, name=\"Vector field\", \n", " plotpoints=11, color=\"limegreen\", axes_scale=axes_scale, visible=\"legendonly\"))\n", " point = figure.add(plotly_points([(xmin, ymin)], showlegend=False))\n", "\n", " def update_point(change):\n", " value = change[\"new\"].strip()\n", " if value[0] != \"(\" or value[-1] != \")\":\n", " return\n", " try:\n", " x0, y0 = [float(num) for num in value[1:-1].split(\",\")]\n", " except:\n", " return\n", " with figure.batch_update():\n", " point.x = [x0]\n", " point.y = [y0]\n", "\n", " def add_vector(widget):\n", " x0, y0 = point.x[0], point.y[0]\n", " vec = f(x0, y0) * scalefactor\n", " with figure.batch_update():\n", " figure.add(plotly_vector(vec, start=(x0, y0), color=\"purple\", line_width=2, \n", " axes_scale=axes_scale, showlegend=False))\n", "\n", " point_input = Text(description=\"Point:\")\n", " point_input.observe(update_point, names=\"value\")\n", " point_input.value = f\"({round(xmin)},{round(ymin)})\"\n", " button = Button(description=\"Add vector\")\n", " button.on_click(add_vector)\n", " label = fr\"{latex(f(x,y)[0])} \\\\ {latex(f(x,y)[1])}\"\n", " label = fr\"$\\qquad {funcname}({x}, {y}) = \\begin{{bmatrix}} {label} \\end{{bmatrix}}$\"\n", " label = HTMLMath(label)\n", "\n", " return VBox((HBox((point_input, button, label)), figure))\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "# Learning goals: \n", "\n", "- Be able to compute the ***gradient vector field*** of a function $f\\!: \\mathbb{R}^n \\to \\mathbb{R}$. \n", "- Be able to explain how the *direction* of the gradient vector field of $f$ relates to the “slope” of the function $f$. \n", "- 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*. \n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "d984844c631b461394a1f7fd2c909ceb", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Interactive function .update at 0x7fab0cfbca60> with 3 widgets\n", " whic…" ] }, "execution_count": 12, "metadata": { }, "output_type": "execute_result" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "9a40564d5b25446a8c685c7da517bcf3", "version_major": 2, "version_minor": 0 }, "text/plain": [ "FigureWidget({\n", " 'data': [{'colorscale': [[0.0, 'lightblue'], [1.0, 'lightblue']],\n", " 'opacity': …" ] }, "execution_count": 12, "metadata": { }, "output_type": "execute_result" } ], "source": [ "gradient3d_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": [ "## The gradient of a function $f\\!: \\mathbb{R}^n \\to \\mathbb{R}$\n", "
\n", "\n", "Consider a function $f$ of two (or more) variables, say $f(x, y)$. \n", "\n", "We know $f$ will have two partial derivatives: \n", "$$ \\frac{\\partial f}{\\partial x} \\qquad \\text{and} \\qquad \\frac{\\partial f}{\\partial y} $$\n", "\n", "\n", " 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: \n", "\n", "$$ \\begin{bmatrix} \\frac{\\partial f}{\\partial x} \\\\ \\frac{\\partial f}{\\partial y} \\end{bmatrix} $$\n", "\n", "the result will be a function from from $\\mathbb{R}^2$ to $\\mathbb{R}^2$... a vector field!\n", "\n", "This vector field is called the gradient of $f$.\n", "\n", "
\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "slideshow": { "slide_type": "skip" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(0, 0): Eigenvalues [10, 40]\n", "(40/3, 0): Eigenvalues [-1240/9, -10]\n", "(0, 20/3): Eigenvalues [-310/9, -40]\n", "(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]\n", "(-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]\n", "(-8, -4): Eigenvalues [-sqrt(4177) + 15, sqrt(4177) + 15]\n", "(5, 5/2): Eigenvalues [-35, 65/4]\n" ] } ], "source": [ "g(x, y) = 5*x^2 + 20*y^2 - 1/4*x^3 - 2*y^3 - 1/2*x^2*y^2\n", "H = g.derivative(2)\n", "for x0, y0 in solve(list(g.derivative()(x, y)), (x, y)):\n", " x0, y0 = x0.rhs(), y0.rhs()\n", " if not (x0 in RR and y0 in RR):\n", " continue\n", " print(f\"({x0}, {y0}): Eigenvalues {H(x0,y0).eigenvalues()}\")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "## An example of the gradient of a function\n", "
\n", "\n", "Define $f\\!: \\mathbb{R}^2 \\to \\mathbb{R}$ by \n", "$$ f(x, y) = 5x^2 + 20y^2 - \\tfrac{1}{4} x^3 - 2y^3 - \\tfrac{1}{2} x^2 y^2 $$\n", "\n", "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$\n", "\n", "And the partial derivative of $f$ with respect to $y$ is $\\quad \\frac{\\partial f}{\\partial y} = 40y - 6y^2 - x^2 y$\n", "\n", "
\n", "\n", "So the gradient of $f$ is a function, which we'll write as $\\mathrm{grad}f$, defined by \n", "\n", "$$ \\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} $$\n", "
\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "7ce4bf96120742aca15f6b928524c3f9", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Text(value='(4,2)', description='Point:'), Button(description='Add vector', styl…" ] }, "execution_count": 14, "metadata": { }, "output_type": "execute_result" } ], "source": [ "options = dict(scalefactor=0.005, axes_scale=(1,1), label=r\"\\textrm{grad}f\")\n", "vectorfield_interactive(g.derivative(), (x, 3.6, 6.4), (y, 1.8, 3.2), **options)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "## Notation and a few details about the gradient\n", "
\n", "\n", "**Notation:** The gradient of $f$ is written as $\\ \\mathrm{grad}f \\ $ or $\\ \\vec{\\nabla}f$. \n", "\n", "So for our previous example, we could write \n", "\n", "$$ \\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} $$\n", "\n", "
\n", "\n", "
\n", "\n", "Dimensions: If $f\\!: \\mathbb{R}^n \\to \\mathbb{R}$, then $\\mathrm{grad} f\\!: \\mathbb{R}^n \\to \\mathbb{R}^n$ is defined by \n", "\n", "$$ \\mathrm{grad}f = \\begin{bmatrix} \\frac{\\partial f}{\\partial x} \\\\ \\frac{\\partial f}{\\partial y} \\\\ \\vdots \\end{bmatrix} $$\n", "\n", "
\n", "
\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "## The gradient, and critical points of $f$\n", "
\n", "\n", "

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

\n", "\n", "

\n", "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! \n", "

\n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "## How to visualize the gradient vector field" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "83814a4455874e2f8c49bfdd8d66637a", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Interactive function .update at 0x7fab0d00b820> with 3 widgets\n", " whic…" ] }, "execution_count": 15, "metadata": { }, "output_type": "execute_result" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "8d30481e8a8a42a0bf95efefb6b98d00", "version_major": 2, "version_minor": 0 }, "text/plain": [ "FigureWidget({\n", " 'data': [{'colorscale': [[0.0, 'lightblue'], [1.0, 'lightblue']],\n", " 'opacity': …" ] }, "execution_count": 15, "metadata": { }, "output_type": "execute_result" } ], "source": [ "gradient3d_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": "-" } }, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "## How to visualize the gradient vector field" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "01300574aa65411a8106b0db5a492f93", "version_major": 2, "version_minor": 0 }, "text/plain": [ "HBox(children=(Button(description='Zoom in', style=ButtonStyle()), FigureWidget({\n", " 'data': [{'fill': 'tosel…" ] }, "execution_count": 16, "metadata": { }, "output_type": "execute_result" } ], "source": [ "gradient2d_interactive(f, (x, 0, 8), (y, 0, 8), critpt_dict, scale=(1,1))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "source": [ "## The direction of the gradient vector field\n", "
\n", "\n", "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.* \n", "\n", "In other words, in short, it points in the “uphill” direction. \n", "\n", "

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

\n", "
    \n", "
  • A local maximum of $f$ will be a stable equilibrium point (sink) of the $\\,\\mathrm{grad}f \\,$ vector field.
    • \n", "
  • A local minimum of $f$ will be an unstable equilibrium point (source) of the $\\,\\mathrm{grad}f \\,$ vector field.
    • \n", "
  • A saddle point of $f$ will be a saddle point of the $\\,\\mathrm{grad}f \\,$ vector field.
    • \n", "
\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}$, its gradient is the vector field $\\mathrm{grad}f\\!: \\mathbb{R}^n \\to \\mathbb{R}^n$ defined by \n", "$$ \\mathrm{grad}f (x, y, \\dotsc) = \\begin{bmatrix} \\frac{\\partial f}{\\partial x} \\\\ \\frac{\\partial f}{\\partial y} \\\\ \\vdots \\end{bmatrix} $$\n", "\n", "- Therefore the critical points of $f$ are the same as the equilibrium points of the gradient vector field of $f$. \n", "\n", "- 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) \n", "\n", "- Therefore we can use the gradient of $f$ to classify the critical points of $f$: \n", " - A local maximum of $f$ will be a stable equilibrium point (sink) of the gradient vector field of $f$. \n", " - A local minimum of $f$ will be an unstable equilibrium point (source) of the gradient vector field of $f$. \n", " - A saddle point of $f$ will be a saddle point of the gradient vector field of $f$. \n", "\n", "
\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 }