{ "cells": [ { "cell_type": "markdown", "id": "33376f75", "metadata": {}, "source": [ "(higher-order-odes-section)=\n", "\n", "# Higher-order ODEs\n", "\n", "The {prf:ref}`Euler method` solves a first-order ODE of the form $y' = f(t, y)$, so what happens when we want to solve a higher order ODE? The solution is we write a higher order ODE as system of first-order ODEs and apply the Euler method to solve the system. For example consider the $N$-th order ODE\n", "\n", "$$ y^{(N)} = f(t, y, y', y'' ,\\ldots ,y^{(N-1)}). $$\n", "\n", "If we introduce functions $y_1, y_2, \\ldots, y_N$ where $y_1=y$, $y_2 =y'$, $y_3 =y''$ and so on then\n", "\n", "\n", "$$ \\begin{align*}\n", " y_1' &= y_2,\\\\\n", " y_2' &= y_3,\\\\\n", " &\\vdots \\\\\n", " y_N' &= f(t, y_1 , y_2 , y_3 , \\ldots, y_{N-1}).\n", "\\end{align*} $$\n", "\n", "This is a system of $N$ first-order ODEs, so we can apply the Euler method to solve the system to give an equivalent solution to the $N$-th order ODE.\n", "\n", "`````{prf:example}\n", ":label: spring-example\n", "\n", "A mass-spring-damper model consists of objects connected via springs and dampers. A simple example of the application of a model is a single object connected to a surface is shown in {numref}`spring-model-figure`.\n", "\n", "```{figure} ../_images/01_mass-spring-damper.svg\n", ":width: 300\n", ":name: spring-model-figure\n", "\n", "The mass-spring-damper model {cite:p}`wikipedia:2008`.\n", "```\n", "\n", "The displacement of the object over time can be modelled by the second-order ODE\n", "\n", "$$ m \\ddot{y} + c \\dot{y} + k y = 0, $$\n", "\n", "where $y$ represents the displacement of e object, $\\dot{y}$ represents the change in the displacement over time (i.e., velocity), $m$ represents the mass of the object, $c$ is the damping coefficient and $k$ is a spring constant based on the length and elasticity of the spring.\n", "\n", "An object of mass 1kg is connected to a dampened spring with $c = 2$ and $k = 4$. The object is displaced by 1m and then released. Use the Euler method to compute the displacement of the object over the first 5 seconds after it was released.\n", "\n", "````{dropdown} Solution (click to show)\n", "We first need to rewrite the governing equation as a system of first-order ODEs. Let $y_1 = y$ and $y_2 = \\dot{y}$ then we have the system of two first-order ODEs\n", "\n", "$$ \\begin{align*}\n", " \\dot{y}_1 &= y_2, \\\\\n", " \\dot{y}_2 &= \\frac{1}{m}(- c y_2 - k y_1).\n", "\\end{align*} $$\n", "\n", "Since $m = 1$, $c = 2$ and $k = 4$ then writing the system in vector form $\\dot{\\mathbf{y}} = \\mathbf{f}(t, \\mathbf{y})$ gives\n", "\n", "$$ \\begin{align*}\n", " \\mathbf{y} &= (y_1, y_2), \\\\\n", " \\mathbf{f}(t, \\mathbf{y}) &= (y_2, \\quad -2y_2 - 4y_1).\n", "\\end{align*} $$\n", "\n", "The initial conditions are $y_1 = 1$ (displacement) and $y_2 = 0$ (velocity of the object). Using the Euler method with a step length of $h = 0.1$\n", "\n", "$$ \\begin{align*}\n", " \\mathbf{y}_1 &= \\mathbf{y}_0 + h \\mathbf{f}(t_0, \\mathbf{y}_0), \\\\\n", " &= (1, 0) + 0.1 (0, \\quad -2(0) - 4(1)) \\\\\n", " &= (1, -0.04), \\\\\n", " t_1 &= t_0 + h = 0 + 0.1 = 0.1, \\\\\n", " \\\\\n", " \\mathbf{y}_2 &= \\mathbf{y}_1 + h \\mathbf{f}(t_1, \\mathbf{y}_1), \\\\\n", " &= (1, -0.04) + 0.1 (-0.04, \\quad -2(-0.04) - 4(1)) \\\\\n", " &= (0.96, -0.72), \\\\\n", " t_2 &= t_1 + h = 0.1 + 0.1 = 0.2, \\\\\n", " \\\\\n", " \\mathbf{y}_3 &= \\mathbf{y}_2 + h \\mathbf{f}(t_2, \\mathbf{y}_2), \\\\\n", " &= (0.96, -0.72) + 0.1 (-0.72, \\quad -2(-0.72) - 4(0.96)) \\\\\n", " &= (0.888, -0.96), \\\\\n", " t_3 &= t_2 + h = 0.2 + 0.1 = 0.3, \\\\\n", " & \\vdots\n", "\\end{align*} $$\n", "\n", "A plot of the displacement of the object in the first 5 seconds of being released is shown in \n", "\n", "```{glue:figure} spring_plot\n", ":name: spring-plot-figure\n", ":figwidth: 600\n", ":alt: Euler method solutions for the mass-spring-damper model.\n", "\n", "Euler method solutions for the mass-spring-damper model with $m = 1$, $c = 2$, $k = 4$ and initial displacement 1m.\n", "```\n", "\n", "````\n", "\n", "\n", "`````" ] }, { "cell_type": "code", "execution_count": 61, "id": "69258ebf", "metadata": { "tags": [ "remove-cell" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGxCAYAAACDV6ltAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABLj0lEQVR4nO3deViU9f7/8ecAgiu4oLjhlpm7Ji6hmZqKqXly6WhZaWX+8pR5zExFT2Zmh9Ni2aZZalbH1Ba1zTTMXBItFzyVe2nhAuIKqIgC8/vj8x0UQWSQ4WZmXo/ruq+ZueeemTdozcvParPb7XZEREREvJiP1QWIiIiIWE2BSERERLyeApGIiIh4PQUiERER8XoKRCIiIuL1FIhERETE6ykQiYiIiNfzs7oAd5CZmcmRI0coV64cNpvN6nJEREQkH+x2OykpKVSvXh0fn7zbgBSI8uHIkSOEhoZaXYaIiIgUwMGDB6lZs2ae1ygQ5UO5cuUA8wsNDAy0uBoRERHJj+TkZEJDQ7O+x/OiQJQPjm6ywMBABSIRERE3k5/hLhpULSIiIl5PgUhERES8ngKRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEImIiIjXUyASERERr6dAJCIiIl5PgUhERES8ntsFonXr1tGnTx+qV6+OzWZj2bJl13zN2rVrCQsLo2TJktSrV4933nnH9YWKiIiI23C7QHT27FlatGjBW2+9la/rDxw4QK9evejYsSOxsbFMnDiRUaNG8fnnn7u4UhEREXEXbre5a8+ePenZs2e+r3/nnXeoVasWM2bMAKBRo0Zs2bKFV155hQEDBrioyvxJSoLTp699XcmSEBLi8nJERES8ltsFImdt3LiRiIiIbOd69OjB3LlzuXjxIiVKlMjxmrS0NNLS0rIeJycnu6S2WbMgMjJ/144eDa++CvnYsFdERESc5HZdZs5KSEgg5IrmlZCQENLT0zl+/Hiur4mKiiIoKCjrCA0NdUltfn6m9edaB8CMGZDPXkIRERFxkscHIgDbFc0qdrs91/MOkZGRJCUlZR0HDx50SV1jx0Jq6rWPl14y148eDd9955JSREREvJrHB6KqVauSkJCQ7VxiYiJ+fn5UqlQp19cEBAQQGBiY7bDS2LEwdChkZsLAgbBnj6XliIiIeByPD0Th4eFER0dnO/fdd9/RunXrXMcPFUc2G8yeDe3bm4HYffrAyZNWVyUiIuI53C4QnTlzhu3bt7N9+3bATKvfvn07cXFxgOnuGjJkSNb1I0aM4K+//mLMmDHs2rWLefPmMXfuXMaOHWtF+QUWEABLl0KtWrBvn2kpunjR6qpEREQ8g9sFoi1btnDzzTdz8803AzBmzBhuvvlmJk+eDEB8fHxWOAKoW7cuy5cvZ82aNbRs2ZLnn3+eN954w/Ip9wVRpQp89RWUKQPffw9PPml1RSIiIp7BZneMMJarSk5OJigoiKSkJMvHEwF88QX06wd2O8ycCf/4h9UViYiIFD/OfH+7XQuRwF13wQsvmPtPPAGrV1tbj4iIiLtTIHJTEybA/fdDRgbcfbcZVyQiIiIFo0Dkpmw2eO89aNcOTp0yM8/ysw2IiIiI5KRA5MZKloRly6BmTbM20ZgxVlckIiLinhSI3FzVqrBokbn/0UfgokW1RUREPJoCkQfo0AE6d4b0dLPnmYiIiDhHgchDjB9vbt9914wpEhERkfxTIPIQPXpAs2Zw5gy8847V1YiIiLgXBSIPYbPBuHHm/uuvw/nz1tYjIiLiThSIPMigQRAaCkePwocfWl2NiIiI+1Ag8iAlSlyaev/KK2bRRhEREbk2BSIP88gjUKGCWbn6iy+srkZERMQ9KBB5mLJl4bHHzP0XXzQbwIqIiEjeFIg80KhREBAAP/8M69dbXY2IiEjxp0DkgapUgYceMvdfesnaWkRERNyBApGHeuopMxX/m2/gt9+srkZERKR4UyDyUPXrw4AB5v7LL1tbi4iISHGnQOTBHAs1fvyxNn0VERHJiwKRB2vTRpu+ioiI5IcCkYdztBJp01cREZGrUyDycHfcoU1fRURErkWByMNp01cREZFrUyDyApdv+vrBB1ZXIyIiUvwoEHmBEiVg9Ghzf84cS0sREREplhSIvMT994OvL2zZYjZ+FRERkUsUiLxElSrQtau5v3ixtbWIiIgUNwpEXuTee83twoVgt1tbi4iISHGiQORF+vYFf3/YuRN+/dXqakRERIoPBSIvUr489Opl7i9aZGkpIiIixYoCkZe55x5zu2iRus1EREQcFIi8TJ8+UKYMHDgAP/9sdTUiIiLFgwKRlyldGu66y9xfuNDaWkRERIoLBSIv5Og2W7wYMjKsrUVERKQ4UCDyQj16QIUKkJAA69ZZXY2IiIj1FIi8kL8/DBhg7qvbTERERIHIazm6zT7/HC5csLYWERERqykQeanOnSEkBE6ehOhoq6sRERGxlgKRl/L1hYEDzX0t0igiIt5OgciLOfY2W7YMzp2ztBQRERFLKRB5sVtugdq14cwZ+OYbq6sRERGxjgKRF7PZsm/lISIi4q0UiLyco9vsm28gKcnaWkRERKyiQOTlmjeHhg0hLQ2++MLqakRERKyhQOTlbLZLrURapFFERLyVWwaimTNnUrduXUqWLElYWBjr16/P8/oFCxbQokULSpcuTbVq1XjooYc4ceJEEVVb/DnGEUVHw7Fj1tYiIiJiBbcLRIsXL2b06NFMmjSJ2NhYOnbsSM+ePYmLi8v1+h9//JEhQ4YwbNgwduzYwaeffsrmzZt55JFHirjy4qtBA2jVymz0+vnnVlcjIiJS9NwuEL366qsMGzaMRx55hEaNGjFjxgxCQ0OZNWtWrtdv2rSJOnXqMGrUKOrWrcutt97Ko48+ypYtW676GWlpaSQnJ2c7PJ26zURExJu5VSC6cOECW7duJSIiItv5iIgIYmJicn1N+/btOXToEMuXL8dut3P06FE+++wzevfufdXPiYqKIigoKOsIDQ0t1J+jOHKsWr1+PRw6ZG0tIiIiRc2tAtHx48fJyMggJCQk2/mQkBASEhJyfU379u1ZsGABgwYNwt/fn6pVq1K+fHnefPPNq35OZGQkSUlJWcfBgwcL9ecojmrVgltvBbsdPvnE6mpERESKllsFIgebzZbtsd1uz3HOYefOnYwaNYrJkyezdetWVqxYwYEDBxgxYsRV3z8gIIDAwMBshzdwDK7WOCIREfE2bhWIgoOD8fX1zdEalJiYmKPVyCEqKooOHTrw9NNP07x5c3r06MHMmTOZN28e8fHxRVG22/jb38ztpk1w/Li1tYiIiBQltwpE/v7+hIWFER0dne18dHQ07du3z/U1586dw8cn+4/p6+sLmJYluSQ0FFq0gMxMWLHC6mpERESKjlsFIoAxY8YwZ84c5s2bx65du3jyySeJi4vL6gKLjIxkyJAhWdf36dOHJUuWMGvWLPbv38+GDRsYNWoUbdu2pXr16lb9GMWWY6z5119bW4eIiEhR8rO6AGcNGjSIEydOMHXqVOLj42natCnLly+ndu3aAMTHx2dbk+jBBx8kJSWFt956i6eeeory5ctz++238+KLL1r1IxRrd94J//63aSG6eBFKlLC6IhEREdez2dVvdE3JyckEBQWRlJTk8QOsMzKgalUzhmjNGujUyeqKRERECsaZ72+36zIT1/L1hV69zH11m4mIiLdQIJIc7rzT3CoQiYiIt1AgkhwiIsDPD3bvht9/t7oaERER11MgkhyCgqBjR3P/m2+srUVERKQoKBBJrhzdZgpEIiLiDRSIJFeOQLRmDaSkWFqKiIiIyykQSa4aNID69c1aRFcsDC4iIuJxFIjkqjTbTEREvIUCkVyVIxAtX272NxMREfFUCkRyVR07QrlycPQobN1qdTUiIiKuo0AkV+Xvb9YkAnWbiYiIZ1MgkjxpHJGIiHgDBSLJU8+eYLPBtm1w5IjV1YiIiLiGApHkKSQE2rY195cvt7YWERERV1Egkmvq3dvcqttMREQ8lQKRXJNjHFF0NJw/b20tIiIirqBAJNfUsiVUrw7nzsHatVZXIyIiUvgUiOSabDbNNhMREc+mQCT5cnkgstutrUVERKSwKRBJvtx+OwQEwJ9/ws6dVlcjIiJSuBSIJF/KlDGhCNRtJiIinkeBSPLN0W32zTfW1iEiIlLYFIgk3xzrEW3YACdPWluLiIhIYVIgknyrXRuaNoXMTFixwupqRERECo8CkThF0+9FRMQTKRCJU3r1MrfR0aalSERExBMoEIlTbrkFypaF48dh+3arqxERESkcCkTilBIlLk2//+47a2sREREpLApE4rSICHOrQCQiIp5CgUic5ghEP/4IZ89aW4uIiEhhUCASp9WvD3XqwMWLsHat1dWIiIhcPwUicZrNpm4zERHxLApEUiAKRCIi4kkUiKRAbr8dfHxg1y44eNDqakRERK6PApEUSIUK0LatuR8dbW0tIiIi18uvIC/asWMHGzZs4PDhw6SmphIcHEzjxo257bbbCAwMLOwapZiKiIBNm0y32cMPW12NiIhIweU7EJ06dYrZs2fz7rvv8tdff2G323O+mZ8fvXr1YtSoUdzuWL1PPFZEBEydalqIMjLA19fqikRERAomX11mb7zxBvXr1+eVV16hZ8+eLFq0iH379pGUlERaWhrx8fFs2LCB//znP5w6dYru3btzxx138Pvvv7u6frFQ27YQGAgnT0JsrNXViIiIFJzNnltTzxXq16/P5MmTuffeeylRosQ13/SPP/7ghRdeoH79+kycOLFQCrVScnIyQUFBJCUlqUvwCv36wbJl8MIL4AF/1CIi4kGc+f7OVyBKT0/Hz8/54UYZGRn4ekA/igLR1c2aBY89Bp06wZo1VlcjIiJyiTPf3/nqMitIGAI8IgxJ3hzrEcXEwJkz1tYiIiJSUNc17f7cuXOcPHkyxyHe44YboF49beMhIiLuzelAdO7cOUaPHk3lypUpV64clStXznGId+ne3dxq1WoREXFXTveFjRw5ko8++og+ffrQqFEj/P39XVFXnmbOnMnLL79MfHw8TZo0YcaMGXTs2PGq16elpTF16lT++9//kpCQQM2aNZk0aRIPa/GcQhERAbNnKxCJiIj7cjoQffXVV0RFRTF27FhX1HNNixcvZvTo0cycOZMOHTowe/Zsevbsyc6dO6lVq1aurxk4cCBHjx5l7ty51K9fn8TERNLT04u4cs/l2MZj926Ii4Or/DGIiIgUW/maZXa5ypUrs2jRIrp27eqqmvLUrl07WrVqxaxZs7LONWrUiL59+xIVFZXj+hUrVnDPPfewf/9+KlasmK/PSEtLIy0tLetxcnIyoaGhmmWWh/btYeNGmDMHhg2zuhoREREXzDK7XP/+/fnOor6RCxcusHXrViIcU5v+T0REBDExMbm+5ssvv6R169a89NJL1KhRgwYNGjB27FhSU1Ov+jlRUVEEBQVlHaGhoYX6c3gixx+Jus1ERMQdOd1lNn36dAYMGMCYMWPo1atXrq0urVq1KpTirnT8+HEyMjIICQnJdj4kJISEhIRcX7N//35+/PFHSpYsydKlSzl+/DiPPfYYJ0+eZN68ebm+JjIykjFjxmQ9drQQydVFRMBzz8GqVdrGQ0RE3I/TgSg1NZX09HRmzJjB66+/nu05u92OzWYjIyOj0ArMjc1my/Vzc5OZmYnNZmPBggUEBQUB8Oqrr3L33Xfz9ttvU6pUqRyvCQgIICAgoPAL92CXb+OxbRu0aWN1RSIiIvnndCAaNmwYmzdvZvTo0UU+yyw4OBhfX98crUGJiYk5Wo0cqlWrRo0aNbLCEJgxR3a7nUOHDnHjjTe6tGZv4ecHXbvC0qWm20yBSERE3InTgeiHH37g1VdfZfjw4a6oJ0/+/v6EhYURHR1Nv379ss5HR0dz11135fqaDh068Omnn3LmzBnKli0LwN69e/Hx8aFmzZpFUre3iIi4FIgmTbK6GhERkfxzelB1uXLlqFOnjgtKyZ8xY8YwZ84c5s2bx65du3jyySeJi4tjxIgRgBn/M2TIkKzrBw8eTKVKlXjooYfYuXMn69at4+mnn+bhhx/OtbtMCu7ybTxSUqytRURExBlOB6IhQ4awaNEiV9SSL4MGDWLGjBlMnTqVli1bsm7dOpYvX07t2rUBiI+PJy4uLuv6smXLEh0dzenTp2ndujX33Xcfffr04Y033rDqR/BY9eqZrTzS07XRq4iIuBen1yFauHAhkyZNokWLFvTu3TvXWWb9+/cvtAKLA+12n3+PPQazZsHIkfDmm1ZXIyIi3syZ72+nA5GPT96NSkUxy6yoKRDl37Jl0K8fNGgAe/ZYXY2IiHgzZ76/CzSoWuRqunQxaxDt3Qt//gkWDjcTERHJN6cDUadOnVxRh3iIoCC45RbYsAGio8GCyYgiIiJOc3pQtci1OGabrVxpbR0iIiL5la9A1KtXL2JjY/P9pmlpabz66qu8/fbbBS5M3Ff37uZ29WqzjYeIiEhxl69AVLVqVdq0aUOHDh2YPXs2e3IZLZuSksKqVat44oknqFGjBm+//TY333xzoRcsxV+bNmYbj1OnwIkcLSIiYpl8BaJ58+axefNmatasyahRo2jcuDFly5albt26NGrUiJCQECpUqECPHj345ptvmDhxIjt37qR9+/aurl+KIT8/M7gazDgiERGR4s7pafeJiYmsXLmSTZs2ceTIEVJTUwkODqZhw4Z07tyZDh06XHWjVXelaffOe+steOIJuP12+P57q6sRERFv5NJ1iLyRApHz9uyBhg3B3990nZUubXVFIiLibZz5/tYsM3GJBg2gZk24cAF+/NHqakRERPKmQCQuYbNdmm22apW1tYiIiFyLApG4TLdu5laBSEREijsFInGZrl3NbWwsHDtmbS0iIiJ5USASlwkJgebNzf3Vq62tRUREJC9OB6J169Zx5syZXJ87c+YM69atu+6ixHM4us20HpGIiBRnTgeiLl26sHPnzlyf27NnD10cK/KJcGlgdXQ0aIEHEREprpwORHktW3Tx4kV8fNQLJ5d07AglSkBcHPzxh9XViIiI5M4vPxclJydz+vTprMcJCQnExcVluyY1NZUPPviAqlWrFmqB4t7KlIH27WHtWtNKVL++1RWJiIjklK9A9NprrzF16lQAbDYb/fr1y/U6u93OxIkTC6868Qjdu5tAtGoV/OMfVlcjIiKSU74CUUREBGXLlsVutzNu3DieeOIJatWqle2agIAAmjVrRqdOnVxSqLivbt3gX/8yM80yMsDX1+qKREREsstXIAoPDyc8PByAs2fPMnz4cKpXr+7SwsRzhIVBUBCcPg1bt0LbtlZXJCIikp3TI6CfffZZhSFxip+f2fUetGq1iIgUT/lqIbrSn3/+ySeffMJff/1FampqtudsNhtz584tlOLEc3TrBkuXmoHVGmYmIiLFjdOB6JtvvqF///5kZGRQpUoVAgICsj1vs9kKrTjxHI71iGJi4OxZM/tMRESkuLDZ81pYKBctW7akYsWKLFq0iCpVqriqrmIlOTmZoKAgkpKSCAwMtLoct2S3Q506Zj2iFSugRw+rKxIREU/nzPe302OI9u3bx/jx470mDEnhsNm0jYeIiBRfTgei2rVrX3UvM5G8OLrNNLBaRESKG6cD0cSJE3nllVc4d+6cK+oRD+aYafa//8HRo9bWIiIicjmnB1X//PPPJCYmUr9+fbp06UKlSpWyPW+z2Xj99dcLrUDxHFWqQIsWJhCtXg333mt1RSIiIobTg6qvtXmrzWYjIyPjuooqbjSouvA8/TS88go8/DBodQYREXEllw6qzszMzPPwtDAkhevygdXORXERERHXcToQiVyPjh3B3x8OHoR9+6yuRkRExChwIFq5ciWRkZEMHz6cuLg4ADZv3syxY8cKrTjxPKVLQ4cO5r6m34uISHHhdCA6d+4c3bt3p2fPnrz00kvMmzeP48ePA/DKK6/w4osvFnqR4lkc3Waafi8iIsWF04Fo0qRJbNmyhc8//5ykpCQuH5MdERHBKn3LyTU41iP64QdIT7e2FhEREShAIPr00095/vnn6devH6VKlcr2XK1atbK6z0SuplUrqFABkpJgyxarqxERESlAIDp27BhNmjTJ/c18fEhNTb3uosSz+fpeWqRRDYoiIlIcOB2IatSowa+//prrc7/88gt169a97qLE82lfMxERKU6cDkT9+/fnhRdeIDY2NuuczWbjr7/+4rXXXuPvf/97oRYonskxjmjjRkhJsbYWERERpwPRs88+S/Xq1Wnbti2tW7fGZrPx0EMP0bRpU6pUqcKECRNcUad4mBtugHr14OJFWLvW6mpERMTbOR2IypUrR0xMDM8//zxly5blhhtuoHTp0kRGRrJu3bocA61FriYiwtx+9521dYiIiDi9l5k30l5mrrF0KfTvDzfdBLt3W12NiIh4GpfuZSZSWLp0MTPO9uwBrdYgIiJWKlAgWrZsGX//+99p27YtzZs3z3a0aNGisGvMYebMmdStW5eSJUsSFhbG+vXr8/W6DRs24OfnR8uWLV1boORL+fLQtq25r9lmIiJiJacD0csvv0z//v1Zt24dJUqUoFKlStmOihUruqLOLIsXL2b06NFMmjSJ2NhYOnbsSM+ePa+5IGRSUhJDhgyha9euLq1PnKNxRCIiUhw4PYaobt26dO3aldmzZ+Pr6+uquq6qXbt2tGrVilmzZmWda9SoEX379iUqKuqqr7vnnnu48cYb8fX1ZdmyZWzfvj3fn6kxRK4TE2M2e61YERITTReaiIhIYXDpGKITJ04wePBgS8LQhQsX2Lp1KxGOZoX/ExERQUxMzFVf9/777/PHH3/w7LPP5utz0tLSSE5OznaIa7RtC4GBcPIkbNtmdTUiIuKtnA5EHTp0YNeuXa6o5ZqOHz9ORkYGISEh2c6HhISQkJCQ62v27dvHhAkTWLBgAX5+fvn6nKioKIKCgrKO0NDQ665dcufnB45eTHWbiYiIVZwORDNmzODtt9/myy+/5MKFC66o6ZpsNlu2x3a7Pcc5gIyMDAYPHsxzzz1HgwYN8v3+kZGRJCUlZR0HDx687prl6hyrVisQiYiIVfLXZHKZ+vXr061bN/r164fNZqN06dLZnrfZbCQlJRVagZcLDg7G19c3R2tQYmJijlYjgJSUFLZs2UJsbCwjR44EIDMzE7vdjp+fH9999x23O3YZvUxAQAABAQEu+RkkJ0cPqGMbj3LlrK1HRES8j9OBaNy4cbz11lu0bNmSRo0a4e/v74q6cuXv709YWBjR0dH069cv63x0dDR33XVXjusDAwNzbEQ7c+ZMVq9ezWeffaaNaIsJxzYe+/ebbTzuvNPqikRExNs4HYjmz5/P+PHj85zR5UpjxozhgQceoHXr1oSHh/Puu+8SFxfHiBEjANPddfjwYT788EN8fHxo2rRpttdXqVKFkiVL5jgv1oqIgHfeMd1mCkQiIlLUnA5EGRkZdHcM+rDAoEGDOHHiBFOnTiU+Pp6mTZuyfPlyateuDUB8fPw11ySS4ufyQCQiIlLUnF6HaODAgbRs2ZKJEye6qqZiR+sQud7p01CpEmRmwl9/Qa1aVlckIiLuzpnvb6dbiJ555hkGDRpEmTJl6N27d64rU7t6tWrxPOXLQ7t2ZmB1dDQMG2Z1RSIi4k2cbiHy8TEz9XOb5u6QkZFxfVUVM2ohKhpTpsBzz8HAgbB4sdXViIiIu3NpC9HkyZPzDEMiBRURYQLRqlWQkaFtPEREpOg43ULkjdRCVDTS0804ouRk+PlnaNPG6opERMSduXQvs8ulpqZy+PBh0tPTr+dtRACzjYdjnUzNNhMRkaJUoED0ww8/EB4eTrly5ahduza//PILAI8//jhLliwp1ALFuzhWrVYgEhGRouR0IFq9ejURERGcP3+esWPHkpmZmfVccHAw8+fPL8z6xMtcuY2HiIhIUXA6EE2ePJlevXoRGxvLtGnTsj3XokULtm/fXli1iRdybONx8aLZxkNERKQoOB2IYmNjefTRR4GcU+8rV65MYmJi4VQmXkvdZiIiUtScDkR+fn5cvHgx1+cSExMpp63K5To5doZRIBIRkaLidCBq06YNH330Ua7PffbZZ4SHh193UeLdbr8dfHxgzx7QtnQiIlIUnA5EEyZMYOnSpfTr148vv/wSm83GTz/9xMiRI/nss88YN26cK+oUL+LYxgPMNh4iIiKu5nQg6tatGx988AHr169nwIAB2O12Hn/8cT7++GPmz5/Prbfe6oo6xctoHJGIiBSlAq9UnZqaSkxMDEePHiU4OJgOHTpQpkyZwq6vWNBK1UUvJgY6dICKFSExUdt4iIiI81y6l5lDqVKl6Nq1a0FfLpKntm0hMBBOnoRt27SNh4iIuJbTXWbvv/8+U6ZMyfW5KVOm8OGHH15vTSLaxkNERIqU04HojTfeoEKFCrk+FxwczBtvvHHdRYnApXFEK1daW4eIiHg+pwPR77//TtOmTXN9rnHjxuzbt++6ixIB6NHD3MbEwOnTlpYiIiIerkCbuyYlJV31vHa+l8JSrx40bAgZGZp+LyIiruV0IGrWrBmLFi3K9bmFCxfSrFmz6y5KxKFXL3O7fLm1dYiIiGdzOhA5FmAcOnQoP/30E4cPH+ann37iwQcf5PPPP+eJJ55wRZ3ipRyB6NtvITPT2lpERMRzOT3tfvDgwezevZuoqCj++9//Zp338fHhX//6F/fdd1+hFijerWNHKFsWjh6F2FgIC7O6IhER8UQFWodo6tSpPPzww0RHR3Ps2DEqV65MREQEtWvXLuz6xMv5+5vNXpcuNd1mCkQiIuIKBV6p2ptopWprzZkDw4fDLbfAxo1WVyMiIu6iSFaqBjh27Bipqak5zteqVet63lYkm549ze1PP8GxY1C5srX1iIiI5ynQtPtp06ZRpUoVqlatSt26dXMcIoWpRg1o0QLsdi3SKCIiruF0IJo3bx7/+c9/GDVqFHa7nYkTJxIZGUnNmjW58cYbmTNnjivqFC+n6fciIuJKTgeit99+OysEAfTr149p06axe/duypUrx/Hjxwu9SBFHIFqxwizUKCIiUpgKtHXHLbfcgo+PeemFCxcAKFWqFE899RTvvvtu4VYoghlQXb48nDplxhKJiIgUJqcDkZ+fGYdts9kIDAzk0KFDWc8FBwdz+PDhwqtO5P/4+V3a20zdZiIiUticDkQ33ngjBw8eBKBNmza89957XLx4kYyMDN59913q1KlT2DWKABpHJCIiruN0IOrVqxfr1q0DIDIyktWrV1O+fHkqVqzI559/zvjx4wu9SBGAO+4wt7GxcOSItbWIiIhnue6FGTdv3syiRYuw2Wz07t2bLl26FFZtxYYWZiw+2raFzZth7lx4+GGrqxERkeKsyBZmBNNt1qZNm+t9G5F86dXLBKLlyxWIRESk8BRoYUYRqzjGEX33HVy8aG0tIiLiOfLVQnT77bfn+w1tNhvff/99gQsSyUvr1mbrjmPHYMMG6NzZ6opERMQT5KuFKDMzE7vdnq8jMzPT1TWLF/PxubS3mWabiYhIYdFu9/mgQdXFy+LFcM890KQJ/Pab1dWIiEhx5cz3t8YQiduJiDAtRTt2wF9/WV2NiIh4ggLNMsvIyOCTTz7hhx9+4MSJE1SqVIkuXbrw97//PWslaxFXqVAB2reHH3803Wb/+IfVFYmIiLtzuoXo+PHjtGvXjvvuu4/58+cTExPD/Pnzue+++2jXrp02d5UioVWrRUSkMDkdiJ588kn27NnDggULSE1NJT4+ntTUVP773/+yb98+nnzySVfUKZKNIxB9/z2cP29tLSIi4v6cDkRfffUV06ZN495778XX1xcAX19fBg8ezNSpU/nqq68KvUiRKzVvDtWrQ2oqrF1rdTUiIuLunA5EdrudJk2a5Ppc06ZNKYpJazNnzqRu3bqULFmSsLAw1q9ff9VrlyxZQvfu3alcuTKBgYGEh4ezcuVKl9cormWzqdtMREQKj9OBqFu3bqxatSrX56Kjo+ns4pXyFi9ezOjRo5k0aRKxsbF07NiRnj17EhcXl+v169ato3v37ixfvpytW7fSpUsX+vTpQ2xsrEvrFNdTIBIRkcLi9DpE27dvp3///gwYMIDBgwdTtWpVEhISWLBgAUuWLGHJkiXUqlUr6/qKFSsWasHt2rWjVatWzJo1K+tco0aN6Nu3L1FRUfl6jyZNmjBo0CAmT56cr+u1DlHxlJwMwcFmC4+9e+HGG62uSEREihOXbu7aqlUrAKZPn86rr76add6Rq8LCwrJdn5GR4exHXNWFCxfYunUrEyZMyHY+IiKCmJiYfL1HZmYmKSkpeQa1tLQ00tLSsh4nJycXrGBxqcBA6NgRVq82rUT//KfVFYmIiLtyOhBNnjwZm83milqu6fjx42RkZBASEpLtfEhICAkJCfl6j+nTp3P27FkGDhx41WuioqJ47rnnrqtWKRq9eplA9PXXCkQiIlJwTgeiKVOmuKAM51wZyOx2e75C2sKFC5kyZQpffPEFVapUuep1kZGRjBkzJutxcnIyoaGhBS9YXKZPHxg7FtasgVOnzKKNIiIiziqUrTvOnz/P7t27C7V7LDfBwcH4+vrmaA1KTEzM0Wp0pcWLFzNs2DA++eQTunXrlue1AQEBBAYGZjukeGrQABo3hvR0+OYbq6sRERF35XQgevPNN3n++eezHm/dupXQ0FCaNGlCgwYNOHjwYKEWeDl/f3/CwsKIjo7Odj46Opr27dtf9XULFy7kwQcf5OOPP6Z3794uq0+s0b+/uV2yxNo6RETEfTkdiObMmUP58uWzHo8fP56KFSvy2muvYbfbmTZtWmHWl8OYMWOYM2cO8+bNY9euXTz55JPExcUxYsQIwHR3DRkyJOv6hQsXMmTIEKZPn84tt9xCQkICCQkJJCUlubROKTqOQLRiBZw7Z20tIiLinpweQxQXF0fDhg0BSElJYd26dSxatIj+/ftToUKFfE9lL6hBgwZx4sQJpk6dSnx8PE2bNmX58uXUrl0bgPj4+GxrEs2ePZv09HQef/xxHn/88azzQ4cOZf78+S6tVYpGy5ZQu7bZ+X7lSujXz+qKRETE3TgdiNLS0ihRogQAGzduJDMzM2tMTp06dfI92+t6PPbYYzz22GO5PndlyFmzZo3L6xFr2Wymlei112DpUgUiERFxntNdZrVq1craKuOLL76gZcuWWYOOjx07pgHIYglHCPrqK7NQo4iIiDOcDkT3338/U6dOJSwsjNmzZ3P//fdnPbdlyxYaNGhQqAWK5Ef79lClCpw+babgi4iIOMPpLrNJkybh5+dHTEwM/fr1Y9SoUVnP/fbbbwwYMKBQCxTJD19fuOsueO89M9use3erKxIREXfi9F5m3kh7mbmHFSugZ0+oWhUOHwafQlllS0RE3JUz39/6yhCPcfvtZn+zhATYtMnqakRExJ3kq8vs4Ycf5plnnqFu3bo8/PDDeV5rs9mYO3duoRQn4gx/f7jzTvj4Y9NtlsdanSIiItnkKxD98MMP/PP/ds5cvXp1nvuGWbXxqwiY6fcff2ym37/8spmSLyIici35CkQHDhzIuv/nn3+6qhaR63bHHVCyJOzfD7/8Ai1aWF2RiIi4A40hEo9Spgz06GHuL11qbS0iIuI+FIjE4zgWadRmryIikl/56jLz8fFxamxQRkZGgQsSuV59+ph1iX79FX7/HerXt7oiEREp7vIViCZPnpwtEL3//vucOXOGPn36ULVqVeLj4/n6668pU6bMNWehibhaxYrQuTN8/73pNnv6aasrEhGR4i5fgWjKlClZ96dPn07VqlVZtWoVZcuWzTqfkpJCt27dKF26dKEXKeKs/v0ViEREJP+cHkM0c+ZMxo0bly0MAZQrV45x48Yxc+bMQitOpKD69jW3GzfCkSOWliIiIm7A6UB0+PBh/Pxyb1jy8/MjISHhuosSuV7Vq8Mtt5j7X3xhbS0iIlL8OR2IGjVqxKuvvsrFixeznb9w4QLTp0+nYcOGhVacyPXo39/caraZiIhci9O73U+bNo2+fftSr149+vfvT9WqVUlISGDJkiUkJCSwbNkyF5Qp4rx+/WDcOPjhBzh50gy2FhERyY3TLUS9e/dmxYoV1KhRg7fffptJkybx1ltvUbNmTb799lt69+7tijpFnFa/PjRrBhkZ8PXXVlcjIiLFmdMtRABdu3ala9eunDt3jlOnTlGhQgXNLpNiqV8/sx7RkiUwZIjV1YiISHF1XStVly5dmho1aigMSbHlGEe0ciWcPWttLSIiUnxp6w7xaM2bQ926cP48rFhhdTUiIlJcKRCJR7PZNNtMRESuTYFIPN6AAeb2iy/UbSYiIrlTIBKPd8stcMMNJgxpVQgREcmNApF4PJsN7r/f3P/oI2trERGR4kmBSLyCIxBFR0N8vLW1iIhI8aNAJF6hfn0ID4fMTFi40OpqRESkuFEgEq/xwAPm9sMPra1DRESKHwUi8RoDB0KJEvC//5nVq0VERBwUiMRrVKoEjq32NLhaREQup0AkXsXRbbZggdn0VUREBAq4uauIu+rdGypUgCNH4IcfoFs3qysSsc6xY7B/PyQmmvuJidnvO26Tk8Fuz/0Ac+vnB+XLm/++rnYbHAy1apmjZk3ThS1SXCgQiVcJCIBBg+Cdd0y3mQKReIPUVNi504yd+/VX+OUXc3v0aOF9xoULcO6c+cdGfthsUL061K5tjlq1zG2dOtCwobnvoz4MKUI2u92R8eVqkpOTCQoKIikpicDAQKvLkesUEwMdOkCZMuYLoUwZqysSKTx2O+zYAStXwk8/mfCzb59ZciI3NWtCSAhUqWKOypVz3g8MNOHEZsv9ALh4EU6fNsepU+Zw3HfcJiZCXJw5LlzI++coVcoEo0aNsh/164O/f+H9vsSzOfP9rRYi8Trh4WYrjz/+MFt53Hef1RWJXJ+TJ2HVKlixAr77Dg4fznlNpUrQrBk0b25umzWDJk2gbNmirzcz81I4+uuv7Ld//AF795pWrdhYc1zOz8+EopYtzXHzzea2SpWi/znEs6iFKB/UQuR5pkyB556DiAjzL2kRd5KRAT//bP7urlgBmzdnbwEqWRI6dYLbbzdhoVkzqFr1UmtOcZeeDgcOwK5d5ti509zu3g0pKbm/plq17AHp5pvNP3zc5WcW13Dm+1uBKB8UiDzP77/DjTeaboBDh8z/TEWKu337YN48+OCDnFvQNGkCPXqYo2NH0+Xkaex20/r122+wffulY+/eSwO8L1e+PISFQevW0KaNua1VSyHJmygQFTIFIs/Uvj1s3AivvAJPPWV1NSK5O3sWPvsM5s6F9esvnS9fHrp3vxSCata0rETLnTljBolfHpL+9z9IS8t5beXKJhg5QlKbNqb1TDyTAlEhUyDyTLNmwWOPQYsW5n+gIsWF3W66webONXvvObqJfHzgjjtg2DC4804NLs7LxYtmcPnmzbBlizl++cV0x10pNBTatr10hIVBuXJFX7MUPgWiQqZA5JlOnDBdZRcvmv9RNmtmdUXi7c6fN11is2aZbiGHevXg4Ydh6FDvbgm6XufPm//Wt2wxQennn83YpCu/BW02M6OtbdtLrUjNm5tlO8S9KBAVMgUiz9Wvn5lp9vTT8NJLVlcj3io1Fd57D1588dI6PiVLwt13m9ag227TmjyukpIC27aZcOQ44uJyXleihGlNvryrrVEjM+tNii8FokKmQOS5liyBAQPMAnFxceDra3VF4k3OnYN33zVBKCHBnAsNNQH9gQfMOCEpekePXmpB+vln06J04kTO60qXNrPZ2rSBVq3McdNNCknFiQJRIVMg8lxpaabb7NQpiI7WytVSNM6eNaulv/zypdWia9WCiRPhwQfVNVPc2O3w558mJDnGJG3dmvsSAKVKme41R0Bq1crMANSfqTUUiAqZApFn+8c/zJfTkCFmOrOIq5w9CzNnmpmNiYnmXJ06JggNHapB0u4kMxP27LkUkGJjzeSMM2dyXluiBDRubILS5QtjVqumJQBczeMD0cyZM3n55ZeJj4+nSZMmzJgxg44dO171+rVr1zJmzBh27NhB9erVGTduHCNGjMj35ykQeTZt5SGuZrfDokVmeQfH+kH16sGkSaZrTJuceobMTLPG2bZt2Y9Tp3K/3rF6uGMF8aZNzbikoKCirduTeXQgWrx4MQ888AAzZ86kQ4cOzJ49mzlz5rBz505q1aqV4/oDBw7QtGlThg8fzqOPPsqGDRt47LHHWLhwIQMGDMjXZyoQeTa73SzS+McfZsPX+++3uiLxJDt3wsiR8MMP5nG9evDMM2bLGAUhz2e3m/GJsbHZN9fNa3+5qlXNPm6O46abzG2tWhpc7yyPDkTt2rWjVatWzJo1K+tco0aN6Nu3L1FRUTmuHz9+PF9++SW7du3KOjdixAj+97//sXHjxlw/Iy0tjbTLVvRKTk4mNDRUgciDObby6NbNjCUSuV5nzsDUqfDaa2btm1KlTIvQ2LEaTyJmZuGuXZcC0q+/mqUWrlyB/HKlSpl/vNWrZ466dc1Rr57pevXE1cmvl8du7nrhwgW2bt3KhAkTsp2PiIggJiYm19ds3LiRiIiIbOd69OjB3LlzuXjxIiVy+SdaVFQUzz33XOEVLsXe0KHmy2vVKrNfUsOGVlck7spuh08/hTFjLm2yetddMGOG+dISARNeHIOuL5ecbMYm7d5tDsf9fftMiPrlF3Pkplq1SyGpZk2oUcPMoK1RwxzVqmkGXF7c6ldz/PhxMjIyCAkJyXY+JCSEBMec1SskJCTken16ejrHjx+nWi6bWEVGRjJmzJisx44WIvFcdetCnz7w5Zfw5pvw9ttWVyTuaM8e0z22apV5XK8evPEG9O5tbV3iPgIDL61zdLn0dDPTbe9es/Ht/v2XbvfvNzPe4uPNcZX2AWw2CAm5FJSqVIHg4EtH5crZHwcGetegb7cKRA62K/6E7HZ7jnPXuj638w4BAQEEqE3b6/zznyYQffABvPCC1oCR/LtwAZ5/3qwndPGi6RKLjITx480CiyLXy88P6tc3x5Xsdjh58lJA+vNP0zp5+REfb0JVQoI5tm7N32cGBuZ9lCtn/o7ndQQEmDXefH3NGKir3S9f3toB5W4ViIKDg/H19c3RGpSYmJijFcihatWquV7v5+dHpUqVXFaruJ8uXcwsj99+M3tIacNXyY9du8wA6dhY87h3b9MqVK+etXWJ97DZzIy1SpXMStq5ycyEY8cuBaQjR+D48ezHsWOX7p89awLUyZPmKAoTJkAuQ4GLjFsFIn9/f8LCwoiOjqZfv35Z56Ojo7nrrrtyfU14eDhfffVVtnPfffcdrVu3znX8kHgvmw1GjYL/9//grbdg9GitXC1XZ7ebrtWnnzZ7ZFWqZNazGjDAu7oZxD34+JjuspCQnOOWcpOaalbnTkkx45qudqSkmL//1zoyMkwoy8i4+n2r1+Fyu1lmjmn377zzDuHh4bz77ru899577Nixg9q1axMZGcnhw4f58MMPgUvT7h999FGGDx/Oxo0bGTFihKbdS67OnTNbJ5w8abb1uCx3i2Q5csRstrpypXncowe8/74ZtCoixYcz399ut6LBoEGDmDFjBlOnTqVly5asW7eO5cuXU7t2bQDi4+OJu2xnvrp167J8+XLWrFlDy5Ytef7553njjTfyHYbEu5QubVqIwHR7iFxpyRKzkN7KlWZ8xJtvwrffKgyJuDu3ayGyglqIvMvBg2bWWUYG/O9/ZgVZkeRkM/B+/nzz+OabYcECs7KwiBRPHt1CJOJqoaHQv7+5//rr1tYixUNMDLRsacKQzWYGf27apDAk4kkUiERy8c9/mtsFC8zMC/FOdrtZUPG228yU5tq1Yc0aMxPG6gGgIlK4FIhEctG+PYSFQVoavPee1dWIFc6ehcGD4cknTffpPfeYLtTbbrO6MhFxBQUikVzYbJdaiWbONIvtiffYtw9uucXsUO/nZ1qJPv5Yu5CLeDIFIpGrGDjQrNlx+DB8/rnV1UhR+eors7jdb7+ZXcdXrzbhWGsLiXg2BSKRqwgIgBEjzH0NrvZ8GRnwzDPwt7+ZGWUdOsC2bdCxo9WViUhRUCASycOIEVCihJlR9PPPVlcjrnLypNlyY9o08/iJJ0zLkNYWEvEeCkQieaha1QymBS3U6KliY80A+pUroVQp+Ogj82etWWQi3kWBSOQaRo0yt598YnaMFs+xZInpGvvzT7MZ68aNcP/9VlclIlZQIBK5htatzTT8ixdh1iyrq5HCYLfDf/5jNmJNTYU77oAtW6BFC6srExGrKBCJ5INjCv4775idm8V9XbgAw4ZBZKR5PHKkmVlWoYK1dYmItRSIRPKhXz+oWdOsWr1okdXVSEGdOAEREWZneh8fszHrm2+atYZExLspEInkQ4kS8Nhj5v706WaKtriXvXshPBzWroVy5eDrr03rkIgIKBCJ5Nujj5qVin/7zaxaLO5jzRqz8vS+fWY/spgY6NnT6qpEpDhRIBLJp4oVzS7nAP/6l8YSuYv33zfdZKdOQbt28NNP0LSp1VWJSHGjQCTihFGjoEYNiIsze5xJ8ZWZaQZOP/ywmSE4aBD88IPZjkVE5EoKRCJOKF0annvO3H/hBTh92tJy5CrS0sx6Qv/5j3n8zDOmm7NUKWvrEpHiS4FIxElDh0KjRma7hxdftLoaudLp02ZdoYULzeyxDz6AqVPNrDIRkavR/yJEnOTnd6nlYcYMOHzY0nLkMgcPwq23mkHU5crB8uUwZIjVVYmIO1AgEimAPn3Mlg/nz8OUKVZXIwC//GJmku3YAdWrw/r10L271VWJiLtQIBIpAJsNXnrJ3J83D3butLYeb/f996Zl6MgRaNzY7EmmbThExBkKRCIF1L499O1rZjNNnGh1Nd7ro4/MmKGUFOjUCX78EWrVsroqEXE3CkQi1+Hf/zaDdb/4AjZssLoa72K3Q1SUGSOUng733AMrV2pPMhEpGAUikevQqJHZKBRg3DjzJS2ul54O//jHpZa5p5+GBQsgIMDaukTEfSkQiVynZ58169vExMCXX1pdjec7e9Zstjt7thnL9cYbZjyXptWLyPXQ/0JErlONGjB6tLkfGWlaL8Q1jh6FLl3MxqwlS8Lnn8MTT1hdlYh4AgUikUIwfrzZ62zXLpg/3+pqPNPevWYg++bNUKkSrF5tWopERAqDApFIIQgKMhu+gulCO3fO2no8TUyMCUP790O9euZxeLjVVYmIJ1EgEikkjz0GtWubtXBef93qajzH0qXQtSucOAFt2pg1hho0sLoqEfE0CkQihSQgAKZNM/eff950n8n1efNNGDDArAh+551mt/oqVayuSkQ8kQKRSCEaPNhsF5GaCvfea3ZdF+dlZpqp9KNGmaUMRowwLUVlylhdmYh4KgUikULk42N2Vw8Ohv/9z8w6E+ekpppg+cor5nFUFMycaTbVFRFxFQUikUJWrRq8/765/9prsGKFtfW4k/h46NwZFi+GEiXMthwTJpj1hkREXEmBSMQF7rwTRo4094cONevnSN62b4e2beHnn80SBt99B/ffb3VVIuItFIhEXOSll6BpU0hMhAcfNONiJHfLlkGHDnDoEDRsCD/9ZFqKRESKigKRiIuUKgULF5oVlVesMDOmJDu7HV58Efr3N2s3de9uptXXr291ZSLibRSIRFyoaVOYPt3cHzfOdAuJkZYGDz1kxgjZ7fD447B8OZQvb3VlIuKNFIhEXOwf/4A+feDCBTMVX6tYw7FjZrHFDz4AX1946y1zaCaZiFhFgUjExWw2mDfPzD7bvRvGjLG6Imv99psZPL1hg9nyZPly0zokImIlBSKRIhAcDB9+aO7Pnm0WGfRGH39s9iD780+44QbYtAkiIqyuSkREgUikyHTrZlZfBnjkETOjylukpsLw4XDffXDmDHTpYmaSNWxodWUiIoYCkUgRmjYNwsLg5Em4+25ISrK6Itfbvdt0kc2ZY7oPJ0+G6GioVMnqykRELlEgEilC/v5mKn758qaF5Pbb4fhxq6tynQ8/NAHwt98gJMQEoeeeMwOpRUSKE7cKRKdOneKBBx4gKCiIoKAgHnjgAU6fPn3V6y9evMj48eNp1qwZZcqUoXr16gwZMoQjR44UXdEiV7jxRrNre+XKsG0bdOoEnvZX8uxZM6V+6FAzq65rV7PkQNeuVlcmIpI7twpEgwcPZvv27axYsYIVK1awfft2Hnjggatef+7cObZt28YzzzzDtm3bWLJkCXv37uVvf/tbEVYtklPLlrBuHdSoATt3QseOZqCxJ9ixw3SRzZ9vNrudOhVWroSqVa2uTETk6mx2u91udRH5sWvXLho3bsymTZto164dAJs2bSI8PJzdu3dz00035et9Nm/eTNu2bfnrr7+oVatWvl6TnJxMUFAQSUlJBAYGFvhnELnSgQNmsPX+/SYcrVrlvgON7XaYOxdGjTKDqKtVM92DnTpZXZmIeCtnvr/dpoVo48aNBAUFZYUhgFtuuYWgoCBiYmLy/T5JSUnYbDbK57EcblpaGsnJydkOEVeoWxfWr4fGjeHwYbjtNvdczfrXX03wGT7chKGICPNzKAyJiLtwm0CUkJBAlSpVcpyvUqUKCQkJ+XqP8+fPM2HCBAYPHpxnUoyKisoapxQUFERoaGiB6xa5lurVYe1aaNXKrODcpYvZz8sdJCebhSZvvtkEu9Kl4eWX4dtvIZf/XEVEii3LA9GUKVOw2Wx5Hlu2bAHAZrPleL3dbs/1/JUuXrzIPffcQ2ZmJjNnzszz2sjISJKSkrKOgwcPFuyHE8mn4GBYvdrs+H76tNnk9Pvvra7q6ux2s8jiTTfBa69BRgYMGAC7dsHYsWbskIiIO7F856CRI0dyzz335HlNnTp1+OWXXzh69GiO544dO0ZISEier7948SIDBw7kwIEDrF69+pr9iAEBAQQEBFy7eJFCFBRkBh/362emp/fubfb6GjjQrN9TXOzYYbbaWLvWPL7xRnjzTejRw9q6RESuh+WBKDg4mODg4GteFx4eTlJSEj///DNt27YF4KeffiIpKYn27dtf9XWOMLRv3z5++OEHKmk1OCnGypSBr76Ce+6BZcvM7VtvwYsvQh5/zYtESopZQ+j11yE9HUqVgn/9C556CvTvBxFxd27TsN2oUSPuuOMOhg8fzqZNm9i0aRPDhw/nzjvvzDbDrGHDhiz9v42i0tPTufvuu9myZQsLFiwgIyODhIQEEhISuHDhglU/ikieAgLgk08gMhJKloQffzRdaX37min6Re3332H8eKhXD6ZPN2GoXz/TPTZxosKQiHgGtwlEAAsWLKBZs2ZEREQQERFB8+bN+eijj7Jds2fPHpL+bz+EQ4cO8eWXX3Lo0CFatmxJtWrVsg5nZqaJFLUSJeDf/4Z9+8y+Zz4+8MUX0KwZDBsGrh7WduGCCWXdupkusZdeMitq33CD2Z1+yRKoXdu1NYiIFCW3WYfISlqHSKy2axdMmgT/1/hJQIBZ72fCBKhYsfA+5/ff4b334P33zYw3MOOX7rgDHn3UjGvys7yjXUQkf5z5/lYgygcFIikuNm403Vfr15vH5cvD4MHQvDk0bQpNmphz+WG3w4kTphVq1y6ziOKqVZeer1bNtEY98ohag0TEPSkQFTIFIilO7HbTbRUZaRZEvFLNmiYYNW16KSSlp5vWn337st9euRWgzWZmizlag0qUKJIfSUTEJRSICpkCkRRHGRlmXNHGjWY3+d9+g0OHnH+fmjWhfn249VbTIlSnTqGXKiJiCWe+vzUaQMRN+fpC//7mcDh92sxEcwSk334z6wb5+5vB0fXrZ7+tV8+sLi0i4u0UiEQ8SPnyZr0iq9csEhFxN2417V5ERETEFRSIRERExOspEImIiIjXUyASERERr6dAJCIiIl5PgUhERES8ngKRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEImIiIjXUyASERERr6dAJCIiIl5PgUhERES8np/VBbgDu90OQHJyssWViIiISH45vrcd3+N5USDKh5SUFABCQ0MtrkRERESclZKSQlBQUJ7X2Oz5iU1eLjMzkyNHjlCuXDlsNluhvndycjKhoaEcPHiQwMDAQn1vuUS/56Kh33PR0O+5aOj3XDRc+Xu22+2kpKRQvXp1fHzyHiWkFqJ88PHxoWbNmi79jMDAQP0HVwT0ey4a+j0XDf2ei4Z+z0XDVb/na7UMOWhQtYiIiHg9BSIRERHxegpEFgsICODZZ58lICDA6lI8mn7PRUO/56Kh33PR0O+5aBSX37MGVYuIiIjXUwuRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEFlo5syZ1K1bl5IlSxIWFsb69eutLsnjrFu3jj59+lC9enVsNhvLli2zuiSPExUVRZs2bShXrhxVqlShb9++7Nmzx+qyPNKsWbNo3rx51gJ24eHhfPvtt1aX5fGioqKw2WyMHj3a6lI8ypQpU7DZbNmOqlWrWlaPApFFFi9ezOjRo5k0aRKxsbF07NiRnj17EhcXZ3VpHuXs2bO0aNGCt956y+pSPNbatWt5/PHH2bRpE9HR0aSnpxMREcHZs2etLs3j1KxZk//85z9s2bKFLVu2cPvtt3PXXXexY8cOq0vzWJs3b+bdd9+lefPmVpfikZo0aUJ8fHzW8euvv1pWi6bdW6Rdu3a0atWKWbNmZZ1r1KgRffv2JSoqysLKPJfNZmPp0qX07dvX6lI82rFjx6hSpQpr167ltttus7ocj1exYkVefvllhg0bZnUpHufMmTO0atWKmTNnMm3aNFq2bMmMGTOsLstjTJkyhWXLlrF9+3arSwHUQmSJCxcusHXrViIiIrKdj4iIICYmxqKqRApHUlISYL6oxXUyMjJYtGgRZ8+eJTw83OpyPNLjjz9O79696datm9WleKx9+/ZRvXp16tatyz333MP+/fstq0Wbu1rg+PHjZGRkEBISku18SEgICQkJFlUlcv3sdjtjxozh1ltvpWnTplaX45F+/fVXwsPDOX/+PGXLlmXp0qU0btzY6rI8zqJFi9i2bRubN2+2uhSP1a5dOz788EMaNGjA0aNHmTZtGu3bt2fHjh1UqlSpyOtRILKQzWbL9thut+c4J+JORo4cyS+//MKPP/5odSke66abbmL79u2cPn2azz//nKFDh7J27VqFokJ08OBB/vnPf/Ldd99RsmRJq8vxWD179sy636xZM8LDw7nhhhv44IMPGDNmTJHXo0BkgeDgYHx9fXO0BiUmJuZoNRJxF0888QRffvkl69ato2bNmlaX47H8/f2pX78+AK1bt2bz5s28/vrrzJ492+LKPMfWrVtJTEwkLCws61xGRgbr1q3jrbfeIi0tDV9fXwsr9ExlypShWbNm7Nu3z5LP1xgiC/j7+xMWFkZ0dHS289HR0bRv396iqkQKxm63M3LkSJYsWcLq1aupW7eu1SV5FbvdTlpamtVleJSuXbvy66+/sn379qyjdevW3HfffWzfvl1hyEXS0tLYtWsX1apVs+Tz1UJkkTFjxvDAAw/QunVrwsPDeffdd4mLi2PEiBFWl+ZRzpw5w++//571+MCBA2zfvp2KFStSq1YtCyvzHI8//jgff/wxX3zxBeXKlctq+QwKCqJUqVIWV+dZJk6cSM+ePQkNDSUlJYVFixaxZs0aVqxYYXVpHqVcuXI5xsCVKVOGSpUqaWxcIRo7dix9+vShVq1aJCYmMm3aNJKTkxk6dKgl9SgQWWTQoEGcOHGCqVOnEh8fT9OmTVm+fDm1a9e2ujSPsmXLFrp06ZL12NEvPXToUObPn29RVZ7FsXRE586ds51///33efDBB4u+IA929OhRHnjgAeLj4wkKCqJ58+asWLGC7t27W12aiNMOHTrEvffey/Hjx6lcuTK33HILmzZtsux7UOsQiYiIiNfTGCIRERHxegpEIiIi4vUUiERERMTrKRCJiIiI11MgEhEREa+nQCQiIiJeT4FIREREvJ4CkYiIiHg9BSIRKRZiYmKYMmUKp0+fzvFc586dc6yEXRx8+OGHVK5cmZSUlHy/Zu/evfj7+7Nt2zYXViYiztJK1SJSLLzyyis8/fTTHDhwgDp16mR7bufOnQA0btzYgspyd+7cORo0aMDo0aMZO3asU6996KGH2L9/P2vXrnVRdSLiLLUQiUix17hx42IVhgA++OADTpw4wSOPPOL0a0eOHMm6deuIiYlxQWUiUhAKRCJiuSlTpvD0008DULduXWw2GzabjTVr1gA5u8z+/PNPbDYbL7/8Mi+++CJ16tShVKlSdO7cmb1793Lx4kUmTJhA9erVCQoKol+/fiQmJub43MWLFxMeHk6ZMmUoW7YsPXr0IDY2Nl81z5o1iz59+lC+fPls5z/99FPatWtHUFAQpUuXpl69ejz88MPZrgkLC6NRo0a88847+f8liYhLKRCJiOUeeeQRnnjiCQCWLFnCxo0b2bhxI61atcrzdW+//TYbNmzg7bffZs6cOezevZs+ffowbNgwjh07xrx583jppZdYtWpVjpacf//739x77700btyYTz75hI8++oiUlBQ6duyY1UV3NYcOHeLXX3+lS5cu2c5v3LiRQYMGUa9ePRYtWsQ333zD5MmTSU9Pz/EenTt35ttvv0WjFkSKBz+rCxARqVmzJrVq1QLg5ptvzjGG6GrKly/PsmXL8PEx/7Y7fvw4o0ePpmHDhnzxxRdZ1+3evZsZM2aQnJxMYGAgBw8e5Nlnn2XkyJG88cYbWdd1796dG2+8keeee47Fixdf9XMdXV1XBraYmBjsdjvvvPMOQUFBWecffPDBHO/RqlUrZs2axZ49e2jYsGG+fl4RcR21EImI2+rVq1dWGAJo1KgRAL179852neN8XFwcACtXriQ9PZ0hQ4aQnp6edZQsWZJOnTplddVdzZEjRwCoUqVKtvNt2rQBYODAgXzyySccPnz4qu/heG1e14hI0VEgEhG3VbFixWyP/f398zx//vx5AI4ePQqYAFOiRIlsx+LFizl+/Hien5uamgpAyZIls52/7bbbWLZsWVbYqlmzJk2bNmXhwoU53sPxWsd7iYi11GUmIl4nODgYgM8++4zatWsX+PUnT56kWrVq2Z676667uOuuu0hLS2PTpk1ERUUxePBg6tSpQ3h4eNZ1J0+ezPZeImItBSIRKRYCAgKAomkx6dGjB35+fvzxxx8MGDDA6dc7xvz88ccfNGnSJNdrAgIC6NSpE+XLl2flypXExsZmC0T79+/Hx8eHm266qWA/hIgUKgUiESkWmjVrBsDrr7/O0KFDKVGiBDfddBPlypUr9M+qU6cOU6dOZdKkSezfv5877riDChUqcPToUX7++WfKlCnDc889d9XXt2vXjlKlSrFp0yb+9re/ZZ2fPHkyhw4domvXrtSsWZPTp0/z+uuvU6JECTp16pTtPTZt2kTLli2pUKFCof98IuI8jSESkWKhc+fOREZG8tVXX3HrrbfSpk0btm7d6rLPi4yM5LPPPmPv3r0MHTqUHj16MG7cOP766y9uu+22PF/r7+/P3XffnW0mG5iglJCQwPjx44mIiOD//b//R6lSpVi9enW2lqQzZ87w/fffc99997nkZxMR52nrDhGRAtiyZQtt2rRh06ZNtGvXzqnXzp07l3/+858cPHhQLUQixYQCkYhIAQ0aNIizZ8/y9ddf5/s16enpNG7cmKFDhzJp0iQXVicizlCXmYhIAU2fPp02bdo4tdv9wYMHuf/++3nqqadcWJmIOEstRCIiIuL11EIkIiIiXk+BSERERLyeApGIiIh4PQUiERER8XoKRCIiIuL1FIhERETE6ykQiYiIiNdTIBIRERGv9/8BWAXjX6STu+AAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/papermill.record/image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGxCAYAAACDV6ltAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABLj0lEQVR4nO3deViU9f7/8ecAgiu4oLjhlpm7Ji6hmZqKqXly6WhZaWX+8pR5zExFT2Zmh9Ni2aZZalbH1Ba1zTTMXBItFzyVe2nhAuIKqIgC8/vj8x0UQWSQ4WZmXo/ruq+ZueeemTdozcvParPb7XZEREREvJiP1QWIiIiIWE2BSERERLyeApGIiIh4PQUiERER8XoKRCIiIuL1FIhERETE6ykQiYiIiNfzs7oAd5CZmcmRI0coV64cNpvN6nJEREQkH+x2OykpKVSvXh0fn7zbgBSI8uHIkSOEhoZaXYaIiIgUwMGDB6lZs2ae1ygQ5UO5cuUA8wsNDAy0uBoRERHJj+TkZEJDQ7O+x/OiQJQPjm6ywMBABSIRERE3k5/hLhpULSIiIl5PgUhERES8ngKRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEImIiIjXUyASERERr6dAJCIiIl5PgUhERES8ntsFonXr1tGnTx+qV6+OzWZj2bJl13zN2rVrCQsLo2TJktSrV4933nnH9YWKiIiI23C7QHT27FlatGjBW2+9la/rDxw4QK9evejYsSOxsbFMnDiRUaNG8fnnn7u4UhEREXEXbre5a8+ePenZs2e+r3/nnXeoVasWM2bMAKBRo0Zs2bKFV155hQEDBrioyvxJSoLTp699XcmSEBLi8nJERES8ltsFImdt3LiRiIiIbOd69OjB3LlzuXjxIiVKlMjxmrS0NNLS0rIeJycnu6S2WbMgMjJ/144eDa++CvnYsFdERESc5HZdZs5KSEgg5IrmlZCQENLT0zl+/Hiur4mKiiIoKCjrCA0NdUltfn6m9edaB8CMGZDPXkIRERFxkscHIgDbFc0qdrs91/MOkZGRJCUlZR0HDx50SV1jx0Jq6rWPl14y148eDd9955JSREREvJrHB6KqVauSkJCQ7VxiYiJ+fn5UqlQp19cEBAQQGBiY7bDS2LEwdChkZsLAgbBnj6XliIiIeByPD0Th4eFER0dnO/fdd9/RunXrXMcPFUc2G8yeDe3bm4HYffrAyZNWVyUiIuI53C4QnTlzhu3bt7N9+3bATKvfvn07cXFxgOnuGjJkSNb1I0aM4K+//mLMmDHs2rWLefPmMXfuXMaOHWtF+QUWEABLl0KtWrBvn2kpunjR6qpEREQ8g9sFoi1btnDzzTdz8803AzBmzBhuvvlmJk+eDEB8fHxWOAKoW7cuy5cvZ82aNbRs2ZLnn3+eN954w/Ip9wVRpQp89RWUKQPffw9PPml1RSIiIp7BZneMMJarSk5OJigoiKSkJMvHEwF88QX06wd2O8ycCf/4h9UViYiIFD/OfH+7XQuRwF13wQsvmPtPPAGrV1tbj4iIiLtTIHJTEybA/fdDRgbcfbcZVyQiIiIFo0Dkpmw2eO89aNcOTp0yM8/ysw2IiIiI5KRA5MZKloRly6BmTbM20ZgxVlckIiLinhSI3FzVqrBokbn/0UfgokW1RUREPJoCkQfo0AE6d4b0dLPnmYiIiDhHgchDjB9vbt9914wpEhERkfxTIPIQPXpAs2Zw5gy8847V1YiIiLgXBSIPYbPBuHHm/uuvw/nz1tYjIiLiThSIPMigQRAaCkePwocfWl2NiIiI+1Ag8iAlSlyaev/KK2bRRhEREbk2BSIP88gjUKGCWbn6iy+srkZERMQ9KBB5mLJl4bHHzP0XXzQbwIqIiEjeFIg80KhREBAAP/8M69dbXY2IiEjxp0DkgapUgYceMvdfesnaWkRERNyBApGHeuopMxX/m2/gt9+srkZERKR4UyDyUPXrw4AB5v7LL1tbi4iISHGnQOTBHAs1fvyxNn0VERHJiwKRB2vTRpu+ioiI5IcCkYdztBJp01cREZGrUyDycHfcoU1fRURErkWByMNp01cREZFrUyDyApdv+vrBB1ZXIyIiUvwoEHmBEiVg9Ghzf84cS0sREREplhSIvMT994OvL2zZYjZ+FRERkUsUiLxElSrQtau5v3ixtbWIiIgUNwpEXuTee83twoVgt1tbi4iISHGiQORF+vYFf3/YuRN+/dXqakRERIoPBSIvUr489Opl7i9aZGkpIiIixYoCkZe55x5zu2iRus1EREQcFIi8TJ8+UKYMHDgAP/9sdTUiIiLFgwKRlyldGu66y9xfuNDaWkRERIoLBSIv5Og2W7wYMjKsrUVERKQ4UCDyQj16QIUKkJAA69ZZXY2IiIj1FIi8kL8/DBhg7qvbTERERIHIazm6zT7/HC5csLYWERERqykQeanOnSEkBE6ehOhoq6sRERGxlgKRl/L1hYEDzX0t0igiIt5OgciLOfY2W7YMzp2ztBQRERFLKRB5sVtugdq14cwZ+OYbq6sRERGxjgKRF7PZsm/lISIi4q0UiLyco9vsm28gKcnaWkRERKyiQOTlmjeHhg0hLQ2++MLqakRERKyhQOTlbLZLrURapFFERLyVWwaimTNnUrduXUqWLElYWBjr16/P8/oFCxbQokULSpcuTbVq1XjooYc4ceJEEVVb/DnGEUVHw7Fj1tYiIiJiBbcLRIsXL2b06NFMmjSJ2NhYOnbsSM+ePYmLi8v1+h9//JEhQ4YwbNgwduzYwaeffsrmzZt55JFHirjy4qtBA2jVymz0+vnnVlcjIiJS9NwuEL366qsMGzaMRx55hEaNGjFjxgxCQ0OZNWtWrtdv2rSJOnXqMGrUKOrWrcutt97Ko48+ypYtW676GWlpaSQnJ2c7PJ26zURExJu5VSC6cOECW7duJSIiItv5iIgIYmJicn1N+/btOXToEMuXL8dut3P06FE+++wzevfufdXPiYqKIigoKOsIDQ0t1J+jOHKsWr1+PRw6ZG0tIiIiRc2tAtHx48fJyMggJCQk2/mQkBASEhJyfU379u1ZsGABgwYNwt/fn6pVq1K+fHnefPPNq35OZGQkSUlJWcfBgwcL9ecojmrVgltvBbsdPvnE6mpERESKllsFIgebzZbtsd1uz3HOYefOnYwaNYrJkyezdetWVqxYwYEDBxgxYsRV3z8gIIDAwMBshzdwDK7WOCIREfE2bhWIgoOD8fX1zdEalJiYmKPVyCEqKooOHTrw9NNP07x5c3r06MHMmTOZN28e8fHxRVG22/jb38ztpk1w/Li1tYiIiBQltwpE/v7+hIWFER0dne18dHQ07du3z/U1586dw8cn+4/p6+sLmJYluSQ0FFq0gMxMWLHC6mpERESKjlsFIoAxY8YwZ84c5s2bx65du3jyySeJi4vL6gKLjIxkyJAhWdf36dOHJUuWMGvWLPbv38+GDRsYNWoUbdu2pXr16lb9GMWWY6z5119bW4eIiEhR8rO6AGcNGjSIEydOMHXqVOLj42natCnLly+ndu3aAMTHx2dbk+jBBx8kJSWFt956i6eeeory5ctz++238+KLL1r1IxRrd94J//63aSG6eBFKlLC6IhEREdez2dVvdE3JyckEBQWRlJTk8QOsMzKgalUzhmjNGujUyeqKRERECsaZ72+36zIT1/L1hV69zH11m4mIiLdQIJIc7rzT3CoQiYiIt1AgkhwiIsDPD3bvht9/t7oaERER11MgkhyCgqBjR3P/m2+srUVERKQoKBBJrhzdZgpEIiLiDRSIJFeOQLRmDaSkWFqKiIiIyykQSa4aNID69c1aRFcsDC4iIuJxFIjkqjTbTEREvIUCkVyVIxAtX272NxMREfFUCkRyVR07QrlycPQobN1qdTUiIiKuo0AkV+Xvb9YkAnWbiYiIZ1MgkjxpHJGIiHgDBSLJU8+eYLPBtm1w5IjV1YiIiLiGApHkKSQE2rY195cvt7YWERERV1Egkmvq3dvcqttMREQ8lQKRXJNjHFF0NJw/b20tIiIirqBAJNfUsiVUrw7nzsHatVZXIyIiUvgUiOSabDbNNhMREc+mQCT5cnkgstutrUVERKSwKRBJvtx+OwQEwJ9/ws6dVlcjIiJSuBSIJF/KlDGhCNRtJiIinkeBSPLN0W32zTfW1iEiIlLYFIgk3xzrEW3YACdPWluLiIhIYVIgknyrXRuaNoXMTFixwupqRERECo8CkThF0+9FRMQTKRCJU3r1MrfR0aalSERExBMoEIlTbrkFypaF48dh+3arqxERESkcCkTilBIlLk2//+47a2sREREpLApE4rSICHOrQCQiIp5CgUic5ghEP/4IZ89aW4uIiEhhUCASp9WvD3XqwMWLsHat1dWIiIhcPwUicZrNpm4zERHxLApEUiAKRCIi4kkUiKRAbr8dfHxg1y44eNDqakRERK6PApEUSIUK0LatuR8dbW0tIiIi18uvIC/asWMHGzZs4PDhw6SmphIcHEzjxo257bbbCAwMLOwapZiKiIBNm0y32cMPW12NiIhIweU7EJ06dYrZs2fz7rvv8tdff2G323O+mZ8fvXr1YtSoUdzuWL1PPFZEBEydalqIMjLA19fqikRERAomX11mb7zxBvXr1+eVV16hZ8+eLFq0iH379pGUlERaWhrx8fFs2LCB//znP5w6dYru3btzxx138Pvvv7u6frFQ27YQGAgnT0JsrNXViIiIFJzNnltTzxXq16/P5MmTuffeeylRosQ13/SPP/7ghRdeoH79+kycOLFQCrVScnIyQUFBJCUlqUvwCv36wbJl8MIL4AF/1CIi4kGc+f7OVyBKT0/Hz8/54UYZGRn4ekA/igLR1c2aBY89Bp06wZo1VlcjIiJyiTPf3/nqMitIGAI8IgxJ3hzrEcXEwJkz1tYiIiJSUNc17f7cuXOcPHkyxyHe44YboF49beMhIiLuzelAdO7cOUaPHk3lypUpV64clStXznGId+ne3dxq1WoREXFXTveFjRw5ko8++og+ffrQqFEj/P39XVFXnmbOnMnLL79MfHw8TZo0YcaMGXTs2PGq16elpTF16lT++9//kpCQQM2aNZk0aRIPa/GcQhERAbNnKxCJiIj7cjoQffXVV0RFRTF27FhX1HNNixcvZvTo0cycOZMOHTowe/Zsevbsyc6dO6lVq1aurxk4cCBHjx5l7ty51K9fn8TERNLT04u4cs/l2MZj926Ii4Or/DGIiIgUW/maZXa5ypUrs2jRIrp27eqqmvLUrl07WrVqxaxZs7LONWrUiL59+xIVFZXj+hUrVnDPPfewf/9+KlasmK/PSEtLIy0tLetxcnIyoaGhmmWWh/btYeNGmDMHhg2zuhoREREXzDK7XP/+/fnOor6RCxcusHXrViIcU5v+T0REBDExMbm+5ssvv6R169a89NJL1KhRgwYNGjB27FhSU1Ov+jlRUVEEBQVlHaGhoYX6c3gixx+Jus1ERMQdOd1lNn36dAYMGMCYMWPo1atXrq0urVq1KpTirnT8+HEyMjIICQnJdj4kJISEhIRcX7N//35+/PFHSpYsydKlSzl+/DiPPfYYJ0+eZN68ebm+JjIykjFjxmQ9drQQydVFRMBzz8GqVdrGQ0RE3I/TgSg1NZX09HRmzJjB66+/nu05u92OzWYjIyOj0ArMjc1my/Vzc5OZmYnNZmPBggUEBQUB8Oqrr3L33Xfz9ttvU6pUqRyvCQgIICAgoPAL92CXb+OxbRu0aWN1RSIiIvnndCAaNmwYmzdvZvTo0UU+yyw4OBhfX98crUGJiYk5Wo0cqlWrRo0aNbLCEJgxR3a7nUOHDnHjjTe6tGZv4ecHXbvC0qWm20yBSERE3InTgeiHH37g1VdfZfjw4a6oJ0/+/v6EhYURHR1Nv379ss5HR0dz11135fqaDh068Omnn3LmzBnKli0LwN69e/Hx8aFmzZpFUre3iIi4FIgmTbK6GhERkfxzelB1uXLlqFOnjgtKyZ8xY8YwZ84c5s2bx65du3jyySeJi4tjxIgRgBn/M2TIkKzrBw8eTKVKlXjooYfYuXMn69at4+mnn+bhhx/OtbtMCu7ybTxSUqytRURExBlOB6IhQ4awaNEiV9SSL4MGDWLGjBlMnTqVli1bsm7dOpYvX07t2rUBiI+PJy4uLuv6smXLEh0dzenTp2ndujX33Xcfffr04Y033rDqR/BY9eqZrTzS07XRq4iIuBen1yFauHAhkyZNokWLFvTu3TvXWWb9+/cvtAKLA+12n3+PPQazZsHIkfDmm1ZXIyIi3syZ72+nA5GPT96NSkUxy6yoKRDl37Jl0K8fNGgAe/ZYXY2IiHgzZ76/CzSoWuRqunQxaxDt3Qt//gkWDjcTERHJN6cDUadOnVxRh3iIoCC45RbYsAGio8GCyYgiIiJOc3pQtci1OGabrVxpbR0iIiL5la9A1KtXL2JjY/P9pmlpabz66qu8/fbbBS5M3Ff37uZ29WqzjYeIiEhxl69AVLVqVdq0aUOHDh2YPXs2e3IZLZuSksKqVat44oknqFGjBm+//TY333xzoRcsxV+bNmYbj1OnwIkcLSIiYpl8BaJ58+axefNmatasyahRo2jcuDFly5albt26NGrUiJCQECpUqECPHj345ptvmDhxIjt37qR9+/aurl+KIT8/M7gazDgiERGR4s7pafeJiYmsXLmSTZs2ceTIEVJTUwkODqZhw4Z07tyZDh06XHWjVXelaffOe+steOIJuP12+P57q6sRERFv5NJ1iLyRApHz9uyBhg3B3990nZUubXVFIiLibZz5/tYsM3GJBg2gZk24cAF+/NHqakRERPKmQCQuYbNdmm22apW1tYiIiFyLApG4TLdu5laBSEREijsFInGZrl3NbWwsHDtmbS0iIiJ5USASlwkJgebNzf3Vq62tRUREJC9OB6J169Zx5syZXJ87c+YM69atu+6ixHM4us20HpGIiBRnTgeiLl26sHPnzlyf27NnD10cK/KJcGlgdXQ0aIEHEREprpwORHktW3Tx4kV8fNQLJ5d07AglSkBcHPzxh9XViIiI5M4vPxclJydz+vTprMcJCQnExcVluyY1NZUPPviAqlWrFmqB4t7KlIH27WHtWtNKVL++1RWJiIjklK9A9NprrzF16lQAbDYb/fr1y/U6u93OxIkTC6868Qjdu5tAtGoV/OMfVlcjIiKSU74CUUREBGXLlsVutzNu3DieeOIJatWqle2agIAAmjVrRqdOnVxSqLivbt3gX/8yM80yMsDX1+qKREREsstXIAoPDyc8PByAs2fPMnz4cKpXr+7SwsRzhIVBUBCcPg1bt0LbtlZXJCIikp3TI6CfffZZhSFxip+f2fUetGq1iIgUT/lqIbrSn3/+ySeffMJff/1FampqtudsNhtz584tlOLEc3TrBkuXmoHVGmYmIiLFjdOB6JtvvqF///5kZGRQpUoVAgICsj1vs9kKrTjxHI71iGJi4OxZM/tMRESkuLDZ81pYKBctW7akYsWKLFq0iCpVqriqrmIlOTmZoKAgkpKSCAwMtLoct2S3Q506Zj2iFSugRw+rKxIREU/nzPe302OI9u3bx/jx470mDEnhsNm0jYeIiBRfTgei2rVrX3UvM5G8OLrNNLBaRESKG6cD0cSJE3nllVc4d+6cK+oRD+aYafa//8HRo9bWIiIicjmnB1X//PPPJCYmUr9+fbp06UKlSpWyPW+z2Xj99dcLrUDxHFWqQIsWJhCtXg333mt1RSIiIobTg6qvtXmrzWYjIyPjuooqbjSouvA8/TS88go8/DBodQYREXEllw6qzszMzPPwtDAkhevygdXORXERERHXcToQiVyPjh3B3x8OHoR9+6yuRkRExChwIFq5ciWRkZEMHz6cuLg4ADZv3syxY8cKrTjxPKVLQ4cO5r6m34uISHHhdCA6d+4c3bt3p2fPnrz00kvMmzeP48ePA/DKK6/w4osvFnqR4lkc3Waafi8iIsWF04Fo0qRJbNmyhc8//5ykpCQuH5MdERHBKn3LyTU41iP64QdIT7e2FhEREShAIPr00095/vnn6devH6VKlcr2XK1atbK6z0SuplUrqFABkpJgyxarqxERESlAIDp27BhNmjTJ/c18fEhNTb3uosSz+fpeWqRRDYoiIlIcOB2IatSowa+//prrc7/88gt169a97qLE82lfMxERKU6cDkT9+/fnhRdeIDY2NuuczWbjr7/+4rXXXuPvf/97oRYonskxjmjjRkhJsbYWERERpwPRs88+S/Xq1Wnbti2tW7fGZrPx0EMP0bRpU6pUqcKECRNcUad4mBtugHr14OJFWLvW6mpERMTbOR2IypUrR0xMDM8//zxly5blhhtuoHTp0kRGRrJu3bocA61FriYiwtx+9521dYiIiDi9l5k30l5mrrF0KfTvDzfdBLt3W12NiIh4GpfuZSZSWLp0MTPO9uwBrdYgIiJWKlAgWrZsGX//+99p27YtzZs3z3a0aNGisGvMYebMmdStW5eSJUsSFhbG+vXr8/W6DRs24OfnR8uWLV1boORL+fLQtq25r9lmIiJiJacD0csvv0z//v1Zt24dJUqUoFKlStmOihUruqLOLIsXL2b06NFMmjSJ2NhYOnbsSM+ePa+5IGRSUhJDhgyha9euLq1PnKNxRCIiUhw4PYaobt26dO3aldmzZ+Pr6+uquq6qXbt2tGrVilmzZmWda9SoEX379iUqKuqqr7vnnnu48cYb8fX1ZdmyZWzfvj3fn6kxRK4TE2M2e61YERITTReaiIhIYXDpGKITJ04wePBgS8LQhQsX2Lp1KxGOZoX/ExERQUxMzFVf9/777/PHH3/w7LPP5utz0tLSSE5OznaIa7RtC4GBcPIkbNtmdTUiIuKtnA5EHTp0YNeuXa6o5ZqOHz9ORkYGISEh2c6HhISQkJCQ62v27dvHhAkTWLBgAX5+fvn6nKioKIKCgrKO0NDQ665dcufnB45eTHWbiYiIVZwORDNmzODtt9/myy+/5MKFC66o6ZpsNlu2x3a7Pcc5gIyMDAYPHsxzzz1HgwYN8v3+kZGRJCUlZR0HDx687prl6hyrVisQiYiIVfLXZHKZ+vXr061bN/r164fNZqN06dLZnrfZbCQlJRVagZcLDg7G19c3R2tQYmJijlYjgJSUFLZs2UJsbCwjR44EIDMzE7vdjp+fH9999x23O3YZvUxAQAABAQEu+RkkJ0cPqGMbj3LlrK1HRES8j9OBaNy4cbz11lu0bNmSRo0a4e/v74q6cuXv709YWBjR0dH069cv63x0dDR33XVXjusDAwNzbEQ7c+ZMVq9ezWeffaaNaIsJxzYe+/ebbTzuvNPqikRExNs4HYjmz5/P+PHj85zR5UpjxozhgQceoHXr1oSHh/Puu+8SFxfHiBEjANPddfjwYT788EN8fHxo2rRpttdXqVKFkiVL5jgv1oqIgHfeMd1mCkQiIlLUnA5EGRkZdHcM+rDAoEGDOHHiBFOnTiU+Pp6mTZuyfPlyateuDUB8fPw11ySS4ufyQCQiIlLUnF6HaODAgbRs2ZKJEye6qqZiR+sQud7p01CpEmRmwl9/Qa1aVlckIiLuzpnvb6dbiJ555hkGDRpEmTJl6N27d64rU7t6tWrxPOXLQ7t2ZmB1dDQMG2Z1RSIi4k2cbiHy8TEz9XOb5u6QkZFxfVUVM2ohKhpTpsBzz8HAgbB4sdXViIiIu3NpC9HkyZPzDEMiBRURYQLRqlWQkaFtPEREpOg43ULkjdRCVDTS0804ouRk+PlnaNPG6opERMSduXQvs8ulpqZy+PBh0tPTr+dtRACzjYdjnUzNNhMRkaJUoED0ww8/EB4eTrly5ahduza//PILAI8//jhLliwp1ALFuzhWrVYgEhGRouR0IFq9ejURERGcP3+esWPHkpmZmfVccHAw8+fPL8z6xMtcuY2HiIhIUXA6EE2ePJlevXoRGxvLtGnTsj3XokULtm/fXli1iRdybONx8aLZxkNERKQoOB2IYmNjefTRR4GcU+8rV65MYmJi4VQmXkvdZiIiUtScDkR+fn5cvHgx1+cSExMpp63K5To5doZRIBIRkaLidCBq06YNH330Ua7PffbZZ4SHh193UeLdbr8dfHxgzx7QtnQiIlIUnA5EEyZMYOnSpfTr148vv/wSm83GTz/9xMiRI/nss88YN26cK+oUL+LYxgPMNh4iIiKu5nQg6tatGx988AHr169nwIAB2O12Hn/8cT7++GPmz5/Prbfe6oo6xctoHJGIiBSlAq9UnZqaSkxMDEePHiU4OJgOHTpQpkyZwq6vWNBK1UUvJgY6dICKFSExUdt4iIiI81y6l5lDqVKl6Nq1a0FfLpKntm0hMBBOnoRt27SNh4iIuJbTXWbvv/8+U6ZMyfW5KVOm8OGHH15vTSLaxkNERIqU04HojTfeoEKFCrk+FxwczBtvvHHdRYnApXFEK1daW4eIiHg+pwPR77//TtOmTXN9rnHjxuzbt++6ixIB6NHD3MbEwOnTlpYiIiIerkCbuyYlJV31vHa+l8JSrx40bAgZGZp+LyIiruV0IGrWrBmLFi3K9bmFCxfSrFmz6y5KxKFXL3O7fLm1dYiIiGdzOhA5FmAcOnQoP/30E4cPH+ann37iwQcf5PPPP+eJJ55wRZ3ipRyB6NtvITPT2lpERMRzOT3tfvDgwezevZuoqCj++9//Zp338fHhX//6F/fdd1+hFijerWNHKFsWjh6F2FgIC7O6IhER8UQFWodo6tSpPPzww0RHR3Ps2DEqV65MREQEtWvXLuz6xMv5+5vNXpcuNd1mCkQiIuIKBV6p2ptopWprzZkDw4fDLbfAxo1WVyMiIu6iSFaqBjh27Bipqak5zteqVet63lYkm549ze1PP8GxY1C5srX1iIiI5ynQtPtp06ZRpUoVqlatSt26dXMcIoWpRg1o0QLsdi3SKCIiruF0IJo3bx7/+c9/GDVqFHa7nYkTJxIZGUnNmjW58cYbmTNnjivqFC+n6fciIuJKTgeit99+OysEAfTr149p06axe/duypUrx/Hjxwu9SBFHIFqxwizUKCIiUpgKtHXHLbfcgo+PeemFCxcAKFWqFE899RTvvvtu4VYoghlQXb48nDplxhKJiIgUJqcDkZ+fGYdts9kIDAzk0KFDWc8FBwdz+PDhwqtO5P/4+V3a20zdZiIiUticDkQ33ngjBw8eBKBNmza89957XLx4kYyMDN59913q1KlT2DWKABpHJCIiruN0IOrVqxfr1q0DIDIyktWrV1O+fHkqVqzI559/zvjx4wu9SBGAO+4wt7GxcOSItbWIiIhnue6FGTdv3syiRYuw2Wz07t2bLl26FFZtxYYWZiw+2raFzZth7lx4+GGrqxERkeKsyBZmBNNt1qZNm+t9G5F86dXLBKLlyxWIRESk8BRoYUYRqzjGEX33HVy8aG0tIiLiOfLVQnT77bfn+w1tNhvff/99gQsSyUvr1mbrjmPHYMMG6NzZ6opERMQT5KuFKDMzE7vdnq8jMzPT1TWLF/PxubS3mWabiYhIYdFu9/mgQdXFy+LFcM890KQJ/Pab1dWIiEhx5cz3t8YQiduJiDAtRTt2wF9/WV2NiIh4ggLNMsvIyOCTTz7hhx9+4MSJE1SqVIkuXbrw97//PWslaxFXqVAB2reHH3803Wb/+IfVFYmIiLtzuoXo+PHjtGvXjvvuu4/58+cTExPD/Pnzue+++2jXrp02d5UioVWrRUSkMDkdiJ588kn27NnDggULSE1NJT4+ntTUVP773/+yb98+nnzySVfUKZKNIxB9/z2cP29tLSIi4v6cDkRfffUV06ZN495778XX1xcAX19fBg8ezNSpU/nqq68KvUiRKzVvDtWrQ2oqrF1rdTUiIuLunA5EdrudJk2a5Ppc06ZNKYpJazNnzqRu3bqULFmSsLAw1q9ff9VrlyxZQvfu3alcuTKBgYGEh4ezcuVKl9cormWzqdtMREQKj9OBqFu3bqxatSrX56Kjo+ns4pXyFi9ezOjRo5k0aRKxsbF07NiRnj17EhcXl+v169ato3v37ixfvpytW7fSpUsX+vTpQ2xsrEvrFNdTIBIRkcLi9DpE27dvp3///gwYMIDBgwdTtWpVEhISWLBgAUuWLGHJkiXUqlUr6/qKFSsWasHt2rWjVatWzJo1K+tco0aN6Nu3L1FRUfl6jyZNmjBo0CAmT56cr+u1DlHxlJwMwcFmC4+9e+HGG62uSEREihOXbu7aqlUrAKZPn86rr76add6Rq8LCwrJdn5GR4exHXNWFCxfYunUrEyZMyHY+IiKCmJiYfL1HZmYmKSkpeQa1tLQ00tLSsh4nJycXrGBxqcBA6NgRVq82rUT//KfVFYmIiLtyOhBNnjwZm83milqu6fjx42RkZBASEpLtfEhICAkJCfl6j+nTp3P27FkGDhx41WuioqJ47rnnrqtWKRq9eplA9PXXCkQiIlJwTgeiKVOmuKAM51wZyOx2e75C2sKFC5kyZQpffPEFVapUuep1kZGRjBkzJutxcnIyoaGhBS9YXKZPHxg7FtasgVOnzKKNIiIiziqUrTvOnz/P7t27C7V7LDfBwcH4+vrmaA1KTEzM0Wp0pcWLFzNs2DA++eQTunXrlue1AQEBBAYGZjukeGrQABo3hvR0+OYbq6sRERF35XQgevPNN3n++eezHm/dupXQ0FCaNGlCgwYNOHjwYKEWeDl/f3/CwsKIjo7Odj46Opr27dtf9XULFy7kwQcf5OOPP6Z3794uq0+s0b+/uV2yxNo6RETEfTkdiObMmUP58uWzHo8fP56KFSvy2muvYbfbmTZtWmHWl8OYMWOYM2cO8+bNY9euXTz55JPExcUxYsQIwHR3DRkyJOv6hQsXMmTIEKZPn84tt9xCQkICCQkJJCUlubROKTqOQLRiBZw7Z20tIiLinpweQxQXF0fDhg0BSElJYd26dSxatIj+/ftToUKFfE9lL6hBgwZx4sQJpk6dSnx8PE2bNmX58uXUrl0bgPj4+GxrEs2ePZv09HQef/xxHn/88azzQ4cOZf78+S6tVYpGy5ZQu7bZ+X7lSujXz+qKRETE3TgdiNLS0ihRogQAGzduJDMzM2tMTp06dfI92+t6PPbYYzz22GO5PndlyFmzZo3L6xFr2Wymlei112DpUgUiERFxntNdZrVq1craKuOLL76gZcuWWYOOjx07pgHIYglHCPrqK7NQo4iIiDOcDkT3338/U6dOJSwsjNmzZ3P//fdnPbdlyxYaNGhQqAWK5Ef79lClCpw+babgi4iIOMPpLrNJkybh5+dHTEwM/fr1Y9SoUVnP/fbbbwwYMKBQCxTJD19fuOsueO89M9use3erKxIREXfi9F5m3kh7mbmHFSugZ0+oWhUOHwafQlllS0RE3JUz39/6yhCPcfvtZn+zhATYtMnqakRExJ3kq8vs4Ycf5plnnqFu3bo8/PDDeV5rs9mYO3duoRQn4gx/f7jzTvj4Y9NtlsdanSIiItnkKxD98MMP/PP/ds5cvXp1nvuGWbXxqwiY6fcff2ym37/8spmSLyIici35CkQHDhzIuv/nn3+6qhaR63bHHVCyJOzfD7/8Ai1aWF2RiIi4A40hEo9Spgz06GHuL11qbS0iIuI+FIjE4zgWadRmryIikl/56jLz8fFxamxQRkZGgQsSuV59+ph1iX79FX7/HerXt7oiEREp7vIViCZPnpwtEL3//vucOXOGPn36ULVqVeLj4/n6668pU6bMNWehibhaxYrQuTN8/73pNnv6aasrEhGR4i5fgWjKlClZ96dPn07VqlVZtWoVZcuWzTqfkpJCt27dKF26dKEXKeKs/v0ViEREJP+cHkM0c+ZMxo0bly0MAZQrV45x48Yxc+bMQitOpKD69jW3GzfCkSOWliIiIm7A6UB0+PBh/Pxyb1jy8/MjISHhuosSuV7Vq8Mtt5j7X3xhbS0iIlL8OR2IGjVqxKuvvsrFixeznb9w4QLTp0+nYcOGhVacyPXo39/caraZiIhci9O73U+bNo2+fftSr149+vfvT9WqVUlISGDJkiUkJCSwbNkyF5Qp4rx+/WDcOPjhBzh50gy2FhERyY3TLUS9e/dmxYoV1KhRg7fffptJkybx1ltvUbNmTb799lt69+7tijpFnFa/PjRrBhkZ8PXXVlcjIiLFmdMtRABdu3ala9eunDt3jlOnTlGhQgXNLpNiqV8/sx7RkiUwZIjV1YiISHF1XStVly5dmho1aigMSbHlGEe0ciWcPWttLSIiUnxp6w7xaM2bQ926cP48rFhhdTUiIlJcKRCJR7PZNNtMRESuTYFIPN6AAeb2iy/UbSYiIrlTIBKPd8stcMMNJgxpVQgREcmNApF4PJsN7r/f3P/oI2trERGR4kmBSLyCIxBFR0N8vLW1iIhI8aNAJF6hfn0ID4fMTFi40OpqRESkuFEgEq/xwAPm9sMPra1DRESKHwUi8RoDB0KJEvC//5nVq0VERBwUiMRrVKoEjq32NLhaREQup0AkXsXRbbZggdn0VUREBAq4uauIu+rdGypUgCNH4IcfoFs3qysSsc6xY7B/PyQmmvuJidnvO26Tk8Fuz/0Ac+vnB+XLm/++rnYbHAy1apmjZk3ThS1SXCgQiVcJCIBBg+Cdd0y3mQKReIPUVNi504yd+/VX+OUXc3v0aOF9xoULcO6c+cdGfthsUL061K5tjlq1zG2dOtCwobnvoz4MKUI2u92R8eVqkpOTCQoKIikpicDAQKvLkesUEwMdOkCZMuYLoUwZqysSKTx2O+zYAStXwk8/mfCzb59ZciI3NWtCSAhUqWKOypVz3g8MNOHEZsv9ALh4EU6fNsepU+Zw3HfcJiZCXJw5LlzI++coVcoEo0aNsh/164O/f+H9vsSzOfP9rRYi8Trh4WYrjz/+MFt53Hef1RWJXJ+TJ2HVKlixAr77Dg4fznlNpUrQrBk0b25umzWDJk2gbNmirzcz81I4+uuv7Ld//AF795pWrdhYc1zOz8+EopYtzXHzzea2SpWi/znEs6iFKB/UQuR5pkyB556DiAjzL2kRd5KRAT//bP7urlgBmzdnbwEqWRI6dYLbbzdhoVkzqFr1UmtOcZeeDgcOwK5d5ti509zu3g0pKbm/plq17AHp5pvNP3zc5WcW13Dm+1uBKB8UiDzP77/DjTeaboBDh8z/TEWKu337YN48+OCDnFvQNGkCPXqYo2NH0+Xkaex20/r122+wffulY+/eSwO8L1e+PISFQevW0KaNua1VSyHJmygQFTIFIs/Uvj1s3AivvAJPPWV1NSK5O3sWPvsM5s6F9esvnS9fHrp3vxSCata0rETLnTljBolfHpL+9z9IS8t5beXKJhg5QlKbNqb1TDyTAlEhUyDyTLNmwWOPQYsW5n+gIsWF3W66webONXvvObqJfHzgjjtg2DC4804NLs7LxYtmcPnmzbBlizl++cV0x10pNBTatr10hIVBuXJFX7MUPgWiQqZA5JlOnDBdZRcvmv9RNmtmdUXi7c6fN11is2aZbiGHevXg4Ydh6FDvbgm6XufPm//Wt2wxQennn83YpCu/BW02M6OtbdtLrUjNm5tlO8S9KBAVMgUiz9Wvn5lp9vTT8NJLVlcj3io1Fd57D1588dI6PiVLwt13m9ag227TmjyukpIC27aZcOQ44uJyXleihGlNvryrrVEjM+tNii8FokKmQOS5liyBAQPMAnFxceDra3VF4k3OnYN33zVBKCHBnAsNNQH9gQfMOCEpekePXmpB+vln06J04kTO60qXNrPZ2rSBVq3McdNNCknFiQJRIVMg8lxpaabb7NQpiI7WytVSNM6eNaulv/zypdWia9WCiRPhwQfVNVPc2O3w558mJDnGJG3dmvsSAKVKme41R0Bq1crMANSfqTUUiAqZApFn+8c/zJfTkCFmOrOIq5w9CzNnmpmNiYnmXJ06JggNHapB0u4kMxP27LkUkGJjzeSMM2dyXluiBDRubILS5QtjVqumJQBczeMD0cyZM3n55ZeJj4+nSZMmzJgxg44dO171+rVr1zJmzBh27NhB9erVGTduHCNGjMj35ykQeTZt5SGuZrfDokVmeQfH+kH16sGkSaZrTJuceobMTLPG2bZt2Y9Tp3K/3rF6uGMF8aZNzbikoKCirduTeXQgWrx4MQ888AAzZ86kQ4cOzJ49mzlz5rBz505q1aqV4/oDBw7QtGlThg8fzqOPPsqGDRt47LHHWLhwIQMGDMjXZyoQeTa73SzS+McfZsPX+++3uiLxJDt3wsiR8MMP5nG9evDMM2bLGAUhz2e3m/GJsbHZN9fNa3+5qlXNPm6O46abzG2tWhpc7yyPDkTt2rWjVatWzJo1K+tco0aN6Nu3L1FRUTmuHz9+PF9++SW7du3KOjdixAj+97//sXHjxlw/Iy0tjbTLVvRKTk4mNDRUgciDObby6NbNjCUSuV5nzsDUqfDaa2btm1KlTIvQ2LEaTyJmZuGuXZcC0q+/mqUWrlyB/HKlSpl/vNWrZ466dc1Rr57pevXE1cmvl8du7nrhwgW2bt3KhAkTsp2PiIggJiYm19ds3LiRiIiIbOd69OjB3LlzuXjxIiVy+SdaVFQUzz33XOEVLsXe0KHmy2vVKrNfUsOGVlck7spuh08/hTFjLm2yetddMGOG+dISARNeHIOuL5ecbMYm7d5tDsf9fftMiPrlF3Pkplq1SyGpZk2oUcPMoK1RwxzVqmkGXF7c6ldz/PhxMjIyCAkJyXY+JCSEBMec1SskJCTken16ejrHjx+nWi6bWEVGRjJmzJisx44WIvFcdetCnz7w5Zfw5pvw9ttWVyTuaM8e0z22apV5XK8evPEG9O5tbV3iPgIDL61zdLn0dDPTbe9es/Ht/v2XbvfvNzPe4uPNcZX2AWw2CAm5FJSqVIHg4EtH5crZHwcGetegb7cKRA62K/6E7HZ7jnPXuj638w4BAQEEqE3b6/zznyYQffABvPCC1oCR/LtwAZ5/3qwndPGi6RKLjITx480CiyLXy88P6tc3x5Xsdjh58lJA+vNP0zp5+REfb0JVQoI5tm7N32cGBuZ9lCtn/o7ndQQEmDXefH3NGKir3S9f3toB5W4ViIKDg/H19c3RGpSYmJijFcihatWquV7v5+dHpUqVXFaruJ8uXcwsj99+M3tIacNXyY9du8wA6dhY87h3b9MqVK+etXWJ97DZzIy1SpXMStq5ycyEY8cuBaQjR+D48ezHsWOX7p89awLUyZPmKAoTJkAuQ4GLjFsFIn9/f8LCwoiOjqZfv35Z56Ojo7nrrrtyfU14eDhfffVVtnPfffcdrVu3znX8kHgvmw1GjYL/9//grbdg9GitXC1XZ7ebrtWnnzZ7ZFWqZNazGjDAu7oZxD34+JjuspCQnOOWcpOaalbnTkkx45qudqSkmL//1zoyMkwoy8i4+n2r1+Fyu1lmjmn377zzDuHh4bz77ru899577Nixg9q1axMZGcnhw4f58MMPgUvT7h999FGGDx/Oxo0bGTFihKbdS67OnTNbJ5w8abb1uCx3i2Q5csRstrpypXncowe8/74ZtCoixYcz399ut6LBoEGDmDFjBlOnTqVly5asW7eO5cuXU7t2bQDi4+OJu2xnvrp167J8+XLWrFlDy5Ytef7553njjTfyHYbEu5QubVqIwHR7iFxpyRKzkN7KlWZ8xJtvwrffKgyJuDu3ayGyglqIvMvBg2bWWUYG/O9/ZgVZkeRkM/B+/nzz+OabYcECs7KwiBRPHt1CJOJqoaHQv7+5//rr1tYixUNMDLRsacKQzWYGf27apDAk4kkUiERy8c9/mtsFC8zMC/FOdrtZUPG228yU5tq1Yc0aMxPG6gGgIlK4FIhEctG+PYSFQVoavPee1dWIFc6ehcGD4cknTffpPfeYLtTbbrO6MhFxBQUikVzYbJdaiWbONIvtiffYtw9uucXsUO/nZ1qJPv5Yu5CLeDIFIpGrGDjQrNlx+DB8/rnV1UhR+eors7jdb7+ZXcdXrzbhWGsLiXg2BSKRqwgIgBEjzH0NrvZ8GRnwzDPwt7+ZGWUdOsC2bdCxo9WViUhRUCASycOIEVCihJlR9PPPVlcjrnLypNlyY9o08/iJJ0zLkNYWEvEeCkQieaha1QymBS3U6KliY80A+pUroVQp+Ogj82etWWQi3kWBSOQaRo0yt598YnaMFs+xZInpGvvzT7MZ68aNcP/9VlclIlZQIBK5htatzTT8ixdh1iyrq5HCYLfDf/5jNmJNTYU77oAtW6BFC6srExGrKBCJ5INjCv4775idm8V9XbgAw4ZBZKR5PHKkmVlWoYK1dYmItRSIRPKhXz+oWdOsWr1okdXVSEGdOAEREWZneh8fszHrm2+atYZExLspEInkQ4kS8Nhj5v706WaKtriXvXshPBzWroVy5eDrr03rkIgIKBCJ5Nujj5qVin/7zaxaLO5jzRqz8vS+fWY/spgY6NnT6qpEpDhRIBLJp4oVzS7nAP/6l8YSuYv33zfdZKdOQbt28NNP0LSp1VWJSHGjQCTihFGjoEYNiIsze5xJ8ZWZaQZOP/ywmSE4aBD88IPZjkVE5EoKRCJOKF0annvO3H/hBTh92tJy5CrS0sx6Qv/5j3n8zDOmm7NUKWvrEpHiS4FIxElDh0KjRma7hxdftLoaudLp02ZdoYULzeyxDz6AqVPNrDIRkavR/yJEnOTnd6nlYcYMOHzY0nLkMgcPwq23mkHU5crB8uUwZIjVVYmIO1AgEimAPn3Mlg/nz8OUKVZXIwC//GJmku3YAdWrw/r10L271VWJiLtQIBIpAJsNXnrJ3J83D3butLYeb/f996Zl6MgRaNzY7EmmbThExBkKRCIF1L499O1rZjNNnGh1Nd7ro4/MmKGUFOjUCX78EWrVsroqEXE3CkQi1+Hf/zaDdb/4AjZssLoa72K3Q1SUGSOUng733AMrV2pPMhEpGAUikevQqJHZKBRg3DjzJS2ul54O//jHpZa5p5+GBQsgIMDaukTEfSkQiVynZ58169vExMCXX1pdjec7e9Zstjt7thnL9cYbZjyXptWLyPXQ/0JErlONGjB6tLkfGWlaL8Q1jh6FLl3MxqwlS8Lnn8MTT1hdlYh4AgUikUIwfrzZ62zXLpg/3+pqPNPevWYg++bNUKkSrF5tWopERAqDApFIIQgKMhu+gulCO3fO2no8TUyMCUP790O9euZxeLjVVYmIJ1EgEikkjz0GtWubtXBef93qajzH0qXQtSucOAFt2pg1hho0sLoqEfE0CkQihSQgAKZNM/eff950n8n1efNNGDDArAh+551mt/oqVayuSkQ8kQKRSCEaPNhsF5GaCvfea3ZdF+dlZpqp9KNGmaUMRowwLUVlylhdmYh4KgUikULk42N2Vw8Ohv/9z8w6E+ekpppg+cor5nFUFMycaTbVFRFxFQUikUJWrRq8/765/9prsGKFtfW4k/h46NwZFi+GEiXMthwTJpj1hkREXEmBSMQF7rwTRo4094cONevnSN62b4e2beHnn80SBt99B/ffb3VVIuItFIhEXOSll6BpU0hMhAcfNONiJHfLlkGHDnDoEDRsCD/9ZFqKRESKigKRiIuUKgULF5oVlVesMDOmJDu7HV58Efr3N2s3de9uptXXr291ZSLibRSIRFyoaVOYPt3cHzfOdAuJkZYGDz1kxgjZ7fD447B8OZQvb3VlIuKNFIhEXOwf/4A+feDCBTMVX6tYw7FjZrHFDz4AX1946y1zaCaZiFhFgUjExWw2mDfPzD7bvRvGjLG6Imv99psZPL1hg9nyZPly0zokImIlBSKRIhAcDB9+aO7Pnm0WGfRGH39s9iD780+44QbYtAkiIqyuSkREgUikyHTrZlZfBnjkETOjylukpsLw4XDffXDmDHTpYmaSNWxodWUiIoYCkUgRmjYNwsLg5Em4+25ISrK6Itfbvdt0kc2ZY7oPJ0+G6GioVMnqykRELlEgEilC/v5mKn758qaF5Pbb4fhxq6tynQ8/NAHwt98gJMQEoeeeMwOpRUSKE7cKRKdOneKBBx4gKCiIoKAgHnjgAU6fPn3V6y9evMj48eNp1qwZZcqUoXr16gwZMoQjR44UXdEiV7jxRrNre+XKsG0bdOoEnvZX8uxZM6V+6FAzq65rV7PkQNeuVlcmIpI7twpEgwcPZvv27axYsYIVK1awfft2Hnjggatef+7cObZt28YzzzzDtm3bWLJkCXv37uVvf/tbEVYtklPLlrBuHdSoATt3QseOZqCxJ9ixw3SRzZ9vNrudOhVWroSqVa2uTETk6mx2u91udRH5sWvXLho3bsymTZto164dAJs2bSI8PJzdu3dz00035et9Nm/eTNu2bfnrr7+oVatWvl6TnJxMUFAQSUlJBAYGFvhnELnSgQNmsPX+/SYcrVrlvgON7XaYOxdGjTKDqKtVM92DnTpZXZmIeCtnvr/dpoVo48aNBAUFZYUhgFtuuYWgoCBiYmLy/T5JSUnYbDbK57EcblpaGsnJydkOEVeoWxfWr4fGjeHwYbjtNvdczfrXX03wGT7chKGICPNzKAyJiLtwm0CUkJBAlSpVcpyvUqUKCQkJ+XqP8+fPM2HCBAYPHpxnUoyKisoapxQUFERoaGiB6xa5lurVYe1aaNXKrODcpYvZz8sdJCebhSZvvtkEu9Kl4eWX4dtvIZf/XEVEii3LA9GUKVOw2Wx5Hlu2bAHAZrPleL3dbs/1/JUuXrzIPffcQ2ZmJjNnzszz2sjISJKSkrKOgwcPFuyHE8mn4GBYvdrs+H76tNnk9Pvvra7q6ux2s8jiTTfBa69BRgYMGAC7dsHYsWbskIiIO7F856CRI0dyzz335HlNnTp1+OWXXzh69GiO544dO0ZISEier7948SIDBw7kwIEDrF69+pr9iAEBAQQEBFy7eJFCFBRkBh/362emp/fubfb6GjjQrN9TXOzYYbbaWLvWPL7xRnjzTejRw9q6RESuh+WBKDg4mODg4GteFx4eTlJSEj///DNt27YF4KeffiIpKYn27dtf9XWOMLRv3z5++OEHKmk1OCnGypSBr76Ce+6BZcvM7VtvwYsvQh5/zYtESopZQ+j11yE9HUqVgn/9C556CvTvBxFxd27TsN2oUSPuuOMOhg8fzqZNm9i0aRPDhw/nzjvvzDbDrGHDhiz9v42i0tPTufvuu9myZQsLFiwgIyODhIQEEhISuHDhglU/ikieAgLgk08gMhJKloQffzRdaX37min6Re3332H8eKhXD6ZPN2GoXz/TPTZxosKQiHgGtwlEAAsWLKBZs2ZEREQQERFB8+bN+eijj7Jds2fPHpL+bz+EQ4cO8eWXX3Lo0CFatmxJtWrVsg5nZqaJFLUSJeDf/4Z9+8y+Zz4+8MUX0KwZDBsGrh7WduGCCWXdupkusZdeMitq33CD2Z1+yRKoXdu1NYiIFCW3WYfISlqHSKy2axdMmgT/1/hJQIBZ72fCBKhYsfA+5/ff4b334P33zYw3MOOX7rgDHn3UjGvys7yjXUQkf5z5/lYgygcFIikuNm403Vfr15vH5cvD4MHQvDk0bQpNmphz+WG3w4kTphVq1y6ziOKqVZeer1bNtEY98ohag0TEPSkQFTIFIilO7HbTbRUZaRZEvFLNmiYYNW16KSSlp5vWn337st9euRWgzWZmizlag0qUKJIfSUTEJRSICpkCkRRHGRlmXNHGjWY3+d9+g0OHnH+fmjWhfn249VbTIlSnTqGXKiJiCWe+vzUaQMRN+fpC//7mcDh92sxEcwSk334z6wb5+5vB0fXrZ7+tV8+sLi0i4u0UiEQ8SPnyZr0iq9csEhFxN2417V5ERETEFRSIRERExOspEImIiIjXUyASERERr6dAJCIiIl5PgUhERES8ngKRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEImIiIjXUyASERERr6dAJCIiIl5PgUhERES8np/VBbgDu90OQHJyssWViIiISH45vrcd3+N5USDKh5SUFABCQ0MtrkRERESclZKSQlBQUJ7X2Oz5iU1eLjMzkyNHjlCuXDlsNluhvndycjKhoaEcPHiQwMDAQn1vuUS/56Kh33PR0O+5aOj3XDRc+Xu22+2kpKRQvXp1fHzyHiWkFqJ88PHxoWbNmi79jMDAQP0HVwT0ey4a+j0XDf2ei4Z+z0XDVb/na7UMOWhQtYiIiHg9BSIRERHxegpEFgsICODZZ58lICDA6lI8mn7PRUO/56Kh33PR0O+5aBSX37MGVYuIiIjXUwuRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEFlo5syZ1K1bl5IlSxIWFsb69eutLsnjrFu3jj59+lC9enVsNhvLli2zuiSPExUVRZs2bShXrhxVqlShb9++7Nmzx+qyPNKsWbNo3rx51gJ24eHhfPvtt1aX5fGioqKw2WyMHj3a6lI8ypQpU7DZbNmOqlWrWlaPApFFFi9ezOjRo5k0aRKxsbF07NiRnj17EhcXZ3VpHuXs2bO0aNGCt956y+pSPNbatWt5/PHH2bRpE9HR0aSnpxMREcHZs2etLs3j1KxZk//85z9s2bKFLVu2cPvtt3PXXXexY8cOq0vzWJs3b+bdd9+lefPmVpfikZo0aUJ8fHzW8euvv1pWi6bdW6Rdu3a0atWKWbNmZZ1r1KgRffv2JSoqysLKPJfNZmPp0qX07dvX6lI82rFjx6hSpQpr167ltttus7ocj1exYkVefvllhg0bZnUpHufMmTO0atWKmTNnMm3aNFq2bMmMGTOsLstjTJkyhWXLlrF9+3arSwHUQmSJCxcusHXrViIiIrKdj4iIICYmxqKqRApHUlISYL6oxXUyMjJYtGgRZ8+eJTw83OpyPNLjjz9O79696datm9WleKx9+/ZRvXp16tatyz333MP+/fstq0Wbu1rg+PHjZGRkEBISku18SEgICQkJFlUlcv3sdjtjxozh1ltvpWnTplaX45F+/fVXwsPDOX/+PGXLlmXp0qU0btzY6rI8zqJFi9i2bRubN2+2uhSP1a5dOz788EMaNGjA0aNHmTZtGu3bt2fHjh1UqlSpyOtRILKQzWbL9thut+c4J+JORo4cyS+//MKPP/5odSke66abbmL79u2cPn2azz//nKFDh7J27VqFokJ08OBB/vnPf/Ldd99RsmRJq8vxWD179sy636xZM8LDw7nhhhv44IMPGDNmTJHXo0BkgeDgYHx9fXO0BiUmJuZoNRJxF0888QRffvkl69ato2bNmlaX47H8/f2pX78+AK1bt2bz5s28/vrrzJ492+LKPMfWrVtJTEwkLCws61xGRgbr1q3jrbfeIi0tDV9fXwsr9ExlypShWbNm7Nu3z5LP1xgiC/j7+xMWFkZ0dHS289HR0bRv396iqkQKxm63M3LkSJYsWcLq1aupW7eu1SV5FbvdTlpamtVleJSuXbvy66+/sn379qyjdevW3HfffWzfvl1hyEXS0tLYtWsX1apVs+Tz1UJkkTFjxvDAAw/QunVrwsPDeffdd4mLi2PEiBFWl+ZRzpw5w++//571+MCBA2zfvp2KFStSq1YtCyvzHI8//jgff/wxX3zxBeXKlctq+QwKCqJUqVIWV+dZJk6cSM+ePQkNDSUlJYVFixaxZs0aVqxYYXVpHqVcuXI5xsCVKVOGSpUqaWxcIRo7dix9+vShVq1aJCYmMm3aNJKTkxk6dKgl9SgQWWTQoEGcOHGCqVOnEh8fT9OmTVm+fDm1a9e2ujSPsmXLFrp06ZL12NEvPXToUObPn29RVZ7FsXRE586ds51///33efDBB4u+IA929OhRHnjgAeLj4wkKCqJ58+asWLGC7t27W12aiNMOHTrEvffey/Hjx6lcuTK33HILmzZtsux7UOsQiYiIiNfTGCIRERHxegpEIiIi4vUUiERERMTrKRCJiIiI11MgEhEREa+nQCQiIiJeT4FIREREvJ4CkYiIiHg9BSIRKRZiYmKYMmUKp0+fzvFc586dc6yEXRx8+OGHVK5cmZSUlHy/Zu/evfj7+7Nt2zYXViYiztJK1SJSLLzyyis8/fTTHDhwgDp16mR7bufOnQA0btzYgspyd+7cORo0aMDo0aMZO3asU6996KGH2L9/P2vXrnVRdSLiLLUQiUix17hx42IVhgA++OADTpw4wSOPPOL0a0eOHMm6deuIiYlxQWUiUhAKRCJiuSlTpvD0008DULduXWw2GzabjTVr1gA5u8z+/PNPbDYbL7/8Mi+++CJ16tShVKlSdO7cmb1793Lx4kUmTJhA9erVCQoKol+/fiQmJub43MWLFxMeHk6ZMmUoW7YsPXr0IDY2Nl81z5o1iz59+lC+fPls5z/99FPatWtHUFAQpUuXpl69ejz88MPZrgkLC6NRo0a88847+f8liYhLKRCJiOUeeeQRnnjiCQCWLFnCxo0b2bhxI61atcrzdW+//TYbNmzg7bffZs6cOezevZs+ffowbNgwjh07xrx583jppZdYtWpVjpacf//739x77700btyYTz75hI8++oiUlBQ6duyY1UV3NYcOHeLXX3+lS5cu2c5v3LiRQYMGUa9ePRYtWsQ333zD5MmTSU9Pz/EenTt35ttvv0WjFkSKBz+rCxARqVmzJrVq1QLg5ptvzjGG6GrKly/PsmXL8PEx/7Y7fvw4o0ePpmHDhnzxxRdZ1+3evZsZM2aQnJxMYGAgBw8e5Nlnn2XkyJG88cYbWdd1796dG2+8keeee47Fixdf9XMdXV1XBraYmBjsdjvvvPMOQUFBWecffPDBHO/RqlUrZs2axZ49e2jYsGG+fl4RcR21EImI2+rVq1dWGAJo1KgRAL179852neN8XFwcACtXriQ9PZ0hQ4aQnp6edZQsWZJOnTplddVdzZEjRwCoUqVKtvNt2rQBYODAgXzyySccPnz4qu/heG1e14hI0VEgEhG3VbFixWyP/f398zx//vx5AI4ePQqYAFOiRIlsx+LFizl+/Hien5uamgpAyZIls52/7bbbWLZsWVbYqlmzJk2bNmXhwoU53sPxWsd7iYi11GUmIl4nODgYgM8++4zatWsX+PUnT56kWrVq2Z676667uOuuu0hLS2PTpk1ERUUxePBg6tSpQ3h4eNZ1J0+ezPZeImItBSIRKRYCAgKAomkx6dGjB35+fvzxxx8MGDDA6dc7xvz88ccfNGnSJNdrAgIC6NSpE+XLl2flypXExsZmC0T79+/Hx8eHm266qWA/hIgUKgUiESkWmjVrBsDrr7/O0KFDKVGiBDfddBPlypUr9M+qU6cOU6dOZdKkSezfv5877riDChUqcPToUX7++WfKlCnDc889d9XXt2vXjlKlSrFp0yb+9re/ZZ2fPHkyhw4domvXrtSsWZPTp0/z+uuvU6JECTp16pTtPTZt2kTLli2pUKFCof98IuI8jSESkWKhc+fOREZG8tVXX3HrrbfSpk0btm7d6rLPi4yM5LPPPmPv3r0MHTqUHj16MG7cOP766y9uu+22PF/r7+/P3XffnW0mG5iglJCQwPjx44mIiOD//b//R6lSpVi9enW2lqQzZ87w/fffc99997nkZxMR52nrDhGRAtiyZQtt2rRh06ZNtGvXzqnXzp07l3/+858cPHhQLUQixYQCkYhIAQ0aNIizZ8/y9ddf5/s16enpNG7cmKFDhzJp0iQXVicizlCXmYhIAU2fPp02bdo4tdv9wYMHuf/++3nqqadcWJmIOEstRCIiIuL11EIkIiIiXk+BSERERLyeApGIiIh4PQUiERER8XoKRCIiIuL1FIhERETE6ykQiYiIiNdTIBIRERGv9/8BWAXjX6STu+AAAAAASUVORK5CYII=", "application/papermill.record/text/plain": "
" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "spring_plot" } }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# Define Euler method function\n", "def euler(f, tspan, y0, h):\n", "\n", " nsteps = int((tspan[1] - tspan[0]) / h)\n", " t = np.zeros(nsteps + 1)\n", " y = np.zeros((nsteps + 1, len(y0)))\n", " t[0] = tspan[0]\n", " y[0,:] = y0\n", " \n", " for n in range(nsteps):\n", " y[n+1,:] = y[n,:] + h * f(t[n], y[n,:])\n", " t[n+1] = t[n] + h\n", "\n", " return t, y\n", "\n", "\n", "# Define Mass-Spring model\n", "def spring(t, y):\n", " return np.array([ y[1], (- c * y[1] - k * y[0]) / m ])\n", " \n", "\n", "# Define IVP parameters\n", "tspan = [0, 5] # boundaries of the t domain\n", "y0 = [1, 0] # initial values\n", "h = 0.1 # step length\n", "m = 1; # mass of object\n", "c = 2; # damping coefficient\n", "k = 4; # spring constant\n", "\n", "# Solve the IVP using the Euler method\n", "t, y = euler(spring, tspan, y0, h)\n", "\n", "# Plot solution\n", "fig, ax = plt.subplots()\n", "plt.plot(t, y[:,0], \"b-\")\n", "plt.xlabel(\"time (s)\", fontsize=12)\n", "plt.ylabel(\"displacement (m)\", fontsize=12)\n", "plt.show()\n", "\n", "from myst_nb import glue\n", "glue(\"spring_plot\", fig, display=False)" ] }, { "cell_type": "markdown", "id": "8ac52498", "metadata": {}, "source": [ "## Code\n", "\n", "The code below sets up and solves the mass-sprint-damper model example from {prf:ref}`spring-example` using the Euler method.\n", "\n", "`````{tab-set}\n", "````{tab-item} Python\n", "```python\n", "# Define mass-spring-damper model\n", "def spring(t, y):\n", " return np.array([ y[1], (- c * y[1] - k * y[0]) / m ])\n", " \n", "\n", "# Define IVP parameters\n", "tspan = [0, 5] # boundaries of the t domain\n", "y0 = [1, 0] # initial values\n", "h = 0.1 # step length\n", "m = 1; # mass of object\n", "c = 2; # damping coefficient\n", "k = 4; # spring constant\n", "\n", "# Solve the IVP using the Euler method\n", "t, y = euler(spring, tspan, y0, h)\n", "print(y)\n", "\n", "# Plot solution\n", "fig, ax = plt.subplots()\n", "plt.plot(t, y[:,0], \"b-\")\n", "plt.xlabel(\"time (s)\", fontsize=12)\n", "plt.ylabel(\"displacement (m)\", fontsize=12)\n", "plt.show()\n", "```\n", "````\n", "\n", "````{tab-item} MATLAB\n", "```matlab\n", "% Define mass-spring-damper model\n", "spring = @(t, y, m, c, k) [ y(2), (-c * y(2) - k * y(1)) / m];\n", "\n", "% Define IVP parameters\n", "tspan = [0, 5]; % boundaries of the t domain\n", "y0 = [1, 0]; % initial values [S, I, R]\n", "h = 0.1; % step length\n", "m = 1; % mass of object\n", "c = 2; % damping coefficient\n", "k = 4; % spring constant\n", "\n", "% Solve the IVP using the Euler method\n", "[t, y] = euler(@(t, y)spring(t, y, m, c, k), tspan, y0, h);\n", "\n", "% Plot solution\n", "plot(t, y(:, 1), 'b-', LineWidth=2)\n", "xlabel('time (s)', Fontsize=14, Interpreter='latex')\n", "ylabel('displacement (m)', Fontsize=14, Interpreter='latex')\n", "axis padded\n", "```\n", "````\n", "`````" ] } ], "metadata": { "celltoolbar": "Tags", "kernelspec": { "display_name": "base", "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.11.5" } }, "nbformat": 4, "nbformat_minor": 5 }