Notebooks we created and walked through.

Author: mikeckennedy

apps/py/ch11_notebooks/Lorenz.ipynb

@@ -0,0 +1,233 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# The Lorenz Differential Equations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Before we start, we import some preliminary libraries. We will also import (below) the accompanying `lorenz.py` file, which contains the actual solver and plotting routine."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "from ipywidgets import interactive, fixed"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We explore the Lorenz system of differential equations:\n",
+    "\n",
+    "$$\n",
+    "\\begin{aligned}\n",
+    "\\dot{x} & = \\sigma(y-x) \\\\\n",
+    "\\dot{y} & = \\rho x - y - xz \\\\\n",
+    "\\dot{z} & = -\\beta z + xy\n",
+    "\\end{aligned}\n",
+    "$$\n",
+    "\n",
+    "Let's change (\\\\(\\sigma\\\\), \\\\(\\beta\\\\), \\\\(\\rho\\\\)) with ipywidgets and examine the trajectories."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "7b36bce3b70f4e85aa5de91273312367",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=10.0, description='sigma', max=50.0), FloatSlider(value=2.666666666666…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "from lorenz import solve_lorenz\n",
+    "w=interactive(solve_lorenz,sigma=(0.0,50.0),rho=(0.0,50.0))\n",
+    "w"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the default set of parameters, we see the trajectories swirling around two points, called attractors. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The object returned by `interactive` is a `Widget` object and it has attributes that contain the current result and arguments:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t, x_t = w.result"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'sigma': 35.1, 'beta': 0.63333, 'rho': 30.4}"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "w.kwargs"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "After interacting with the system, we can take the result and perform further computations. In this case, we compute the average positions in \\\\(x\\\\), \\\\(y\\\\) and \\\\(z\\\\)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xyz_avg = x_t.mean(axis=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(30, 3)"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "xyz_avg.shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Creating histograms of the average positions (across different trajectories) show that, on average, the trajectories swirl about the attractors."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib import pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEKCAYAAAALoA6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAARrElEQVR4nO3de5BkZX3G8e8jKxcJAXVH7rCQEKJiFDIavEQNEEU0oIkmoEZQUhvLGDVllVlilVZZJt5vlIlkgyiWiije8BZFlJCkAN1F7qtyW2UR2MELIiQi+MsffZY0w85Mz3TP9L7s91PVtefy9nl/c6b3mXfe06cnVYUkqT0PGncBkqSFMcAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS4tkSRvSfKaOdp8K8mjl6omtc0A11glOS/JT5NsN+5aFlOSCeAlwL9O274hycF9m94JvGkpa1O7DHCNTZIVwB8CBRy9CMdfNupjDuEE4MtV9T+bNiRZDuwKXNXX7mzgj5LstrTlqUUGuMbpJcCFwIeB4zdtTPL3Sc7qb5jkfUlO7pb3SPLpJFNJrk/yqr5267vnXwbckWRZklVJrk1ye5Krkjyvr/0hSb7T7ftUkjOTvLlv/4x9TZfk7Uk+17f+jiTnJtkWeBbwH337fhu4gd7/wR8n+XGSZVX1v8Ba4JnzPZnaClWVDx9jeQDXAK8Afh/4FbBrt31f4E5gp259G+Am4FB6gbcWeAOwLbA/cB3wzK7teuASYG9gh27bC4A9uuf+BXAHsHv3/B8ArwYeDPwpcBfw5u55s/a1ma/n4cBtwMHAy4HLgZ27fVPA46e1/xvgzM0c52Tg3eP+/vjY8h+OwDUWSZ5CL6g/WVVrgWuBFwJU1Q+Ai4FNI+XDgDur6kLg8cBEVb2pqu6qquuAfwOO7Tv8yVV1Q3XTFVX1qar6UVX9uqrOBK4GnkDvB8Kyrv2vquozwLf6jjNIX/eqqh8D7wFOB04Cjqqq27rduwC3T3vKY+n9sJnu9q69NCsDXONyPPC1qrq1W/84fdMo3fpx3fILu3Xohf4eSX626QH8A7255E1u6O8oyUuSXNLX/iBgOb1R+Y1VVTM8d5C+pvsO8BjgpKrqP9ZPgZ2mtX0ccOlmjrET8LNZ+pCA3uhDWlJJdgD+HNgmyc3d5u2AXZI8tqouBT4FvCvJXvRG4k/s2t0AXF9VB8zSxb2BnGRfeqPmw4ELquqeJJcAoTcts2eS9IX43vR+Gxi0r/6v6zHAB+iNwF/G///QAbgM+B3g213bB9H7QbK5EfgjgY8O0qe2bo7ANQ7PBe4BHkVvFPo4eqH1n/QubFJVU8B5wIfohei67rnfAm7vLlTukGSbJAclefwMfe1IL9CnAJK8lF5wAlzQ1fHK7mLnMfSmVjYZuK8kewJfoDf3/QrgMUme3tfky8DT+tZ36B73+T+YZHt61wTOmeHrke5lgGscjgc+VFU/rKqbNz2A9wMv6nv738eBI+gbyVbVPcBz6IX+9cCtwKnAzpvrqKquAt5FL6xvoTe98d/dvrvoXbg8kd6UxYuBLwK/nE9fSX6TXkC/u6rOrqo7gXcA/9jX7CPAUd1vH1TVHcApwFVJNvS1+xPgvKr60RznUCL3nf6Ttm5JLgJOqaoPLcKx/wnYWFXvnaP/E6vqilH3rwceA1xbtSRPA75Hb3T9Inqj4v2r6qaxFiYNwIuY2todCHyS3lz5dcDzDW+1whG4JDXKi5iS1KglnUJZvnx5rVixYim7lKTmrV279taqmpi+fUkDfMWKFaxZs2Ypu5Sk5iX5wea2O4UiSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGjVngCc5LcnGJPf7cJ0kr01S3R9nlSQtoUFG4B8Gjpy+McnewDOAH464JknSAOYM8Ko6H/jJZna9B3gdfX/9RJK0dBZ0J2b3l0turKpLk8zVdiWwEmCfffZZSHfaiqxY9aWx9Lv+rc8eS7/SMOZ9ETPJQ+j9Ydc3DNK+qlZX1WRVTU5M3O9WfknSAi3kXSi/BewHXJpkPbAXcHGS3UZZmCRpdvOeQqmqy4FHbFrvQnyyqm4dYV2SpDkM8jbCM+j9QdgDk2xIcuLilyVJmsucI/CqOm6O/StGVo0kaWDeiSlJjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYN8lfpT0uyMckVfdvekeS7SS5L8tkkuyxumZKk6QYZgX8YOHLatnOAg6rq94DvAyeNuC5J0hzmDPCqOh/4ybRtX6uqu7vVC4G9FqE2SdIsRjEH/jLgKzPtTLIyyZoka6ampkbQnSQJhgzwJK8H7gY+NlObqlpdVZNVNTkxMTFMd5KkPssW+sQkJwDPAQ6vqhpZRZKkgSwowJMcCbwOeFpV3TnakiRJgxjkbYRnABcABybZkORE4P3ATsA5SS5Jcsoi1ylJmmbOEXhVHbeZzR9chFokSfPgnZiS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjZozwJOclmRjkiv6tj0syTlJru7+fejililJmm6QEfiHgSOnbVsFnFtVBwDnduuSpCU0Z4BX1fnAT6ZtPgY4vVs+HXjuiOuSJM1hoXPgu1bVTd3yzcCuMzVMsjLJmiRrpqamFtidJGm6oS9iVlUBNcv+1VU1WVWTExMTw3YnSeosNMBvSbI7QPfvxtGVJEkaxEID/Gzg+G75eODzoylHkjSoQd5GeAZwAXBgkg1JTgTeCvxxkquBI7p1SdISWjZXg6o6boZdh4+4FknSPHgnpiQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGjVUgCf5uyRXJrkiyRlJth9VYZKk2S04wJPsCbwKmKyqg4BtgGNHVZgkaXbDTqEsA3ZIsgx4CPCj4UuSJA1iwQFeVTcC7wR+CNwE3FZVX5veLsnKJGuSrJmamlp4pZKk+xhmCuWhwDHAfsAewI5JXjy9XVWtrqrJqpqcmJhYeKWSpPsYZgrlCOD6qpqqql8BnwGeNJqyJElzGSbAfwgcmuQhSQIcDqwbTVmSpLkMMwd+EXAWcDFweXes1SOqS5I0h2XDPLmq3gi8cUS1SJLmwTsxJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUUPdibmUVqz60tj6Xv/WZ4+tby0NX19qkSNwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUUMFeJJdkpyV5LtJ1iV54qgKkyTNbtgPs3of8O9V9fwk2wIPGUFNkqQBLDjAk+wMPBU4AaCq7gLuGk1ZkqS5DDOFsh8wBXwoyXeSnJpkx+mNkqxMsibJmqmpqSG6kyT1GybAlwGHAB+oqoOBO4BV0xtV1eqqmqyqyYmJiSG6kyT1GybANwAbquqibv0seoEuSVoCCw7wqroZuCHJgd2mw4GrRlKVJGlOw74L5W+Bj3XvQLkOeOnwJUmSBjFUgFfVJcDkiGqRJM2Dd2JKUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRQwd4km2SfCfJF0dRkCRpMKMYgb8aWDeC40iS5mGoAE+yF/Bs4NTRlCNJGtSwI/D3Aq8Dfj2CWiRJ87DgAE/yHGBjVa2do93KJGuSrJmamlpod5KkaYYZgT8ZODrJeuATwGFJPjq9UVWtrqrJqpqcmJgYojtJUr8FB3hVnVRVe1XVCuBY4BtV9eKRVSZJmpXvA5ekRi0bxUGq6jzgvFEcS5I0GEfgktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY1acIAn2TvJN5NcleTKJK8eZWGSpNktG+K5dwOvraqLk+wErE1yTlVdNaLaJEmzWPAIvKpuqqqLu+XbgXXAnqMqTJI0u2FG4PdKsgI4GLhoM/tWAisB9tlnn1F0J0kLsmLVl8bW9/q3Pnvkxxz6ImaS3wA+Dbymqn4+fX9Vra6qyaqanJiYGLY7SVJnqABP8mB64f2xqvrMaEqSJA1imHehBPggsK6q3j26kiRJgxhmBP5k4C+Bw5Jc0j2OGlFdkqQ5LPgiZlX9F5AR1iJJmgfvxJSkRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElq1FABnuTIJN9Lck2SVaMqSpI0twUHeJJtgH8GngU8CjguyaNGVZgkaXbDjMCfAFxTVddV1V3AJ4BjRlOWJGkuy4Z47p7ADX3rG4A/mN4oyUpgZbf6iyTfG6LPYS0Hbp3vk/K2Rahk/hZU+xbC2mexiK8vz/vSm7HuIb/P+25u4zABPpCqWg2sXux+BpFkTVVNjruOhbD28bD28Wi19qWue5gplBuBvfvW9+q2SZKWwDAB/m3ggCT7JdkWOBY4ezRlSZLmsuAplKq6O8krga8C2wCnVdWVI6tscWwRUzkLZO3jYe3j0WrtS1p3qmop+5MkjYh3YkpSowxwSWrUAzrAk7wgyZVJfp1kxrf2JFmf5PIklyRZs5Q1zmQetW9xH2eQ5GFJzklydffvQ2dod093zi9JMtYL4HOdxyTbJTmz239RkhVLX+X9DVD3CUmm+s7zX42jzs1JclqSjUmumGF/kpzcfW2XJTlkqWucyQC1Pz3JbX3n/Q2LUkhVPWAfwCOBA4HzgMlZ2q0Hlo+73vnWTu/i8bXA/sC2wKXAo7aA2t8OrOqWVwFvm6HdL8Zd66DnEXgFcEq3fCxwZiN1nwC8f9y1zlD/U4FDgCtm2H8U8BUgwKHAReOueR61Px344mLX8YAegVfVuqoa552fCzZg7VvqxxkcA5zeLZ8OPHeMtQxikPPY/zWdBRyeJEtY4+Zsqd//gVTV+cBPZmlyDPCR6rkQ2CXJ7ktT3ewGqH1JPKADfB4K+FqStd2t/63Y3McZ7DmmWvrtWlU3dcs3A7vO0G77JGuSXJhknCE/yHm8t01V3Q3cBjx8Saqb2aDf/z/rpiDOSrL3ZvZvqbbU1/egnpjk0iRfSfLoxehg0W+lX2xJvg7stpldr6+qzw94mKdU1Y1JHgGck+S73U/YRTWi2sdittr7V6qqksz0XtV9u/O+P/CNJJdX1bWjrnUr9wXgjKr6ZZK/pvdbxGFjrmlrcDG91/cvkhwFfA44YNSdNB/gVXXECI5xY/fvxiSfpfer6aIH+AhqH9vHGcxWe5JbkuxeVTd1v/JunOEYm877dUnOAw6mN6e71AY5j5vabEiyDNgZ+PHSlDejOeuuqv4aT6V3faIVzX5cR1X9vG/5y0n+JcnyqhrpB3Rt9VMoSXZMstOmZeAZwGavLG+BttSPMzgbOL5bPh64328TSR6aZLtueTnwZOCqJavwvgY5j/1f0/OBb1R3tWqM5qx72pzx0cC6JaxvWGcDL+nejXIocFvf1NwWLclum66RJHkCvawd/Q/8cV/NXcwH8Dx682a/BG4Bvtpt3wP4cre8P72r95cCV9Kbvmii9m79KOD79EauW0rtDwfOBa4Gvg48rNs+CZzaLT8JuLw775cDJ4655vudR+BNwNHd8vbAp4BrgG8B+4/7PA9Y91u61/WlwDeB3x13zX21nwHcBPyqe62fCLwceHm3P/T+aMy13WtkxneSbYG1v7LvvF8IPGkx6vBWeklq1FY/hSJJrTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqP+D67Oa2aj0uSKAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.hist(xyz_avg[:,0])\n",
+    "plt.title('Average $x(t)$');"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEKCAYAAAALoA6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAShUlEQVR4nO3de7SldV3H8fcnJgSVJeKcuI8jpaRSJp3MS3nDlMBEy1pQBhSt6bLUbFU21ipXrkoyK3VZ2USILhUvpEVeSlKJanFxQEYug3JxlMHBOUgZWongtz/2M7g5nnP2Pnvvs/f8nPdrrb3Oc/nt5/edh83nPOf3PM9+UlVIktrzbbMuQJI0GgNckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeDSlCR5dZKXDWhzRZLHTqsmtc0A10wluTjJfyZ5wKxrWUtJ5oDTgb9etHxnksf3LXot8Kpp1qZ2GeCamSQbgR8GCnjeGmx/3aS3OYYzgQ9W1f/uWZBkPXAocH1fuwuBZyQ5bLrlqUUGuGbpdOAy4DzgjD0Lk/xWkgv6GyZ5fZI3dNNHJPm7JAtJPpPkpX3tdnTv/yTwlSTrkmxOcnOSu5Jcn+QFfe2PT/KJbt17krwryR/0rV+2r0X1PTjJvUkO71t2XJJdSQ4CfhT417513wXcSu//wS8m+WKSdVX1f8CVwHNG2aHatxjgmqXTgbd3r+ckObRb/k7gpC74SLIf8FPAO5J8G/CPwDbgSOAE4GVJ+gPvNOBk4OCquge4md6R/kOA3wfeluTwJPsD76P3C+QQ4HygP9yH6QuAqvoycANwfN/is4E/qqq7gO8BPtXX/ibgN4ALqurBVfWwrlaA7cDjhtqD2qcZ4JqJJD8EPBx4d1VdSS9kfxqgqj4LXMU3wvSZwP9U1WXADwBzVfWqqrq7qm4B/gY4tW/zb6iqW/cMV1TVe6rq81X19ap6F3Aj8ATgicC6rv3Xquq9wBV92xmmr34fpwvwJE8FHsM3xrwPBu5a1P5xwNVLbOeurr20IgNcs3IG8OGquqObfwd9wyjd/Gnd9E9389AL/SOS/NeeF/Db9MaS97i1v6Mkpye5uq/9ccB64Ajgtrr/dyr3v3eYvvrdF+DAa4Dfraq7u/n/BA5a1P776B3dL3YQ8F/L9CHdZ286yaN9RJID6Q2J7Jfk9m7xA4CDkzyuqrYB7wH+NMlR9I7En9S1uxX4TFU9coUu7gvkJA+nd9R8AnBpVd2b5GogwC7gyCTpC/Gj6f01MGxf/T4OvDzJTwAH8I1fOgCfBB7VtdkzPHMcSx+BPxp425B9ah/mEbhm4fnAvfSGGL6vez0a+Dd64+JU1QJwMfBmeiG6vXvvFcBd3YnKA5Ps150s/IFl+noQvUBfAEjyc/SCE+DSro4Xdyc7T6E3tLLHavvaBhwG/CnwikVH9h8EntY3f2D3ut//g0kOAL4fuGiZPqT7GOCahTOAN1fV56rq9j0v4I3Az/Rd/vcO4Fn0HclW1b3Ac+mF/meAO4Bz6J2g/CZVdT29QL0U+AK9k4n/0a27G/hx4Cx6QxYvAt4PfHXEvr4KXAPsqKoPLVr9VnonZg/s2n4FeBNwfZKdfe1+DLi4qj6/5J6T+sRHqknfkORy4E1V9eYR3rs/cBPwU90J18Xr/wjYXVWvG9D/WVV17Wr7177HANc+LcnT6F3edwfwM/SOio+pql0jbOsPu/eeNrCxNAGexNS+7ljg3fTGym8BXrja8E5yPPAxeicqXzCguTQxHoFLUqM8iSlJjZrqEMr69etr48aN0+xSkpp35ZVX3lFVc4uXTzXAN27cyNatW6fZpSQ1L8lnl1ruEIokNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElq1MAAT3Jukt1Jrl20/CVJbkhyXZLXrF2JkqSlDHMEfh5wYv+CJM8ATgEeV1WPBV47+dIkSSsZGOBVdQlw56LFvwyc3X3/MVW1ew1qkyStYNQ7MR8F/HD39Zn/B/xGVX18qYZJNgGbADZs2DBid9pXbNz8gZn0u+Psk2fSrzSOUU9irgMOofdU798E3p0kSzWsqi1VNV9V83Nz33QrvyRpRKMG+E7gvdVzBfB1ek/5liRNyagB/vfAMwCSPArYn94TTSRJUzJwDDzJ+cDTgfXdw1dfCZwLnNtdWng3cEb5ZAhJmqqBAb7C8/1eNOFaJEmr4J2YktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGDQzwJOcm2d09fWfxul9PUkl8HqYkTdkwR+DnAScuXpjkaODZwOcmXJMkaQgDA7yqLgHuXGLVnwMvB3wWpiTNwEhj4ElOAW6rqm0TrkeSNKSBDzVeLMkDgd+mN3wyTPtNwCaADRs2rLY7SdIyRjkC/07gEcC2JDuAo4Crkhy2VOOq2lJV81U1Pzc3N3qlkqT7WfUReFVdA3zHnvkuxOer6o4J1iVJGmCYywjPBy4Fjk2yM8lZa1+WJGmQgUfgVXXagPUbJ1aNJGlo3okpSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjRrmkWrnJtmd5Nq+ZX+S5IYkn0zyviQHr22ZkqTFhjkCPw84cdGyi4Djqup7gU8Dr5hwXZKkAQYGeFVdAty5aNmHq+qebvYy4Kg1qE2StIJJjIH/PPCh5VYm2ZRka5KtCwsLE+hOkgRjBniS3wHuAd6+XJuq2lJV81U1Pzc3N053kqQ+60Z9Y5IzgecCJ1RVTawiSdJQRgrwJCcCLweeVlX/M9mSJEnDGOYywvOBS4Fjk+xMchbwRuAg4KIkVyd50xrXKUlaZOAReFWdtsTiv12DWiRJq+CdmJLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSoYR6pdm6S3Umu7Vt2SJKLktzY/Xzo2pYpSVpsmCPw84ATFy3bDHykqh4JfKSblyRN0cAAr6pLgDsXLT4FeEs3/Rbg+ROuS5I0wKhj4IdW1a5u+nbg0OUaJtmUZGuSrQsLCyN2J0labOyTmFVVQK2wfktVzVfV/Nzc3LjdSZI6owb4F5IcDtD93D25kiRJwxg1wC8EzuimzwD+YTLlSJKGNcxlhOcDlwLHJtmZ5CzgbOBHktwIPKublyRN0bpBDarqtGVWnTDhWiRJq+CdmJLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSosQI8ya8luS7JtUnOT3LApAqTJK1s5ABPciTwUmC+qo4D9gNOnVRhkqSVjTuEsg44MMk64IHA58cvSZI0jJEDvKpuA14LfA7YBXypqj68uF2STUm2Jtm6sLAweqWSpPsZZwjlocApwCOAI4AHJXnR4nZVtaWq5qtqfm5ubvRKJUn3M84QyrOAz1TVQlV9DXgv8OTJlCVJGmScAP8c8MQkD0wS4ARg+2TKkiQNMs4Y+OXABcBVwDXdtrZMqC5J0gDrxnlzVb0SeOWEapEkrYJ3YkpSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1KixrgOfpo2bPzCzvnecffLM+pak5XgELkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjRorwJMcnOSCJDck2Z7kSZMqTJK0snFvpX898E9V9cIk+wMPnEBNkqQhjBzgSR4CPBU4E6Cq7gbunkxZkqRBxhlCeQSwALw5ySeSnJPkQYsbJdmUZGuSrQsLC2N0J0nqN06ArwOOB/6qqh4PfAXYvLhRVW2pqvmqmp+bmxujO0lSv3ECfCews6ou7+YvoBfokqQpGDnAq+p24NYkx3aLTgCun0hVkqSBxr0K5SXA27srUG4Bfm78kiRJwxgrwKvqamB+QrVIklbBOzElqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWrUuDfySN8SNm7+wMz63nH2yTPrW23zCFySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckho1doAn2S/JJ5K8fxIFSZKGM4kj8F8Ftk9gO5KkVRgrwJMcBZwMnDOZciRJwxr3CPx1wMuBry/XIMmmJFuTbF1YWBizO0nSHiMHeJLnArur6sqV2lXVlqqar6r5ubm5UbuTJC0yzhH4U4DnJdkBvBN4ZpK3TaQqSdJAIwd4Vb2iqo6qqo3AqcBHq+pFE6tMkrQirwOXpEZN5Ik8VXUxcPEktiVJGo5H4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSocZ5Kf3SSjyW5Psl1SX51koVJklY2ziPV7gF+vaquSnIQcGWSi6rq+gnVJklawThPpd9VVVd103cB24EjJ1WYJGllE3mocZKNwOOBy5dYtwnYBLBhw4ZJdCdJI9m4+QMz63vH2SdPfJtjn8RM8mDg74CXVdV/L15fVVuqar6q5ufm5sbtTpLUGSvAk3w7vfB+e1W9dzIlSZKGMc5VKAH+FtheVX82uZIkScMY5wj8KcDPAs9McnX3OmlCdUmSBhj5JGZV/TuQCdYiSVoF78SUpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRo37UOMTk3wqyU1JNk+qKEnSYOM81Hg/4C+AHwUeA5yW5DGTKkyStLJxjsCfANxUVbdU1d3AO4FTJlOWJGmQkR9qDBwJ3No3vxP4wcWNkmwCNnWzX07yqTH6HGQ9cMekN5o/nvQWl7Um9U+JtY9oAp8v9/1srKr2Mf87P3ypheME+FCqaguwZa37AUiytarmp9HXWmi5fmufnZbrt/bxjDOEchtwdN/8Ud0ySdIUjBPgHwcemeQRSfYHTgUunExZkqRBRh5Cqap7krwY+GdgP+DcqrpuYpWNZipDNWuo5fqtfXZart/ax5CqmnUNkqQReCemJDXKAJekRjUd4El+Msl1Sb6eZNnLeZLsSHJNkquTbJ1mjStZRf173VcWJDkkyUVJbux+PnSZdvd2+/3qJDM9yT1oPyZ5QJJ3desvT7Jx+lUubYjaz0yy0Levf2EWdS4lyblJdie5dpn1SfKG7t/2ySTHT7vG5QxR+9OTfKlvv//eVAusqmZfwKOBY4GLgfkV2u0A1s+63lHqp3eC+GbgGGB/YBvwmL2g9tcAm7vpzcAfL9Puy7Ouddj9CPwK8KZu+lTgXbOuexW1nwm8cda1LlP/U4HjgWuXWX8S8CEgwBOBy2dd8ypqfzrw/lnV1/QReFVtr6q1vLNzTQ1Z/976lQWnAG/ppt8CPH+GtQxjmP3Y/2+6ADghSaZY43L21s/AUKrqEuDOFZqcAry1ei4DDk5y+HSqW9kQtc9U0wG+CgV8OMmV3a39LVnqKwuOnFEt/Q6tql3d9O3Aocu0OyDJ1iSXJZllyA+zH+9rU1X3AF8CHjaV6lY27GfgJ7ohiAuSHL3E+r3V3voZH9aTkmxL8qEkj51mx2t+K/24kvwLcNgSq36nqv5hyM38UFXdluQ7gIuS3ND9Zl1zE6p/JlaqvX+mqirJctejPrzb98cAH01yTVXdPOlaxT8C51fVV5P8Ir2/JJ4545r2BVfR+4x/OclJwN8Dj5xW53t9gFfVsyawjdu6n7uTvI/en6RTCfAJ1D+zryxYqfYkX0hyeFXt6v7c3b3MNvbs+1uSXAw8nt547rQNsx/3tNmZZB3wEOCL0ylvRQNrr6r+Os+hd46iFc1+LUdV/Xff9AeT/GWS9VU1lS/o+pYfQknyoCQH7ZkGng0seUZ5L7W3fmXBhcAZ3fQZwDf9NZHkoUke0E2vB54CXD+1Cu9vmP3Y/296IfDR6s5UzdjA2heNGT8P2D7F+sZ1IXB6dzXKE4Ev9Q3P7dWSHLbnPEmSJ9DL1On90p/1Wd4xzxC/gN542VeBLwD/3C0/AvhgN30MvbP224Dr6A1dzLz2Yevv5k8CPk3vyHWvqJ/e2PBHgBuBfwEO6ZbPA+d0008Grun2/TXAWTOu+Zv2I/Aq4Hnd9AHAe4CbgCuAY2a9n1dR+6u7z/c24GPAd8+65r7azwd2AV/rPu9nAb8E/FK3PvQeDnNz9zlZ9oqyvbD2F/ft98uAJ0+zPm+ll6RGfcsPoUjStyoDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXq/wEXGCNkf4/zeAAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.hist(xyz_avg[:,1])\n",
+    "plt.title('Average $y(t)$');"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "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.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

apps/py/ch11_notebooks/lorenz.py

@@ -0,0 +1,46 @@
+from matplotlib import pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import numpy as np
+from scipy import integrate
+
+def solve_lorenz(sigma=10.0, beta=8./3, rho=28.0):
+    """Plot a solution to the Lorenz differential equations."""
+
+    max_time = 4.0
+    N = 30
+
+    fig = plt.figure()
+    ax = fig.add_axes([0, 0, 1, 1], projection='3d')
+    ax.axis('off')
+
+    # prepare the axes limits
+    ax.set_xlim((-25, 25))
+    ax.set_ylim((-35, 35))
+    ax.set_zlim((5, 55))
+    
+    def lorenz_deriv(x_y_z, t0, sigma=sigma, beta=beta, rho=rho):
+        """Compute the time-derivative of a Lorenz system."""
+        x, y, z = x_y_z
+        return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]
+
+    # Choose random starting points, uniformly distributed from -15 to 15
+    np.random.seed(1)
+    x0 = -15 + 30 * np.random.random((N, 3))
+
+    # Solve for the trajectories
+    t = np.linspace(0, max_time, int(250*max_time))
+    x_t = np.asarray([integrate.odeint(lorenz_deriv, x0i, t)
+                      for x0i in x0])
+    
+    # choose a different color for each trajectory
+    colors = plt.cm.viridis(np.linspace(0, 1, N))
+
+    for i in range(N):
+        x, y, z = x_t[i,:,:].T
+        lines = ax.plot(x, y, z, '-', c=colors[i])
+        plt.setp(lines, linewidth=2)
+    angle = 104
+    ax.view_init(30, angle)
+    plt.show()
+
+    return t, x_t