{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Euler's Method\n",
"\n",
"### Modules - Ordinary Differential Equations\n",
"\n",
"By Magnus A. Gjennestad, Vegard Hagen, Aksel Kvaal, Morten Vassvik, Trygve B. Wiig and Peter Berg\n",
"\n",
"Last edited: March 11th 2018 \n",
"___\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Question:__\n",
"\n",
"How can we solve a first-order differential equation of the form\n",
"$$\n",
"\\frac{d}{dt}x(t)=g(x(t),t),\n",
"$$\n",
"with the initial condition $x(t_0)=x_0$, if we cannot solve it analytically?\n",
"\n",
"__Example 1:__\n",
"\n",
"We want to solve the ordinary differential equation (ODE)\n",
"\\begin{equation}\n",
"\\frac{d}{dt}x(t)=\\cos(x(t))+\\sin(t)\\qquad\\qquad\\qquad\\qquad\n",
"\\label{eq:1}\n",
"\\end{equation}\n",
"with $x(0)=0$, i.e. we need to find the right function $x(t)$ that fulfills the ODE\n",
"and the initial condition (IC)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given the initial condition $x(0)=0$, we want to know $x(t)$ for $t>0$.\n",
"We will now find an approximate numerical solution of the exact solution by computing\n",
"the values of the function only at discrete values of $t$.\n",
"\n",
"To do so, we define a discrete set of $t$-values, called grid points, by\n",
"\n",
"$$\n",
"t_n=t_0+n\\cdot h~~~~~\\mathrm{with}~~~~n=0,1,2,3,...,N.\n",
"$$\n",
"\n",
"The distance between two adjacent grid points is $h$. The largest value is $t_N=t_0+N*h$. \n",
"Depending on the problem, $t_N$ might be given and $h$ is then determined by how\n",
"many grid points $N$ we choose\n",
"\n",
"$$\n",
"h=\\frac{t_N-t_0}{N}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The key is now to approximate the derivative of $x(t)$ at a point $t_n$ by\n",
"\n",
"\\begin{equation}\n",
"\\frac{dx}{dt}_{t=t_n}\\approx \\frac{x(t_{n+1})-x(t_n)}{h},~~~~~h>0.\n",
"\\label{eq:2}\n",
"\\end{equation}\n",
"\n",
"We know that this relation is exact in the limit $h\\to 0$, since $x(t)$ is differentiable\n",
"according to equation \\eqref{eq:1}. For $h>0$, however, equation \\eqref{eq:2} is only\n",
"an approximation that takes into account the current value of $x(t)$ and the value \n",
"at the next (forward) grid point. Hence, this method is called a \n",
"__forward difference__ approximation. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In equation \\eqref{eq:2}, we approximate the slope of the tangent line at $t_n$ (\"the derivative\") by the slope of the chord that connects the point $(t_n,x(t_n))$ with the point $(t_{n+1},x(t_{n+1}))$. This is illustrated in the figure below: blue - graph;\n",
"dotted - tangent line; green - chord.\n",
"\n",
"![picture](https://www.numfys.net/media/notebooks/images/eulers_method.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Substituting the approximation \\eqref{eq:2} into \\eqref{eq:1}, we obtain\n",
"\n",
"\\begin{equation*}\n",
"\\frac{x(t_{n+1})-x(t_n)}{h} \\approx \\cos(x(t_n))+\\sin(t_n).\n",
"\\end{equation*}\n",
"\n",
"Rearranging the equation, using the notation $x_n=x(t_n)$ and writing this as \n",
"an equality (rather than an approximation) yields\n",
"\n",
"$$\n",
"x_{n+1} = x_n + h\\left[ \\cos(x_n)+\\sin(t_n)\\right].\n",
"$$\n",
"\n",
"This describes an iterative method to compute the values of the function successively\n",
"at all grid points $t_n$ (with $t_n>0$), starting at $t_0=0$ and $x_0=0$ in our case.\n",
"This is called __Euler's method__."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For example, the value of $x$ at the next grid point, $t_1=h$, \n",
"after the starting point is\n",
"\n",
"\\begin{eqnarray*}\n",
"x_{1} &=& x_0 + h\\left[ \\cos(x_0)+\\sin(t_0)\\right] \\\\\n",
" &=& 0 + h\\left[ \\cos(0) +\\sin(0) \\right] \\\\\n",
" &=& h.\n",
"\\end{eqnarray*}\n",
"\n",
"Similarly, we find at $t_2=2h$\n",
"\n",
"\\begin{eqnarray*}\n",
"x_{2} &=& x_1 + h\\left[ \\cos(x_1)+\\sin(t_1)\\right] \\\\\n",
" &=& h + h\\left[ \\cos(h) +\\sin(h) \\right] .\n",
"\\end{eqnarray*}\n",
"\n",
"It is now a matter of what value to choose for $h$. To look at this, we will write some code which uses Euler's method to calculate $x(t)$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# Set common figure parameters\n",
"newparams = {'figure.figsize': (16, 6), 'axes.grid': True,\n",
" 'lines.linewidth': 1.5, 'lines.markersize': 10,\n",
" 'font.size': 14}\n",
"plt.rcParams.update(newparams)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x_N = 1.742849\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA78AAAF/CAYAAACBl34HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XdY1XXDx/H3lw2CiIoDFbc4cCFijrRpNjUrLXPbUNu77u6627u8G2bDmZo2bGppaWpuFPcG3AMBB0M4jHN+zx9Sj09P5Qj4cc75vK6L6wrO79jHrjjnfH7fZSzLQkRERERERMST+dgdQERERERERKSsqfyKiIiIiIiIx1P5FREREREREY+n8isiIiIiIiIeT+VXREREREREPJ7Kr4iIiIiIiHg8lV8RERERERHxeCq/IiIiIiIi4vFUfkVERERERMTjqfyKiIiIiIiIx/OzO0BZq169utWgQQO7Y4iIiIiIiEgZSEpKyrQsK/JM13l8+W3QoAFr1qyxO4aIiIiIiIiUAWPM3rO5TtOeRURERERExOOp/IqIiIiIiIjHU/kVERERERERj6fyKyIiIiIiIh5P5VdEREREREQ8nsqviIiIiIiIeDyVXxEREREREfF4Kr8iIiIiIiLi8VR+RURERERExOOp/IqIiIiIiIjHU/kVERERERERj+dndwAREREREfFcRU4X2flF5BYUk+M49WVZFlbJ4waoFOhHaJAfYYF+VA72J8jf187I4qFUfkVERERE5B/JKywmJT2XHWk5JKfncuB4HgdPODh8Ip+M3AIs68x/xunCg/2pHR5EVJVg6kYE07RGKM1qhtGsZhgRlQLK5i8hHk/lV0REREREzprTZbHzSA5Je4+zdu9x1u47zt5jeb8X3AA/H+pGBFOnSjAxMZHUDg8mIsSfsCB/woL8CA30w9fHAGCMwWVZ5BWeGhHOLSjmRF4Rh7PyOXzCweEsB4m7j5FbUPz7v79W5SA61I8grn4EHepH0LJ2ZQL8tJpTzkzlV0RERERE/lZaloNFO9JZtCODZamZ5DhOldHqoYHERVfh+vZ1iakVRkytMKKrhvxebkuDZVkcznKw80gOO4/ksPlgNkl7jzNn02EAKgX40qVJdS6KieSimBrUqRJcav9u8SwqvyIiIiIi8v+kZuQye8Nhftx8mO1pOcCpUderYmvTqVFVOtSPILpqCMaUXtH9M8YYoqoEE1UlmItiavz+8yPZDpL2HmdZSiaLdmTw89YjADSvFca1baO4pk1t6lerVKbZxL0Y61wn4LuZ+Ph4a82aNXbHEBERERGp8A6eyOebdQeZvfEw2w5nYwx0bFCVS5vX4KKYGjSrGVrmZfd8WJZFakYuC7dnMHdLGkl7jwPQuk44fdrXoW/7Olor7MGMMUmWZcWf8TqVXxERERER71XkdPHL9nRmJO5j8c4MLAvioqtwTZsorm5Tm5qVg+yOeM4OnshnzsZDfL/hMJsOZhHg58NVsbW4JSGahIZVK2SBl/On8ltC5VdERERE5P9Lz3EwdcVeZq7eT0ZOAbUqB9Evvi43xdejXtUQu+OVmm2Hs5mZuI+v1h0kx1FMs5qh3H5hI3q3q6ONsjyEym8JlV8RERERkf+VfCSH8Ut28/W6gxS5XFwSU4NbEqK5KCYSP1/PLYP5hU6+33iIiUt3sz0th5qVAxnWtSEDOkVTOcjf7njyD6j8llD5FRERERGBpL3Hee+XZBbuyCDI34ebOtRjeLeGNKzuXZtCWZbFr8mZfPRrKstSjhIW5MftFzZiWNcGhKkEuyWV3xIqvyIiIiLizTbsP8GY+TtZtCODqpUCGNK5AYM616eqNoBi88Es3lmQzE9bjxAR4s/IHo0Z3LkBwQG+dkeTc6DyW0LlV0RERES80ZZDWYz5eSfzt6VTJcSfO7s3ZnDn+lQK1Gmnf7Rh/wne+nkni3dmUD00kEeuaMaNHeqV6nnFUnZUfkuo/IqIiIiIN0nLcvD6vB18te4AYYF+3NG9EUO6aErv2Vi95xgv/7CNtftO0CqqMv+5thUJDavaHUvOQOW3hMqviIiIiHiDvMJiPvp1Fx8u3oXTZTGsawNGX9yE8GCV3nNhWRbfbTjEqz9u51CWg6tb1+aJq5pTN8JzdsD2NGdbfjXnQURERETEjf1W1l7+YTtp2Q6ual2Lx3u1ILqaytr5MMbQu10derasxUe/7mLc4hQW7kjnoZ4xDO3SQFOh3ZhGfkVERERE3NSujFye+nYzy1KO0rpOOE9f25KODTRNtzQdOJ7H099u4Zft6bSuE87LfVsTWyfc7lhyGk17LqHyKyIiIiKexlHk5IPFqby/MJVAPx8e7RXDgE71NSpZRizLYs6mwzzz3VaO5xUyoltDHry8GUH+2hW6ItC0ZxFxK5ZlceB4PnuOnmTv0Tz2Hcvj4PF8jp0s5Hjeqa+TBU6KXS5cLih2ufDz9aFSgC8hAX6EBPhSPTSQWuFB1KwcRFSVIBpHhtK0ZiiRoYEYow8DIiLiGVbtOsoTX21iV+ZJrm0bxVNXt6BG5SC7Y3k0YwzXtIniwqaRvPLjdj76dRcLt6fzVr92tK6rUWB3oZFfEbHFsZOFrNx1lPX7T7D5YBabD2aR7Sj+/fEAPx/qVgmmWmgAVUICiAjxJzTQH39fg4+PwdcYilwu8gudnCxwcrKgmIzcAtKyHKTnOChy/u9rW3iwPzG1wmgfXYUO0RHE1Y+gemigHX9tERGR85Zf6OS1eduZtGwP0VVDeKFPLN2bRdodyyv9ujODR7/cSGZuAfde2pTRFzXGz9fH7lheS9OeS6j8ilQMRU4XK3cdZeH2DFbsOsq2w9kABPj60Lx2GK2iwmkVVZkmNUKpXy2EmmFB+Jzn1C2XyyIjt4CU9FySj+SQnJ7L5kPZbD2U9XspbhxZiR7NanBRTCQJDatq2pKIiFRoSXuP8fAXG9mdeZIhnevz2JXNCQnQJE47ZeUV8fR3m/l2/SHa1qvCmH5taRQZancsr6TyW0LlV8Q+jiIni3ZkMG9LGgu2HSHbUUygnw/xDSLo3KganRtXp3WdcAL8yudOqaPIyeaDWazZe5zlqUdZuesohcUugvx96N40kmvaRnFZixr6MCEiIhWGo8jJWz/v5OMlu6hTJZjXbmxDl8bV7Y4lp5m98RD//mYzhcUuXrq+NX3a17E7ktdR+S2h8itSvizLYsuhbD5fs59v1x8iK7+IKiH+XNq8Jle0qkn3ZpEVZpQ1v9DJyt1HWbQ9nblb0jiSXUCwvy+XtqhB37g69GhWQxuHiIiIbZKP5HDPjHVsT8thQKdo/nVVC0IDdYO2Ijqclc+9M9axes9x+sfX45nrWhEcUDE+73gDld8SKr8i5cNR5OTb9QeZvHwv2w5nE+DnwxWtanFjh7p0aVwN/wq+Dsblsli95xjfbzzED5vSOHaykKjwIPp3jKZ/x3rUCtdGIiIiUj4sy2Lm6v08+/0WKgX48Ua/tlwcU8PuWHIGxU4XY+bv5P1FqTStEcrYAXE0rRlmdyyvoPJbQuVXpGxl5BQwbeVepq3cy9GThTSvFcatF9TnujZRhIf42x3vvBQ5XczfeoRPE/exJDkTXx9Dz5Y1ubNHY9rVq2J3PBER8WBZeUU88fVGftiUxoVNq/Nmv7bUCNMNWHfy684MHvhsPXmFTl7uq2nQ5UHlt4TKr0jZSM9x8MGiXUxftZeCYheXNq/BiG4N6dy4mkcdK7T36Ek+TdzHjFX7yHYU06lhVUb2aMxFMZEe9fcUERH7Je09xr0z1nMk28EjV8Rw+4WNznvzR7HXkWwH98xYR+LuY9zWrSGPX9lcu0GXIZXfEiq/IqUrM7eADxenMnXlXoqcFte3r8OoixrT2MN3N8wtKGZm4j4mLN3N4SwHzWuF8VDPGC5rUUMlWERE/hHLspiyfA8vzNlGVJVg3rmlvWYaeYAip4sXZm9lyoq9dG1SjXdviaNqpQC7Y3kkld8SKr8ipcNR5GTC0t2MXZiCo8hJn/Z1uPeSpjSoXsnuaOWqyOniu/WHeG9hCrszT9K2XhUe6RlD1yaeNeItIiLlI6+wmMdnbeK7DYe4rEUN3uzXjvBg91w2JH/uizX7efKbzUSGBvLR4A60igq3O5LHUfktofIr8s9YlsXsjYd55cftHDyRT69WtXi0V4zXn2NX7HQxa+0B3p6fzKEsB50bVePJq1sQW0dvaCIicnZ2ZeQyatpadqbn8HDPGEb1aKxpzh5qw/4T3Dk1iRP5hbx5UzuublPb7kgeReW3hMqvyPnbkZbDk19vYs3e47SoXZmnr2lJ58bV7I5VoRQUO5mxah/v/JLC8bxCbu5Yj4d6xlA9NNDuaCIiUoH9tCWNhz7fgJ+v4e2b29O9WaTdkaSMZeQUMHJaEkl7j/Nor1M3OzRrrHSo/JZQ+RU5d44iJ+/9ksIHi1MJC/LjsV7NuSm+ns68/RtZ+UW8syCZKcv3EOzvy72XNmVIlwYE+GlzCxER+V+WZfHOghTGzN9Jm7rhvH9rHHUjQuyOJeXEUeTk0S838t2GQ/SPr8cL18dW+OMg3YHKbwmVX5Fzs3LXUf711SZ2ZZ6kb/s6/Pualtqc4RykZuTywuytLNyRQUzNMF6+oTVx0RF2xxIRkQrAUeTkkS838v2GQ/RtX4eX+rYmyN/X7lhSzizLYszPO3nnlxS6NqnG+7d20Drvf0jlt4TKr8jZcRQ5eeXH7Uxevod6VYN56frWXNhUU7DO14JtR3jqm80cznYw6IL6PHJFDGFBemMTEfFWR7Id3PHJGjYezOLRK5ozskcjTXn1crOSDvD4VxupX60Sk4Z2pF5VzQA4Xyq/JVR+Rc5s88EsHvhsPcnpuQzt0oDHejUnOEB3ov+p3IJi3pi3gykr9lAzLIjn+8RyecuadscSEZFytulAFrd9spocRzFv39xe7wXyu5W7jnLn1CT8fX2YMryjdoI+Tyq/JVR+Rf6a02Xx8ZJdvPnTDiJCAnj9prb00IYbpW79/hM8Pmsj29NyuCGuLv+5riWVNQosIuIVfth0mAc/X0+1SoGMHxJPi9qV7Y4kFUxKei6DJ6wix1HMR4PjtbnoeVD5LaHyK/LnjuYWcP9n61mSnEmvVrV4uW9rIrS2t8wUOV28uyCZ9xamUDs8mDduaqs3NxERDzd+yS5emLONDvUj+HBQB50EIH/pcFY+gycksvdYHu/c3I5esToK6VycbfnV1mIiXihp73GueXcpq3Yf4+W+rRk3ME7Ft4z5+/rwYM8YvhzVhQA/HwaMX8kLs7fiKHLaHU1EREqZ02Xx7PdbeGHONq6MrcX02zqp+Mrfqh0ezBcjOxMbVZnR09fy6ap9dkfySOVefo0xTxhjVhtjso0xGcaY740xsWfxvNbGmMXGmHxjzEFjzNNGuwSInBPLspi0bDf9P1yBn6/hq1FduCUhWhtulKO46Ajm3NuNgZ3qM37pbvqMXUZqRq7dsUREpJQ4ipzc/elaJi3bw/CuDRk7IE47OstZqRISwPTbLqBHs0j+9fUm3lmQjKfP0i1vdoz8XgS8D3QBLgGKgfnGmKp/9QRjTGXgZ+AI0BG4D3gEeLCsw4p4CkeRk/tmrufZ77dyUUwNZt99IbF1tKmCHUIC/Hi+TyyThnYkPaeAa99dyjfrDtodS0RE/qETeYUMHL+KuVvS+PfVLXj62pb4+OgGs5y94ABfPhocT9+4Orz1805emLNNBbgU2b7m1xgTCmQBfSzL+v4vrhkFvArUtCwrv+Rn/wZGAXWtv/lLaM2vCKRnO7i95HiFh3vGMKpHY70ZVxBpWQ7unbGOxD3H6B9fj2eua6WdtkVE3ND+Y3kMnZTI/uP5jOnXjqvbaM2mnD+Xy+K52VuZvHwPt3aK5vnesfrs9jfOds2vX3mEOYMwTo1AH/+bazoDS34rviXmAc8DDYDdZZZOxM1tPpjFbVPWkO0o4oOBHbiiVS27I8lpaoUH8entnRgzfydjF6ayfv8Jxt4aR5MaoXZHExGRs7QjLYdBE1ZRUOxi2ohOJDT8ywmNImfFx8fwn2tbEuTvyweLU3EUuXjtxjb4qgD/IxVhw6u3gfXAir+5phanpjyf7shpj/0fxpg7jDFrjDFrMjIySieliBuau/kwN32wAh8DX47souJbQfn5+vDIFc2ZMjyBzNwC+oxdxk9b0uyOJSIiZ2HdvuP0+3AFxsAXIzur+EqpMcbwWK8YHrisGbPWHuC+mesocrrsjuXWbC2/xpi3gG7ADZZlldqWp5ZlfWRZVrxlWfGRkTqzVLzTxKW7GTltLc1rh/HN3V1pGaVzBSu6Hs0i+f6ebjSKrMQdU5MY8/NOXC6t8xERqaiWpWRy6/hVhAf78+XILjSrGWZ3JPEwxhjuu6wpT1zZnNkbDzNq2loKinVSxPmyrfwaY8YAtwCXWJa16wyXpwE1//Czmqc9JiIlLMvi1bnbeW72Vq5oVZMZt19AjbAgu2PJWYqqEsznd3bmhri6vL0gmTumJpHjKLI7loiI/MG8LWkMm7SaehEhfDmyM/WqhtgdSTzYnT0a8+x1rZi/7Qi3TVmjoxLPky3l1xjzNv9bfLefxVNWABcaY07/BH85cAjYU/oJRdxTsdPFo19uZNyiVAZ0iub9WzvoeAU3FOTvyxs3teGZa1uycEc6vccuIyVdxyGJiFQUs5IOMHr6WlpGVeazOy+gRmXdZJayN6RLA169oTVLUzK5/RMV4PNhxzm/Y4FhwADguDGmVslX6GnXvGyMWXDa0z4F8oDJxphYY0xf4HHgrb/b6VnEm+QXOrlzahJfJB3g/sua8mKfWG2K4MaMMQzt2pBpIzpxIq+I699fxrKUTLtjiYh4vUnLdvPQFxu4oFFVpt/WiSohAXZHEi/Sv2M0r97QhqUpmdwxNUkF+BzZMfI7mlM7PC8ADp/29fBp19QGGv/2jWVZWZwa6Y0C1gBjgTeBt8onskjFlltQzJBJiSzckc4LfWK5/7JmGKPi6wk6N67Gd3d3JSo8mCETE5mZuM/uSCIiXuv9RSk8+/2pZUUTh3akUmBFODhFvE2/+Hq80rc1v+7MYOS0JK0BPge2n/Nb1nTOr3i6bEcRQycmsuFAFv/t345r20bZHUnKQLajiLumr2VJciYjezTm0StidN6fiEg5endBMm/+vJPr2kbxVr+2+PlWhENTxJvNSNzHE19t4pLmNRg3MI5AP+9d6na25/zqt1bEjWXlFTFw/Co2Hcxi7IA4FV8PVjnIn0lDOzKgUzQfLE7l7hlrNdVJRKQcWJbFmJ938ubPO+nbvg5j+rdT8ZUK4ZaEaF68PpZftqdz1/S1FBbrGKQz0W+uiJs6drKQWz5eyfbDOXwwsAO9YnWGr6fz8/XhxT6xPHlVC37cnEb/j1aSmVtgdywREY9lWRZv/LSDtxckc1OHurx+U1vtpyEVyq2d6vN871bM35bO3Z+u1TnAZ6DyK+KGjp0s5JaPVpKakcv4IfFc2uKPJ4GJpzLGcHv3RnwwsAM70rK5cdxy9h3NszuWiIjHsSyLV37cztiFqdyScGqTIRVfqYgGdW7As9e14qetR7h3xjqKVYD/ksqviJvJyi9i0IRV7Dl6kklDO9K9WaTdkcQGV7SqxfTbLuB4XhF9xy1ny6EsuyOJiHgMy7J4Yc42Pvx1F4MuqM+LfWK1z4JUaEO6NOCpa1ry4+Y0Hp21EZfLs/d1Ol8qvyJuJLegmKGTEkk+kstHg+Pp0qS63ZHERh3qR/DlyM74+xr6f7iS5ak6CklE5J+yLItnv9/KhKW7GdqlAc/1bqXiK25hRLeGPHh5M75ae5Bnvt+Cp29sfD5UfkXcRH6hkxGTV7PxQBbvDWhPD434CtC0ZhizRnWhdngQQyeu5odNh+2OJCLitn4b8Z28fA+3dWvIf65tqaMDxa3cc0kT7ujeiE9W7OW1eTvsjlPhqPyKuIGCYid3Tksicc8x3urXlp6ttLmV/K+oKsF8MbIzreuGc9ena5m6Yo/dkURE3I5lWbw2b8fvI75PXt1CxVfcjjGGJ65szq2dohm3KJWxC1PsjlShqPyKVHBOl8V9M9bz684MXu3bht7t6tgdSSqgKiEBTBvRiUtiavDUt1v0Zicico7+Oz+ZcYtSubVTtEZ8xa0ZY3i+dyx92kXx+rwdTFm+x+5IFYbKr0gFZlkW//luM3O3pPH0NS3p17Ge3ZGkAgsO8OXDQR1+f7N7be52rfcRETkLYxem8PaCZPrF1+X53rEqvuL2fHwMr9/Ulstb1uQ/323hizX77Y5UIaj8ilRgYxemMG3lPkb2aMzwbg3tjiNuwM/Xh7f6teOWhGjeX5TKs99v1Y6PIiJ/4+Nfd/H6vB1c374OL/dto82txGP4+/rw3oD2XNi0Oo/N2qh9QVD5FamwPl+9nzd+2knf9nV4rFeM3XHEjfj4GF66PpbhXRsyefkenvhqE04VYBGR/2fyst28+MM2rm5dm9dv1Dm+4nkC/U7NCouLjuC+metYmuzdJ0Oo/IpUQL9sP8ITX2+ie7NIXr2xjaZfyTkzxvDUNS2495ImfLZmPw98tp4iHXovIvK7T1ft45nvt9KzZU3+e3M7/Hz1sVg8U0iAHxOGdqRxZCh3TF3DxgMn7I5kG/2Wi1Qw6/efYPT0tbSsXZlxt8bhrzdjOU/GGB7sGcNjvZrz3YZDjJ6+FkeR0+5YIiK2+2rtAf719SYujonk3QHt9V4rHi882J8pwxOoWimAoZNWsysj1+5IttBvukgFcvBEPrdNWUNkWCATh3akUqCf3ZHEA4y6qDHP9W7Fz1uPMHJakgqwiHi1n7ak8ciXG+nSuBrjBnYg0M/X7kgi5aJm5SCmjuiEAQZNSORItsPuSOVO5VekgsgtKGbE5NUUFDuZNLQjkWGBdkcSDzK4cwNe6duaRTsyVIBFxGstT8nk7k/XEVsnnI8GxxPkr+Ir3qVh9UpMHpbAibxCBk9IJCuvyO5I5UrlV6QCcLos7p+5juT0XMYOiKNJjTC7I4kHujkhWgVYRLzW+v0nuO2TNTSoHsKUYR0J1ewq8VKt6566+bM78yQjpqwmv9B7Pg+o/IpUAK/8uI3529J55tqWdG8WaXcc8WCnF+A7p6oAi4h32JGWw9BJiVQPDWTqiE5UCQmwO5KIrbo2qc6Y/u1I2necuz9dS7GXbIqp8itisxmJ+/h4yW6GdmnAoM4N7I4jXuC3Arx4pwqwiHi+fUfzGDRhFQG+Pkwb0YmalYPsjiRSIVzdpjbP9Y5lwfZ0Hv9qE5bl+cciqvyK2GjlrqM89c1mejSL5N9Xt7A7jniRmxOiefUGFWAR8Wzp2Q4GTlhFQbGLqSM6EV0txO5IIhXKoAvqc/9lTfky6QCvzN1ud5wyp/IrYpNDJ/K5a/pa6lcL4d0B7XW+oJS7/h1VgEXEc53IK2TQhEQycwuYPKwjMbW0n4bIn7nv0qYMuqA+Hy7excSlu+2OU6b0aVvEBo4iJ6OmJVFQ7OLDQfFUDvK3O5J4qdML8Ojpayks9o41PyLi2U4WFDN00mp2Z57k48HxtI+OsDuSSIVljOGZ61pxRauaPD9nK3M2HrY7UplR+RUpZ5Zl8fS3m9lwIIs3+7WlSY1QuyOJl+vfMZqXrm/NL9vTuXfGOq/Z9EJEPFNhsYuR05LYdDCLdwe0p2uT6nZHEqnwfH0Mb9/cng7RETzw2XpW7Tpqd6QyofIrUs6mr9rH52sOcM8lTbiiVS2744gAMKBTNE9d05K5W9J46IsNOF2ev+mFiHgel8vi0S83sCQ5k5f7ttb7rMg5CPL3ZfyQeOpVDeb2T9aw80iO3ZFKncqvSDlK2nucZ7/fwkUxkdx/WTO744j8HyO6NeSRK2L4dv0hnvx6Ey4VYBFxMy//uI1v1h/ikSti6Bdfz+44Im6nSkgAU4YnEOTvy5CJiaRlOeyOVKpUfkXKSXqOg1HTkoiqEszb/dvj62PsjiTy/9x1cRPuvaQJM1fv59nvt3jFsQci4hk+/nUXHy/ZzZDO9Rl9UWO744i4rboRIUwa1pEcRzFDJyWS7SiyO1KpUfkVKQdOl8X9M9eT7Sjig4EdCA/RBldScT1weTNuv7AhU1bs5ZUft6sAi0iF9826g7z4wzaual2Lp69thTG6wSzyT7SKCueDgR1ISc/lzk+SKCj2jBMhVH5FysG7vySzPPUoz/WOpUXtynbHEflbxhj+dVULBl4QzYe/7uLtBcl2RxIR+Uu/7szg4S82cEGjqrzVr51mVomUkm5Nq/P6TW1YsesoH/+6y+44pcLP7gAinm55SiZvL0imb1wdbupQ1+44ImfFGMNz18VSUOTiv/OTCfL3ZWQPTSMUkYpl04EsRk1LokmNUD4aHE+Qv6/dkUQ8yvXt61IpwI/uzSLtjlIqVH5FylB6joN7Z66ncWQoL/SJ1TQscSs+PoZXbmiDo9jFKz9uJyTAl8GdG9gdS0QEgD2ZJxk6KfH3DXoqB2lJkUhZ6OlBu6ar/IqUkd/W+eYWFDH9tk6EBOjXTdyPr4/hrX5tcRQ5efrbLVQO8qdP+zp2xxIRL5eRU8DgiYm4LItPRiRQs3KQ3ZFExA1oza9IGTl9nW9MrTC744icN39fH969pT2dG1XjoS82MH/rEbsjiYgXyy0oZtjkRDJyCpg4tCONI0PtjiQibkLlV6QMrEg9qnW+4lGC/H35eEg8sVGVGf3pWlakHrU7koh4ocJiF6OmJbHtcA7v3xpH++gIuyOJiBtR+RUpZSfyCnngs/U0rFaJ53trna94jtBAPyYPS6B+1RBu/2QNGw+csDuSiHgRy7J4fNZGliRn8krf1lzcvIbdkUTEzaj8ipQiy7J48uvNZOYW8PbN7akUqHW+4lkiKgUwdUQnqoT4M2RiIinpOXZHEhEvMWZ+Ml+tO8iDlzfjpvh6dscRETek8itSimatPcicTYd5qGcMreuG2x1HpEzUCg9i2ohO+Pr4MHB8IgeO59kdSUQ83Oer9/POgmT6xdflnkua2B1HRNyUyq9IKdl79CT/+XYznRpW5Y7ujeyOI1KmGlSvxNQRCeQVFjNw/CoycgrsjiQiHurXnRk88fUmLmxanReSWdoOAAAgAElEQVSvb63lRCJy3lR+RUpBsdPF/Z+tx9fHMKZ/O3x99MYsnq9F7cpMGtaRI9mnjhzJyi+yO5KIeJith7IZPX0tTWuE8v6tcfj76qOriJw/vYKIlIJ3f0lh3b4TvNS3NVFVgu2OI1JuOtSvygeDOpCSnsPwyavJKyy2O5KIeIjDWfkMn7yasKBTm+2FBfnbHUlE3JzKr8g/lLT3OO/+cupYo2vaRNkdR6Tc9WgWyds3t2fdvuOMmraWwmKX3ZFExM1lO4oYNmk1uQXFTBzakVrhQXZHEhEPoPIr8g/kFRbz0OfriaoSzLPXtbI7johtrmpdm5f7tmbxzgwe+mIDLpdldyQRcVNFThejp60lJT2XcQPjaFG7st2RRMRD6BwWkX/gtbk72HM0j5l3XKDpWOL1+neM5nheEa/8uJ2qIf48c10rbUwjIufEsiye+GoTS1Myef3GNlzYNNLuSCLiQVR+Rc7TitSjTF6+h6FdGnBBo2p2xxGpEO7s3ojMnALGL91N9dBA7rm0qd2RRMSNvLMghS+TDnDfpU11lq+IlDqVX5HzcLKgmEe+3ECDaiE81qu53XFEKgxjDP+6qgXHThby5s87qRYayIBO0XbHEhE38GXSAcbM38kNcXW5/zLdOBOR0qfyK3IeXvphGwdP5PPFnZ0JDvC1O45IheLjY3j1xjYczyvk399somolf3rF1rY7lohUYMtSMnl81ka6NanOy311lq+IlA1teCVyjpYkZzB91T5u69aQ+AZV7Y4jUiH5+/ow9tY42tWrwr0z1rM8NdPuSCJSQW1Py2bk1CSa1Ajl/YFxBPjp46mIlA29uoicgxxHEY99uZHGkZV4qGeM3XFEKrSQAD8mDu1I/Woh3PFJEpsPZtkdSUQqmLQsB8MmrSYk0JeJQztSWZtHikgZUvkVOQcvzN5GWraDN25qS5C/pjuLnEmVkAA+GZFA5SA/hk5KZE/mSbsjiUgFkVtQzLDJq8lxFDNpaAJRVYLtjiQiHk7lV+QsLU3O5LM1+7mje2PaR0fYHUfEbdQOD+aTEZ1wuiwGT0wkPdthdyQRsVmR08Xo6WvZeSSHsbfG0TJKZ/mKSNlT+RU5C/mFTv719SYaVq+kHShFzkOTGqFMGpZAZm4BQyatJttRZHckEbGJZVk89c1mft2ZwUvXx9Kjmc7yFZHyofIrchb+O38n+47l8XLf1pruLHKe2tWrwgcDO5CSnsNtU9bgKHLaHUlEbPD+olRmrt7PPZc0oX9HHYUmIuVH5VfkDDYdyOLjJbu4JaEeFzSqZnccEbfWvVkkb9zUlsTdx7h3xjqKnS67I4lIOfp2/UFen7eD69vX4cHLm9kdR0S8jMqvyN8ocrp4bNZGqoUG8viVLeyOI+IRererw3+ubclPW4/w7282Y1mW3ZFEpBwk7j7GI19sJKFhVV65QWf5ikj587M7gEhFNn7JbrYezuaDgXGEB+v4BZHSMqxrQ47mFvLewhSqhQbwyBXN7Y4kImUoNSOXO6auoW7VYD4a1IFAPy0hEpHyZ8vIrzGmuzHmO2PMQWOMZYwZeobrG5Rc98evXuUUWbzQ7syT/Hf+Tq5oVZNesbXtjiPicR7q2YxbEuoxdmEqE5futjuOiJSRo7kFDJu0Gl9jmDw0gSohAXZHEhEvZdfIbyiwGfik5Ots9QI2nPb9sdIMJfIby7J44quNBPj68FzvWLvjiHgkYwwv9GnN8ZNFPDd7K1UrBdCnfR27Y4lIKXIUObn9kzUcyXYw844LiK4WYnckEfFitoz8Wpb1g2VZ/7Is60vgXHY7OWpZVtppX4VllVG82xdrDrBy1zGeuKoFNSsH2R1HxGP5+hj+e3M7LmhUlYe/2MCiHel2RxKRUuJyWTz4+XrW7T/B2ze3o310hN2RRMTLuduGV18ZY9KNMcuMMTfaHUY807GThbz04zbi60dwc8d6dscR8XhB/r58NDieZjXDGDVtLev2Hbc7koiUglfnbueHTWk8eVULLR8SkQrBXcpvLvAw0A+4ClgAfGaMGWhrKvFIr83dTo6jmBeuj8XHRztRipSHykH+TB7ekciwQIZNXk1Keo7dkUTkH5i2ci8f/rqLwZ3rM6JbQ7vjiIgAblJ+LcvKtCzrTcuyVlqWtcayrKeBD4FH/+x6Y8wdxpg1xpg1GRkZ5RtW3FrS3uPMXL2fEd0a0rxWZbvjiHiVGmFBTB2RgJ+PD4MmJHLoRL7dkUTkPCzcns7T327m0uY1ePqaljrSSEQqDLcov39hFdD0zx6wLOsjy7LiLcuKj4yMLOdY4q6KnS7+/c1malUO4r5L//R/LREpY/WrVWLK8I7kOooZNGEVx05qawcRd7L5YBZ3fbqWFrUr884t7fHzdeePmiLiadz5FakdcNjuEOI5pqzYy7bD2fzn2pZUCtQR2CJ2aRUVzvgh8Rw4ns/QSYnkFhTbHUlEzsKhE/mMmLKaKsH+TBzaUe+lIlLh2HXOb6gxpp0xpl1JhuiS76NLHn/ZGLPgtOuHGGMGGGNaGGNijDEPA3cB79qRXzxPWpaDt37aQY9mkfSKrWV3HBGv16lRNcYOiGPLoWzu+GQNBcVOuyOJyN/IcRQxfPJq8gqcTBzWUScliEiFZNfIbzywruQrGHi25J+fK3m8NtD4D8/5N7AGWA3cDAy3LGtMuaQVj/fCnK0UuSye691Ka5NEKojLWtbktRvasDz1KPfNWI/TZdkdSUT+RJHTxejpa0lJz+X9gXHaM0NEKixb5qNYlrUI+MuGYVnW0D98PwWYUrapxFstSc5g9sbDPHBZM+pXq2R3HBE5zQ0d6nIiv4jnZ2/lya838XLf1rpBJVKBWJbFU99sZklyJq/e0JoLm2qvFRGpuLQYQ7xaQbGTp7/dQoNqIdzZo5HdcUTkT4zo1pDjJwt5b2EKEZUCeKxXc7sjiUiJcYtTmbl6P3df3IT+HaPtjiMi8rdUfsWrTVi6m92ZJ5kyPIEgf1+744jIX3ioZzOO5RUyblEqESH+3NH9jytjRKS8fbfhEK/N3UHvdlE81LOZ3XFERM5I5Ve81pFsB+/9ksLlLWvSo5mmaYlUZMYYnu8dS1Z+ES/9sJ0qIQH0i69ndywRr7V6zzEe/nwDCQ2q8tqNbbQcQUTcgsqveK1XftxOsdPi31e3sDuKiJwFXx/DmH7tyM4v4vFZGwkP9ueKVtqdXaS87crI5fZP1lA3IpgPB3Ug0E8zp0TEPbjzOb8i5y1p73G+XneQ2y5sqE2uRNxIgJ8PHwzsQJu6VbhnxjqWp2baHUnEqxzNLWDY5NX4GMOkYR2JqBRgdyQRkbOm8itex+WyeOa7LdSsHMhdFzexO46InKNKgX5MGtqR+lVDuOOTJDYdyLI7kohXyCssZviUNaRlOfh4cLxuHouI21H5Fa/zRdJ+Nh3M4okrW1ApUDP/RdxRRKUApo7oRHiwP0MmJZKakWt3JBGPVux0cc+n69h04ATv3tKeDvUj7I4kInLOVH7Fq2TlF/Ha3B10qB9B73ZRdscRkX+gVngQU0ckYIDBExI5nJVvdyQRj2RZFk9/t4UF29N5tncsPbXWXkTclMqveJV3FiRzLK+QZ65tpZ0pRTxAo8hQpgxPIDu/iEETEjl+stDuSCIe5/1FqXy6ah+jLmrMoAvq2x1HROS8qfyK10hJz2HK8j30j69H67rhdscRkVISWyecj4fEs+9YHkMnJZJbUGx3JBGPMSvpAK/P20GfdlE80jPG7jgiIv+Iyq94BcuyeG72NoIDfHn4Cr15i3iaCxpVY+yAODYfymb45NXkFzrtjiTi9pYkZ/DYrI10bVKN125si4+PZkyJiHtT+RWvsGhHBr/uzOC+S5tSPTTQ7jgiUgYub1mTMf3bsXrPMUZOS6KgWAVY5HxtOZTFqGlraVIjlHEDOxDgp4+MIuL+9EomHq/Y6eLFH7bRoFoIgzs3sDuOiJSh69pG8Urf1izemcG9M9ZR7HTZHUnE7Rw8kc+wSasJC/Jj0rCOVA7ytzuSiEipUPkVjzdz9X5S0nN5/MrmunMt4gX6d4zm6WtaMm/LER7+YgMul2V3JBG3kZVXxNCJieQXOZk8LIHa4cF2RxIRKTU65FQ8Wo6jiDE/7yShQVWu0NEMIl5jeLeG5Bc5eX3eDkIC/XixT6x2eBc5A0eRk9unrmHv0TymDE8gplaY3ZFEREqVyq94tHGLUjl6spCJQ1vog6+Il7nr4iacLCjm/UWphPj78uTVeh0Q+StOl8VDX2wgcfcx3rmlPZ0bV7M7kohIqVP5FY918EQ+E5bupne7KNrWq2J3HBGxwSNXxJBX6GT80t1UCvTjgcub2R1JpMKxLItnv9/CnI2HefKqFlzXNsruSCIiZULlVzzW63O3Y3Hqw6+IeCdjDE9f05K8wmLeXpBMSIAvd/ZobHcskQrlnQUpfLJiL3d2b8Tt3RvZHUdEpMyo/IpH2njgBN+sP8SoixpTNyLE7jgiYiMfH8PLfduQV+jk5R+3ExLgyyDt/C4CwLSVexkzfyc3xNXl8Sub2x1HRKRMqfyKx7EsixfmbKNapQBGX6QRHhEBXx/DmP7tcBQ5eerbLQT5+3JTfD27Y4nY6odNh3nq281c2rwGr9zQWmviRcTj6dwX8Tg/bT1C4u5j3H95M8J0NqGIlPD39eG9AXFc2LQ6j87ayKykA3ZHErHNspRM7p+5ng7REbw3IA5/X30kFBHPp1c68SiFxS5e+XE7TWqEcktHjeqIyP8V5O/Lx4Pj6dq4Og9/uYFv1h20O5JIudt0IIs7PllDw+qVmDCkI8EBvnZHEhEpFyq/4lFmJO5jd+ZJ/nVVc/x0F1tE/sRvBbhzo2o8+Pl6vl2vAizeY3fmSYZOSqRKSABThicQHqIZUiLiPdQOxGOcLCjm3V+S6dSwKhfH1LA7johUYMEBvowfEk9Cw6o88Nl6Zm88ZHckkTJ3JNvBoAmrsIBPRiRQKzzI7kgiIuVK5Vc8xvglu8nMLeTxK5tr0w4ROaOQAD8mDu1IfP2q3DdzPT9sOmx3JJEycyKvkMETEjl2spBJQzvSODLU7kgiIuVO5Vc8QmZuAR/9mkqvVrVoHx1hdxwRcRMhAX5MGtaR9vWqcO+MdczdrAIsnifHUcSQiYnszjzJR4PiaVuvit2RRERsofIrHuG9X1LIL3Ly8BUxdkcRETdTKdCPycMTaFM3nLs/Xce8LWl2RxIpNfmFTkZMWcPmQ9mMvTWObk2r2x1JRMQ2Kr/i9vYfy2P6qr3071iPJjU0jUtEzl1ooB9ThicQWyecuz9dy08qwOIBCotdjJyWxOo9x3irX1sub1nT7kgiIrY65/JrjAk0xjQ0xrQ0xkSWRSiRc/HmTzvwMYb7Lm1mdxQRcWNhQf58MiKBllHhjJ6+ljkbNQVa3Fex08V9M9exeGcGL1/fmt7t6tgdSUTEdmdVfo0xYcaYUcaYX4EsIAXYDKQZY/YZYz42xnQsy6Aif2broWy+3XCIYV0batdKEfnHKgf5M21EAu2jq3DPjLV8ve6A3ZFEzpnLZfHorI38uDmNp65pyc0J0XZHEhGpEM5Yfo0xDwJ7gOHAz0BvoB3QDOgMPAP4AT8bY+YaY5qWVViRP3pt3nbCAv0Y1aOx3VFExEOEBfkzZXgCFzSqxoOfb+Cz1fvsjiRy1izL4pnvt/DV2oM8eHkzRnRraHckEZEKw+8srrkA6GFZ1ua/eDwRmGiMGcWpgtwDSC6lfCJ/aUXqURbtyOCJK5sTHuJvdxwR8SC/HYN059QkHpu1iYJiF4M7N7A7lsjfsiyLV+fu4JMVe7mzeyPuuaSJ3ZFERCqUM478WpbV77fia4zZbIwJ/4vrHJZlvW9Z1vjSDinyR5Zl8crc7dQOD2JIlwZ2xxERDxTk78tHgztwecuaPP3tFsYv2WV3JJG/ZFkWb/60kw8WpzLwgmideS8i8ifOdcOrlkDgH39ojAk3xowtnUgiZzZvSxob9p/g/suaEuTva3ccEfFQgX6+vH9rHFe3rs0Lc7bx3i+a2CQV05j5yby3MIVbEurx3HWxKr4iIn/ibKY9Y4z5kVPTmy2gHpD+h0tCgDuBu0o1ncifKHa6eH3eDprUCOWGuLp2xxERD+fv68PbN7cj0M+HN37aeepM8Z4xKhdSYfx3/k7eWZBM//h6vNinNT4++n9TROTPnFX5BTZxai2vARKNMTnABmAdsBFoDuhMCCkX36w/RGrGST4YGIefr46qFpGy5+frwxs3tSXQ34exC1M5kVfEc71j8VXJEJu9PT+Z/85P5qYOdXm5r4qviMjfOavya1nWowDGmAJO7fAcxakdn9sBV5f8OY+WUUaR3xUWu3h7wU5i61Tmila17I4jIl7Ex8fw0vWtqRISwLhFqWTlF/FWv3YE+OkmnNjj3QXJjJm/kxvi6vLKDW1UfEVEzuBsR35/E2pZVhGwFphdBnlE/tbna/az/1g+zw3TeiYRKX/GGB7r1ZyIEH9e+mE7WflFfDioAyEB5/p2KvLPjF2Ywps/76Rv+zq8dmMbzUIQETkL53S7uqT4itjCUeTk3V+S6VA/gouaRdodR0S82B3dG/PajW1YlpLJreNXcSKv0O5I4iUsy+LdBcm8Pm8HfdpF8fpNbVV8RUTO0hnLrzHmrE9HN6fU+2eRRP7c9FX7OJJdwEM9m2nUV0Rs1y++Hu/f2oEtB7Pp9+EK0rIcdkcSD2dZFq/N2/H7iO+b/dqp+IqInIOzGfldYYyZYIzp/FcXGGMijDGjgK1A71JLJ1LiZEEx4xal0KVxNbo0rm53HBERAHrF1mLy8I4cOuHgxg+Wsysj1+5I4qFcLotnv9/KuEWpDOgUzRsa8RUROWdnU36bA8eAOcaYTGPMPGPMJGPMOGPMTGPMRk4dfTQQuN+yrPfKMrB4pykr9pCZW8hDPZvZHUVE5P/o0rg6M26/gPxCJ33HLWfNnmN2RxIP43RZ/OvrTUxevocR3RryYp9YbW4lInIezlh+Lcs6YVnWI0AdYCSwDagCNASKgSlAe8uyulqWNa8sw4p3ynYU8eHiXVwcE0mH+lXtjiMi8v+0rhvOV6O7EBESwIDxq/hhk07/k9JR7HTx4Ofrmbl6P/de0oR/X91CS39ERM7TWW9PaVlWvjFmsWVZX5ZlIJE/mrBkN1n5RTzUM8buKCIif6l+tUrMGtWF2z9Zw12fruXJq1owoltDFRU5bwXFTu6bsZ65W9J4tFcMoy9qYnckERG3dq6HEy43xjQqkyQif+L4yUImLN3NlbG1iK0TbnccEZG/VbVSANNv68SVsbV4Yc42nvluC06XZXcscUO5BcUMn7yauVvSePqaliq+IiKl4FzL7w+cKsBxp//QGNPdGLOs9GKJnPLhr7s4WVjMA5drra+IuIcgf1/euyWO27o1ZMqKvYyclkR+odPuWOJGMnMLuPmjFazcdYw3b2rL8G5nffCGiIj8jXM95/c+4A3gF2NMT2NMO2PMXGAhsK8sAor3Ss9xMHn5bnq3jaJZzTC744iInDUfH8O/r2nJf65tyfxtR3QUkpy1fUfzuHHcclLTTzJ+cDw3dKhrdyQREY9xriO/WJb1BvAyMBtYDeQAbSzLuqWUs4mXG7colSKnxX2XadRXRNzTsK4N+XhQPLsycrnuvaWs33/C7khSgW05lEXfccs5kV/E9Ns7cXHzGnZHEhHxKOdUfo0x9YwxHwLPcar4FgBzLMvaUhbhxHsdzspn+sp93BhXl4bVK9kdR0TkvF3WsiZfje5KgJ8P/T5cwbfrD9odSSqg5amZ9P9wJQG+hi9HdiYuOsLuSCIiHudcR36TgfbANZZldQWuA/5rjHmy1JOJV3t/YSouy+LuS7TBh4i4v5haYXx3dzfa16vCfTPX8+rc7bi0EZaU+DLpAEMmJlI7PIhZo7vQpIaW+oiIlIVzLb+3WpaVYFnWzwCWZf0C9ABGG2PeL/V04pUOZ+Xz2er93BRfj3pVQ+yOIyJSKqpWCmDqiE7ckhDNuEWp3DF1DbkFxXbHEhu5XBZvzNvBw19sIKFhVb4c1YXa4cF2xxIR8VjnuuHVrD/52QagC3BRKWUSLzdu0alR37submx3FBGRUhXg58NL18fy7HWtWLgjg97vLSUlPcfuWGIDR5GTe2eu472FKfSPr8fkYQmEB/vbHUtExKOd84ZXf8ayrL1A19L4s8S7Hc7KZ2bifm6Kr0vdCI36iojnMcYwpEsDpo5IICu/iOveW8b3Gw7ZHUvKUWZuAQM+XsnsjYd54srmvHJDa/x9S+UjmYiI/I1Se6W1LOt4af1Z4r1+G/UdfZHW+oqIZ+vSuDqz77mQFrUrc8+MdTz7/RYKi112x5IytuVQFn3GLmPr4Ww+GBjHnT0aY4yxO5aIiFew5TajMaa7MeY7Y8xBY4xljBl6Fs9pbYxZbIzJL3ne00bvFh4lLcvx+6iv1vqKiDeoFR7EzDsuYHjXhkxatodbPl6p84A92LfrD3LDuOUUOy0+u6MzvWJr2x1JRMSr2DXHJhTYDNwH5J/pYmNMZeBn4AjQseR5jwAPlmFGKWfjFqVo1FdEvI6/rw9PX9uS9wa0Z9vhbK5+ZwmLdqTbHUtKUbHTxYtztnLfzPW0qVOF7+/pRtt6VeyOJSLidWwpv5Zl/WBZ1r8sy/oSOJs5XrcCIcAQy7I2lzzvVeBBjf56hrQsBzMS93NjB436ioh3uqZNFN/d3ZXIsECGTlrNC7O3UlDstDuW/EPHThYyZFIiHy/ZzeDO9Zl2WyciwwLtjiUi4pXcZXeFzsASy7JOHyWeB0QBDWxJJKXqt1Hfuy7WqK+IeK8mNcL45q6uDOlcn/FLd9P3/eWkZuTaHUvO09p9x7n23aWs3nOc125sw3O9Ywnwc5ePXiIinsddXoFrcWrK8+mOnPbY/2GMucMYs8YYsyYjI6PMw8k/k5blYMZqjfqKiAAE+fvybO9YPh4cz6ET+VzzzlI+W70Py7LsjiZnyeWy+HBxKv0+WIEx8MWdnekXX8/uWCIiXs9dyu85sSzrI8uy4i3Lio+MjLQ7jpzBB4tTcbk06isicrrLW9bkx/u60z66Co/N2sTIaUlk5BTYHUvO4GhuAcOnrOblH7dzecuazLn3Qq3vFRGpINyl/KYBNf/ws5qnPSZu6ki2g08T93FDnEZ9RUT+qFZ4EFNHdOKJK5uzcEcGPccs5vsNhzQKXEGtSD3KVe8sYXnKUZ7v3Yr3b40jPNjf7lgiIlLCXcrvCuBCY0zQaT+7HDgE7LElkZSKcYs06isi8nd8fQx39mjMD/d2I7paJe6ZsY7R09eSmatR4IrCUeTkhdlbGTB+JcH+vnw1uguDOjfQ+b0iIhWMXef8hhpj2hlj2pVkiC75Prrk8ZeNMQtOe8qnQB4w2RgTa4zpCzwOvGXp9rfb+m3Ut29cHaKradRXROTvNKkRxqyRnXm0VwwLtqXTc8yvfLv+oEaBbbbpQBbXvruU8Ut3c2unaObceyGxdcLtjiUiIn/CrpHfeGBdyVcw8GzJPz9X8nhtoPFvF1uWlcWpkd4oYA0wFngTeKv8IktpG7coFafL4u6Lm9odRUTELfj5+jD6oibMvrcb9SKCuW/megZPTGR35km7o3mdIqeLt+cnc/37y8h2FDFleAIv9GlNpUA/u6OJiMhfMJ5+xzg+Pt5as2aN3THkD9KzHVz42kJ6t4vitRvb2h1HRMTtOF0W01ft5fW5Oyhwuhh9UWNG9mhMkL+v3dE83rp9x3niq01sT8uhd7sonrsulvAQre0VEbGLMSbJsqz4M12n25Nii3GLUynWqK+IyHnz9TEM7tyAXq1q8fycbfx3fjLfrj/Ef65tyUUxNeyO55GyHUW8MW8HU1fupWZYEB8O6sAVrf7fiYsiIlJBucuGV+JBMnIK+HTVPq5vr7W+IiL/VI3KQbx7S3umjkjAsiyGTlrN4ImJ7EjLsTuax7AsizkbD3P5W4uZunIvQzo3YP5DPVR8RUTcjEZ+pdxNWLqbopIpeiIiUjoubBrJTw/04JMVe3hnQTJXvv0rNydE88BlzYgMC7Q7ntvafDCL52ZvJXH3MVrWrsxHg+J1bq+IiJtS+ZVydSKvkKkr9nB1mygaRYbaHUdExKME+Plw24WNuCGuLm8vSGbayr18t/4QI7o1ZHi3hjpz9hykZzt446cdfJF0gIiQAF68Ppb+8fXw89WkORERd6XyK+Vq0rI9nCx0ctfFGvUVESkrEZUCeOa6VgzuXJ/X5u7g7QXJTFq2m9svbMTQrg0IC1IJ/isn8gr5eMkuJi3bQ5HTxe0XNuLuS5pQWf/NRETcnsqvlJscRxGTlu2mZ8uaNK9V2e44IiIer1FkKB8M6sDmg1n8d34yb/68kwklJXhgp/raofg0OY7/ae/Oo7yu73uPP9+zwIDs+7AMiyCoFUVxw6hIJIva3JjcamIatVmMa2yz9aZNb9P2xtyaHK2JUpfkJGpbTZvmNjEhiYkCEjUScMENZVhFdtmXAWbmc/+YMSUIyjIzn9/yfJwzB37f7+c3vM7he36/3+v3+X6/nz189zdL+O7sJWzd1chF42v5/HvGMrLfUbmjSZLaiOVXHeb+3y5jS0Mj108ZnTuKJJWVPxrSk+9cMZH5KzZx669e5Ru/fIU7ZtRz6anD+MRZIxnWp3xvPrh+2y7ue3IZ9z25lE079vDe4wfyF1OP8UtaSSpBll91iJ27m/ju7CWce0x/xg/1RiGSlMP4ob343p+dxksrt/Cd2Yu5/8ll3PvEUt5/Qi1XThrBxOG9iYjcMTvE4nXbuGf2Ev7z6RXsbmzm/GMHcuO7x3DC0J65o0mS2onlVx3igTnLeWP7bm5w1leSsjtucA9uufQkvvi+sXz/8aX821PL+dn8VYwe0I2PnlbHhyYMofdRnXLHbHN7mpp55OU1PDDnNavT0M4AABbWSURBVB5buI7qygo+fPIQPvmuUYwe4E0YJanURUopd4Z2NXHixDR37tzcMcrarsYmzrl5BiP7HcWDV52ZO44kaR87djfy0+dW8cDvlvPM8k10qqpg6nED+cCJgzn3mP7UVFfmjnhE6tdu40dPr+Df565g/bZdDOpRwyWnDuPjZwx3GShJKgERMS+lNPGdxjnzq3b3w3krWLNlF7dcclLuKJKk/ejaqYpLTh3GJacO4+VVW3hwznIemr+Kn81fRffOVUw9fiAXja9l0tH9iqYIL12/nZ89v4qHnlvJgtVbqQiYMm4AHzm1jslj+7tkkSSVIWd+1a72NDVz3jdn0q9bZ/7ftZPK5loySSp2jU3NPLHoDR56biW/eHE1Wxsa6VxVwRmj+nLe2P5MHjuA4X27Fszr+q7GJn63ZCOzXl3LrFfX8eqabQCcXNeLi8YP5sLxtQzsUZM5pSSpPRzszK/lV+3qh/NW8IX/eI7vXjGRdx87MHccSdJh2NXYxJOL3mDmK+uY9eo6lqzfDsCA7p2ZOKI3J9f1ZuKIPowb1L3DZoY37djNM8s38fTyjS0/yzaxc08TnSorOG1kHyaP7c/7T6hlSK8uHZJHkpSPpz0ru6bmxLQZ9Rxb24Mp4wbkjiNJOkydqyqZPHYAk8e2vJYve2M7jy1cz9ylG5i7dCPTn18NQEXAiH5HMXZgd44Z2J0R/boyuGcXBvfqwqCeNVQf4qnGO3c3sWrzTlZvbmDFpp0sXLOVhWu3sXDNNl7ftBOAyopg3KDuXDJxKOcc058zj+5L105+vJEkvZXvDmo3059fxeL125n2sZML5rQ4SdKRG973KD7e9yg+fsZwAFZvbuDp5RtZsGoLr6zZyoLVW/nFi6vZ++SyCOjeuYoeXarpXlNN95oqqiuDIHjzLWLH7ia272pka0MjWxv2sKWh8Q/+3U5VFYzu342JI3rzsUF1TBjWmxOH9bTsSpIOiu8WahfNzYk7ZtQzekA33nf8oNxxJEntaFDPGi44oZYLTqj9/baGPU2s3LSTlZsaWLlpJ69v2snmnXvYsrOl1G5p2EPDnmZSSrzZkbt2qqTvUV3pVlNF985VDOhRQ23PGmp7dmFwrxqG9u5KZYVfpkqSDo/lV+3ikQVrWbB6K7deeiIVflCRpLJTU13JqP7dGNXf9XMlSYXB+/yrzaWUuP3RhdT16cofjx+cO44kSZIkWX7V9mYvXM9zKzZz7eSjXUdRkiRJUkGwmajN3f5oPbU9a/jQyUNzR5EkSZIkwPKrNvbU4jeYs3QDnzlnFJ2qPLwkSZIkFQbbidrU7TPq6detEx85rS53FEmSJEn6Pcuv2syzr21i9sL1fPrsUdRUV+aOI0mSJEm/Z/lVm7n90Xp6da3mY2cMzx1FkiRJkv6A5Vdt4qWVW/j1y2v4xFkj6dbZ5aMlSZIkFRbLr9rEHTPr6d65iismjcgdRZIkSZLewvKrI1a/dhvTn1/F5ZOG07NLde44kiRJkvQWll8dsWkz66mpquQTZ43MHUWSJEmS9svyqyOy/I0d/PjZlVx2eh19u3XOHUeSJEmS9svyqyNy52OLqIzgqnNG5Y4iSZIkSQdk+dVhW7V5Jz+cu4JLTh3KwB41ueNIkiRJ0gFZfnXY7n5sMc0p8Zlzjs4dRZIkSZLeluVXh2X9tl08MGc5F08YwrA+XXPHkSRJkqS3ZfnVYfnO7CXsbmzmmsnO+kqSJEkqfJZfHbJNO3Zz/5NLuXD8YEb175Y7jiRJkiS9I8uvDtn3n1jK9t1NXHees76SJEmSioPlV4dk265Gvvf4UqYeN5Bxg3rkjiNJkiRJB8Xyq0PyL79dxuade7j+vNG5o0iSJEnSQbP86qA17GniO7MXc/aYfpw4rFfuOJIkSZJ00Cy/OmgPzlnO+m27uWHKmNxRJEmSJOmQWH51UHY3NnPXY4s5bUQfThvZJ3ccSZIkSTokll8dlB89vYJVmxu4forX+kqSJEkqPpZfvaPGpmamzVzE+KE9OXtMv9xxJEmSJOmQWX71jn46fxXLN+zg+vNGExG540iSJEnSIbP86m01NyfumFHP2IHdOf/YgbnjSJIkSdJhsfzqbT380moWrt3GdVNGU1HhrK8kSZKk4mT51QGllPj2o/WM7HcUF55QmzuOJEmSJB02y68OaOar63hx5RaumXw0lc76SpIkSSpill/tV0qJ2x+tZ0ivLlw8YUjuOJIkSZJ0RCy/2q/fLt7AvGUbufrcUVRXephIkiRJKm62Gu3X7TMW0r97Z/5k4rDcUSRJkiTpiFl+9RZPL9/I4/VvcNXZo6iprswdR5IkSZKOmOVXb3HHo/X06lrNZafX5Y4iSZIkSW0iW/mNiGsjYklENETEvIg4+23GTo6ItJ+fcR2ZuRy8uHIzjyxYyyfPGslRnatyx5EkSZKkNpGl/EbEpcBtwE3ABOAJ4OcR8U5TjccDtXv9LGzPnOVo2oxFdO9cxeWTRuSOIkmSJEltJtfM7+eA76eU7kkpvZxSugFYBVzzDs9bm1JavddPU/tHLR/1a7cy/YVVXD5pOD27VOeOI0mSJEltpsPLb0R0Ak4BHt5n18PApHd4+tyIWBURj0TEee0SsIxNm7mImqpKPnHWyNxRJEmSJKlN5Zj57QdUAmv22b4GGHSA57w5K/xh4EPAK8Ajb3edsA7N0vXb+fGzK7ns9Dr6duucO44kSZIktamiuKNRSukVWgrvm56MiBHAF4HZ+46PiKuAqwDq6rxj8cG4Y0Y9VRXBZ84ZlTuKJEmSJLW5HDO/64EmYOA+2wcCqw/h9zwFjNnfjpTS3SmliSmlif379z+8lGXktQ07+NEzr/PR0+oY0KMmdxxJkiRJanMdXn5TSruBecDUfXZNpeWuzwfrJFpOh9YRumNGPZURXH3u0bmjSJIkSVK7yHXa8y3A/RExB3gcuBoYDNwJEBH3AaSULm99/OfAUuBFoBPwp8AHabkGWEdgxcYd/HDeCi47vY5BPZ31lSRJklSaspTflNIPIqIv8BVa1ut9AbggpbSsdci+F+p2Ar4BDAV20lKCL0wpTe+gyCXrn2cuIgJnfSVJkiSVtGw3vEopTQOmHWDf5H0e3wzc3AGxysrKTTv597mvccnEYQzu1SV3HEmSJElqNzlueKUCceesRQBcM9lZX0mSJEmlzfJbplZvbuDBOa/xP08ZytDeXXPHkSRJkqR2ZfktU3fOWkRzSlw7eXTuKJIkSZLU7iy/ZWjtlgYemLOcD508hGF9nPWVJEmSVPosv2XorscW09icuO48Z30lSZIklQfLb5lZt3UX//rUMv7HSYMZ3veo3HEkSZIkqUNYfsvMPbMXs7uxmeud9ZUkSZJURiy/ZWT9tl3c/+QyPnDiYEb175Y7jiRJkiR1GMtvGbln9mIaGpu4fsqY3FEkSZIkqUNZfsvEuq27uO+JZfzx+MGMHuCsryRJkqTyYvktE3fOWsSuxiZuPN9ZX0mSJEnlx/JbBtZsaeBffruMiycM5Wiv9ZUkSZJUhiy/ZWDajHoamxM3vttZX0mSJEnlyfJb4l7ftJMH5rzGn5wylLq+XXPHkSRJkqQsLL8l7vZH60kkrp/iur6SJEmSypflt4S9tmEH/zH3NT5yah1DezvrK0mSJKl8WX5L2LceWUhFRXDdec76SpIkSSpvlt8StWT9dn70zOt87PQ6BvWsyR1HkiRJkrKy/Jaobz2ykOrK4JrJR+eOIkmSJEnZWX5LUP3arfzXs69zxZkjGNDdWV9JkiRJsvyWoFt/vZAu1ZVcdc6o3FEkSZIkqSBYfkvMSyu38LP5q/izs0bQt1vn3HEkSZIkqSBYfkvMN365gB41VVx1ttf6SpIkSdKbLL8lZM6SDcx4ZR3XTB5Nz67VueNIkiRJUsGw/JaIlBI3/2IBA7p35spJI3LHkSRJkqSCYvktEY8uWMvcZRu58fwxdOlUmTuOJEmSJBUUy28JaGpO3PyLVxjRtyuXTByWO44kSZIkFRzLbwn4yXOv88qarXz+PWOprvS/VJIkSZL2ZVMqcrsbm7nlV69y/OAeXHhCbe44kiRJklSQLL9F7sHfLee1DTv50vvGUVERueNIkiRJUkGy/Bax7bsa+dYj9Zw+sg/njOmXO44kSZIkFSzLbxH73uNLWL9tF1963zginPWVJEmSpAOx/BapN7bt4q5Zi5l63EBOGd47dxxJkiRJKmiW3yL1rUcWsmNPE3/5vnG5o0iSJElSwbP8FqFF67bxr08t57LT6hg9oFvuOJIkSZJU8Cy/Regff76AmupKbjx/TO4okiRJklQULL9F5qnFb/DwS2u4ZvLR9OvWOXccSZIkSSoKlt8i0tycuGn6y9T2rOGT7xqZO44kSZIkFQ3LbxF5aP5KnluxmS+8Zyw11ZW540iSJElS0bD8FomGPU3c/ItXOK62BxdPGJI7jiRJkiQVFctvkbj3iaW8vmknX7nwWCoqInccSZIkSSoqlt8isGH7bm6fUc+UcQOYNLpf7jiSJEmSVHQsv0Xgmw+/wo7dTXz5/eNyR5EkSZKkomT5LXAvrtzMA3OWc/mZwxkzsHvuOJIkSZJUlCy/BSylxN/95CV6d+3En59/TO44kiRJklS0LL8F7KfzVzFn6Qa++N6x9OxSnTuOJEmSJBUty2+B2rm7ia9Pf5njB/fgkonDcseRJEmSpKJm+S1Q/zxrESs3N/DVDxxPpUsbSZIkSdIRsfwWoBUbd3DXrEV84MTBnDqiT+44kiRJklT0LL8F6KbpL1MRwZcvcGkjSZIkSWoLlt8C85uF65n+/GqunXw0tT275I4jSZIkSSXB8ltAGvY08Tc/foERfbvy6XNG5Y4jSZIkSSWjKncA/be7Zi1myfrt3P/J06iprswdR5IkSZJKhjO/BWLp+u3cMbOei8bXcvaY/rnjSJIkSVJJyVZ+I+LaiFgSEQ0RMS8izn6H8ee2jmuIiMURcXVHZW1vKSX+909epFNlBX9z0XG540iSJElSyclSfiPiUuA24CZgAvAE8POIqDvA+JHA9NZxE4CvA9+OiA93TOL2Nf351Tz26jo+/55jGNijJnccSZIkSSo5uWZ+Pwd8P6V0T0rp5ZTSDcAq4JoDjL8aWJlSuqF1/D3AvcAXOihvu9nasIe//+mLHD+4Bx8/Y3juOJIkSZJUkjq8/EZEJ+AU4OF9dj0MTDrA087cz/hfAhMjorptE3asW3+1kLVbd/G1i0+gqtJLsCVJkiSpPeRoW/2ASmDNPtvXAIMO8JxBBxhf1fr7/kBEXBURcyNi7rp1644wbvvq2qmSK84cwUnDeuWOIkmSJEklqySXOkop3Q3cDTBx4sSUOc7b+sJ7x5JSQUeUJEmSpKKXY+Z3PdAEDNxn+0Bg9QGes/oA4xtbf19Ri4jcESRJkiSppHV4+U0p7QbmAVP32TWVlrs578+TBxg/N6W0p20TSpIkSZJKTa47LN0CXBkRn4qIYyPiNmAwcCdARNwXEfftNf5OYEhE/FPr+E8BVwLf7OjgkiRJkqTik+Wa35TSDyKiL/AVoBZ4AbggpbSsdUjdPuOXRMQFwK20LIe0EvhsSuk/OzC2JEmSJKlIZbvhVUppGjDtAPsm72fbLODkdo4lSZIkSSpBLiwrSZIkSSp5ll9JkiRJUsmz/EqSJEmSSp7lV5IkSZJU8iy/kiRJkqSSZ/mVJEmSJJU8y68kSZIkqeRZfiVJkiRJJS9SSrkztKuIWAcsy53jHfQD1ucOobLncahC4bGoQuBxqELgcahCUejH4vCUUv93GlTy5bcYRMTclNLE3DlU3jwOVSg8FlUIPA5VCDwOVShK5Vj0tGdJkiRJUsmz/EqSJEmSSp7ltzDcnTuAhMehCofHogqBx6EKgcehCkVJHIte8ytJkiRJKnnO/EqSJEmSSp7lV5IkSZJU8iy/mUXEtRGxJCIaImJeRJydO5PKR0R8OSJ+FxFbImJdRDwUEX+UO5fKW+txmSLi9txZVH4iojYi7m19TWyIiJci4tzcuVQ+IqIyIv5hr8+HSyLi/0REVe5sKl0RcU5E/CQiXm99D75yn/0REV+NiJURsTMiZkbE8ZniHjbLb0YRcSlwG3ATMAF4Avh5RNRlDaZyMhmYBkwCpgCNwK8jok/OUCpfEXEGcBUwP3cWlZ+I6AU8DgRwIXAscAOwNmculZ2/BK4DPguMA25sffzlnKFU8roBL9ByvO3cz/4vAZ+n5TXxVFpeF38VEd07LGEb8IZXGUXEU8D8lNKn99q2EPhhSskXOHW4iOgGbAY+mFJ6KHcelZeI6Ak8DXwK+FvghZTS9XlTqZxExE3AuSmls3JnUfmKiJ8Cb6SUrthr271A35TSRfmSqVxExDbg+pTS91sfB7ASuD2l9LXWbV1oKcBfSCndlSvroXLmN5OI6AScAjy8z66HaZmFk3LoTsvrwsbcQVSW7qbly78ZuYOobH0QeCoifhARayPi2Yi4vvWDn9RRfgOcFxHjACLiOFrOzpqeNZXK2UhgEHv1lpTSTuAxiqy3eO1APv2ASmDNPtvXAOd3fBwJaDkN/1ngydxBVF4i4tPAaOBPc2dRWRsFXAvcCvxf4CTg2637vAZdHeUfafky+qWIaKLl8/rXUkrT8sZSGRvU+uf+esuQDs5yRCy/kgCIiFuAdwHvSik15c6j8hERY2m598G7Ukp7cudRWasA5u516dEzETGGlustLb/qKJcClwOXAS/S8iXMbRGxJKX03azJpCLnac/5rAeagIH7bB8IrO74OCpnEXEr8FFgSkppce48Kjtn0nI2zIsR0RgRjcC5wLWtjzvnjacysgp4aZ9tLwPeiFId6RvAN1NKD6aUnk8p3Q/cgje8Uj5vdpOi7y2W30xSSruBecDUfXZNpeWuz1KHiIjb+O/iuyB3HpWl/wJOoGV2482fucCDrX/fnS+ayszjwNh9th0DLMuQReWrKy0TJHtrws/tymcJLSX3970lImqAsymy3uJpz3ndAtwfEXNoecO9GhgM3Jk1lcpGRNwBfJyWm7xsjIg3r+nYllLali+ZyklKaROwae9tEbEd2JBSeiFPKpWpW4EnIuKvgR/QsgzhZ4G/yppK5eYh4H9FxBJaTnueAHwOuC9rKpW01hU/Rrc+rADqIuIkWt6Ll0fEPwF/FRELgFeBrwDbgH/LEvgwudRRZhFxLS3rZtXSsrbWX6SUHsubSuUiIg70AvB3KaWvdmQWaW8RMROXOlIGEXEhLdegjwWW03Kt77eTH5jUQVrXTf0H4GJgAC2n4z8I/H1KqSFnNpWuiJgM7G+1hXtTSle23vX+b4HPAL2Bp4Driu1LasuvJEmSJKnkee2AJEmSJKnkWX4lSZIkSSXP8itJkiRJKnmWX0mSJElSybP8SpIkSZJKnuVXkiRJklTyLL+SJEmSpJJn+ZUkqYRExDci4pe5c0iSVGgsv5IklZbTgDm5Q0iSVGgipZQ7gyRJOkIR0QnYBlTvtfnllNJxmSJJklRQnPmVJKk0NAJntv79dKAWOCtfHEmSCktV7gCSJOnIpZSaI6IW2Ar8LnlqlyRJf8CZX0mSSscE4DmLryRJb2X5lSSpdJwEPJM7hCRJhcjyK0lS6TgRmJ87hCRJhcjyK0lS6agCxkXE4IjolTuMJEmFxPIrSVLp+GvgI8AK4OuZs0iSVFBc51eSJEmSVPKc+ZUkSZIklTzLryRJkiSp5Fl+JUmSJEklz/IrSZIkSSp5ll9JkiRJUsmz/EqSJEmSSp7lV5IkSZJU8iy/kiRJkqSSZ/mVJEmSJJW8/w/Mju9DOIo6WwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"N = 10000 # number of steps\n",
"h = 0.001 # step size\n",
"\n",
"# initial values\n",
"t_0 = 0\n",
"x_0 = 0\n",
"\n",
"t = np.zeros(N+1)\n",
"x = np.zeros(N+1)\n",
"t[0] = t_0\n",
"x[0] = x_0\n",
"\n",
"t_old = t_0\n",
"x_old = x_0\n",
"\n",
"for n in range(N):\n",
" x_new = x_old + h*(np.cos(x_old)+np.sin(t_old)) # Euler's method\n",
" \n",
" t[n+1] = t_old+h\n",
" x[n+1] = x_new\n",
" \n",
" t_old = t_old+h\n",
" x_old = x_new\n",
"\n",
"print(r'x_N = %f' % x_old)\n",
"\n",
"# Plot x(t)\n",
"plt.figure()\n",
"plt.plot(t, x)\n",
"plt.ylabel(r'$x(t)$')\n",
"plt.xlabel(r'$t$')\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the code above, we have chosen $h=0.001$ and $N=10000$,\n",
"and so $t_N=10$.\n",
"In the plot of $x(t)$, the discrete points have been connected by straight lines.\n",
"\n",
"What happens to $x_N$ when we decrease $h$ by a factor of $10$?\n",
"(Remember to increase $N$ simultaneously by a factor of $10$ so as \n",
"to obtain the same value for $t_N$.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Accuracy:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the value of $x_N$ depends on the step size $h$. \n",
"In theory, a higher accuracy of the numerical solution in comparison \n",
"to the exact solution can be achieved by decreasing $h$ since our approximation\n",
"of the derivative $\\frac{d}{dt}x(t)$ becomes more accurate.\n",
"\n",
"However, we cannot decrease $h$ indefinitely since, eventually, we are hitting \n",
"the limits set by the machine precision. Also, lowering $h$ requires more steps \n",
"and, hence, more computational time. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For Euler's method, it turns out that the global error (error at a given $t$) is proportional to the step \n",
"size $h$ while the local error (error per step) is proportional to $h^2$. This is\n",
"called a first-order method. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now summarize __Euler's method__.\n",
"\n",
"Given the ODE\n",
"$$\n",
"\\frac{d}{dt}x(t)=g(x(t),t)~~~\\mathrm{with}~~~x(t_0)=x_0,\n",
"$$\n",
"we can approximate the solution numerically in the following way: \n",
"1. Choose a step size $h$.\n",
"2. Define grid points: $t_n=t_0+n*h~~\\mathrm{with}~~n=0,1,2,3,...,N$.\n",
"3. Compute iteratively the values of the function at these grid points:\n",
"$x_{n+1}=x_n+h*g(x_n,t_n)$. Start with $n=0$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Instability:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Apart from its fairly poor accuracy, the main problem with Euler's method is that it can\n",
"be unstable, i.e. the numerical solution can start to deviate from the exact solution\n",
"in dramatic ways. Usually, this happens when the numerical solution grows large in\n",
"magnitude while the exact solution remains small.\n",
"\n",
"A popular example to demonstrate this feature is the ODE\n",
"$$\n",
"\\frac{dx}{dt}=-x~~~\\mathrm{with}~~~x(0)=1.\n",
"$$\n",
"The exact solution is simply $x(t)=e^{-t}$. It fulfills the ODE and the initial condition.\n",
"\n",
"On the other hand, our Euler method reads\n",
"$$\n",
"x_{n+1}=x_n+h\\cdot(-x_n)=(1-h)x_n.\n",
"$$\n",
"\n",
"Clearly, if $h>1$, $x(t_n)$ will oscillate between negative and positive numbers and\n",
"grow without bounds in magnitude as $t_n$ increases. We know that this is incorrect\n",
"since we know the exact solution in this case.\n",
"\n",
"On the other hand, when $0towards the center\n",
"of a planet of mass $M$. Let us assume that the atmosphere exerts a force\n",
"$$\n",
"F_\\mathrm{drag}=Dv^2\n",
"$$\n",
"onto the particle which is proportional to the square of the velocity. Here, $D$ is the \n",
"drag coefficient. Note that the $x$-axis is pointing away from the planet. \n",
"Hence, we only consider $v\\leq 0$.\n",
"\n",
"The particle motion is described by the following governing equation \n",
"($G$: gravitational constant)\n",
"$$\n",
"m\\frac{d^2 x}{dt^2}=Dv^2-\\frac{GmM}{x^2}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Dividing each side by $m$ gives\n",
"$$\n",
"\\frac{d^2 x}{dt^2}=\\frac{D}{m}v^2-\\frac{GM}{x^2}.\n",
"$$\n",
"Following our recipe above, we re-cast this as two first-order ODEs\n",
"\\begin{eqnarray*}\n",
"\\frac{dx}{dt} &=& v, \\\\\n",
"\\frac{dv}{dt} &=& \\frac{D}{m}v^2-\\frac{GM}{x^2}.\n",
"\\end{eqnarray*} \n",
"We choose $D=0.0025\\,\\mathrm{kg}\\,\\mathrm{m}^{-1}$, $m=1\\,\\mathrm{kg}$ and \n",
"$M=M_\\mathrm{Earth}$, i.e. the mass of the Earth."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Accordingly, our algorithm now reads\n",
"\\begin{eqnarray*}\n",
"x_{n+1} &=& x_n+h*v_n, \\\\\n",
"v_{n+1} &=& v_n+h*\\left[ \\frac{D}{m}v_n^2-\\frac{GM}{x_n^2} \\right] .\n",
"\\end{eqnarray*} \n",
"Let us specify the following initial conditions and step size:\n",
"$$\n",
"t_0=0,~~x(t_0)=x_0=7000.0\\,\\mathrm{km},~~v(t_0)=v_0\n",
"=0\\,\\mathrm{m/s},~~h=0.001\\,\\mathrm{s}.\n",
"$$\n",
"We could now iterate the above equations until the particle hits the ground, i.e. until\n",
"$x=R_\\mathrm{Earth}$, where $R_\\mathrm{Earth}$ is the radius of Earth. This occurs in finite time\n",
"both in reality and in our code."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Moreover, the particle would also reach $x=0$ in finite time, given the above \n",
"equations, while the speed grows to infinity. However, the code would crash well before $x$ approaches zero due to the speed reaching very large values.\n",
"\n",
"__Note__: The governing equation actually changes when $|x|"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA8cAAAF/CAYAAABkCbtrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmcnnV97//XZ2YyWSaZJclkmUlCCISQBQhJUIMgIKJYPZae9mhVrF1OAa2tPa3aVvvrT4+0p/3VutRqBWvVxnK0tT9bRFn0KJUlCGEJSQhrAgkhy0z2STKTWb7nj/ueMIQsk8nMfd3L6/l4XI/c97Xc9zvhIvCe6/p+r0gpIUmSJElSJavKOoAkSZIkSVmzHEuSJEmSKp7lWJIkSZJU8SzHkiRJkqSKZzmWJEmSJFU8y7EkSZIkqeJZjiVJkiRJFc9yLEmSJEmqeJZjSZIkSVLFsxxLkiRJkipeTdYBsjZ58uQ0e/bsrGNIkiRJkobZww8/3J5Sah7MvhVfjmfPns2qVauyjiFJkiRJGmYR8cJg9/W2akmSJElSxbMcS5IkSZIqnuVYkiRJklTxLMeSJEmSpIpnOZYkSZIkVTzLsSRJkiSp4lmOJUmSJEkVz3IsSZIkSap4ZVWOI+KDEbExIjoj4uGIuDTrTJIkSZKk4lc25Tgi3gV8AfgL4ELgfuD2iJiVaTBJkiRJUtErm3IM/AHwjZTSV1NK61NKvwtsBT6QcS5JkiRJUpEri3IcEbXAUuCuozbdBVxc+ETDY2P7AX72dFvWMSRJkiSp7JVFOQYmA9XA9qPWbwemHb1zRFwXEasiYlVbW/GWz5t/9hy/9+1HSSllHUWSJEmSylq5lONTklK6OaW0LKW0rLm5Oes4x7WgpYE9B7vZsudQ1lEkSZIkqayVSzluB3qBqUetnwpsK3yc4bGopR6AtVv2ZZxEkiRJkspbWZTjlNJh4GHgqqM2XUVu1uqSNH96PdVVwRMv7c06iiRJkiSVtZqsAwyjzwIrIuJB4D7gBqAF+EqmqU7DmFHVnNVcx9qXvHIsSZIkSSOpbMpxSuk7ETEJ+FNgOrAW+IWU0gvZJjs9i1oauPfZ9qxjSJIkSVJZK4vbqvullL6cUpqdUhqdUlqaUvpZ1plO18LWBnbs72LH/s6so0iSJElS2SqrclyO+iflWuet1ZIkSZI0YizHRW5Bfzne4qRckiRJkjRSLMdFbsKYUcyeNM7HOUmSJEnSCLIcl4CFrQ2s2+qVY0mSJEkaKZbjErCopYHNuw6x92B31lEkSZIkqSxZjkvAwiOTcnn1WJIkSZJGguW4BPSX47WWY0mSJEkaEZbjEjBp/GhaGsb4OCdJkiRJGiGW4xKxoKWBtT7OSZIkSZJGhOW4RCxqrWdD+wEOdPVkHUWSJEmSyo7luEQsamkgJXhym7dWS5IkSdJwsxyXiEWtDQCs3WI5liRJkqThZjkuEVPrRzOprtZxx5IkSZI0AizHJSIiWNjawFpnrJYkSZKkYWc5LiGLWup5Zvt+unp6s44iSZIkSWXFclxCFrY00NOXeHpbR9ZRJEmSJKmsWI5LyKLWegDWvuS4Y0mSJEkaTpbjEjJr4jgmjKlhneVYkiRJkoaV5biERAQLW+p9nJMkSZIkDTPLcYlZ2NLA+q376OntyzqKJEmSJJUNy3GJWdRaT1dPH8+1Hcg6iiRJkiSVDctxiVnU0gDguGNJkiRJGkaW4xIzp3k8Y0ZVOe5YkiRJkoaR5bjEVFcF86fX+zgnSZIkSRpGluMStKilgSde2kdfX8o6iiRJkiSVBctxCVrUWk9HVw+bdh3MOookSZIklQXLcQlamJ+Uy1urJUmSJGl4WI5L0Nyp4xlVHU7KJUmSJEnDxHJcgkbXVHPO1Ak+zkmSJEmShonluEQtbKln3Uv7SMlJuSRJkiTpdFmOS9Si1gZ2HTjM1r2dWUeRJEmSpJJnOS5RRybl2uKt1ZIkSZJ0uizHJWr+9AlUBax7yUm5JEmSJOl0WY5L1LjaGuY0j3dSLkmSJEkaBpbjEraopd7HOUmSJEnSMLAcl7BFrQ1s29dJe0dX1lEkSZIkqaRZjkvYgpZ6wHHHkiRJknS6LMclzBmrJUmSJGl4WI5LWMPYUcyaOM5JuSRJkiTpNFmOS9yi1npvq5YkSZKk02Q5LnELWxp4YedB9h7qzjqKJEmSJJWsoi/HEXFdRPw0IvZERIqI2cfYpykiVkTE3vyyIiIaC5+28BbmJ+V6wqvHkiRJkjRkRV+OgXHAXcAnT7DPLcAS4Or8sgRYMeLJikD/pFyOO5YkSZKkoavJOsDJpJQ+DxARy461PSLmkyvEl6SUVubXXQ/cExHzUkpPFSxsBponjGZa/RjHHUuSJEnSaSiFK8cnsxzoAO4fsO4+4ABwcSaJCmxhS72Pc5IkSZKk01AO5Xga0JZSSv0r8q935Le9Sn4c86qIWNXW1lagmCNnYWsDz7V1cPBwT9ZRJEmSJKkkZVKOI+LG/ORaJ1ouH6nvTyndnFJallJa1tzcPFJfUzCLWurpS7B+6/6so0iSJElSScpqzPHngW+dZJ9Ng/ysbUBzRET/1eOICGBKflvZO29GblKutVv2svSMpozTSJIkSVLpyaQcp5TagfZh+riVwHhyY4/7xx0vB+p45TjksjWtfgyTx4/m8RcddyxJkiRJQ1H0s1VHxDRyY4fPya9akH+G8aaU0q6U0vqIuAO4KSKuy+9zE3Bbuc9U3S8iOK+1njVb9mQdRZIkSZJKUilMyHUD8Cjwz/n3P8i/f8eAfd4DrAbuzC+rgfcVMGPmzpvRyLM7nJRLkiRJkoai6K8cp5Q+CXzyJPvsBq4tRJ5idX5rA30JnnhpH8tmT8w6jiRJkiSVlFK4cqxB6J+Uy3HHkiRJknTqLMdlYmr9GKbWj2bNFsuxJEmSJJ0qy3EZOa+1kcdfdFIuSZIkSTpVluMycv6MBja0H2B/Z3fWUSRJkiSppFiOy8h5MxpICda9tC/rKJIkSZJUUizHZeS81tykXGuclEuSJEmSTonluIxMHj+a1saxPO6kXJIkSZJ0SizHZWZRaz1rnJRLkiRJkk6J5bjMnD+jked3HmTvISflkiRJkqTBshyXmf5xx+u8tVqSJEmSBs1yXGb6y7HjjiVJkiRp8CzHZaaprpaZE8c6Y7UkSZIknQLLcRk6v7WRx7c4KZckSZIkDZbluAydN6OBzbsOsfvA4ayjSJIkSVJJsByXofPz447XOO5YkiRJkgbFclyGFlqOJUmSJOmUWI7LUMPYUZw5uY7HX3TcsSRJkiQNhuW4TC1qbWDtln1Zx5AkSZKkkmA5LlPntzawZc8h2ju6so4iSZIkSUXPclymzpvhuGNJkiRJGizLcZla2FJPBKx50XIsSZIkSSdjOS5TE8aMYs7kOh63HEuSJEnSSVmOy9j5MxpZ/eIeUkpZR5EkSZKkomY5LmMXzGigbX8X2/Z1Zh1FkiRJkoqa5biMLZ7VBMBjm3zesSRJkiSdiOW4jM2fPoHa6ioee9FyLEmSJEknYjkuY6NrqpnfUu+VY0mSJEk6CctxmVs8o4E1W/bS2+ekXJIkSZJ0PJbjMnfBzEYOHu7l2R0dWUeRJEmSpKJlOS5zi2c2AvDY5t0ZJ5EkSZKk4mU5LnOzJ9VRP6aGxzbvzTqKJEmSJBUty3GZq6oKLpjZyGObnZRLkiRJko7HclwBFs9s5Ont+zl4uCfrKJIkSZJUlCzHFWDxzEZ6+xJrt+zLOookSZIkFSXLcQU4f0ZuUq7V3lotSZIkScdkOa4AzRNG09o4lsdetBxLkiRJ0rFYjivE4lmNPLbJcixJkiRJx2I5rhCLZzSyZc8h2vZ3ZR1FkiRJkoqO5bhCLJ7luGNJkiRJOh7LcYVY1NJAdVWw2nHHkiRJkvQqluMKMba2mnlTJ/CYV44lSZIk6VWKuhxHxMSI+GJEPBkRhyJic0T8fURMOmq/pohYERF788uKiGjMKnexumBmI6s376GvL2UdRZIkSZKKSlGXY6AFaAU+BpwHXAu8AfjfR+13C7AEuDq/LAFWFC5maVg8s4F9nT08v/NA1lEkSZIkqajUZB3gRFJKa4H/OmDVsxHxUeC2iKhPKe2LiPnkCvElKaWVABFxPXBPRMxLKT1V+OTFafHMJgAe27yHOc3jM04jSZIkScWj2K8cH0s90AUczL9fDnQA9w/Y5z7gAHBxYaMVt7OnjKeutppHfd6xJEmSJL1CSZXj/DjiTwNfTSn15FdPA9pSSkcG0uZf78hvO9bnXBcRqyJiVVtb20jHLhrVVcHiWY08sml31lEkSZIkqahkUo4j4saISCdZLj/qmPHA94Et5MYgD1lK6eaU0rKU0rLm5ubT+aiSs2RWE09u28/Bwz0n31mSJEmSKkRWY44/D3zrJPts6n+RL8Y/zL99e0qpc8B+24DmiIj+q8cREcCU/DYNsOSMJnr7Eqs372X5WZNOfoAkSZIkVYBMynFKqR1oH8y+ETEBuB0I4OqUUsdRu6wExpMbe9w/7ng5UMcrxyELWJKflOuRTbstx5IkSZKUV9SzVeeL8V3kJuG6BqiLiLr85l0ppcMppfURcQdwU0Rcl992E3CbM1W/WsO4UZzVXMcjLzjuWJIkSZL6FfuEXEuB1wELgKeBrQOWgTNRvwdYDdyZX1YD7yto0hKy9IwmHtm0mwFzmEmSJElSRSvqK8cppbvJ3U59sv12A9eOeKAysWRWE/+y6kU2th/weceSJEmSRPFfOdYIWHpGbtzxw95aLUmSJEmA5bgindU8nvoxNTyyaU/WUSRJkiSpKFiOK1BVVXDhrCYn5ZIkSZKkPMtxhVoyq4mnd+xnX2d31lEkSZIkKXOW4wq19IwmUoLHvLVakiRJkizHleqCmQ1EwCObvLVakiRJkizHFWrCmFHMmzrBGaslSZIkCctxRVtyRhOPbdpDX1/KOookSZIkZcpyXMGWzmpif1cPz+zoyDqKJEmSJGXKclzBlpzRBDjuWJIkSZIsxxVs9qRxTKyr9XnHkiRJkiqe5biCRQRLZjWxynIsSZIkqcJZjivcRbOb2Nh+gB37O7OOIkmSJEmZsRxXuIvOnAjAque9eixJkiSpclmOK9yilgbGjKriwY27so4iSZIkSZmxHFe42poqLpzZxEPPW44lSZIkVS7LsbjozIms37qP/Z3dWUeRJEmSpExYjsVrZk+kL8Ejm/ZkHUWSJEmSMmE5FhfOaqS6KnjIcceSJEmSKpTlWNSNrmFRSz0POu5YkiRJUoWyHAuAi2ZP5LHNe+jq6c06iiRJkiQVnOVYQG5SrsM9fax5cW/WUSRJkiSp4CzHAnJXjgFvrZYkSZJUkSzHAmBiXS1nTxnvpFySJEmSKpLlWEdcNHsiq17YTW9fyjqKJEmSJBWU5VhHvObMJvZ39vDUtv1ZR5EkSZKkgjrlchwRoyPizIhYEBHNIxFK2egfd/yQ444lSZIkVZhBleOImBARH4iInwF7gWeBtcC2iNgUEV+NiItGMqhG3oymcbQ0jOFBxx1LkiRJqjAnLccR8QfA88BvAj8CfhFYDJwDLAc+CdQAP4qIOyJi7kiF1ch77ZxJ/HzjTlJy3LEkSZKkylEziH1eB1yWUlp7nO0PAv8YETcAvwVcBjwzTPlUYMvnTOJ7j27h2R0dzJ06Ies4kiRJklQQJy3HKaV3DuaDUkpdwJdPO5EytfysSQCs3LDTcixJkiSpYpzShFwRsTYiGkYqjLI3c+I4WhvHsvK5nVlHkSRJkqSCOdXZqhcAo49eGRENEfGl4YmkrC0/axIPbNhJn887liRJklQhBjtb9e0R8SkgATOPscs44PrhDKbsLJ8zid0Hu3lqu887liRJklQZBjMhF8AachNtBfBgROwHVgOPAo8D5wJbRyShCu51/eOOn9vJ/On1GaeRJEmSpJE3qHKcUvoYQER0kXt8Uwu5xzktBt6W/5yPjVBGFVhr41hmTRzHyg07+c1Lzsw6jiRJkiSNuMFeOe43PqXUDTwC3DYCeVQkls+ZxO1rt9Lbl6iuiqzjSJIkSdKIOumY44g4cukwX4xPtG9ExLHGJKvELD9rEvs6e1i/dV/WUSRJkiRpxA1mQq6VEfG1iFh+vB0ioikiPgA8AfzisKVTZvqfd/zABh/pJEmSJKn8Dea26nOBTwA/iIg+4GHgJaATaCL3eKf5wIPA76eU7hyhrCqgqfVjmDO5jpXP7eS/Xzon6ziSJEmSNKJOeuU4pbQnpfRRoBW4AVgPNAJnAj3AN4ELU0qvtxiXl9fOmcSDG3fR09uXdRRJkiRJGlGDes4xQErpEPB0Sun3U0q/lFK6OqV0bUrpb1JKa0cqYER8NSKei4hDEdEWEf8REfOP2qcpIlZExN78siIiGkcqU6VYftYk9nf1sO4lxx1LkiRJKm+DLsd5D0XENyJixoikObZVwK+Tu3X7LeSetfzjiBg1YJ9bgCXA1fllCbCigBnL0vI5uXHH9z/nuGNJkiRJ5e1Uy/FicrdUPx0R/18hrs6mlG5KKd2TUno+pfQI8KfknrM8ByB/Fflq4LqU0sqU0krgeuDtETFvpPOVs+YJozl32gTufbYt6yiSJEmSNKJOqRynlNanlK4B3gQsB56LiI9ExOgRSXeUiKgDfgPYBDyfX70c6ADuH7DrfcAB4OJC5Cpnl86dzEMbd3PocG/WUSRJkiRpxJzqlWMAUkr3p5QuBX6T3C3PT0fErw1nsIEi4oMR0UGuBL8VuDKl1JXfPA1oSymlAfkSsCO/7Vifd11ErIqIVW1tXhU9kUvmNnO4t48Hn9+VdRRJkiRJGjFDKscD/Az4HeBF4OuDPSgiboyIdJLl8gGH/DNwIXAZ8DTwrxExbqihU0o3p5SWpZSWNTc3D/VjKsJrZk+ktrqKe5/xhwiSJEmSytdgnnN8RER8kNxzjfuXZqAP2Aj8+yl81OeBb51kn039L1JKe4G9wDMR8QCwG/hlcpNubQOaIyL6rx5HRABT8tt0GsbWVrNsdhP3PNOedRRJkiRJGjGnVI6BPwPWAI+RK6ZrgHX5xzwNWkqpHRhq24r80j/OeSUwntzY4/5xx8uBOl45DllDdOncZv7qjifZsb+TKRPGZB1HkiRJkobdKZXjlNIxx/COlIg4m9wV4h8DbcAM4I+BLuC2fKb1EXEHcFNEXJc/9CbgtpTSU4XMW64unTuZv7oD7nu2nV+6sJBP8ZIkSZKkwjjdMccjrQu4HLgdeBb4DrAfWJ5SGnjL9HuA1cCd+WU18L6CJi1jC6bXM7Gu1lurJUmSJJWtU72tuqBSSpvJzU59sv12A9eOfKLKVFUVXHzWJO59pp2UErkh3ZIkSZJUPor9yrGKxBvmNrNjfxdPb+/IOookSZIkDTvLsQblkrmTAbjHRzpJkiRJKkOWYw1KS+NY5jTXce+zjjuWJEmSVH4sxxq0N8xt5ucbdtHZ3Zt1FEmSJEkaVpZjDdpl85o51N3Lgxt3ZR1FkiRJkoaV5ViDtnzOJEbXVPGTJ3dkHUWSJEmShpXlWIM2ZlQ1F581ibufshxLkiRJKi+WY52SN547hed3HmRj+4Gso0iSJEnSsLEc65RcPm8KAD/11mpJkiRJZcRyrFMyc+I4zp4ynp96a7UkSZKkMmI51im7Yl7ukU4HD/dkHUWSJEmShoXlWKfsinlTONzbx33P7sw6iiRJkiQNC8uxTtmy2ROpq6321mpJkiRJZcNyrFNWW1PFJXMnc/eTO0gpZR1HkiRJkk6b5VhDcsW8Kby0t5Ont3dkHUWSJEmSTpvlWEPS/0inH6/fnnESSZIkSTp9lmMNybSGMVwwo4G7nrAcS5IkSSp9lmMN2VULprJ68x627+vMOookSZIknRbLsYbszQunAfAjrx5LkiRJKnGWYw3Z3CnjmT1pnLdWS5IkSSp5lmMNWURw1YKprHyunf2d3VnHkSRJkqQhsxzrtLx54TS6exN3P9WWdRRJkiRJGjLLsU7LkllNTKqr9dZqSZIkSSXNcqzTUl0VXDl/Cnc/uYPDPX1Zx5EkSZKkIbEc67S9ecE09nf18MCGnVlHkSRJkqQhsRzrtF0ydzJjR1Vz57ptWUeRJEmSpCGxHOu0jRlVzeXzmrlz3XZ6+1LWcSRJkiTplFmONSzedv502ju6eHDjrqyjSJIkSdIpsxxrWLzx3CmMGVXFD9dszTqKJEmSJJ0yy7GGxbjaGt547hRuX7vVW6slSZIklRzLsYbN285rob3jMD/f6KzVkiRJkkqL5VjD5opzmxk7qpofPO6t1ZIkSZJKi+VYw2ZcbQ1vnD+FO9dto6e3L+s4kiRJkjRolmMNq7edN532jsPOWi1JkiSppFiONayumDeFsaOquc1ZqyVJkiSVEMuxhtXY2mqunD+FO9Zuo9tbqyVJkiSVCMuxht0vLm5l14HD3PNMW9ZRJEmSJGlQLMcadped00zjuFF879GXso4iSZIkSYNiOdawq62p4u3nT+euddvY39mddRxJkiRJOinLsUbEL104g66ePu5Yuy3rKJIkSZJ0UpZjjYglsxo5Y9I4/v2xLVlHkSRJkqSTKplyHDm3R0SKiF85altTRKyIiL35ZUVENGaVVRARXLO4lfuf28nWvYeyjiNJkiRJJ1Qy5Rj4Q+B4zwa6BVgCXJ1flgArCpRLx3HNha2kBLc+5sRckiRJkopbSZTjiLgI+DDwG8fYNp9cIb4upbQypbQSuB54e0TMK2xSDXTm5DoWz2zke496a7UkSZKk4lb05TgiJpC7MnxdSmnHMXZZDnQA9w9Ydx9wALh45BPqRP7rklae3LaftVv2Zh1FkiRJko6r6Msx8BXgjpTS7cfZPg1oSyml/hX51zvy214lIq6LiFURsaqtrW3YA+tl77ighdqaKv5l1easo0iSJEnScWVSjiPixvzEWidaLo+I9wEXAB8dzu9PKd2cUlqWUlrW3Nw8nB+tozSOq+Wti6bxvUe30Nndm3UcSZIkSTqmrK4cfx6Yf5LlQeBKYAHQERE9EdGTP/47EXFv/vU2oDkiov/D86+n5LcpY++6aCb7O3t85rEkSZKkolWTxZemlNqB9pPtFxGfAD5z1Oo1wEeA/8i/XwmMJzf2uH/c8XKgjleOQ1ZGXnfmJM6YNI5vP7SJay5szTqOJEmSJL1KUY85TiltSSmtHbjkN21OKW3I77MeuAO4KSKWR8Ry4CbgtpTSUxlF1wBVVcE7l83kgQ272Nh+IOs4kiRJkvQqRV2OT8F7gNXAnfllNfC+TBPpFX5l6QyqAifmkiRJklSUMrmt+nSklOIY63YD12YQR4M0tX4Mbzx3Ct99+EX+8KpzqKkul5/LSJIkSSoHNhQVzLsumkXb/i5+vH571lEkSZIk6RUsxyqYK+Y109o4ln9a+ULWUSRJkiTpFSzHKpia6ire+7pZ3P/cTp7Zvj/rOJIkSZJ0hOVYBfWuZTOprany6rEkSZKkomI5VkFNGj+a/3J+C//2yIvs6+zOOo4kSZIkAZZjZeD9F5/BwcO9/P8Pv5h1FEmSJEkCLMfKwPkzGrlwViP/tPIF+vpS1nEkSZIkyXKsbLx/+Ww2tB/gnmfbs44iSZIkSZZjZeOt501jyoTR/MM9G7KOIkmSJEmWY2VjdE01v/H6M7nnmXbWbtmbdRxJkiRJFc5yrMy857WzqKut5qtePZYkSZKUMcuxMtMwdhTvfs0sbnt8Ky/uPph1HEmSJEkVzHKsTP3mJWcSwD/e+3zWUSRJkiRVMMuxMtXSOJb/ckEL335oE3sPdmcdR5IkSVKFshwrc9e9YQ4HD/fyjfufzzqKJEmSpAplOVbm5k+v503zp/CP921kX6dXjyVJkiQVnuVYReHDV57D3kPdfPO+57OOIkmSJKkCWY5VFM6b0cCV507hH+7dyH6vHkuSJEkqMMuxisaH3zSXvYe6+aeVL2QdRZIkSVKFsRyraJw/o5E3njuFr96zgY6unqzjSJIkSaoglmMVlQ9fOZc9B7v5+r0bs44iSZIkqYJYjlVULpjZyJsXTOWmn21gZ0dX1nEkSZIkVQjLsYrOx66ex8HDPXzxJ89mHUWSJElShbAcq+icPWUC77poJv/88xfYtPNg1nEkSZIkVQDLsYrS77/pHKqrgs/c9VTWUSRJkiRVAMuxitLU+jH81iVncuvql1jz4t6s40iSJEkqc5ZjFa3rLzuLiXW1fPq2J0gpZR1HkiRJUhmzHKto1Y8ZxcfeMo8Hn9/FratfyjqOJEmSpDJmOVZRe+eymZw/o4G/+OF6DnT1ZB1HkiRJUpmyHKuoVVUFn3rHQrbv6/LRTpIkSZJGjOVYRe/CWU38ytIZfO3eDWxo68g6jiRJkqQyZDlWSfijq89lTE01f/rva52cS5IkSdKwsxyrJDRPGM2f/MJ87n9uJ/+yanPWcSRJkiSVGcuxSsavXjST1545kRt/sJ4d+zqzjiNJkiSpjFiOVTKqqoK//OXzOdzTx5/9x7qs40iSJEkqI5ZjlZQzJ9fx+286hzvWbeOHa7ZmHUeSJElSmbAcq+T89qVncv6MBj7+vTVs9/ZqSZIkScPAcqySU1NdxefetZjO7l4+8q+r6etz9mpJkiRJp8dyrJJ0VvN4/vRtC7jnmXa+ufL5rONIkiRJKnGWY5Ws9752Fm+aP4X/dfuTPLVtf9ZxJEmSJJUwy7FKVkRu9ur6MaP44D8/TEdXT9aRJEmSJJWooi/HEXF3RKSjlm8ftU9TRKyIiL35ZUVENGaVWYUzefxo/vbdi9nYfoA//rfHScnxx5IkSZJOXdGX47yvA9MHLNcftf0WYAlwdX5ZAqwoZEBl5+KzJvPRt5zLbY9v5Rv3P591HEmSJEklqCbrAIN0MKW07VgbImI+uUJ8SUppZX7d9cA9ETEvpfRUAXMqIzdcNodHNu3mz3+wnvNaG1g2e2LWkSRJkiSVkFK5cvyrEdEeEesi4jMRMWHAtuVAB3D/gHX3AQeAiwsZUtmJCD7z3y6gtWksN3zrYTbvOph1JEmSJEklpBTK8S3Ae4ErgE8Dvwz824Dt04C2NGCwaf71jvy2V4mI6yJiVUSsamtrG7FQnwuIAAASPElEQVTgKqyGsaP42vsvoqunj9/65kPs6+zOOpIkSZKkEpFJOY6IG48xydbRy+UAKaWbU0p3ppTWpJS+DbwLuCoilgz1+/OfuSyltKy5uXmYflcqBmdPGc9Xrl3KhrYDfOiWR+np7cs6kiRJkqQSkNWV488D80+yPHicY1cBvcDc/PttQHNERP8O+ddT8ttUYV5/9mQ+fc0ifvZ0G3926zpnsJYkSZJ0UplMyJVSagfah3j4eUA1sDX/fiUwntzY4/5xx8uBOl45DlkV5N2vmcWmXQf5+7ufo2ncKD76lnOzjiRJkiSpiBX1bNURcRa58cY/JFemFwB/AzxKbtItUkrrI+IO4KaIuC5/6E3Abc5UXdk+9pZ57DnYzZd++hwTxozihsvOyjqSJEmSpCJV1OUYOAxcCXyY3NXhzcAPgE+llHoH7Pce4IvAnfn3twIfKmBOFaGI4MZrFrG/s5u/vP1JJoyp4b2vPSPrWJIkSZKKUFGX45TSZuCyQey3G7h25BOp1FRXBZ9952IOdPXwie+tpa8v8b7ls7OOJUmSJKnIlMKjnKTTUltTxd9fu5Q3zZ/C//Mf6/iHezZkHUmSJElSkbEcqyKMGVXNl9+7lLcumsaNP1jPl376rLNYS5IkSTrCcqyKUVtTxRfffSHXLG7hr+98ik/euo7ePguyJEmSpCIfcywNt5rqKj77zsVMqR/DzT/bwJY9nfztuxczrtZ/FSRJkqRK5pVjVZyqquDjvzCfT71jIT95cjvvvvkBduzrzDqWJEmSpAxZjlWx3n/xbG563zKe3t7B2754Lw89vyvrSJIkSZIyYjlWRbtqwVS+9zsXU1dbzbtvfoCv37fRibokSZKkCmQ5VsU7d1o9t/7uJVw+bwqf+v4TfOiWR9lz8HDWsSRJkiQVkOVYAurHjOLm9y3lj64+lzvXbePqz9/Dfc+2Zx1LkiRJUoFYjqW8qqrgA5efxfc++HrqRlfz3n/4OZ++7QkOHe7NOpokSZKkEWY5lo5y3owGbvvdS/m15WfwtXs3ctXn/pOfPrUj61iSJEmSRpDlWDqGsbXV/M9fXMR3rnsdo2uq+I2vP8SHbnnERz5JkiRJZcpyLJ3Aa+dM4ocfvpQ/uOoc7lq3ncs/czdf+PEzHDzck3U0SZIkScPIciydxOiaan7vyrnc9T/ewGXnNPO5Hz/NFZ+5m39ZtZnePh/7JEmSJJUDy7E0SLMn1/H31y7luzcsZ3rDWD723ce56rP/yXcffpHu3r6s40mSJEk6DZZj6RQtmz2R733wYr5y7VLGjKrmI/+6mjf+zd387wc30dntzNaSJElSKYqUKvu20GXLlqVVq1ZlHUMlKqXET57cwd/+n2dY/eJeJtbV8u7XzOTa153B9IaxWceTJEmSKlpEPJxSWjaofS3HlmOdvpQSK5/byTfuf54fr99ORHD1wmm866KZvP7syVRXRdYRJUmSpIpzKuW4ZqTDSJUgIrj47MlcfPZkNu86yLceeIFvP7SZH6zZyrT6MVxzYSu/srSVs6dMyDqqJEmSpGPwyrFXjjVCunp6+T/rd/Ddh1/kP59uo7cvsWB6PW9dNI2rF03j7CnjifCKsiRJkjRSvK36FFiOVQg79ndy62MvcfvabTz8wm4A5kyu480Lp3H5vGaWzGqitsb58SRJkqThZDk+BZZjFdr2fZ3c9cR27ly7jQc27KSnLzGutprXzZnEJWdP5tK5k72qLEmSJA0Dy/EpsBwrS/s6u3nguZ3c80w79z7bzsb2AwA0jhvFkllNLD0jt1wwo5GxtdUZp5UkSZJKixNySSWifswo3rxwGm9eOA2AzbsOsvK5nTz8wm5WvbCLnzy5A4CaqmDu1AksmF7PwpZ6FrTUM396PQ1jR2UZX5IkSSobXjn2yrGK2O4Dh3lk024efmE3617axxNb99G2v+vI9hlNYzl7ynjmTB7PnOa63DJ5PFPrR3tbtiRJkiqeV46lMtFUV8uV86dy5fypR9bt2N/JE/mivH7rfp7b0cHPN+ziUHfvkX3qaquZNamO1saxtDaOobVpLK2N4/K/jmXy+FrLsyRJkjSA5VgqMVMmjGHKvDFcPm/KkXV9fYlt+zrZ0HaADe0dbGg7wKZdB9m86yAPbNhJR1fPKz6jtrqKyeNraZ4wmsnjR9M84eVl8vjc0jB2FI3jRtEwdhRjRjneWZIkSeXNciyVgaqqoKVxLC2NY7lk7uRXbd97qJstuw+xZc8htuw+yNZ9nbTvP0xbRxdb93by+Ja97Ozoou84oyxG11QdKcqNY2upH5t7PWFMDXWjqxlXW0NdbTXjRtdQV5tbVze6hnG11fn3NYytraa2uopR1eFVa0mSJBUdy7FUARryZXZBS/1x9+ntS+w6cJj2ji7aO7rYe6ibvYe62XOwm335X/ce6mbPocNs2XOIJ17aS0dXDwcO99J7vFZ9DBG5sj26pjr366iXX9fWVB21rZpRVUFNdVBTXUVNVVBTVZV7XzVgXf/7qlz5rs7v0/96VFVQVRVURVAVvPJ1BBFQHf37QETk3ue3VUVQVZXbJwYcV101YHt+PQFBHPm9Bhz5YUAcWZf/4UC8vA/5/fr3yW0OBv4cof/Yl7e/fMyR7f7gQZIkaUgsx5IAqK6KI7dWn4qUEod7+zjY1cuBwz0cyP/a//7g4R46unrp6u6lq6fv5V97+ujq6aWre8Drnj66uvvYc6j7yH7dvX309CZ6+hI9ff2vX16n4xtUOT/J8cfddsIjT3bsiY47/tYTfuNJfiYw5O8cgd/HyY8d2pFDzXry7zzeZ576Qaf6PUP5UU+hfkA0pD+zIv1zzn3PEI4ZwhcN6Z9Omf1+ivk8GIqi/nMr4j/rof27UJx/1tcsbuX9F88ewjcVH8uxpNMSEfkrvdU01dUW9LtTSvT29RfnRE9vX/7XRHdvX35bH929ib6USCl3hbwvJfpS7vi+/Lr+130p0Zvy7/s4sm/u15ePy30O+c8d8DkvhyPlfjmS9RXv8+sG7E4iHbU9t46jPufl/Qd+Xjqy7pjff9T2gZ997D/cIW16RcZjbxva5574uBMnGupDGU74+ziN7ztR3pH48znZP7Gh/PkM6ZiTnjnD8R1DOKYAv5f8QYU45ITn7fB+zxCOGdL3FOb3M7R/PkPIVqT/vg39e4ZwTIH+AZXbOVrMv5/amqohfFNxshxLKlkR/bdcZ51EkiRJpa58ar4kSZIkSUNkOZYkSZIkVTzLsSRJkiSp4lmOJUmSJEkVz3IsSZIkSap4lmNJkiRJUsWzHEuSJEmSKp7lWJIkSZJU8SzHkiRJkqSKVxLlOCJeExE/ioiOiNgfEfdHxOQB25siYkVE7M0vKyKiMcvMkiRJkqTSUfTlOCJeC9wF3A28DlgKfAboHrDbLcAS4Or8sgRYUdCgkiRJkqSSVZN1gEH4HPCllNKfD1j3dP+LiJhPrhBfklJamV93PXBPRMxLKT1V0LSSJEmSpJJT1FeOI2IKsBzYGhH3RsSOiLgnIq4csNtyoAO4f8C6+4ADwMWFSytJkiRJKlXFfuV4Tv7XTwEfBR4F/htwZ0QsTSmtBqYBbSml1H9QSilFxI78tleJiOuA6/JvOyKimK8uTwbasw6hiud5qGLhuahi4HmoYuB5qGJR7OfiGYPdMZNyHBE3Ap84yW5XAIfzr29KKf1j/vWjEXEFcAPwgaF8f0rpZuDmoRxbaBGxKqW0LOscqmyehyoWnosqBp6HKgaehyoW5XQuZnXl+PPAt06yzyZgav71E0dtewKYlX+9DWiOiOi/ehwRAUzJb5MkSZIk6YQyKccppXYGcek9Ip4HXgLmHbXpHGBN/vVKYDy5scf9446XA3W8chyyJEmSJEnHVNRjjvNjh/8a+FREPE5uzPE7yT3S6UP5fdZHxB3ATfmxxAA3AbeVyUzVJXH7t8qe56GKheeiioHnoYqB56GKRdmcizFgHquiFRF/BPwOMAlYB3w8pfTjAdubgC8C78ivuhX4UEppT6GzSpIkSZJKT0mUY0mSJEmSRlJRP+dYkiRJkqRCsBwXsYj4YERsjIjOiHg4Ii7NOpPKV0T8SUQ8FBH7IqItIr4fEYuO2ici4pMR8VJEHIqIuyNiYVaZVf7y52WKiL8bsM7zUCMuIqZHxDfzfx92RsQTEXHZgO2ehxpxEVEdEZ8e8P+DGyPixoioGbCP56KGVUS8ISJujYgt+f8G//pR2096zkVEU0SsiIi9+WVFRDQW9DcyBJbjIhUR7wK+APwFcCG5mbdvj4hZJzxQGrrLgS8DFwNvBHqAH0fExAH7fAz4Q+B3gYuAHcCPImJCYaOqEkTE64DrgMeP2uR5qBGV/x+4+4AA3gbMJ3e+7Riwm+ehCqF/3p3fA84FPpx//ycD9vFc1HAbD6wld74dOsb2wZxztwBLgKvzyxJgxQhmHhaOOS5SEfFz4PGU0m8PWPcM8N2U0p8c/0hpeETEeGAvcE1K6fv554e/BPxdSunP8/uMJfcX4kdSSjdll1blJiIagEeA/w78v8DalNKHPA9VCBHxF8BlKaXXH2e756EKIiJuA3amlN4/YN03gUkppbd7LmqkRUQHuYmOv5F/f9JzLiLmA08Al6SU7svvcwlwD3BuMT9RyCvHRSgiaoGlwF1HbbqL3FU9qRAmkPs7Ynf+/ZnANAaclymlQ8DP8LzU8LuZ3A8Df3rUes9DFcI1wM8j4jsRsSMiHouI/h/OgOehCude4IqIOBcgIhaQu7vrh/ntnosqtMGcc8uBDnJ3vva7DzhAkZ+XRf2c4wo2GagGth+1fjvwpsLHUYX6AvAYsDL/flr+12Odl62FCqXyFxG/DZwNXHuMzZ6HKoQ5wAeBzwF/CSwm98hIgL/D81CF81fkflj9RET0kvt/9z9PKX05v91zUYU2mHNuGtCWBtyinFJKEbFjwPFFyXIs6VUi4rPAJeRuh+nNOo8qR0TMIzfXwiUppe6s86hiVQGrBgxjejQi5pIb6/l3xz9MGnbvAn4NeA+wjtwPar4QERtTSl/LNJlUhrytuji1A73A1KPWTwW2FT6OKklEfA54N/DGlNKGAZv6zz3PS42k5eTunlkXET0R0QNcBnww/3pnfj/PQ42kreTGyw20HuifFNO/D1Uofw18JqX07ZTSmpTSCuCzvDwhl+eiCm0w59w2oHnAUJT+scpTKPLz0nJchFJKh4GHgauO2nQVr7x3XxpWEfEFXi7GTx61eSO5v9CuGrD/GOBSPC81fP4dOI/c1ZH+ZRXw7fzrp/E81Mi7D5h31LpzgBfyr/37UIUyjtwFk4F6efn/4T0XVWiDOedWkpvxevmA45YDdRT5eelt1cXrs8CKiHiQ3H+kbwBagK9kmkplKyK+BLyP3EQ0uyOif0xIR0qpIz9W5PPAxyPiSXIl5U/JTbhwSyahVXZSSnuAPQPXRcQBYFdKaW3+veehRtrngPsj4hPAd8g9UvH3gI/DkbFznocqhO8DfxwRG8ndVn0h8AfAP4HnokZG/oklZ+ffVgGzImIxuf8WbzrZOZdSWh8RdwA3RcR1+c+5CbitmGeqBh/lVNQi4oPkniM2ndyzxv5HSuln2aZSuYqI4/1l8KmU0ifz+wS5x+pcDzQBPwd+p7+0SCMhIu4m/yin/HvPQ424iHgbufHv84BN5MYaf7F/ghnPQxVC/rmxnwZ+idwtqVvJ3UnzP1NKnfl9PBc1rCLicuDop0UAfDOl9OuDOeciooncRIbvyK+6ldwjofZQxCzHkiRJkqSK55hjSZIkSVLFsxxLkiRJkiqe5ViSJEmSVPEsx5IkSZKkimc5liRJkiRVPMuxJEmSJKniWY4lSZIkSRXPcixJUoWIiL+OiDuzziFJUjGyHEuSVDleAzyYdQhJkopRpJSyziBJkkZQRNQCHcCoAavXp5QWZBRJkqSi45VjSZLKXw+wPP/6tcB04PXZxZEkqfjUZB1AkiSNrJRSX0RMB/YDDyVvG5Mk6VW8cixJUmW4EFhtMZYk6dgsx5IkVYbFwKNZh5AkqVhZjiVJqgwXAI9nHUKSpGJlOZYkqTLUAOdGREtENGYdRpKkYmM5liSpMnwC+FXgReB/ZZxFkqSi43OOJUmSJEkVzyvHkiRJkqSKZzmWJEmSJFU8y7EkSZIkqeJZjiVJkiRJFc9yLEmSJEmqeJZjSZIkSVLFsxxLkiRJkiqe5ViSJEmSVPEsx5IkSZKkivd/AfiYWed6pZ6dAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"G = 6.67300e-11 # gravitational constant\n",
"M = 5.97219e24 # mass of Earth\n",
"D = 0.0025 # drag coefficient \n",
"m = 1.0 # mass of particle\n",
"\n",
"N = 100000 # number of steps\n",
"h = 0.001 # step size\n",
"\n",
"# initial values\n",
"t_0 = 0.0\n",
"x_0 = 7.0e6\n",
"v_0 = 0.0\n",
"\n",
"t = np.zeros(N+1)\n",
"x = np.zeros(N+1)\n",
"v = np.zeros(N+1)\n",
"t[0] = t_0\n",
"x[0] = x_0\n",
"v[0] = v_0\n",
"\n",
"for n in range(N):\n",
" # Euler's method\n",
" x_new = x[n] + h*v[n]\n",
" v_new = v[n] + h*((D/m)*v[n]**2-G*m*M/x[n]**2)\n",
" \n",
" t[n+1] = t[n] + h\n",
" x[n+1] = x_new\n",
" v[n+1] = v_new\n",
"\n",
"plt.figure()\n",
"plt.plot(t,x) # plotting the position vs. time: x(t)\n",
"plt.xlabel(r'$t$')\n",
"plt.ylabel(r'$x(t)$')\n",
"plt.grid()\n",
"\n",
"plt.figure()\n",
"plt.plot(t,v) # plotting the velocity vs. time: v(t)\n",
"plt.xlabel(r'$t$')\n",
"plt.ylabel(r'$v(t)$')\n",
"plt.grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What do you observe? What happens if you choose $x_0=10000.0\\,\\mathrm{km}$ as your \n",
"initial height?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Summary:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have learned how to use Euler’s method to solve first-order and second-order ODEs numerically. It is the simplest method in this context and very easy to implement.\n",
"However:\n",
"1. There are more precise methods for solving ODEs.\n",
"2. There are more stable methods for solving ODEs.\n",
"\n",
"Notwithstanding these issues, Euler’s method is often useful:\n",
"it is easy to use (the coding is less error prone); it can provide a helpful first impression of the solution; modern-day computer power makes computational expense less of an issue.\n",
"Simply put, sometimes it is sufficient."
]
}
],
"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.5.4"
}
},
"nbformat": 4,
"nbformat_minor": 1
}