Reordered material to more reasonable order.
This commit is contained in:
Roger Labbe committed
1 parent
148db5f889
commit
95cbe78f46
1 file changed
+477
-170
+477
-170
@@ -109,7 +109,7 @@
|
||||
" color: #1d3b84;\n",
|
||||
" font-size: 16pt;\n",
|
||||
" margin-bottom: 0em;\n",
|
||||
" margin-top: 1.5em;\n",
|
||||
" margin-top: 0.5em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" }\n",
|
||||
@@ -133,6 +133,11 @@
|
||||
" overflow-y: scroll;\n",
|
||||
" max-height: 50000px;\n",
|
||||
" }\n",
|
||||
" div.output_wrapper{\n",
|
||||
" margin-top:0.2em;\n",
|
||||
" margin-bottom:0.2em;\n",
|
||||
"}\n",
|
||||
"\n",
|
||||
" code{\n",
|
||||
" font-size: 70%;\n",
|
||||
" }\n",
|
||||
@@ -438,6 +443,154 @@
|
||||
"However, we took advantage of something I left mostly unstated in the last chapter. We were able to provide a definition for $\\mathbf{F}$ because we were able to take advantage of the exact solution that Newtonian equations provide us. However, if you have an engineering background you will realize what a small class of problems that covers. If you don't, I will explain it next, and provide you with several ways to compute $\\mathbf{F}$ for arbitrary systems."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Walking Through the Kalman Filter Equations"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"I promised that you would not have to understand how to derive Kalman filter equations, and that is true. However, I do think it is worth walking through the equations one by one and becoming familiar with the variables. If this is your first time through the material feel free to skip ahead to the next section. However, you will eventually want to work through this material, so why not now? You will need to have passing familiarity with these equations to read material written about the Kalman filter, as they all presuppose that you are familiar with them. I will reiterate them here for easy reference.\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"\\text{Predict Step}\\\\\n",
|
||||
"\\mathbf{x} &= \\mathbf{F x} + \\mathbf{B u}\\;\\;\\;\\;&(1) \\\\\n",
|
||||
"\\mathbf{P} &= \\mathbf{FP{F}}^\\mathsf{T} + \\mathbf{Q}\\;\\;\\;\\;&(2) \\\\\n",
|
||||
"\\\\\n",
|
||||
"\\text{Update Step}\\\\\n",
|
||||
"\\textbf{y} &= \\mathbf{z} - \\mathbf{H}\\mathbf{x}\\;\\;\\;&(3) \\\\\n",
|
||||
"\\mathbf{S} &= \\mathbf{HPH}^\\mathsf{T} + \\mathbf{R} \\;\\;\\;&(4) \\\\\n",
|
||||
"\\mathbf{K} &= \\mathbf{PH}^\\mathsf{T}\\mathbf{S}^{-1}\\;\\;\\;&(5) \\\\\n",
|
||||
"\\mathbf{x} &= \\mathbf{x} +\\mathbf{K}\\mathbf{y} \\;\\;\\;&(6)\\\\\n",
|
||||
"\\mathbf{P} &= (\\mathbf{I}-\\mathbf{K}\\mathbf{H})\\mathbf{P}\\;\\;\\;&(7)\n",
|
||||
"\\end{aligned}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"I will start with the update step, as that is what we started with in the one dimensional Kalman filter case. The first equation is\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\mathbf{y} = \\mathbf{z} - \\mathbf{H x}\\tag{3}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"On the right we have $\\mathbf{Hx}$. That should be recognizable as the measurement function. Multiplying $\\mathbf{H}$ with $\\mathbf{x}$ puts $\\mathbf{x}$ into *measurement space*; in other words, the same basis and units as the sensor's measurements. The variable $\\mathbf{z}$ is just the measurement; it is typical, but not universal to use $\\mathbf{z}$ to denote measurements in the literature ($\\mathbf{y}$ is also sometimes used). Do you remember this chart?"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAEjCAYAAAD5ZS3PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYVGX/BvD7DPs6AgKCCrgg5oK470oKpkJi5mtqpmJu\nqYWavUq9uaS/LC3LrNc2A1NTMdfUUFHAEEjBHXcB93FBQUFZ5/n9wcvkOMMAOuz357q49DznOed8\nz5mFew7POSMJIQSIiIiIiOiFySq7ACIiIiKimoLhmoiIiIhITxiuiYiIiIj0hOGaiIiIiEhPGK6J\niIiIiPSE4ZqIiIiISE8YromqITc3NzRq1Kiyy6gStB2L0NBQyGQyrF69upKqqt6ioqIgk8mwYMEC\ntXZvb2/IZLXz18b8+fMhk8lw8ODByi6lzGrz40ZUGfhqI9IzmUxW7C+y1NRUeHh4QCaTISgoCC9y\nm3lJkp572erAzc1NdSy1/cycORNA4XEo7lg82y6TyarMh5KxY8dq7JOFhQWaN2+OoKAg3Lp1q7JL\n1Dh+uo51aRR96Hk2tNdERYE2Ojq62D5Fgb0ijoc+3i+q0uuHqCozrOwCiGoibb/ITpw4gQEDBuD2\n7dtYvHgxZs+eXQmVVT/Tp09HnTp1NNq7dOkCADhw4ECZ1lfVPpQMHjwYXl5eAIA7d+4gPDwcK1as\nwIYNG/D333/Dzc2tcgt8yq+//oonT5688Hqq2mNQXkr7YaQ6HY/qVCtRZWG4JqoAkZGRGDx4MLKz\nsxEaGoq33nqrskuqFiRJwvTp0+Hi4lJsn+p+Jm3w4MEYPXq0ajo/Px+vvPIKIiMjsWjRIvz888+V\nWJ26hg0b6mU9teWLgWvLfhKROg4LISpnYWFh6N+/P4QQ+OOPPzSCdV5eHr799lsMHDgQrq6uMDU1\nha2tLXx8fLBr165Sb+fpP7knJCSgf//+kMvlsLW1xb/+9S/cuHEDAHD+/HkMHToUdevWhbm5Ofr0\n6YNTp05prO/ChQuYM2cOOnToAHt7e5iamsLNzQ0TJkzAtWvXNPoXjdMNDAxEamoqhg8fjrp168LM\nzAwdO3Ys076URWnGnxfVBhQOzXl6KEZgYKBa30uXLmH8+PGqx8LBwQFDhgzBsWPHNNZb9Gf91atX\nY/fu3ejVqxesra1ha2v73PtjaGiIyZMnAwAOHz6stp8ymQy5ubmYP38+3N3dYWJighkzZjxX7QBw\n+/ZtvP3223B0dIS5uTnatm2LX3/9tdjadI3d3bdvHwYNGgRHR0eYmpqiYcOG8Pf3x86dOwEUDoMZ\nN24cAGDBggVqj0FpxjFv27YNo0aNQrNmzWBpaQlLS0u0b98ey5cvh1Kp1OhfNOwmOjoav//+Ozp1\n6gQLCwvY2dlhxIgRuHnzptbtJCYmon///rCysoJcLoevry/i4+NLrE+fnn6sP/zwQ7i5ucHU1BTu\n7u5YtGgR8vLytC63YcMGtG/fHubm5nB0dMTo0aOL3c+yvO+U1+uHqKbimWuicrRixQpMnz4d9vb2\n2LVrF9q3b6/RJy0tDdOnT0f37t3xyiuvwN7eHjdv3sQff/yBV199Fd9//z0mTpxY6m0eOXIES5Ys\ngY+PDyZNmoT4+Hhs3rwZp06dwqZNm9CzZ0+0b98e48aNQ1JSEv7880/4+PggOTkZFhYWqvVs2bIF\nP/zwA/r06YMePXrA2NgYp0+fxi+//II//vgDiYmJqF+/vsb2r1y5gs6dO6NJkyYYM2YM0tLSsHHj\nRgQEBCAiIgLe3t5lOoalOftX0p+qGzVqhHnz5mHBggWQy+VqgbRoSAZQOMQkICAAubm58Pf3h7u7\nO65fv44tW7bgzz//xPbt29GvXz+N9W/atAl79uyBv78/pk6ditu3b5dhDzUVhUVt+1UUVAYMGIC6\ndeuqPliUtfZ79+6hW7duSElJQffu3dGrVy/cvHkT77zzDnx8fIqtTVtN8+bNw8KFC2FpaYnBgwfD\nxcUFt27dQnx8PH755Rf4+/vjtddeQ0ZGBrZv3w5vb2+154Grq2uJxyQ4OBgGBgbo2rUr6tevj4yM\nDOzfvx8zZszA4cOHsW7dOq3L/fe//8WOHTsQEBCAl19+GfHx8di4cSNOnDiB48ePw9jYWNU3NjYW\nPj4+yMvLw5AhQ+Du7o4TJ07A29sbffr0KbFGfRs6dCiOHj2KoUOHwsjICFu3bsXcuXORmJiIrVu3\nqvX96quv8P7776NOnToYPXo0bGxssGfPHnTv3h1yuVxj3WV53ynv1w9RjSOISK8kSRIymUzMmTNH\nSJIkmjZtKi5fvlxs/5ycHHHjxg2N9oyMDNGqVStha2srnjx5ojbP1dVVNGrUSK0tJCRESJIkJEkS\nW7ZsUZvXv39/IUmSkMvlYtmyZWrzJk6cKCRJEsuXL1drv3HjhsjNzdWoa+/evcLAwEBMnjxZrT0y\nMlK1/U8++URt3p49e4QkSWLgwIHFHAVNrq6uQpIkMX36dDFv3jy1n88++6xUx2L16tVq7ZIkafQt\nkp6eLuzs7ETdunXF2bNn1eadPXtWWFlZCWdnZ5GTk6NqnzdvnpAkSRgYGIg9e/aUet+EEGLMmDFa\na8zNzRXe3t5CkiQxYcIEjePRpk0bkZaW9sK1T5gwQUiSJN577z21/seOHRPGxsZCkiSxYMECtXm9\ne/cWMplMra3osW3UqJG4fv26xn4+3Vb0uDy73tJITk7WaFMqlarjGB8frzavqF0ul4vTp0+rzRs5\ncqSQJEmEhYWprcvDw0PIZDKN1893332nel1HR0eXqt7evXsLSZJ09i96/jx7PIoeaw8PD5Genq5q\nf/LkiejUqZOQJEmsX79e1Z6SkiKMjIyEjY2NSElJUdunYcOGqWp/2vO87+j79UNUUzFcE+lZUcCU\nJEkYGxuLixcvPve6vvzySyFJkjh48KBau65A+fLLL2usZ82aNaqg/6yDBw8KSZLEuHHjSl1X69at\nRePGjdXaisJ1o0aNhFKp1FjGxcVF2Nvbl3obRQFD24+NjY1aP32E62+++UZIkiRWrFihdf6MGTOE\nJEli9+7dqraicDRkyJBS71eRovA3ePBg1YeGKVOmiMaNGwtJkoSjo6NITU3VOB47dux44dpzc3OF\nubm5sLKyUgtvRcaNG1fqcO3v7y8kSRK///57ifv8IuG6OImJiUKSJLFw4UK19qLj+/HHH2ssU/Rc\n/eCDD1RtMTExQpIk0aNHD43+SqVSuLu7lxiWn6aPcL127VqNZfbu3SskSRK+vr6qtkWLFglJksR/\n/vMfjf4pKSnCwMBA43HTpbj3HX2/fohqKg4LISon/fv3R3h4OEaOHIk9e/bAxsam2L5JSUlYunQp\nDh48CIVCgezsbLX5xY2b1KZt27YabfXq1QMAeHp6asxzcnICAFy/fl1j3tq1axEaGooTJ04gPT0d\nBQUFqnkmJiZat+/l5aV16EDDhg3x999/l24n/keSJKSkpOi8oFFfDh06BKDwri7z58/XmH/+/HkA\nwNmzZzFgwAC1eZ06dXru7W7fvh3bt28HAJiamsLV1RXvvvsu5syZo3psikiSpHVbZa393LlzePLk\nCbp166Z1yECvXr0QEhJSqvrj4+MhSZLGMdG3tLQ0LF26FLt370ZycjIeP36sNr/omoJndejQQaOt\nQYMGAIAHDx6o2o4ePQoA6N27t0Z/SZLQvXt3XLp06bnrLytJkrTW0rNnTwDA8ePHVW26andzc0PD\nhg1x9epVjXn6fN95kdcPUU3DcE1UDiRJwo4dOzBs2DBs27YN3t7eiIiIgL29vUbf+Ph49OnTB0ql\nEn379sXgwYNhbW0NmUyGY8eOYfv27cjJySn1trWFJUNDwxLnPXuR1IwZM7B8+XI4OztjwIABqF+/\nPszMzAAAISEhWn9ZA9B627yi7Wi78KyqSEtLAwCsWrWq2D6SJCErK0ujvejDy/MIDQ1Vu1tISRwd\nHTXaylp7RkZGsevS1a5Neno6rK2tYW5uXuplyio9PR0dO3ZEamoqOnfujLFjx8LW1haGhoZ48OAB\nli9fXuxrRNvzseg5//SHRX0eEwCqCwB1PeeL5hV3kai2bZqamsLa2lpVL1C62p99ver7fedFXj9E\nNQ3DNVE5MTQ0xKZNmzBmzBj89ttv6NWrFyIiIjQuAly0aBGys7MRFRWFXr16qc1bvHix6qxmRbpz\n5w6++eYbtG7dGrGxsWoXOgIo9uKx6qzog8fRo0fVLtIqjcq+929Zay/qX9yFl2W5ILNOnTq4f/8+\nsrKyNJ4n+vLzzz8jNTUV8+fPx9y5c9XmxcXFYfny5S+8DX0ek6fXVxQ6tbl37x6A4j+QKhQKjdsf\nZmdn4+HDh6hbt67Gtm7fvo3WrVuXqnZ9v++8yOuHqKbhrfiIypGBgQHWrFmD8ePH4/z58+jVqxdS\nU1PV+ly6dAl2dnYav+AA6Px2t/KUnJwMIQT69eunEZiuX7+O5OTkSqnrRUmSpHa28mndunUDgGr5\n9dZlrf2ll16CmZkZTp48ifT0dI35ZXnede3aFUII/PnnnyX2NTAwAIBiH4PiFA3HeP311zXm6es1\nUnQnn6ioKI15SqUSMTExZVpfUcDUtVxsbCwAoE2bNhrzhBBa963oMX56+Jeu2lNSUrTeOvN53ndq\n6uuHSN8YronKmSRJ+PHHHxEUFISUlBT07NkTFy5cUM1v1KgR0tLSNO41vWrVKuzdu7eiy1XVBAB/\n/fWX2p+1MzMzMWHChDKHo6rCzs4Od+/e1RhbCgCBgYGwsbHBwoULtd7XWAiBmJiYYu8xXJnKWruh\noSFGjRqFzMxMjTPBx44dw9q1a0u97XfffRcA8MEHH2gdt//0WOiis61Xrlwp9fqBfLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x518ca90>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import mkf_internal\n",
|
||||
"mkf_internal.show_residual_chart()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The blue point labeled \"prediction\" is the output of $\\mathbf{Hx}$, and the dot labeled \"measurement\" is $\\mathbf{z}$. Therefore, $\\mathbf{y} = \\mathbf{z} - \\mathbf{Hx}$ is how we compute the residual, drawn in red, where $\\mathbf{z}$ is the residual.\n",
|
||||
"\n",
|
||||
"The next two lines are the formidable:\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"\\mathbf{S} &= \\textbf{HPH}^\\mathsf{T} + \\textbf{R} \\;\\;\\;&(4) \\\\\n",
|
||||
"\\textbf{K} &= \\textbf{PH}^\\mathsf{T}\\mathbf{S}^{-1}\\;\\;\\;&(5) \\\\\n",
|
||||
"\\end{aligned}\n",
|
||||
"$$\n",
|
||||
"Unfortunately it is a fair amount of linear algebra to derive this. The derivation can be quite elegant, and I urge you to look it up if you have the mathematical education to follow it. But $\\mathbf{K}$ is just the *Kalman gain* - the ratio of how much measurement vs prediction we should use to create the new estimate. $\\mathbf{R}$ is the *measurement noise*, and $\\mathbf{P}$ is our *uncertainty covariance matrix* from the prediction step.\n",
|
||||
"\n",
|
||||
"So let's work through this expression by expression. Start with $\\mathbf{HPH}^\\mathsf{T}$. The linear equation $\\mathbf{ABA}^T$ can be thought of as changing $\\mathbf{B}$ to the coordinate system of $\\mathbf{A}$. So $\\mathbf{HPH}^\\mathsf{T}$ is taking the covariance $\\mathbf{P}$ and putting it in measurement ($\\mathbf{H}$) space. \n",
|
||||
"\n",
|
||||
"In English, consider the problem of reading a temperature with a thermometer that provides readings in volts. Our state is in terms of temperature, but we are now doing calculations in *measurement space* - volts. So we need to convert $\\mathbf{P}$ from applying to temperatures to volts. The linear algebra form $\\textbf{H}\\textbf{P}\\textbf{H}^\\mathsf{T}$ takes $\\mathbf{P}$ to the basis used by $\\mathbf{H}$, namely volts. \n",
|
||||
"\n",
|
||||
"Then, once in measurement space, we can add the measurement noise $\\mathbf{R}$ to it. Hence, the expression for the uncertainty once we include the measurement is:\n",
|
||||
"\n",
|
||||
"$$\\mathbf{S} = \\mathbf{HP}\\mathbf{H}^\\mathsf{T} + \\mathbf{R}$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The next equation is\n",
|
||||
"$$\\textbf{K} = \\textbf{P}\\textbf{H}^T\\mathbf{S}^{-1}\\\\\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"$\\mathbf{K}$ is the *Kalman gain* - the ratio that chooses how far along the residual to select between the measurement and prediction in the graph above.\n",
|
||||
"\n",
|
||||
"We can think of the inverse of a matrix as linear algebra's way of computing $\\frac{1}{x}$. So we can read the equation for $\\textbf{K}$ as\n",
|
||||
"\n",
|
||||
"$$ \\textbf{K} = \\frac{\\textbf{P}\\textbf{H}^\\mathsf{T}}{\\mathbf{S}} $$\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\textbf{K} = \\frac{uncertainty_{prediction}}{uncertainty_{measurement}}\\textbf{H}^\\mathsf{T}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"In other words, the *Kalman gain* equation is doing nothing more than computing a ratio based on how much we trust the prediction vs the measurement. If we are confident in our measurements and unconfident in our predictions $\\mathbf{K}$ will favor the measurement, and vice versa. The equation is complicated because we are doing this in multiple dimensions via matrices, but the concept is simple - scale by a ratio.\n",
|
||||
"\n",
|
||||
"Without going into the derivation of $\\mathbf{K}$, I'll say that this equation is the result of finding a value of $\\mathbf{K}$ that optimizes the *mean-square estimation error*. It does this by finding the minimal values for $\\mathbf{P}$ along its diagonal. Recall that the diagonal of $\\mathbf{P}$ is just the variance for each state variable. So, this equation for $\\mathbf{K}$ ensures that the Kalman filter output is optimal. To put this in concrete terms, for our dog tracking problem this means that the estimates for both position and velocity will be optimal - a value of $\\mathbf{K}$ that made the position extremely accurate but the velocity very inaccurate would be rejected in favor of a $\\mathbf{K}$ that made both position and velocity just somewhat accurate."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Our next line is:\n",
|
||||
" $$\\mathbf{x}=\\mathbf{x}' +\\mathbf{Ky}\\tag{5}$$\n",
|
||||
"\n",
|
||||
"This just multiplies the residual by the Kalman gain, and adds it to the state variable. In other words, this is the computation of our new estimate.\n",
|
||||
"\n",
|
||||
"Finally, we have:\n",
|
||||
"\n",
|
||||
"$$\\mathbf{P}=(\\mathbf{I}-\\mathbf{KH})\\mathbf{P}\\tag{6}$$\n",
|
||||
"\n",
|
||||
"$I$ is the identity matrix, and is the way we represent $1$ in multiple dimensions. $H$ is our measurement function, and is a constant. So, simplified, this is simply $P = (1-cK)P$. $K$ is our ratio of how much prediction vs measurement we use. So, if $K$ is large then $(1-cK)$ is small, and P will be made smaller than it was. If $K$ is small, then $(1-cK)$ is large, and P will be made larger than it was. So we adjust the size of our uncertainty by some factor of the *Kalman gain*. I would like to draw your attention back to the g-h filter, which included this Python code:\n",
|
||||
"\n",
|
||||
" # update filter \n",
|
||||
" w = w * (1-scale_factor) + z * scale_factor\n",
|
||||
"\n",
|
||||
"This multidimensional Kalman filter equation is partially implementing this calculation for the variance instead of the state variable."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now we have the measurement steps. The first equation is\n",
|
||||
"\n",
|
||||
"$$\\mathbf{x} = \\mathbf{Fx} + \\mathbf{Bu}\\tag{1}$$\n",
|
||||
"\n",
|
||||
"This is just our state transition equation which we have already discussed. $\\mathbf{Fx}$ multiplies $\\mathbf{x}$ with the state transition matrix to compute the next state. $B$ and $u$ add in the contribution of the control input $\\mathbf{u}$, if any.\n",
|
||||
"\n",
|
||||
"The final equation is:\n",
|
||||
"$$\\mathbf{P} = \\mathbf{FPF}^\\mathsf{T} + \\mathbf{Q}\\tag{2}$$\n",
|
||||
"\n",
|
||||
"$\\mathbf{FPF}^\\mathsf{T}$ is the way we put $\\mathbf{P}$ into the process space using linear algebra so that we can add the process noise $\\mathbf{Q}$ to it."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
@@ -451,7 +604,10 @@
|
||||
"source": [
|
||||
"Modeling dynamic systems is properly the topic of at least one undergraduate course in mathematics. To an extent there is no substitute for a few semesters of ordinary and partial differential equations. If you are a hobbyist, or trying to solve one very specific filtering problem at work you probably do not have the time and/or inclination to devote a year or more to that education.\n",
|
||||
"\n",
|
||||
"However, I can present enough of the theory to allow us to create the system equations for many different Kalman filters, and give you enough background to at least follow the mathematics in the literature. My goal is to get you to the stage where you can read a Kalman filtering book or paper and understand it well enough to implement the algorithms. The background math is deep, but we end up using a few simple techniques over and over again in practice. \n",
|
||||
"However, I can present enough of the theory to allow us to create the system equations for many different Kalman filters, and give you enough background to at least follow the mathematics in the literature. My goal is to get you to the stage where you can read a Kalman filtering book or paper and understand it well enough to implement the algorithms. The background math is deep, but we end up using a few simple techniques over and over again in practice.\n",
|
||||
"\n",
|
||||
"I struggle a bit with the proper way to present this material. If you have not encountered this math before I fear reading this section will not be very profitable for you. In the **Extended Kalman Filter** chapter I take a more ad-hoc way of presenting this information where I expose a problem that the KF needs to solve, then provide the math without a lot of supporting theory. This gives you the motivation behind the mathematics at the cost of not knowing why the math I give you is correct. On the other hand, the following section gives you the math, but somewhat divorced from the specifics of the problem we are trying to solve. Only you know what kind of learner your are. If you like the presentation of the book so far (practical first, then the math) you may want to wait until you read the **Extended Kalman Filter** before \n",
|
||||
"In particular, if your intent is to work with Extended Kalman filters (a very prelevant form of nonlinear Kalman filtering) you will need to understand this math at least at the level I present it. If that is not your intent this section may still prove to be beneficial if you need to simulate a nonlinear system in order to test your filter.\n",
|
||||
"\n",
|
||||
"Let's lay out the problem and discuss what the solution will be. We model *dynamic systems* with a set of first order *differential equations*. This should not be a surprise as calculus is the math of things that vary. For example, we say that velocity is the derivative of distance with respect to time\n",
|
||||
"\n",
|
||||
@@ -535,7 +691,11 @@
|
||||
"\n",
|
||||
"$$x(t) = \\Phi(t-t_0)x(t_0)$$\n",
|
||||
"\n",
|
||||
"In other words, we just want to compute the value of $x$ at time $t$ by multiplying its previous value by some matrix $\\Phi$. This is not trivial to do because the original equations do not include time "
|
||||
"which we can equivalently write as\n",
|
||||
"\n",
|
||||
"$$x(t_k) = \\Phi(\\Delta t)x(t_{k-1})$$\n",
|
||||
"\n",
|
||||
"In other words, we just want to compute the value of $x$ at time $t$ by multiplying its previous value by some matrix $\\Phi$. This is not trivial to do because the original equations do not include time."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -551,9 +711,104 @@
|
||||
"\n",
|
||||
"$$ \\Phi(t) = e^{\\mathbf{F}t} = \\mathbf{I} + \\mathbf{F}t + \\frac{(\\mathbf{F}t)^2}{2!} + \\frac{(\\mathbf{F}t)^3}{3!} + ... $$\n",
|
||||
"\n",
|
||||
"This is much easy to compute, and thus is the typical approach used in Kalman filter design when the filter is reasonably small. If you are wondering where $e$ came from, I point you to the Wikipedia article on the matrix exponential [1]. Here the important point is to recognize the very simple and regular form this equation takes. \n",
|
||||
"This is easy to compute, and thus is the typical approach used in Kalman filter design when the filter is reasonably small. If you are wondering where $e$ came from, I point you to the Wikipedia article on the matrix exponential [1]. Here the important point is to recognize the very simple and regular form this equation takes. \n",
|
||||
"\n",
|
||||
"Before applying it to a Kalman filter, lets do the Taylor expansion of a math function since this is much easier to visualize than a matrix. I choose sin(x). The Taylor series for a real or complex function f(x) at "
|
||||
"Before applying it to a Kalman filter, lets do the Taylor expansion of a math function since this is much easier to visualize than a matrix. I choose sin(x). The Taylor series for a real or complex function f(x) about x=a is\n",
|
||||
"\n",
|
||||
"$$f(x) = \\sum_{n=0}^\\infty\\frac{f^{(n)}(a)}{x!}(x-a)^n$$\n",
|
||||
"\n",
|
||||
"where $f^{n}$ is the nth derivative of f. To compute the Taylor series for $f(x)=sin(x)$ at $x=0$ Let's first work out the terms for f.\n",
|
||||
"\n",
|
||||
"$$\\begin{aligned}\n",
|
||||
"f^{0}(x) &= sin(x) ,\\ \\ &f^{0}(0) &= 0 \\\\\n",
|
||||
"f^{1}(x) &= cos(x),\\ \\ &f^{1}(0) &= 1 \\\\\n",
|
||||
"f^{2}(x) &= -sin(x),\\ \\ &f^{2}(0) &= 0 \\\\\n",
|
||||
"f^{3}(x) &= -cos(x),\\ \\ &f^{3}(0) &= -1 \\\\\n",
|
||||
"f^{4}(x) &= sin(x),\\ \\ &f^{4}(0) &= 0 \\\\\n",
|
||||
"f^{5}(x) &= cos(x),\\ \\ &f^{5}(0) &= 1\n",
|
||||
"\\end{aligned}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"Now we can substitute these values into the equation.\n",
|
||||
"\n",
|
||||
"$$f(x) = \\frac{0}{0!}(x)^0 + \\frac{1}{1!}(x)^1 + \\frac{0}{2!}(x)^2 + \\frac{-1}{3!}(x)^3 + \\frac{0}{4!}(x)^4 + \\frac{-1}{5!}(x)^5 + ... $$\n",
|
||||
"\n",
|
||||
"And let's test this with some code:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 34,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"estimate of sin(.3) is 0.30452025\n",
|
||||
"exact value of sin(.3) is 0.295520206661\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"x = .3\n",
|
||||
"estimate = x + x**3/6 + x**5/120\n",
|
||||
"exact = np.sin(.3)\n",
|
||||
"\n",
|
||||
"print('estimate of sin(.3) is', estimate)\n",
|
||||
"print('exact value of sin(.3) is', exact)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This is not bad for only three terms. If you are curious, go ahead and implement this as a Python function to compute the series for an arbitrary number of terms. But I will forge ahead to the matrix form of the equation. \n",
|
||||
"\n",
|
||||
"Let's consider tracking an object moving in a vacuum. In one dimesion the differential equation for motion with zero acceleration is\n",
|
||||
"\n",
|
||||
"$$ v = \\dot{x}\\\\a=\\ddot{x} =0,$$\n",
|
||||
"\n",
|
||||
"which we can put in state space matrix form as\n",
|
||||
"\n",
|
||||
"$$\\begin{bmatrix}\\dot{x} \\\\ \\ddot{x}\\end{bmatrix} =\\begin{bmatrix}0&1\\\\0&0\\end{bmatrix} \\begin{bmatrix}x \\\\ \\dot{x}\\end{bmatrix}$$\n",
|
||||
"\n",
|
||||
"This is a first order differential equation, so we can set $\\mathbf{F}=\\begin{bmatrix}0&1\\\\0&0\\end{bmatrix}$ and solve the following equation.\n",
|
||||
"\n",
|
||||
"$$\\Phi(t) = e^{\\mathbf{F}t} = \\mathbf{I} + \\mathbf{F}t + \\frac{(\\mathbf{F}t)^2}{2!} + \\frac{(\\mathbf{F}t)^3}{3!} + ... $$\n",
|
||||
"\n",
|
||||
"If you perform the multiplication you will find that $\\mathbf{F}^2=\\begin{bmatrix}0&0\\\\0&0\\end{bmatrix}$, which means that all higher powers of $\\mathbf{F}$ are also $\\mathbf{0}$. This makes the computation very easy.\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"\\Phi(t) &=\\mathbf{I} + \\mathbf{F}t + \\mathbf{0} \\\\\n",
|
||||
"&= \\begin{bmatrix}1&0\\\\0&1\\end{bmatrix} + \\begin{bmatrix}0&1\\\\0&0\\end{bmatrix}t\\\\\n",
|
||||
"&= \\begin{bmatrix}1&t\\\\0&1\\end{bmatrix}\n",
|
||||
"\\end{aligned}$$\n",
|
||||
"\n",
|
||||
"We plug this into $x(t_k) = \\Phi(\\Delta t)x(t_{k-1})$ to get\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"x(t_k) &= \\Phi(\\Delta t)x(t_{k-1}) \\\\\n",
|
||||
"x^- &= \\Phi(\\Delta t)x \\\\\n",
|
||||
"x^- &=\\begin{bmatrix}1&\\Delta t\\\\0&1\\end{bmatrix}x \\\\\n",
|
||||
"\\begin{bmatrix}x \\\\ \\dot{x}\\end{bmatrix}^- &=\\begin{bmatrix}1&\\Delta t\\\\0&1\\end{bmatrix}\\begin{bmatrix}x \\\\ \\dot{x}\\end{bmatrix}\n",
|
||||
"\\end{aligned}$$\n",
|
||||
"\n",
|
||||
"This should look very familiar to you! This is the equation we used in the **Multivariate Kalman Filter** chapter to track a moving object.\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"{\\begin{bmatrix}x\\\\\\dot{x}\\end{bmatrix}}^- =\\begin{bmatrix}1&t \\\\ 0&1\\end{bmatrix} \\begin{bmatrix}x \\\\ \\dot{x}\\end{bmatrix}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"We derived this equation in that chapter by using techniques that are much easier to understand. The advantage of the Taylor series expansion is that we can use it for any arbitrary set of differential equations which are time invariant. \n",
|
||||
"\n",
|
||||
"Time invariant? This just means equations that do not depend on time. For example, $f(x) = sin(t)t$ is not time invariant; the value will be different at different times due to the multiplication by t. On the other hand $f(x) = sin(t)*3$ is time invariant. $x(t)$ is dependent on time, but that is allowed. This is because to be time invariant it must be true that if $f(t) = x(t)$ then $f(t+\\Delta t) = x(t+\\Delta t)$. This will be true for the second equation, but not the first.\n",
|
||||
"\n",
|
||||
"However, we often use a Taylor expansion even when the equations are not time invariant. The answer will still be reasonably accurate so long as the time step is short and the system is nearly constant over that time step."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -589,7 +844,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -696,154 +951,6 @@
|
||||
"If you practice this a bit you will become adept at it. Just isolate the highest term, define a new variable and its derivatives, and then substitute."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Walking Through the Kalman Filter Equations"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"I promised that you would not have to understand how to derive Kalman filter equations, and that is true. However, I do think it is worth walking through the equations one by one and becoming familiar with the variables. If this is your first time through the material feel free to skip ahead to the next section. However, you will eventually want to work through this material, so why not now? You will need to have passing familiarity with these equations to read material written about the Kalman filter, as they all presuppose that you are familiar with them. I will reiterate them here for easy reference.\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"\\text{Predict Step}\\\\\n",
|
||||
"\\mathbf{x} &= \\mathbf{F x} + \\mathbf{B u}\\;\\;\\;\\;&(1) \\\\\n",
|
||||
"\\mathbf{P} &= \\mathbf{FP{F}}^\\mathsf{T} + \\mathbf{Q}\\;\\;\\;\\;&(2) \\\\\n",
|
||||
"\\\\\n",
|
||||
"\\text{Update Step}\\\\\n",
|
||||
"\\textbf{y} &= \\mathbf{z} - \\mathbf{H}\\mathbf{x}\\;\\;\\;&(3) \\\\\n",
|
||||
"\\mathbf{S} &= \\mathbf{HPH}^\\mathsf{T} + \\mathbf{R} \\;\\;\\;&(4) \\\\\n",
|
||||
"\\mathbf{K} &= \\mathbf{PH}^\\mathsf{T}\\mathbf{S}^{-1}\\;\\;\\;&(5) \\\\\n",
|
||||
"\\mathbf{x} &= \\mathbf{x} +\\mathbf{K}\\mathbf{y} \\;\\;\\;&(6)\\\\\n",
|
||||
"\\mathbf{P} &= (\\mathbf{I}-\\mathbf{K}\\mathbf{H})\\mathbf{P}\\;\\;\\;&(7)\n",
|
||||
"\\end{aligned}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"I will start with the update step, as that is what we started with in the one dimensional Kalman filter case. The first equation is\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\mathbf{y} = \\mathbf{z} - \\mathbf{H x}\\tag{3}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"On the right we have $\\mathbf{Hx}$. That should be recognizable as the measurement function. Multiplying $\\mathbf{H}$ with $\\mathbf{x}$ puts $\\mathbf{x}$ into *measurement space*; in other words, the same basis and units as the sensor's measurements. The variable $\\mathbf{z}$ is just the measurement; it is typical, but not universal to use $\\mathbf{z}$ to denote measurements in the literature ($\\mathbf{y}$ is also sometimes used). Do you remember this chart?"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAEjCAYAAAD5ZS3PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYVGX/BvD7DPs6AgKCCrgg5oK470oKpkJi5mtqpmJu\nqYWavUq9uaS/LC3LrNc2A1NTMdfUUFHAEEjBHXcB93FBQUFZ5/n9wcvkOMMAOuz357q49DznOed8\nz5mFew7POSMJIQSIiIiIiOiFySq7ACIiIiKimoLhmoiIiIhITxiuiYiIiIj0hOGaiIiIiEhPGK6J\niIiIiPSE4ZqIiIiISE8YromqITc3NzRq1Kiyy6gStB2L0NBQyGQyrF69upKqqt6ioqIgk8mwYMEC\ntXZvb2/IZLXz18b8+fMhk8lw8ODByi6lzGrz40ZUGfhqI9IzmUxW7C+y1NRUeHh4QCaTISgoCC9y\nm3lJkp572erAzc1NdSy1/cycORNA4XEo7lg82y6TyarMh5KxY8dq7JOFhQWaN2+OoKAg3Lp1q7JL\n1Dh+uo51aRR96Hk2tNdERYE2Ojq62D5Fgb0ijoc+3i+q0uuHqCozrOwCiGoibb/ITpw4gQEDBuD2\n7dtYvHgxZs+eXQmVVT/Tp09HnTp1NNq7dOkCADhw4ECZ1lfVPpQMHjwYXl5eAIA7d+4gPDwcK1as\nwIYNG/D333/Dzc2tcgt8yq+//oonT5688Hqq2mNQXkr7YaQ6HY/qVCtRZWG4JqoAkZGRGDx4MLKz\nsxEaGoq33nqrskuqFiRJwvTp0+Hi4lJsn+p+Jm3w4MEYPXq0ajo/Px+vvPIKIiMjsWjRIvz888+V\nWJ26hg0b6mU9teWLgWvLfhKROg4LISpnYWFh6N+/P4QQ+OOPPzSCdV5eHr799lsMHDgQrq6uMDU1\nha2tLXx8fLBr165Sb+fpP7knJCSgf//+kMvlsLW1xb/+9S/cuHEDAHD+/HkMHToUdevWhbm5Ofr0\n6YNTp05prO/ChQuYM2cOOnToAHt7e5iamsLNzQ0TJkzAtWvXNPoXjdMNDAxEamoqhg8fjrp168LM\nzAwdO3Ys076URWnGnxfVBhQOzXl6KEZgYKBa30uXLmH8+PGqx8LBwQFDhgzBsWPHNNZb9Gf91atX\nY/fu3ejVqxesra1ha2v73PtjaGiIyZMnAwAOHz6stp8ymQy5ubmYP38+3N3dYWJighkzZjxX7QBw\n+/ZtvP3223B0dIS5uTnatm2LX3/9tdjadI3d3bdvHwYNGgRHR0eYmpqiYcOG8Pf3x86dOwEUDoMZ\nN24cAGDBggVqj0FpxjFv27YNo0aNQrNmzWBpaQlLS0u0b98ey5cvh1Kp1OhfNOwmOjoav//+Ozp1\n6gQLCwvY2dlhxIgRuHnzptbtJCYmon///rCysoJcLoevry/i4+NLrE+fnn6sP/zwQ7i5ucHU1BTu\n7u5YtGgR8vLytC63YcMGtG/fHubm5nB0dMTo0aOL3c+yvO+U1+uHqKbimWuicrRixQpMnz4d9vb2\n2LVrF9q3b6/RJy0tDdOnT0f37t3xyiuvwN7eHjdv3sQff/yBV199Fd9//z0mTpxY6m0eOXIES5Ys\ngY+PDyZNmoT4+Hhs3rwZp06dwqZNm9CzZ0+0b98e48aNQ1JSEv7880/4+PggOTkZFhYWqvVs2bIF\nP/zwA/r06YMePXrA2NgYp0+fxi+//II//vgDiYmJqF+/vsb2r1y5gs6dO6NJkyYYM2YM0tLSsHHj\nRgQEBCAiIgLe3t5lOoalOftX0p+qGzVqhHnz5mHBggWQy+VqgbRoSAZQOMQkICAAubm58Pf3h7u7\nO65fv44tW7bgzz//xPbt29GvXz+N9W/atAl79uyBv78/pk6ditu3b5dhDzUVhUVt+1UUVAYMGIC6\ndeuqPliUtfZ79+6hW7duSElJQffu3dGrVy/cvHkT77zzDnx8fIqtTVtN8+bNw8KFC2FpaYnBgwfD\nxcUFt27dQnx8PH755Rf4+/vjtddeQ0ZGBrZv3w5vb2+154Grq2uJxyQ4OBgGBgbo2rUr6tevj4yM\nDOzfvx8zZszA4cOHsW7dOq3L/fe//8WOHTsQEBCAl19+GfHx8di4cSNOnDiB48ePw9jYWNU3NjYW\nPj4+yMvLw5AhQ+Du7o4TJ07A29sbffr0KbFGfRs6dCiOHj2KoUOHwsjICFu3bsXcuXORmJiIrVu3\nqvX96quv8P7776NOnToYPXo0bGxssGfPHnTv3h1yuVxj3WV53ynv1w9RjSOISK8kSRIymUzMmTNH\nSJIkmjZtKi5fvlxs/5ycHHHjxg2N9oyMDNGqVStha2srnjx5ojbP1dVVNGrUSK0tJCRESJIkJEkS\nW7ZsUZvXv39/IUmSkMvlYtmyZWrzJk6cKCRJEsuXL1drv3HjhsjNzdWoa+/evcLAwEBMnjxZrT0y\nMlK1/U8++URt3p49e4QkSWLgwIHFHAVNrq6uQpIkMX36dDFv3jy1n88++6xUx2L16tVq7ZIkafQt\nkp6eLuzs7ETdunXF2bNn1eadPXtWWFlZCWdnZ5GTk6NqnzdvnpAkSRgYGIg9e/aUet+EEGLMmDFa\na8zNzRXe3t5CkiQxYcIEjePRpk0bkZaW9sK1T5gwQUiSJN577z21/seOHRPGxsZCkiSxYMECtXm9\ne/cWMplMra3osW3UqJG4fv26xn4+3Vb0uDy73tJITk7WaFMqlarjGB8frzavqF0ul4vTp0+rzRs5\ncqSQJEmEhYWprcvDw0PIZDKN1893332nel1HR0eXqt7evXsLSZJ09i96/jx7PIoeaw8PD5Genq5q\nf/LkiejUqZOQJEmsX79e1Z6SkiKMjIyEjY2NSElJUdunYcOGqWp/2vO87+j79UNUUzFcE+lZUcCU\nJEkYGxuLixcvPve6vvzySyFJkjh48KBau65A+fLLL2usZ82aNaqg/6yDBw8KSZLEuHHjSl1X69at\nRePGjdXaisJ1o0aNhFKp1FjGxcVF2Nvbl3obRQFD24+NjY1aP32E62+++UZIkiRWrFihdf6MGTOE\nJEli9+7dqraicDRkyJBS71eRovA3ePBg1YeGKVOmiMaNGwtJkoSjo6NITU3VOB47dux44dpzc3OF\nubm5sLKyUgtvRcaNG1fqcO3v7y8kSRK///57ifv8IuG6OImJiUKSJLFw4UK19qLj+/HHH2ssU/Rc\n/eCDD1RtMTExQpIk0aNHD43+SqVSuLu7lxiWn6aPcL127VqNZfbu3SskSRK+vr6qtkWLFglJksR/\n/vMfjf4pKSnCwMBA43HTpbj3HX2/fohqKg4LISon/fv3R3h4OEaOHIk9e/bAxsam2L5JSUlYunQp\nDh48CIVCgezsbLX5xY2b1KZt27YabfXq1QMAeHp6asxzcnICAFy/fl1j3tq1axEaGooTJ04gPT0d\nBQUFqnkmJiZat+/l5aV16EDDhg3x999/l24n/keSJKSkpOi8oFFfDh06BKDwri7z58/XmH/+/HkA\nwNmzZzFgwAC1eZ06dXru7W7fvh3bt28HAJiamsLV1RXvvvsu5syZo3psikiSpHVbZa393LlzePLk\nCbp166Z1yECvXr0QEhJSqvrj4+MhSZLGMdG3tLQ0LF26FLt370ZycjIeP36sNr/omoJndejQQaOt\nQYMGAIAHDx6o2o4ePQoA6N27t0Z/SZLQvXt3XLp06bnrLytJkrTW0rNnTwDA8ePHVW26andzc0PD\nhg1x9epVjXn6fN95kdcPUU3DcE1UDiRJwo4dOzBs2DBs27YN3t7eiIiIgL29vUbf+Ph49OnTB0ql\nEn379sXgwYNhbW0NmUyGY8eOYfv27cjJySn1trWFJUNDwxLnPXuR1IwZM7B8+XI4OztjwIABqF+/\nPszMzAAAISEhWn9ZA9B627yi7Wi78KyqSEtLAwCsWrWq2D6SJCErK0ujvejDy/MIDQ1Vu1tISRwd\nHTXaylp7RkZGsevS1a5Neno6rK2tYW5uXuplyio9PR0dO3ZEamoqOnfujLFjx8LW1haGhoZ48OAB\nli9fXuxrRNvzseg5//SHRX0eEwCqCwB1PeeL5hV3kai2bZqamsLa2lpVL1C62p99ver7fedFXj9E\nNQ3DNVE5MTQ0xKZNmzBmzBj89ttv6NWrFyIiIjQuAly0aBGys7MRFRWFXr16qc1bvHix6qxmRbpz\n5w6++eYbtG7dGrGxsWoXOgIo9uKx6qzog8fRo0fVLtIqjcq+929Zay/qX9yFl2W5ILNOnTq4f/8+\nsrKyNJ4n+vLzzz8jNTUV8+fPx9y5c9XmxcXFYfny5S+8DX0ek6fXVxQ6tbl37x6A4j+QKhQKjdsf\nZmdn4+HDh6hbt67Gtm7fvo3WrVuXqnZ9v++8yOuHqKbhrfiIypGBgQHWrFmD8ePH4/z58+jVqxdS\nU1PV+ly6dAl2dnYav+AA6Px2t/KUnJwMIQT69eunEZiuX7+O5OTkSqnrRUmSpHa28mndunUDgGr5\n9dZlrf2ll16CmZkZTp48ifT0dI35ZXnede3aFUII/PnnnyX2NTAwAIBiH4PiFA3HeP311zXm6es1\nUnQnn6ioKI15SqUSMTExZVpfUcDUtVxsbCwAoE2bNhrzhBBa963oMX56+Jeu2lNSUrTeOvN53ndq\n6uuHSN8YronKmSRJ+PHHHxEUFISUlBT07NkTFy5cUM1v1KgR0tLSNO41vWrVKuzdu7eiy1XVBAB/\n/fWX2p+1MzMzMWHChDKHo6rCzs4Od+/e1RhbCgCBgYGwsbHBwoULtd7XWAiBmJiYYu8xXJnKWruh\noSFGjRqFzMxMjTPBx44dw9q1a0u97XfffRcA8MEHH2gdt//0WOiis61Xrlwp9fqBfLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x5534278>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import mkf_internal\n",
|
||||
"mkf_internal.show_residual_chart()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The blue point labeled \"prediction\" is the output of $\\mathbf{Hx}$, and the dot labeled \"measurement\" is $\\mathbf{z}$. Therefore, $\\mathbf{y} = \\mathbf{z} - \\mathbf{Hx}$ is how we compute the residual, drawn in red, where $\\mathbf{z}$ is the residual.\n",
|
||||
"\n",
|
||||
"The next two lines are the formidable:\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"\\mathbf{S} &= \\textbf{HPH}^\\mathsf{T} + \\textbf{R} \\;\\;\\;&(4) \\\\\n",
|
||||
"\\textbf{K} &= \\textbf{PH}^\\mathsf{T}\\mathbf{S}^{-1}\\;\\;\\;&(5) \\\\\n",
|
||||
"\\end{aligned}\n",
|
||||
"$$\n",
|
||||
"Unfortunately it is a fair amount of linear algebra to derive this. The derivation can be quite elegant, and I urge you to look it up if you have the mathematical education to follow it. But $\\mathbf{K}$ is just the *Kalman gain* - the ratio of how much measurement vs prediction we should use to create the new estimate. $\\mathbf{R}$ is the *measurement noise*, and $\\mathbf{P}$ is our *uncertainty covariance matrix* from the prediction step.\n",
|
||||
"\n",
|
||||
"So let's work through this expression by expression. Start with $\\mathbf{HPH}^\\mathsf{T}$. The linear equation $\\mathbf{ABA}^T$ can be thought of as changing the basis of $\\mathbf{B}$ to $\\mathbf{A}$. So $\\mathbf{HPH}^\\mathsf{T}$ is taking the covariance $\\mathbf{P}$ and putting it in measurement ($\\mathbf{H}$) space. \n",
|
||||
"\n",
|
||||
"In English, consider the problem of reading a temperature with a thermometer that provides readings in volts. Our state is in terms of temperature, but we are now doing calculations in *measurement space* - volts. So we need to convert $\\mathbf{P}$ from applying to temperatures to volts. The linear algebra form $\\textbf{H}\\textbf{P}\\textbf{H}^\\mathsf{T}$ takes $\\mathbf{P}$ to the basis used by $\\mathbf{H}$, namely volts. \n",
|
||||
"\n",
|
||||
"Then, once in measurement space, we can add the measurement noise $\\mathbf{R}$ to it. Hence, the expression for the uncertainty once we include the measurement is:\n",
|
||||
"\n",
|
||||
"$$\\mathbf{S} = \\mathbf{HP}\\mathbf{H}^\\mathsf{T} + \\mathbf{R}$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The next equation is\n",
|
||||
"$$\\textbf{K} = \\textbf{P}\\textbf{H}^T\\mathbf{S}^{-1}\\\\\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"$\\mathbf{K}$ is the *Kalman gain* - the ratio that chooses how far along the residual to select between the measurement and prediction in the graph above.\n",
|
||||
"\n",
|
||||
"We can think of the inverse of a matrix as linear algebra's way of computing $\\frac{1}{x}$. So we can read the equation for $\\textbf{K}$ as\n",
|
||||
"\n",
|
||||
"$$ \\textbf{K} = \\frac{\\textbf{P}\\textbf{H}^\\mathsf{T}}{\\mathbf{S}} $$\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\textbf{K} = \\frac{uncertainty_{prediction}}{uncertainty_{measurement}}\\textbf{H}^\\mathsf{T}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"In other words, the *Kalman gain* equation is doing nothing more than computing a ratio based on how much we trust the prediction vs the measurement. If we are confident in our measurements and unconfident in our predictions $\\mathbf{K}$ will favor the measurement, and vice versa. The equation is complicated because we are doing this in multiple dimensions via matrices, but the concept is simple - scale by a ratio.\n",
|
||||
"\n",
|
||||
"Without going into the derivation of $\\mathbf{K}$, I'll say that this equation is the result of finding a value of $\\mathbf{K}$ that optimizes the *mean-square estimation error*. It does this by finding the minimal values for $\\mathbf{P}$ along its diagonal. Recall that the diagonal of $\\mathbf{P}$ is just the variance for each state variable. So, this equation for $\\mathbf{K}$ ensures that the Kalman filter output is optimal. To put this in concrete terms, for our dog tracking problem this means that the estimates for both position and velocity will be optimal - a value of $\\mathbf{K}$ that made the position extremely accurate but the velocity very inaccurate would be rejected in favor of a $\\mathbf{K}$ that made both position and velocity just somewhat accurate."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Our next line is:\n",
|
||||
" $$\\mathbf{x}=\\mathbf{x}' +\\mathbf{Ky}\\tag{5}$$\n",
|
||||
"\n",
|
||||
"This just multiplies the residual by the Kalman gain, and adds it to the state variable. In other words, this is the computation of our new estimate.\n",
|
||||
"\n",
|
||||
"Finally, we have:\n",
|
||||
"\n",
|
||||
"$$\\mathbf{P}=(\\mathbf{I}-\\mathbf{KH})\\mathbf{P}\\tag{6}$$\n",
|
||||
"\n",
|
||||
"$I$ is the identity matrix, and is the way we represent $1$ in multiple dimensions. $H$ is our measurement function, and is a constant. So, simplified, this is simply $P = (1-cK)P$. $K$ is our ratio of how much prediction vs measurement we use. So, if $K$ is large then $(1-cK)$ is small, and P will be made smaller than it was. If $K$ is small, then $(1-cK)$ is large, and P will be made larger than it was. So we adjust the size of our uncertainty by some factor of the *Kalman gain*. I would like to draw your attention back to the g-h filter, which included this Python code:\n",
|
||||
"\n",
|
||||
" # update filter \n",
|
||||
" w = w * (1-scale_factor) + z * scale_factor\n",
|
||||
"\n",
|
||||
"This multidimensional Kalman filter equation is partially implementing this calculation for the variance instead of the state variable."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now we have the measurement steps. The first equation is\n",
|
||||
"\n",
|
||||
"$$\\mathbf{x} = \\mathbf{Fx} + \\mathbf{Bu}\\tag{1}$$\n",
|
||||
"\n",
|
||||
"This is just our state transition equation which we have already discussed. $\\mathbf{Fx}$ multiplies $\\mathbf{x}$ with the state transition matrix to compute the next state. $B$ and $u$ add in the contribution of the control input $\\mathbf{u}$, if any.\n",
|
||||
"\n",
|
||||
"The final equation is:\n",
|
||||
"$$\\mathbf{P} = \\mathbf{FPF}^\\mathsf{T} + \\mathbf{Q}\\tag{2}$$\n",
|
||||
"\n",
|
||||
"$\\mathbf{FPF}^\\mathsf{T}$ is the way we put $\\mathbf{P}$ into the process space using linear algebra so that we can add in the process noise $\\mathbf{Q}$ to it."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
@@ -893,11 +1000,18 @@
|
||||
"\n",
|
||||
"for the fundamental matrix, and\n",
|
||||
"\n",
|
||||
"$$\\mathbf{Q_c} = \\begin{bmatrix}0&0&0\\\\0&0&0\\\\0&0&1\\end{bmatrix} \\Phi_s$$\n",
|
||||
"\n",
|
||||
"$$\\mathbf{Q_c} = \\begin{bmatrix}0&0&0\\\\0&0&0\\\\0&0&1\\end{bmatrix} \\Phi_s$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"for the continuous process noise matrix, where $\\Phi_s$ is the spectral density of the white noise.\n",
|
||||
"\n",
|
||||
"We could carry out these computations ourselves, but I prefer using SymPy to solve the equation."
|
||||
"We could carry out these computations ourselves, but I prefer using SymPy to solve the equation.\n",
|
||||
"\n",
|
||||
"$$\\mathbf{Q_c} = \\begin{bmatrix}0&0&0\\\\0&0&0\\\\0&0&1\\end{bmatrix} \\Phi_s$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1429,7 +1543,200 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"** material about Euler here**\n"
|
||||
"## Euler's Method\n",
|
||||
"\n",
|
||||
"Let's say we have the initial condition problem of \n",
|
||||
"\n",
|
||||
"$$ y' = y, \\\\ y(0) = 1$$\n",
|
||||
"\n",
|
||||
"I happen to know the exact answer is $y=e^t$, but of course in general we will not know the exact solution. In general all we know is the derivative of the equation, which is equal to the slope. We also know the initial value: at $t=0$, $y=1$. If we know these two pieces of information we can predict that value for y(1): it is the slope of the function at $t=0$. I've plotted this below"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAEWCAYAAADIP0muAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlU1Pe9P/7nzMAM+7DvyKIIuIAgGAEX1Ghjszb5Zalp\njH577m3vL+lNyMm5qeeYxDTetjFNvm3vqfbcc38Sbpu0SZrFpE2MWDcUjSAgyObCvgw7M6wDzLx/\nfwwMEhhgYBi25+McDzKf93w+L5JPJk/fvj7vt0QIIUBERERERFYjnesCiIiIiIiWGoZwIiIiIiIr\nYwgnIiIiIrIyhnAiIiIiIitjCCciIiIisjKGcCIiIiIiK2MIJyIiIiKysglD+B/+8AfExMRAqVRC\nqVQiKSkJX3311YQnLCwsxNatW+Hg4IDAwEC8+eabFi2YiIiIiGihs5noYFBQEI4cOYLw8HDo9Xq8\n9957eOSRR5CdnY2YmJgx4zUaDXbu3ImUlBTk5OSgpKQE+/fvh6OjI1566aVZ+yGIiIiIiBYSibk7\nZnp4eODXv/41/uVf/mXMsWPHjuHAgQNobGyEQqEAAPznf/4njh07htraWstUTERERES0wE25J1yn\n0+Gvf/0r+vr6sGXLlnHHXL58GZs3bzYGcADYtWsX6uvrUVVVNfNqiYiIiIgWgQnbUQBDj3diYiK0\nWi3s7e3x0UcfISIiYtyxKpUKy5YtG/Waj4+P8VhwcLAFSiYiIiIiWtgmDeGRkZEoKCiAWq3Gxx9/\njKeeegpnz55FfHz8mLESicSsi6vVarPGExERERHNJ0qlclrvmzSE29raIiwsDAAQGxuL7Oxs/OEP\nf0BaWtqYsb6+vlCpVKNea2xsNB4jIiIiIqJprBOu0+mg1+vHPZaYmIjMzExotVrjaxkZGQgICGAr\nChERERHRkAlnwn/+85/jgQceQGBgIDo7O/HBBx/g/PnzOHnyJADgwIEDyM7OxunTpwEAe/bswRtv\nvIF9+/bh4MGDKCsrw1tvvYVDhw5NWsh0p/KJrCUnJwcAxm3FIppveL/SQsL7lRYaS7RUTxjCGxsb\n8aMf/QgqlQpKpRIxMTE4efIkdu7cCcDwsGV5eblxvIuLCzIyMvDcc88hPj4e7u7uePnll5Gamjrj\nQomIiIiIFosJQ/h4fd+THV+zZg3Onz8/s6qIiIiIiBYxs3vCiYiIiIhoZhjCiYiIiIisjCGciIiI\niMjKGMKJiIiIiKyMIZyIiIiIyMoYwomIiIiIrIwhnIiIiIjIyhjCiYiIiIisjCGciIiIiMjKGMKJ\niIiIiKyMIZyIiIiIyMoYwomIiIiIrIwhnIiIiIjIyhjCiYiIiIisjCGciIiIiMjKGMKJiIiIiKZg\nYLAfBXe+tci5bCxyFiIiIiKiRUoIgbxbl/DFpf9Fm6YJb+5Ln/E5GcKJiIiIiEyoVN3Epxf+P1Q2\nlFn0vAzhRERERETf0aZpxpdZf8K1sguzcn6GcCIiIiKiIX39vTid8wnO5n6BAV3/mONSqcwi12EI\nJyIiIqIlT6/X4UrxGfzj8vvo7OkYd8yasA14eNOzFrkeQzgRERERLWll1dfxWWYa6lsqxz0e4BmC\nRzbvR8SyGACAWq2e8TUZwomIiIhoSWpsq8XnF99DUUXOuMddHNxwf9LTuCdqm8XaUIYxhBMRERHR\nktLdq8HX336Ii4Unodfrxhy3lcmxff3D2LH+UdjJ7WelBoZwIiIiIloSBnUDuHD9K3xz9SP0arvH\nHRMfuRUPJv0Ibs5es1oLQzgRERERLWpCCBTc+RYnLr6HFrVq3DFhflH4wZb9CPZdaZWaGMKJiIiI\naNGqbryNzzLTcKeuaNzjHi4+eGjTs1i3IhESicRqdTGEExEREdGi09HVir9n/RnZJecgIMYct5M7\n4HsbHseWmPthayO3en0M4URERES0aGgH+vDPa5/hzLXP0T+oHXNcIpEiec0u7N74Qzg7KOegQgPp\nZAN+9atfISEhAUqlEt7e3njooYdQVDT+dP6wyspKSKXSMb9OnTplscKJiIiIiIbphR7fFp/B4fT/\nFye//XDcAL4qOA4/f/q3eGL7T+c0gANTmAk/f/48nn/+eSQkJECv1+O1117Dvffei+LiYri5uU34\n3m+++QYxMTHG7ycbT0RERERkrlu1N/BZ5nHUNpWPe9zPYxke2bwfUcGxVq7MtElD+MmTJ0d9/6c/\n/QlKpRJZWVm4//77J3yvu7s7vL29Z1YhEREREdE4mjsacOJiOgruXBn3uJO9Evcn7sHG1fdCZuHN\ndmbK7J5wjUYDvV4/pVntRx99FH19fQgPD0dqaioee+yxaRVJRERERDSsp68LJ69+hMzrX0GnHxxz\nXCazwbZ1D2FnwmOwVzjOQYWTMzuEv/DCC4iNjUViYqLJMc7OznjnnXeQnJwMGxsbnDhxAk8++STS\n09Px9NNPz6hgIiIiIlqadLpBXCw8ia+//RA9fZ3jjokNT8ZDyXvhofSxcnXmkQghxq7ZYsJLL72E\njz76CBcvXkRISIhZF3r++eeRmZmJ69evG19Tq9XG39+6dcus8xERERHR0iCEQF37beRUnoamt3Xc\nMZ5O/ogP3Qlvl6BZryc8PNz4e6Vyeg94TnkmPDU1FR999BHOnj1rdgAHgISEBBw/ftzs9xERERHR\n0tXW3Yicigyo1JXjHneQuyAuZDtCPVdbdbOdmZpSCH/hhRfw8ccf4+zZs1i5cnpbeebn58Pf39/k\n8fj4+Gmdl8hacnJyAPBepYWB9ystJLxfaTya7nb84/IHuFJ0etzNdhS2dtgZ/xhS4h6C3EZh1dru\n7uaYrklD+HPPPYc///nP+Pzzz6FUKqFSqQAY+r4dHQ2N7gcOHEB2djZOnz4NAEhPT4dcLse6desg\nlUrx5Zdf4ujRozhy5MiMCyYiIiKixat/UItzuV8gI+cTaAf6xhyXQIKNq+/F/Yl74OK4cJe/njSE\nHzt2DBKJBDt27Bj1+qFDh/Daa68BAFQqFcrLR9ZllEgkOHz4MKqqqiCTyRAREYG0tDTs2bPHwuUT\nERER0WKgF3rklmXiy0t/QntXy7hjVgZF4web9yPAK9TK1VnepCFcr9dPepK0tLRR3+/duxd79+6d\nflVEREREtGSU15fgswvHUdU4/kId3m4BeGTTPqwOjV9Qfd8TMXuJQiIiIiIiS2hVN+KLS/+LvFuX\nxj3uYOeM3fc8iU1r74NMtrhi6+L6aYiIiIho3uvVduNU9t9wLv9L6HTjbLYjtcGWmO/jexuegIOd\n0xxUOPsYwomIiIjIKnR6HS7fyMBXV/6Crt7xVxiJXr4RDyXvhbeb6VX1FgOGcCIiIiKadSVVefg8\nMw0NrdXjHg/0DsMPNv8fhAeusXJlc4MhnIiIiIhmTUNrDT7PTENJVe64x5WO7ngg6UdIiEqBVCK1\ncnVzhyGciIiIiCyus0eNr6/8BVk3TkEvxq62Z2sjx471P8CO9T+AwtZuDiqcWwzhRERERGQxA4P9\nOJ//d5zK/hv6+nvGHbMhahseSPoRXJ08rFzd/MEQTkREREQzJoRA/u0sfHHxf9GqaRx3zPKA1fjB\n5v1Y5rPCytXNPwzhRERERDQjVaqb+PTCcVQ0lI573FPpi4c37UP08nsWzWY7M8UQTkRERETT0qZp\nxt+z/oycsvPjHreXO+B79zyJLTHfh43M1srVzW8M4URERERkFm1/L05f+xRnrp3AgK5/zHGpRIpN\n0btx3z1PwsneZQ4qnP8YwomIiIhoSvR6Hb4tPoN/XP4Amp72ccesDo3HI5v2wcc90MrVLSwM4URE\nREQ0qZs1BfjswnHUtVSOe9zfMwQ/2LwfEctirFvYAsUQTkREREQmNbbX4UTme7hRkT3ucWcHV9yf\n+DQ2rtoOqVRm5eoWLoZwIiIiIhqju1eDk1c/QmbB19DrdWOO28hssT3uYdwb/xjs5PZzUOHCxhBO\nREREREaDugFkXv8aJ69+iF5t97hj1kdswYNJP4K7i7eVq1s8GMKJiIiICEIIFJZ/ixOZ6WhWN4w7\nJsQvAj/Y/H8Q6hdh5eoWH4ZwIiIioiWupukOPstMw+3aG+Med3f2wkObnkVseDI327EQhnAiIiKi\nJaq5owGnrn6MqyVnISDGHFfI7bEr4XGkrHsAtjbyOahw8WIIJyIiIlpiyutLcSb3cxTe+Xbc8C2R\nSJG0Zhe+v/EpODu4zkGFix9DOBEREdESoNfrUHDnW5zJPYFKVZnJcZHL1uGRzfvh7xlsxeqWHoZw\nIiIiokVMO9CHb4v/ibN5X6BV3WhynK97EB7ZvB+rQuKsWN3SxRBOREREtAipu9uQef0rXCw4iR5t\nl8lxAV6h2B73MOJWboaMm+1YDUM4ERER0SLS0FqNM7knkFN2HjrdoMlxq4LjsC3uYawMiuaKJ3OA\nIZyIiIhogRNC4GZNAc7knkBJVa7JcTKZDRIitmJb3MPw81hmxQrpuxjCiYiIiBYonW4Qubcu4kzu\nCdQ1V5gc56Bwwqbo3dgS8324OLpZsUIyhSGciIiIaIHp1XYj68YpnMv/O9RdrSbHeSh9sC32Ydyz\najsUtnZWrJAmwxBOREREtEC0aZpwLv/vuHzjFLQDfSbHhfhFYEfcI1gbtgFSPmw5LzGEExEREc1z\n1Y23cSb3BPJvXYJe6McdI4EE0cvvwba4RxDmH2nlCslcE4bwX/3qV/j0009x8+ZNKBQKbNy4Eb/6\n1a+wevXqCU9aWFiI559/HtnZ2XB3d8dPfvITvPrqqxYtnIiIiGgx0ws9iiuu4Uzu57hdV2RynK2N\nHBtX3YuU2Afh5epnxQppJiYM4efPn8fzzz+PhIQE6PV6vPbaa7j33ntRXFwMN7fxmLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x81bc6a0>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"t = np.linspace(-1, 1, 10)\n",
|
||||
"plt.plot(t, np.exp(t))\n",
|
||||
"t = np.linspace(0, 1, 2)\n",
|
||||
"plt.plot(t,t+1)\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"If the curve is relatively straight than the line formed by the slope will not be very far from the curve. Here the step size of 1 is rather large; you can see that the slope is very close to the line at, say, 0.1. But let's continue with a step size of 1 for a moment. We can see that at $t=1$ the estimated value of $y$ is 2. Now we can compute the value at $t=2$ by taking the slope of the curve at $t=1$ and adding it to our initial estimate. The slope is computed with $y'=y$, so the slope is 2."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAusAAAEtCAYAAACxoC4JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9P/7XzCQzmUySSci+kQRIgIQtJOyCslbAUlrr\ngkuR3t+1tepV/Lb30Xqpdat6rbVVr9iqFWwRqhZFrSiChB2ELGxhSSAkZCfrJJlMJpOZz++PJDMM\nWciEzJyZyev5eOTB5JzPmXnH48FXPn7O+8iEEAJEREREROR25FIXQEREREREvWNYJyIiIiJyUwzr\nRERERERuimGdiIiIiMhNMawTEREREbkphnUiIiIiIjfFsE5ERERE5KYkDesdHR148sknMWrUKKjV\naowaNQq//e1vYTabpSyLiIiIiMgt+Ej54S+88AL++te/4u9//zsmTpyIEydO4IEHHoBKpcK6deuk\nLI2IiIiISHKShvVjx45hxYoVWL58OQBg5MiRuO2223D06FEpyyIiIiIicguSLoNZunQpdu/ejfPn\nzwMAzpw5g6ysLCxbtkzKsoiIiIiI3IKkM+u/+MUvUFZWhvHjx8PHxwcdHR1Yt24dfv7zn0tZFhER\nERGRW5A0rL/++uvYsGED/vnPfyItLQ15eXl47LHHkJiYiJ/+9KfWcTqdTsIqiYiIiIhujFarHdRx\nMiGEGOJaBiwyMhLr1q3Do48+at32+9//Hhs3bkRhYaF1G8M6EREREXmywYZ1SdesCyEgl9uXIJfL\nIeHvD0REREREbkPSZTArV67ESy+9hKSkJKSmpiIvLw9/+tOfsHr16j6PGexvJSS97OxsAEBmZqbE\nldCN4rn0HjyX3oPn0nvwXHqPoVgdImlY/9Of/oSgoCA8/PDDqK6uRnR0NB588EE89dRTUpZFRERE\nROQWJA3rGo0Gr7zyCl555RUpyyAiIiIickuSrlknIiIiIqK+MawTEREREbkphnUiIiIiIjfFsE5E\nRERE5KYY1omIiIiI3BTDOhERERGRm2JYJyIiIiJyUwzrRERERERuimGdiIiIiMhNMawTEREREbkp\nhnUiIiIiIjfFsE5ERERENIQ6zCaUVBUMyXsxrBMRERERDRGDUY+/fPYcXv/XuiF5P58heRciIiIi\nomFO11KPtz57FhW1xUP2ngzrREREREQ3qKq+FG9texYNzTVD+r4M60REREREN+Bi+Rm888ULaDW2\nDPl7M6wTEREREQ3SiQuH8f7Xr6LDbLLbPiFp2pC8P8M6EREREdEg7DvxJbbueRcCwm777AlLcMf8\nn6Gl+cZn2hnWiYiIiIgcYBEWfHHwH/g259Me+5bPugdLpt0BmUw2JJ8laevGxMREyOXyHl+33Xab\nlGUREREREfWqw2zCph2v9Qjqcpkc9yx6FN+bfueQBXVA4pn1nJwcmM1m6/cVFRXIyMjAXXfdJWFV\nREREREQ9GYyt+NuXL6Gg9KTddqWvH3667FdITcwY8s+UNKyHhobaff/OO+9Aq9XizjvvlKgiIiIi\nIqKedC31+Mtnz6L8mh7qgWotfvaD32Jk5BinfK7brFkXQuBvf/sb7rvvPqhUKqnLISIiIiIC0HcP\n9fDgGDy08imEaaOc9tluE9Z37tyJ4uJi/Od//qfUpRARERERAQCKKs7i7c9/36OHekJkMh5csQ6B\n/lqnfr5MCCGuP8z57rjjDpSWluLIkSM99ul0OuvrwsJCV5ZFRERERMPU5bpz2F+wDWZLh932uJBk\nzB37Q/gqlP0en5ycbH2t1Q4u1LvFzPqVK1fw+eefY/369VKXQkRERESEc5XHcLRoR4/tyZHpmDF6\nKeQy1zRVdIuwvnHjRvj5+WHVqlXXHZuZmemCisgZsrOzAfAcegOeS+/Bc+k9eC69B8+ltCzCgn8f\n3NRrUF82c5VDrRmvXh0yWJKHdSEE3n33Xdx9993w9/eXuhwiIiIiGqY6zCZs3vV/yD631267XCbH\nXQt/gVlpi1xek+Rhfc+ePbh48SI2b94sdSlERERENEwZjK1478v/xfnSE3bblT4q/HT5fzulh/pA\nSB7W58+fb/dgJCIiIiIiV9Lp6/GXz55Dec0lu+3O7qE+EJKHdSIiIiIiqVTXl+Gtbc+g/toe6tpo\n/HzlUwgPjpaosk4M60REREQ0LBVVnMXbX7yA1rZmu+2u6qE+EAzrRERERDTsnLhwBH//+lWYzO12\n29MSM/HAsl9C5esnUWX2GNaJiIiIaFjZd2I7tu55BwL2zwadPWEx7pj/cyjkCokq64lhnYiIiIiG\nBSEEvji0Cbuyt/bYt3TmKtzqQA91V2FYJyIiIiKv12E2YcuuN3Hs3B677XKZHHcteAizJiyWprDr\nYFgnIiIiIq/W1m7A3758Cecv9+yhvmbZr5CW5L5Pi2VYJyIiIiKvpdPX46+fPY+ymiK77QFqLX62\nYh0SopIlqmxgGNaJiIiIyCu5ew/1gWBYJyIiIiKvU1RxDm9/8fsePdRHRibjZyv+B4H+wRJV5hiG\ndSIiIiLyKicvHsH7X/XsoZ6amIE1y37lNj3UB4JhnYiIiIi8xv4T2/Gvve9CCIvd9llpi3HnAvfq\noT4QDOtERERE5PGEEPj3oU3Y2VsP9Rl349YZd7ldD/WBYFgnIiIiIo9m6jDhw93rcfRslt12d++h\nPhAM60RERETksep01Xhv+8sovXLRbrsn9FAfCIZ1IiIiIvJIJy9+hw92vg6DUW+33VN6qA8EwzoR\nEREReRSzuQNfHPoHdud+1mOfJ/VQHwiGdSIiIiLyGA3Ntdj41Su4VHmux75Jo2finsWPwF8VIEFl\nzsGwTkREREQe4WxJHv6+40/QG5rstsvlCvxgzmrckv59j+z40h+51AVUVlZi9erViIiIgFqtRlpa\nGvbt2yd1WURERETkJiwWM748vBl/2fZsj6AeHBCKx378e8yfusLrgjog8cx6Y2Mj5syZg3nz5mH7\n9u0IDw9HUVERIiIipCyLiIiIiNxEk74Rf//6jygoO9Vj3/iEqbj/e48jQB0kQWWuIWlYf/nllxEb\nG4uNGzdatyUkJEhXEBERERG5jQvl+dj41Sto0jfYbZfJ5Fg+cxUWTbsdcpnkC0WcStKfbtu2bZg+\nfTruuusuREZGIj09HW+++aaUJRERERGRxCzCgp3HtuKNrb/tEdQD/YPx8A+fwZLpd3h9UAcknlkv\nKirC+vXr8cQTT+DJJ59EXl4eHn30UQDAww8/LGVpRERERCQBvaEJ//jmNZwpzumxb0zcBDxw6/9D\nkCZEgsqkIRNCCKk+XKlUYvr06Thw4IB12//8z//g008/xZkzZ6zbdDqd9XVhYaFLayQiIiIi16hp\nLse+81uhNzb12Dcxbg4mj7zZo2bTk5NtD2XSarWDeg9JZ9ZjYmKQmppqt23cuHG4fPmyRBURERER\nkasJIXCu8hhyinfBIix2+1Q+atyU8gPEhoyRqDppSRrW58yZg3Pn7BvaFxQUIDExsc9jMjMznVwV\nOUt2djYAnkNvwHPpPXguvQfPpfcYbufSYNRj867/w4lLh3vsS4weizVLf4mQwHAJKrtxV68OGSxJ\nw/ratWsxe/ZsvPDCC7jzzjuRl5eHN954Ay+++KKUZRERERGRC5TVFGHDl39Aja6yx7756SuwYs5P\noFAM72d4SvrTZ2ZmYtu2bXjyySfx3HPPISEhAc8//zweeughKcsiIiIiIicSQuBw/k78a8876DCb\n7Paplf64Z/F/YfKYmRJV514k/1Vl2bJlWLZsmdRlEBEREZELGNsN+DDrL8g+t7fHvrjwUViz7FcI\nD46WoDL3JHlYJyIiIqLhobKuFO9t/19U15f12Ddn4q340byfwtdHKUFl7othnYiIiIic7ti5Pfjw\n27fQ3mG026709cPdCx5C5ribJarMvTGsExEREZHTmDrasXXvuzh0+pse+6JDR2LNsl8hakS8BJV5\nBoZ1IiIiInKKmsZKvLf9ZZTXXOqxb/r4+bhj/s+g8vWToDLPwbBOREREREPueOEhbN71f2hrb7Xb\n7qtQ4sfzH8TM1IWQyWQSVec5GNaJiIiIaMh0mE347MD72Hv83z32hQfH4KfLfoXY8CQJKvNMDOtE\nRERENCTqm2qw4as/oKSqoMe+KcmzsWrhI1Cr/CWozHMxrBMRERHRDcu/lI1/fPMaWtua7bYr5D74\n4bw1mDtpGZe9DALDOhERERENmtlixvbDm7Eze2uPfSMCw7Fm2a+QEJUiQWXegWGdiIiIiAZFp6/H\n+1/9ERfK83vsS0vKxH1LHoPGL1CCyrwHwzoREREROayg9CTe/+qPaDbo7LbLZXLcNvs+LMhYCblM\nLlF13oNhnYiIiIgGrN1kxFffbcHu3M8hhMVun1YzAg8s/X8YHZsmUXXeh2GdiIiIiAaksOw0/rnr\nTdToKnvsGxs/GT+5dS0C/YMlqMx7MawTERERUb8MRj0+P/B3HDy9o8c+GWS4dcZd+N70OyCXKySo\nzrsxrBMRERFRn04XHcOHWX+BrqWux76QgDDcs/hRjB05WYLKhgeGdSIiIiLqoblVh0/2voucgv29\n7p87aRm+P+d++CnVLq5seGFYJyIiIiIrIQRyzu/D1r3vQn/NA44AICI4BqsWPcybSF2EYZ2IiIiI\nAAANzbX4KOsvyL+U3WOfXCbHwowf4tYZd8HXRylBdcMTwzoRERHRMGcRFhw69Q0+O/g+jO2GHvtj\nw5Nwz6JHER8xSoLqhjdJw/rTTz+NZ5991m5bVFQUKioqJKqIiIiIaHi50lCBLd++iYu9PIXUR+GL\npTPuxoKpP4BCwTleKUj+T33cuHHYs2eP9XuFgi1/iIiIiJzNbDEjK/czfHXknzCZ2Line truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x81bc2e8>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"t = np.linspace(-1, 2, 20)\n",
|
||||
"plt.plot(t, np.exp(t))\n",
|
||||
"t = np.linspace(0, 1, 2)\n",
|
||||
"plt.plot([1,2, 4])\n",
|
||||
"plt.xlabel('t')\n",
|
||||
"plt.ylabel('y')\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Here we see the next estimate for y is 4. The errors are quickly getting large, and you might be unimpressed. But 1 is a very large step size. Let's put this algorithm in code, and verify that it works by trying to generate the result above."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def euler(t, tmax, y, dx, step=1.):\n",
|
||||
" ys = []\n",
|
||||
" while t < tmax:\n",
|
||||
" y = y + step*dx(t, y)\n",
|
||||
" ys.append(y)\n",
|
||||
" t +=step \n",
|
||||
" return ys"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"2.0\n",
|
||||
"4.0\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def dx(t, y): return y\n",
|
||||
"\n",
|
||||
"print(euler(0, 1, 1, dx, step=1.)[-1])\n",
|
||||
"print(euler(0, 2, 1, dx, step=1.)[-1])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This looks correct. So now lets plot the result of a much smaller step size."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAEjCAYAAADuR70GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYVOfdPvD7DDuIiCj7IggIKLiAikRBXKOJWxI1sTYx\nNk2TmCamS1rT9w0kMab2tf5Ma9psTUIWs5jFmM1gFBXEBVxQxAUEBQQGWWSRnXl+f1gmHgYEYYYz\nA/fnurzqfM+ZM7fEPnw9POd5JCGEABERERERGYxK6QBERERERP0dm24iIiIiIgNj001EREREZGBs\nuomIiIiIDIxNNxERERGRgbHpJiIiIiIyMDbdREZs1apVUKlUyM/PVzoKERER9QKbbrql8+fP4+mn\nn0ZYWBiGDBkCKysruLm5Yf78+XjjjTdQV1endMQ+9d5770GlUuGFF17os8+UJKnPPouIqDdUKlWX\nvw4cOKB0zA7t27cPKpUKDz/8sNJRqJ8yVzoAGa+XXnoJ8fHxEEIgMjISq1atgr29Pa5evYoDBw7g\n8ccfx6ZNm5Cdna101D7Xl40w968iIlMiSRLi4uI6Pe7j49OHaW4fb3SQobDppg5t2LABcXFx8PLy\nwmeffYbJkyfrnJOUlIT/+Z//USCd8tgIExF17vnnn1c6Qo9xfCdD4fQS0nHp0iXEx8fD0tIS3333\nXYcNNwDExsYiKSlJVvvyyy8RGxsLBwcH2NjYICQkBHFxcbh+/brO+6dPnw6VSoVLly5h69atCAkJ\ngY2NDXx9fbFx40bteR999BEiIiIwaNAguLi44KmnnkJjY6PO9VQqFXx9fVFVVYUnnngC7u7usLGx\nwZgxY/Cvf/1L5/y2HyV2NlWkLV+bVatWYfXq1QCAF1544ZY/Lv3iiy8wa9YsDB06FNbW1ggMDMS6\ndetQU1PT4Wf99NNPmDZtGuzs7ODk5IQlS5bg3LlzHZ5LRNQfXLt2Db6+vrCyssLRo0dlx4QQuOuu\nu6BSqWTjd1FREV588UXccccdcHV1hZWVFTw8PLBixQpkZWV1+llpaWm4//774eHhoZ0mOXPmTLz/\n/vsAgPj4eMyYMQMAkJCQIBvfExISDPCnp4GId7pJx7vvvouWlhbcf//9GDNmzC3PtbS01P7++eef\nx/r16+Hk5IQVK1ZgyJAhSExMxEsvvYSdO3ciOTkZgwYN0rnGH/7wByQnJ2PBggWYPXs2tm/fjnXr\n1qGxsRFWVlZ4+eWXsWTJEsTExOCbb77B1q1b0dLS0mEj3dTUhFmzZqGmpgYrV65EQ0MDtm/fjief\nfBIXLlzAli1bdN5zqx8l3nxsyZIlqKqqwtdff43p06dj+vTp2mM3/7j0iSeewOuvvw4vLy/ce++9\ncHR0xKFDh7Bx40Z8//33OHjwoOzr8Pnnn2P58uWwsrLCsmXL4OHhgZSUFERFRSEsLKzzLz4RkQkb\nMmQIPv30U0ybNg3Lly/HiRMnMGTIEADApk2b8MMPP+Dee+/FE088oX3PgQMHsHHjRsyYMQMTJkzA\noEGDcOHCBXzxxRfYuXMnUlJSMG7cONnnvP3223jsscdgbm6OBQsWYNSoUSgrK8OxY8fw6quv4sEH\nH0RsbCwuX76MhIQEjBs3DosXL9a+f/z48X3zBaH+TxC1M2PGDCFJkvjPf/7T7fccOnRISJIkvLy8\nRHFxsezYQw89JCRJEk8++aSsHhMTIyRJEgEBAaK0tFRbv3TpkrC0tBS2trbCxcVF5OTkaI9VVVWJ\n4cOHCysrK9l7hBBCkiQhSZKYNm2aaGpq0tbLysqEr6+vkCRJpKamautJSUlCkiTxwgsvdPhniomJ\nESqVSlZ79913b/meDz74QEiSJO69917R0NAgO7Z+/XohSZL4/e9/r63V1NSIoUOHCgsLC3H06FHZ\n+X/84x+FJElCpVKJy5cvd/h5RETGpG0cjo+PF3FxcTq/4uPjdd6zefNmIUmSWLJkiRBCiMOHDwsL\nCwvh5+cnqqqqZOeWlpaK2tpanWtkZGSIQYMGiTvvvFNWP3PmjDA3NxeOjo4iMzNT532FhYXa3+/b\nt09IkiQefvjhHv3ZibrCppt0BAcHC0mSxI8//tjt9zzyyCNCkiTx73//W+eYWq0WNjY2wt7eXjQ3\nN2vrbU13QkKCzntmzpzZaXP78MMPC0mSxIEDB2T1tgY1JSVF5z1vvvmmkCRJ/PrXv9bWDNF0T5gw\nQVhYWIjKykqdYy0tLWLYsGHC1dVVW/vwww+FJEli5cqVOudXV1cLBwcHNt1EZDLamu7OfrUfU9ss\nXLhQSJIkXnzxRTFixAhhZWWlcyOiKwsWLBDW1taipaVFW3vyySeFJEli06ZNXb6/7XsCm24yFE4v\nIb04fvw4AGjnxN3M2dkZoaGhSEtLw4ULFxASEiI73tGP7lxdXTs95u7uDgAoLCzUOWZubo6oqCid\nekxMDADg5MmTXf1Reqyurg4nTpyAk5NTh9NYgBvTcYqLi1FZWQlHR0ft160t383s7e0xfvx47N+/\n32CZiYj0TZIktLa23tZ73nvvPYwfP1676snf//53TJw4scNzv/vuO7z++utIT09HeXk5WlpaZJ9d\nVlYGFxcXAMDhw4cBAPPmzevJH4VIr9h0kw43NzecO3cOBQUF3X5PVVUVJEnSNssdXbPtvPYcHBx0\naubm5l0ea25u1jk2bNiwDudoOzs7d/r5+lJZWQkAKC8vx4svvtjpeZIkoba2Fo6Ojto8bd8g2uus\nTkTUnzg6OiI2NhYJCQmwt7fHr371qw7Pe/XVV/HMM89g6NChmD17Nry9vWFrawtJkvDVV18hIyND\n9qD9tWvXIEkSPDw8+uqPQtQpNt2kY9q0aUhKSsKePXs6Hfjaa2uOi4uLMXjwYJ3jxcXFsvMMpays\nDEIIncZbrVbrfH7byiQ33yW52bVr127rs9uuHRYW1u076m3vacvXXmd1IqL+5Msvv0RCQgKGDx+O\nq1ev4rHHHsO2bdtk57S0tCA+Ph5ubm44fvy4zk2JgwcP6lx3yJAhEEKgsLDQ4N9/iLrCJQNJx8MP\nPwwLCwt88cUXOHPmzC3PbWpqAgCEh4dDCKGzhCAAlJaWIjMzE4MGDcKoUaMMkrlNS0tLhwNv2xSN\nm6erODo6AkCHW6xXVVXhwoULOnUzMzMA6PBHp4MGDcKYMWNw9uxZlJeXdytveHg4gBvLF7ZXU1OD\nEydOcKMGIurX8vLysHr1ajg6OiItLQ3z5s3DJ598grfeekt2XllZGaqqqhAVFaXTcNfW1uL48eM6\n4+WUKVMAAD/88EOXOW41vhPpA5tu0uHj44P4+Hg0Nzfjrrvu0lk/tc3evXu1c5Hb1q/esGGD7O6s\nEAJ/+tOfUF9fj4ceekg7qLXRd0MphMBzzz2n/ccAcGOgfuWVVyBJkmx73+DgYDg4OGDHjh2yzC0t\nLVi7di0aGhp0rj9s2DAAwOXLlzv8/N///vdobm7Gww8/rJ1ucrOamhrZ13PRokVwdHTEp59+qvN1\nfvHFF1FdXd3NPzkRkelpbm7G8uXLUVNTg3feeQfe3t54//334e7ujrVr1yIzM1N7rrOzM2xtbZGe\nni7b+6G5uRlPP/10hzc7Hn/8cZibm+Pll1/u8CbSzc8GOTk5Aeh8fCfqLU4voQ6tW7cOLS0teOGF\nFxAZGYkpU6YgIiIC9vb2KC0tRXJyMs6fP4/AwEAAQGRkJNatW4dXXnkFY8aMwdKlSzF48GDs3r0b\nJ06cQFhYGF555RWdzxF63vnLzc0N9fX1CA0NxcKFC9HQ0IDPP/8carUaTz/9NCIjI7Xnmpub45ln\nnkF8fDzGjx+PxYsXQ5IkJCUlQZIkjB07FhkZGbLrR0VFwc7ODp988gksLCzg7e0NSZLw4IMPwtvb\nGw899BCOHz+Of/7znxg5ciTmzp0LHx8fXLt2DZcuXUJycjLmzp2LL7/8EgBgZ2eHN998E8uXL0dM\nTAyWL18ONzc3HDx4EJmZmYiOjtbZeIeIyJgJIfDCCy90Or7PmzdPu+nan/70J6Snp+O3v/0tFi1a\nBOBG87tt2zbMnDkTy5cvR1paGmxtbaFSqfDUU0/hr3/9q3aMb2pqQlJSEq5du9bhhm3BwcH417/+\nhcceewzh4eFYsGABAgMDUVFRgRMnTqCpqUn7QHtQUBC8vLyQnJyMlStXIiAgAGZmZli0aBFCQ0MN\n+BWjAaM7S5wUFRWJBx98UAwfPlxYW1uLkJAQsX//ftk5cXFxwt3dXdjY2Ijp06eLM2fO6HWZFVLG\nuXPnxFNPPSVCQ0OFg4ODsLCwEK6urmLu3Lni9ddfF/X19bLzt2/fLmJiYsTgwYOFlZWVCA4OFv/7\nv//b4bqq06dP73Q5vFWrVgmVSqXz90wIIeLj44VKpdJZalCSJOHr6yuqqqrE448/Ltzd3YWVlZUY\nPXq0eO211zr9M27atEkEBAQIS0tL4e7uLp544glRUVGhzdfe7t27xdSpU4W9vb12Caz2OXft2iUW\nLVokXF1dhaWlpXBxcRHh4eHi2WefFSdPnuz0mra2tmLo0KFi8eLF4vz589qvA5cMJEPg2E761jYm\n3mrJwFdffVUIIcTOnTuFJEliwoQJsr0V2rzwwgtCpVLJlvBraWkRmzdvFiEhIcLGxka4ubmJBx98\nUOTn599yvDx8+LC47777tGOyu7u7mDVrlvjwww9l5x0/flzMnj1bDBkyRKhUqg6/1xD1lCTErW81\nXrt2DRMmTEB0dDSefPJJDB8+HLm5uXBzc0NQUBAAYOPGjXj55ZeRkJCAwMBAvPjii0hJScH58+c7\n3IGQyBBUKhVGjBiB3NxcpaMQGT2O7UREfavLpvu5555DcnIykpOTOzwuhIC7uzueeLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0xa0594a8>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"ys = euler(0, 4, 1, dx, step=0.00001)\n",
|
||||
"plt.subplot(1,2,1)\n",
|
||||
"plt.title('Computed')\n",
|
||||
"plt.plot(np.linspace(0, 4, len(ys)),ys)\n",
|
||||
"plt.subplot(1,2,2)\n",
|
||||
"t = np.linspace(0, 4, 20)\n",
|
||||
"plt.title('Exact')\n",
|
||||
"plt.plot(t, np.exp(t))\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"exact answer= 54.5981500331\n",
|
||||
"euler answer= 54.59705808834125\n",
|
||||
"difference = 0.00109194480299\n",
|
||||
"iterations = 400000\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print('exact answer=', np.exp(4))\n",
|
||||
"print('euler answer=', ys[-1])\n",
|
||||
"print('difference =', np.exp(4) - ys[-1])\n",
|
||||
"print('iterations =', len(ys))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Here we see that the error is reasonably small, but it took a very large number of iterations to get three digits of precision. In practice Euler's method is too slow for most problems, and we use higher level methods.\n",
|
||||
"\n",
|
||||
"Before we go on, let's formally derive Euler's method, as it is the basis for the more advanced Runge Kutta methods used in the next section. In fact, Euler's method is the simplest form of Runge Kutta.\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Here are the first 3 terms of the Euler expansion of $y$. An infinite expansion would give an exact answer, so $O(h^4)$ denotes the error due to the finite expansion.\n",
|
||||
"\n",
|
||||
"$$y(t_0 + h) = y(t_0) + h y'(t_0) + \\frac{1}{2!}h^2 y''(t_0) + \\frac{1}{3!}h^3 y'''(t_0) + O(h^4).$$\n",
|
||||
"\n",
|
||||
"Here we can see that Euler's method is just using the first two terms of the Taylor expansion. Each subsequent term is smaller than the previous terms, so we are assured that the estimate will not be too far off from the correct value. "
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1444,7 +1751,7 @@
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"\n",
|
||||
"Runge Kutta integration is the workhorse of numerical integration. As mentioned earlier there are a vast number of methods and literature on the subject. In practice, using the Runge Kutta algorithm that I present here will solve most any problem you will face. It offers a very good balance of speed, precision, and stability, and it is used in vast amounts of scientific software. \n",
|
||||
"Runge Kutta integration is the workhorse of numerical integration. As mentioned earlier there are a vast number of methods in the literature. In practice, using the Runge Kutta algorithm that I present here will solve most any problem you will face. It offers a very good balance of speed, precision, and stability, and it is the 'go to' numerical integration method unless you have a very good reason to choose something different. If you have the knowledge to make that decision you have no need to be reading this section!\n",
|
||||
"\n",
|
||||
"Let's just dive in. We start with some differential equation\n",
|
||||
"\n",
|
||||
@@ -1460,7 +1767,8 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"$$t(t+\\Delta t) = y(t) + \\frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) + O(\\Delta t^4)$$\n",
|
||||
"Deriving these equations is outside the scope of this book, but the Runge Kutta RK4 method is defined with these equations.\n",
|
||||
"$$y(t+\\Delta t) = y(t) + \\frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) + O(\\Delta t^4)$$\n",
|
||||
"\n",
|
||||
"$$\\begin{aligned}\n",
|
||||
"k_1 &= f(y,t)\\Delta t \\\\\n",
|
||||
@@ -1473,7 +1781,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -1511,7 +1819,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -1520,7 +1828,7 @@
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuMAAAEWCAYAAADMypmTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X1c1fX9//HHOVzIpUdADkfBFA0FqZhp5kU/ykRXy2l+\ntzLNVsvNtTkXYxtNs3RlmOXXyLTW182L6Zy15bd20fqSZaXTmphdKOJF4gXgOSoiyPW5+v2hUWeI\nCAIHDs/77cYNeX/e73Nen94STz6+P++Pwe12uxERERERkXZn9HYBIiIiIiJdlcK4iIiIiIiXKIyL\niIiIiHiJwriIiIiIiJcojIuIiIiIeInCuIiIiIiIlyiMi4iIiIh4SZNh3OFwMHfuXPr3709wcDD9\n+/fnsccew+l0evRbsGABsbGxhISEMGbMGPLy8jyO19bWMnv2bKKjowkLC2PSpEkUFRW17tmIiIiI\niHQiTYbxrKwsXn75ZV544QX279/P888/z4svvsiiRYvq+yxevJilS5eyfPlydu7cidlsZty4cVRU\nVNT3SU9PZ9OmTWzcuJGtW7dSXl7OhAkTcLlcbXNmIiIiIiIdnKGpJ3B++9vfpmfPnqxevbq+7f77\n7+fMmTP87W9/w+1207t3b372s58xZ84cAGpqajCbzSxZsoSZM2dSVlaG2WxmzZo1TJ06FYDCwkL6\n9u3LP//5T8aPH9+GpygiIiIi0jE1eWX89ttv591332X//v0A5OXlsWXLFu644w4ACgoKsNlsHoE6\nKCiI1NRUtm/fDsCuXbuw2+0efeLi4khKSqrvIyIiIiLS1fg31eEnP/kJhYWFJCUl4e/vj8PhYN68\neTz00EMAWK1WAGJiYjzGmc1miouL6/v4+fkRFRXl0ScmJgabzdYqJyIiIiIi0tk0GcaXLVvG6tWr\n2bhxI8nJyezevZuHH36Yfv368eCDD15yrMFgaFYxZWVlzeovIiIiItKRmEymZvVvcpnKU089xdy5\nc7n77rtJTk5m+vTpZGRk1N/AabFYABpc4bbZbPXHLBYLTqeTkpISjz5Wq7W+j4iIiIhIV9NkGHe7\n3RiNnt2MRiNf3vcZHx+PxWIhJyen/nhNTQ3btm1j1KhRAAwdOpSAgACPPoWFheTn59f3ERERERHp\nappcpnLnnXfy9NNPEx8fz+DBg9m9ezfPPfcc999/P3B+KUp6ejpZWVkkJiaSkJDAwoULCQ8PZ9q0\nacD5y/UzZswgMzMTs9lMZGQkGRkZpKSkkJaWdtH3be4lfulYcnNzARg2bJiXK5Erpbn0HZpL36G5\n9A2aR99xJUutmwzjzz33HN27d2fWrFnYbDZ69erFzJkzefzxx+v7ZGZmUl1dzaxZsygtLWXEiBHk\n5OQQGhpa3yc7Oxt/f3+mTJlCdXU1aWlprF+/vtnrykVEREREfEWT+4y3p6//VqEr452bftv3HZpL\n36G59B2aS9+gefQdV5Jhm1wzLiIiIiIibUNhXERERETESxTGRURERES8RGFcRERERMRLFMZFRERE\nRLxEYVxERERExEsUxkVEREREvERhXERERETESxTGRURERES8RGFcRERERMRLFMZFRERERLxEYVxE\nRERExEsUxkVEREREvERhXERERETESxTGRURERES8RGFcRERERMRLFMZFRERERLxEYVxERERExEsU\nxkVEREREvKTJMN6vXz+MRmODjwkTJgDgdrtZsGABsbGxhISEMGbMGPLy8jxeo7a2ltmzZxMdHU1Y\nWBiTJk2iqKiobc5IRERERKSTaDKM79q1C6vVWv/x8ccfYzAYmDJlCgDPPPMMS5cuZfny5ezcuROz\n2cy4ceOoqKiof4309HQ2bdrExo0b2bp1K+Xl5UyYMAGXy9V2ZyYiIiIi0sE1GcajoqIwm831H//4\nxz8wmUzcfffduN1usrOzmTNnDpMnTyY5OZm1a9dy7tw5NmzYAEBZWRmrVq1iyZIljB07liFDhrBu\n3To+++wzNm/e3OYnKCIiIiLSUTVrzbjb7eb3v/8906dPp1u3bhQUFGCz2Rg/fnx9n6CgIFJTU9m+\nfTtw/sq63W736BMXF0dSUlJ9HxERERGBolMVTXcSn9KsMP72229z5MgRfvjDHwJgtVoBiImJ8ehn\nNpvrj1mtVvz8/IiKivLoExMTg81ma3HhIiIiIr4k78gZBt3jZPKcz6mqsXu7HGkn/s3pvHLlSoYP\nH861117bZF+DwdDiogByc3OvaLx0DJpH36G59B2aS9+hufQNubm52J0upj9roqpmIG98cA0Dpxxh\n8fcLGRQX7O3y5DIkJCS0eOxlXxk/efIkf/3rX+uvigNYLBaABle4bTZb/TGLxYLT6aSkpMSjj9Vq\nre8jIiIi0pU9u8lOwYmB9V8Xn+7HXz+6sgub0jlc9pXxNWvWEBQUxNSpU+vb4uPjsVgs5OTkMHTo\nUABqamrYtm0bS5YsAWDo0KEEBASQk5NTP7awsJD8/HxGjRrV6PsNGzasRSckHcOXV2s0j52f5tJ3\naC59h+bSN3w5j0fO9eSNf8V5HBt41Rf8ceEo/P30SJjOoKysrMVjLyuMu91ufve733HPPfcQEhJS\n324wGEhPTycrK4vExEQSEhJYuHAh4eHhTJs2DQCTycSMGTPIzMzEbDYTGRlJRkYGKSkppKWltbhw\nERERkc6u5JydHz4ditvtV98W3K2cNxZHKYh3EZcVxt977z2++OKL+u0Kvy4zM5Pq6mpmzZpFaWkp\nI0aMICcnh9DQ0Po+2dnZ+Pv7M2XKFKqrq0lLS2P9+vVXvK5cREREpLNyudw8sjqUsoqeHu1LZtsY\ndNXARkaJr7msMD5mzBicTmejx+fPn8/8+fMbPR4YGMiyZctYtmxZ8ysUERER8UErc6r47AvPpUYT\nRu/hx5Ob3ihDfIf+/UNERESknb2/+wRrc270aOvds5ANCxK9VJF4i8K4iIiISDs6U17DlMfcOJzd\n6tsC/Gv481MBhIUEeLEy8QaFcREREZF29F9zvuBkaS+PtkfuK2DkNdryuStSGBcRERFpJ4//Lo8P\nPhns0Tbqun0seDDJSxWJtymMi4iIiLSD93ef4Ok/xHu0RZlO8MbT8RiN2mGuq1IYFxEREWljX60T\nD6pv8/erJev+Q0SZ9Mj7rkxhXERERKSNTb7IOvHv3/YRQwaENDJCugqFcREREZE29PjKvWy9yDrx\nGWkK4qIwLiIiItJm3t1VzNPr+nu09exh1TpxqacwLiIiItIGbGeqmPKYscE68Y1PuLVOXOopjIuI\niIi0MpfLzR2/PEZJWYxH+6/u/YJbh/b2UlXSESmMi4iIiLSy2Uv38vH+QR5tNw/J46kfJXupIumo\nFMZFREREWtEr7xTw29cTPdp69Szk9cVXe6ki6cgUxkVERERayaHCMn64qAdut199W7fASt5YHIQp\ntJsXK5OOSmFcREREpBXU1jn51i9OU1Hdw6P9mVlFDEuM9lJV0tEpjIuIiIi0gnt/k8ehQs9tDL87\nZg+zvzuokREiCuMiIiIiV+zZDflseu8aj7aBfQ6z7vHBjYwQOU9hXEREROQKvLe7mHkv9/VoCwsp\n5c3/7km3QL9GRomcpzAuIiIi0kLFJRV8d64fdsdXD/YxGh2sfrSc/rEmL1YmncVlhfETJ05w//33\nYzabCQ4OJjk5mQ8++MCjz4IFC4iNjSUkJIQxY8aQl5fncby2tpbZs2cTHR1NWFgYkyZNoqioqPXO\nRERERKQdOZwuvplezJlys0f7L6bu5zu39PNOUdLpNBnGz549y+jRozEYDLz55pvk5+ezfPlyzOav\n/uItXryYpUuXsnz5cnbu3InZbGbcuHFUVFTU90lPT2fTpk1s3LiRrVu3Ul5ezoQJE3C5XG1zZiIi\nIiJt6L4n9rL3cIJH2y3X57H4J9c0MkKkIf+mOjzzzDPExsayZs2a+ra+fb9aF+V2u8nOzmbOnDlM\nnjwZgLVr12I2m9mwYQMzZ86krKyMVatWsWbNGsaOHQvAunXr6Nu3L5s3b2b8+PGtfFoiIiIibef5\nV/fzymbP0H1VzDH+uniglyqSzqrJK+Ovv/46w4cPZ8qUKcTExDBkyBBWrFhRf7ygoACbzeYRqIOC\ngkhNTWX79u0A7Nq1C7vd7tEnLi6OpKSk+j4iIiIincG2T6088mKcR1toUBlvLg0nLCTAS1VJZ9Vk\nGD98+DAvvvgiV199NTk5OTz88MP8+te/rg/kVqsVgJiYGI9xZrO5/pjVasXPz4+oqCiPPjExMdhs\ntlY5EREREZG2VnSqgslzoM4eUt9mNDhYOaeUwf0ivViZdFZNLlNxuVwMHz6cp556CoCUlBQOHjzI\nihUrmDVr1iXHGgyGFheWm5vb4rHScWgefYfm0ndoLn2H5rJ91TlcfP+5cErKEj3a703bwdU9QsjN\nLWnR62oeO7+EhISmOzWiySvjvXv3ZvBgzw3rExMTOXbsGAAWiwWgwRVum81Wf8xiseB0Oikp8fxL\narVa6/uIiIiIdGTz/+jmYKFnEB+euJtZ3wr2UkXiC5q8Mj569Gjy8/M92g4cOEC/fv0AiI+Px2Kx\nkJOTw9ChQwGoqalh27ZtLFmyBIChQ4cSEBBATk4OU6dOBaCwsJD8/HxGjRp10fcdNmxYi09KvO/L\n3/I1j52f5tJ3aC59h+ay/T2+ci/vfOx5cbKv5SibX7i2xevENY++o6ysrMVjmwzjP//5zxk1ahRZ\nWVncfffd7N69mxdeeIFFixYB55eipKenk5WVRWJiIgkJCSxcuJDw8HCmTZsGgMlkYsaMGWRmZmI2\nm4mMjCQjI4OUlBTS0tJaXLyIiIhIW3vtvSMs+oPnMoTwkDO89ZxJN2zKFWsyjA8bNLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x91fbef0>"
|
||||
"<matplotlib.figure.Figure at 0xaa29160>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
@@ -1541,8 +1849,7 @@
|
||||
"y = 1.\n",
|
||||
"dt = .1\n",
|
||||
"\n",
|
||||
"ys = []\n",
|
||||
"ts = []\n",
|
||||
"ys, ts = [], []\n",
|
||||
"\n",
|
||||
"def func(y,t):\n",
|
||||
" return t*math.sqrt(y)\n",
|
||||
@@ -1583,7 +1890,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"execution_count": 22,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -1592,7 +1899,7 @@
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAp8AAAGNCAYAAABXHhoFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VOeZPv77nDNdM2qAOqIJJDoIITBgOjbYuAC2sbEd\n1zhOvN6UzSbxJrv2btrP+WYdO5s41XGJe8EVbFMsU01HVAEC1JBQ76Op55zfH4oJGFRm5syMZub+\nXJcuG8153/NoGGYePW8TVFVVQUREREQUAmK4AyAiIiKi2MHkk4iIiIhChsknEREREYUMk08iIiIi\nChkmn0REREQUMkw+iYiIiChkmHwSERERUchomnxu3boVN954I7KysiCKIl588cVLHr/33nshiuIl\nX7NmzdIyBCIiIiIawDRNPu12OyZNmoRnnnkGZrMZgiBc8rggCFiyZAlqa2svfK1fv17LEIiIiIho\nANNp2dmyZcuwbNkyAN1Vzq9SVRUGgwEpKSla3paIiIiIIkRI53wKgoDt27cjNTUVubm5eOihh9DQ\n0BDKEIiIiIgojDStfPZl6dKlWLVqFUaMGIGysjL85Cc/wcKFC7F//34YDIYL17W1tYUyLCIiIiIK\nkoSEhEv+HNLkc/Xq1Rf+f/z48Zg2bRqGDRuGdevWYcWKFaEMhYiIiIjCIKxbLaWnpyMrKwunT58O\nZxhEREREFCIhrXx+VUNDA6qrq5Gent7jNV8t1Uazffv2AQAKCgrCHMnAxOend3x+esfnp2d8bnrH\n56d3fH56FsvPTW9TKDVNPu12O0pLSwEAiqKgoqICxcXFGDRoEJKTk/H444/jlltuQVpaGsrLy/HY\nY48hNTWVQ+5EREREMULTYfe9e/ciPz8f+fn5cDqdePzxx5Gfn4/HH38ckiTh6NGjuOmmm5Cbm4t7\n770XY8eOxRdffIG4uDgtwyAiIiKiAUrTyuf8+fOhKEqPj3/yySda3o6IiIiIIgzPdiciIiKikGHy\nSUREREQhw+STiIiIiEKGyScRERERhQyTTyIiIiIKGSafRERERBQyTD6JiIiIKGSYfBIRERFRyDD5\nJCIiIqKQYfJJRERERCHD5JOIiIiIQobJJxERERGFDJNPIiIiIgoZJp9EREREFDJMPomIiIgoZJh8\nEhEREVHI6MIdABFRoGRZRru9Da0dLWhrb0VbZ/eXx+OBoipQFBlnz56FIIo423gSkihBp9Mh3pqA\nBGsiEm2JSLAlIcGWCJ3Et0UiomDiuywRRYwuhx0V58tRUVOGipqzqGuqQ1tHC9o726Cqaq9tW9va\nAACna0p6vc5qsSEpPhmDk1MwPGMEstOHY1jmSMTHxWv1YxARxTQmn0Q0IHllL85UlqKs+gwqas6i\nsqYcdU21Qb9vZ1cHOrs6UFVbgYPH9174fnLiYAzLGIHhGSMwPHMkRg/Lg0FvCHo8RETRhsknEQ0Y\nXQ47jpYewqGTB3DkVDG6nF3hDumC5tZGNLc2XkhIDXojxudMxOS8aZg0ZgrirQlhjpCIKDIw+SSi\nsGpoqcehE/tx6MRBnCovgazI4Q6pX9weFw6W7MPBkn0QBAEjs3IwOW8aJuflI2NIJgRBCHeIREQD\nEpNPIgo5r+zFwZJ9+HzPJpwsOx7ucAKmqirOVJXiTFUp1m58HcMyRmB+4RIUTrwKRoMx3OEREQ0o\nTD6JKGSa25qwdd9n2La/CG0dreEOJ2gqasrw4nt/xlufvoLZU+dh3vRFSBucHu6wiIgGBCafRBRU\nqqqi5OxRFO3ZhEMn9kNRlHCHFDJdDjs27lyPjTvXY9yoiZhfuBhT8qZBFLnFMhHFLiafRBQ0J8tK\nsHbj6zhTVRruUMLu+JkjOH7mCDJSsrBy8WpMzsvnvFAiiklMPolIc5Xny/Hupjdx5FRxuEMZcGrq\nz+F3r/4vcrLHYOWS2zFmeF64QyIiCikmn0SkmYbmOry3+S3sPrwz3KEMeKcrT+FXz/0PJuVOxcrF\nq5GVlh3ukIiIQoLJJxEFzOly4N1Nb2LL3s3wyt5whxNRDp88iCOnijFz8hzceu0a7hdKRFGPyScR\nBaTkzFG88N5f0NTaEO5QIpaqqviieBuOnCrGnTfch+kTZoY7JCKioGHySUR+cboceOvT17Bl76Zw\nhwKg+8ShRFsi4m2JSIpPQoI1Eba4eEiSDpIo4syZs1BVFSNHjoRX9qCzqxPtna1o7WhFa0cL2jpa\n4XQ5wvozdHZ14E9v/Bb7ju7GncvvZRWUiKISk08i8lk4q506SYeM1KEYnjECwzJGIDtjOFKT02A2\nWXpdPZ5k2AcAKCgo6PEap8uBhuZ6VJwvv3CefGVtBTwet9Y/Rq/2H9uNk2XHWQUloqjE5JOI+s3l\nduLNT14NabVTJ+mQO2IsJuXmY+TQHGSlZkOv0wflXiajGUPTh2Fo+jDMyZ8HAJBlGecbqnH23Gkc\nOVWMY6ePwO1xBeX+F7u4Cnr3jffDarEF/Z5ERKHA5JOI+qWptRG/f/UpVJ4vD/q94ixWTBozFZPz\n8jF+1ESYTZag37MnkiQhKy0bWWnZmFuwEG6PGyVnj+LQyYM4fOIAWjtagnr//cd2o/J8GR6543tc\nEU9EUYHJJxH1qbTiJJ597TfosLcH7R6SKGHquOmYV7AQY4aPhSRJQbtXIAx6Aybn5mNybj6U5feh\nvPoMtu4rwu4jO4M2PN/QXI9f/uUJPHjLtzB1bM/TBoiIIgGTTyLq1bZ9RXjlo+eDtoVSUnwy5hYs\nxNXTFiAxPiko9wgWURQxcuhojBw6Grdcewe+KN6Gz/dsQl1Treb3crmd+P2rT+HmRbfi+nk383Qk\nIopYTD6J6Iq8shdvffIKNu/6NCj95wzLxZKrlmFyXj50UuS/FVktNiyZdR0WzVyKkrNHsemLT4Jy\nwtN7m99CdX0V7r35IRgNJs37JyIKtsh/xycizXU57PjD68+g5OxRzfsemjYMK5bchomjp0Rl9U4U\nRYzPmYTxOZNwuvIU3t7wGk5XnNT0HnuP7EJdYy3+9a5/j7hqMRERk08iukRnVweefulJlFef1bTf\nwUlDcNPCWzBj0myIoqhp3wNVTvYY/PCB/8LhUwexduMbqK6r0qzvyvPl+NVz/4N/u+/HGJQ4WLN+\niYiCjcknEV3Qbm/Hb174JapqKzTr02Ky4KZFt2Le9EVRMbzuK0EQMDk3HxNHT8HuwzuwduMbaGlv\n1qTv+uY6/Oq5n+L79/0HhiSnatInEVGwxUb5gYj61G5vx/8+/3NNE89JuVPx34/+CotmXhuTiefF\nRFHEVVOuxhP/8v9h1tS5mvXb1NqAXz33U9QHYZETEVEwMPkkInR2deA3L/xSs2Fhi8mC+1Y+jEfv\n/D6S4pM16TNaxJmtuH/lw3j0Lu2em5b2Zvz6hV+gsSX0J04REfmKySdRjOty2PHM33+lWcXzy2rn\n7Klzo3JBkVYm5+ZrWgVtbm3EUy/8QrMhfSKiYGHySRTDvLIXz77+NMrOnQm4L0mUcOcN97Ha6YMv\nq6DfvP07MOiNAfdX31yHp196Eg5nlwbREREFB5NPohj2xsd/x4mzxwLux2qx4Xv3/gcWFC5htdMP\n08YX4rGvP4FBiUMC7qu6rgp/fftZKIqiQWRERNpj8kkUo7bs3Yyi3RsD7icrLRs/fvinyB0xVoOo\nYtfQ9GH48cM/xZjhgT+Ph04ewAefva1BVERE2mPySRSDTpaV4NWPXgi4n6njpuNHDz6OIUkpgQdF\niI+Lx/fufQzzpi8OuK+PtryHPUe+0CAqIiJtMfkkijENLfX4w+tPQ1bkgPpZMus6fHP1t2EymjWK\njABAJ+lw9433Y/WyuwLu64V3/4yKmjINoiIi0g6TT6IY4nQ58Oyrv0FnV0dA/SyfvwK3Lb0zZk4q\nCocls67D3Tc+EFAfbo8Lv3v1KbR3tmkUFRFR4PjJQRRDXl//94C3VLp50a24edGtXFgUAvOmL8J9\nKx8O6LluaWvCc+88C1VVNYyMiMh/TD6JYsSRU8XYfuDzgPpYPn8Fls9foU1A1C+zp84NuAJ67PQR\nbNtfpFFERESBYfJJFAO6HHa89P5fA+rj2jnLcdPCWzSKiHwxt2Ah7rj+noD6eOvTV9DU2qhRRERE\n/mPySRQD3vjk5YBOvpk2fgZWLbmdQ+1htGjmtVgy6zq/2zucDrz0/l84/E5EYcfkkyjKHTlVjB0H\ntvjdfmjaMNy/8htcXDQA3HLNHRifM9Hv9hx+J6KBgJ8mRFEs0OF2W1w8HlnzPRgNJg2jIn9JkoSH\nbnsUqYPS/O6Dw+9EFG5MPomi2NsbX/d7uF0n6fCtO76LwUmBH/lI2okzW/HImn+D2eTf/qoOpwMv\nf/g3jaMiIuo/Jp9EUaqmvhrbAxhiveP6ezB6WK6GEZFWMlIy8fVb/8XvObhHThWj5MxRjaMiIuof\nJp9EUerdzW9CURS/2uaPK8TcgoUaR0RamjRmKhbNvNbv9ms3vcHFR0QUFkw+iaLQ2apSHDy+16+2\nVosNd91wH1e2R4AVi1cjJTnVr7Zl585g//E9GkdERNQ3Jp9EUUZVVbyz8Q2/26+5/h7EWxM0jIiC\nxWgw4t4V3/D7F4V3N70JWZY1joqIqHdMPomizNHSQzhZdtyvtvnjCjF94lUaR0TBNGZ4nt/D73WN\n57HjoP/bcBER+YPJJ1EUURQFa/2senK4PXKtWLza7+2XPixaC4/XrXFEREQ90zT53Lp1K2688UZk\nZWVBFEW8+OKLl13zxBNPIDMzExaLBQsWLMDx4/5VaIjocqfPnURVbYVfbVcuWc3h9ghlNBj9Pn6z\npb0Zh08f1DgiIqKeaZp82u12TJo0Cc888wzMZvNlFZQnn3wSTz31FH73u99h7969SLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0xa293b00>"
|
||||
"<matplotlib.figure.Figure at 0xa971630>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
@@ -1806,7 +2113,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"execution_count": 23,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -1939,7 +2246,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"execution_count": 24,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -1983,7 +2290,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"execution_count": 25,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
|
||||
Reference in new issue
Block a user