Created
May 2, 2020 12:49
-
-
Save prl900/0358908e3c9f3ffac649db59f16730b4 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Numerical methods for solving differential equations" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import numpy as np\n", | |
"\n", | |
"# x'=sin(t)^2*x\n", | |
"def dxdt(t,x):\n", | |
" return np.square(np.sin(t))*x" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"### Eulers method: \n", | |
"\n", | |
"$ x_{t+1} = x_t + \\Delta t x'_t $" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 7, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"(10,)\n", | |
"0.25\n", | |
"[0. 0.25 0.5 0.75 1. 1.25 1.5 1.75 2. 2.25 2.5 2.75 3. 3.25\n", | |
" 3.5 3.75 4. 4.25 4.5 4.75]\n", | |
"(20,)\n", | |
"(10,) (20,)\n" | |
] | |
} | |
], | |
"source": [ | |
"Δt = 0.5\n", | |
"t = np.arange(0,5,Δt)\n", | |
"x = np.zeros(t.shape)\n", | |
"print(x.shape)\n", | |
"\n", | |
"#ci\n", | |
"x[0] = 3\n", | |
"\n", | |
"for i in range(1,x.shape[0]):\n", | |
" x[i] = x[i-1] + Δt*dxdt(t[i-1],x[i-1])\n", | |
" \n", | |
"\n", | |
"\n", | |
"Δt = 0.25\n", | |
"print(Δt)\n", | |
"th = np.arange(0,5,Δt)\n", | |
"print(th)\n", | |
"xh = np.zeros(th.shape)\n", | |
"print(xh.shape)\n", | |
"\n", | |
"#ci\n", | |
"xh[0] = 3\n", | |
"\n", | |
"for i in range(1,xh.shape[0]):\n", | |
" xh[i] = xh[i-1] + Δt*dxdt(th[i-1],xh[i-1])\n", | |
" \n", | |
"print(x.shape, xh.shape)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 8, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"[<matplotlib.lines.Line2D at 0x7fca96e900b8>,\n", | |
" <matplotlib.lines.Line2D at 0x7fca96e90208>]" | |
] | |
}, | |
"execution_count": 8, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD5CAYAAAAp8/5SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xUVf7/8dcnCYHQew2hCSqotFBUXMSCyKrYBRTsqKu76jbXXV1d6351LftTVwVRQAV7W8WCiouFFhClC1ITSoBAQkifnN8fd4CACSkzk5lM3s/HYx4zc++5956M4Z3rmXM/15xziIhI9IoJdwdERCS0FPQiIlFOQS8iEuUU9CIiUU5BLyIS5RT0IiJRLq68BmbWEZgGtAEcMNE5928zaw68DnQGNgCXOud2l7L9lcBd/rcPOOemlnfMli1bus6dO1fwRxARkUWLFu10zrUqbZ2VN4/ezNoB7Zxzi82sEbAIOB+4Cshwzv3TzP4CNHPO3XHYts2BFCAZ74/EIqB/aX8QSkpOTnYpKSkV+uFERATMbJFzLrm0deUO3TjntjrnFvtf7wVWAh2AUcD+s/OpeOF/uLOAWc65DH+4zwJGVP5HEBGRqqrUGL2ZdQb6AvOBNs65rf5V2/CGdg7XAdhc4n2qf5mIiFSTCge9mTUE3gZuc85llVznvPGfgGopmNkEM0sxs5QdO3YEsisRESmhQkFvZnXwQv5V59w7/sXb/eP3+8fx00vZNA3oWOJ9on/ZLzjnJjrnkp1zya1alfp9goiIVEG5QW9mBkwGVjrnHi+x6gPgSv/rK4H3S9n8U2C4mTUzs2bAcP8yERGpJhU5oz8ZGAecZmZL/I+RwD+BM81sDXCG/z1mlmxmLwA45zKA+4GF/sd9/mUiIlJNyp1eGQ6aXikiUjkBTa8UEZFqsPYLmP88+AqDvmsFvYhIuBX74LO7vKAPgXJLIIiISIj9+Aakr4CLX4LYOkHfvc7oRUTCqSgfZj8E7fpAz9IKDAROZ/QiIuGU8iJkboLz/g0xoTn31hm9iEi45GXBnEehy1DodlrIDqOgFxEJl7nPQM4uOOOekB5GQS8iEg7ZO2Du09BzFHToH9JDKehFRMLh639BYS6cdnfID6WgFxGpbrs3wMLJ0PcKaNk95IdT0IuIVLfZD0FMLJz6l2o5nIJeRKQ6bVvmXSA16AZo3L5aDqmgFxGpTl/cB/Uaw5Dbq+2QCnoRkeqy8TtY86kX8gnNqu2wCnoRkergHHx+LzRqBwNvqNZDK+hFRKrDT5/A5vkw9A6Ir1+th1bQi4iEWrEPPv8HNO/mTamsZipqJiISaj++DjtWwiVTQlKGuDw6oxcRCaVqKENcHp3Ri4iE0sLJkLkZznsKzMLSBZ3Ri4iESl6WV9Om66nQbVjYulHuGb2ZvQicA6Q7547zL3sdONrfpCmwxznXp5RtNwB7AR9QVNYdykVEotLcp70yxKeHtgxxeSoydDMFeBqYtn+Bc+6y/a/N7DEg8wjbD3PO7axqB0VEaqTsdPjuaW9cvkO/sHal3KB3zs0xs86lrTMzAy4FQndrFBGRmmjOv6Aor1rKEJcn0DH6U4Dtzrk1Zax3wGdmtsjMJgR4LBGRmiFjvXcv2H7joOVR4e5NwLNuxgAzjrB+iHMuzcxaA7PMbJVzbk5pDf1/CCYAJCUlBdgtEZEw+uphrwzx0OopQ1yeKp/Rm1kccCHwelltnHNp/ud04F1g4BHaTnTOJTvnklu1alXVbomIhNeBMsQ3QuN24e4NENjQzRnAKudcamkrzayBmTXa/xoYDiwL4HgiIpHvi3/4yxDfFu6eHFBu0JvZDGAucLSZpZrZtf5Vozls2MbM2pvZTP/bNsA3ZvYDsAD4yDn3SfC6LiISYTZ8C2s+q/YyxOWpyKybMWUsv6qUZVuAkf7X64DeAfZPRKRmCGMZ4vLoylgRkWBY/TGkLvDuA1vNZYjLo6AXEQlUsc+7RWCLo6BP9ZchLo+KmomIBOqH1/xliKdCbOTFqs7oRUQCUZjnzZtv3xd6jgp3b0oVeX96RERqkhR/GeJRT4etDHF5dEYvIlJVeVleTZuuw7xSxBFKQS8iUlXfPQW5GXBGeMsQl0dBLyJSFdnpMPcZ6HWBNz4fwRT0IiJVMedRrwzxsLvC3ZNyKehFRCorYz2kvAT9xkdEGeLyKOhFRCpr9kMQEwdD7wh3TypEQS8iUhnblsLSN2Fw5JQhLo+CXkSkMr64D+o1gZMjpwxxeRT0IiIVdUgZ4qbh7k2FKehFRCrCOfj8Hq8M8aDIKkNcHgW9iEhFrJ4JqQu9MsR1EsLdm0pR0IuIlOdAGeLuIStDvC+/iK2ZuSHZt4JeRKQ8P7wGO1bB6XeHpAyxc44/vvkD5z/zLfvyi4K+fwW9iMiRFOZ58+bb94NjzwvJIZ7+ci0fL9vG9ad0pUHd4P8hUZliEZEjWfgCZKXC+f8JSRniz5Zv47FZP3FB3w5cO6RL0PcPOqMXESlbXiZ8/Zi/DPHQoO9+zfa93P76Ek5IbMLDFx6PhaiefblBb2Yvmlm6mS0rsexeM0szsyX+x8gyth1hZqvNbK2Z/SWYHRcRCbkDZYjvDfquM3MKuX5aCgnxcTw/rj/16sQG/Rj7VeSMfgowopTlTzjn+vgfMw9faWaxwDPA2UBPYIyZ9QyksyIi1eZAGeILoX2foO66yFfMLTMWk7Ynl+fH9aNdk9BO1yw36J1zc4CMKux7ILDWObfOOVcAvAZE5g0VRUQO979HwFcApwW/DPEjn67m6zU7uW/UcfTv1Dzo+z9cIGP0t5jZj/6hnWalrO8AbC7xPtW/TEQksmWsg0X+MsQtugV11+99n8bEOesYf2InxgxMCuq+y1LVoH8W6Ab0AbYCjwXaETObYGYpZpayY8eOQHcnIlJ1sx+CmDrwqz8HdbdLUzO54+0fGdSlOXefU30j2VUKeufcdueczzlXDEzCG6Y5XBrQscT7RP+ysvY50TmX7JxLbtWqVVW6JSISuK0/+ssQ3xTUMsQ79uYz4eUUWjasy38u70ed2Oqb9FilI5lZyZ/+AmBZKc0WAt3NrIuZxQOjgQ+qcjwRkWrzxX1QrymcfGvQdllQVMxNryxid04BE8f3p0XDukHbd0WUe8GUmc0ATgVamlkqcA9wqpn1ARywAbjB37Y98IJzbqRzrsjMbgE+BWKBF51zy0PyU4iIBMOGb2DtLDjzvqCWIb7ng+WkbNzNU2P60qt9k6Dtt6LKDXrn3JhSFk8uo+0WYGSJ9zOBX0y9FBGJOM7BrHugUXsYOCFou31l3kZmLNjEb07txrm92wdtv5WhK2NFRABWfQRpKUEtQzx/3S7u/WA5w45uxR+GHx2UfVaFgl5ExFdUogzx5UHZZdqeXH7z6mKSWtTn32P6EhsTmvIGFaGiZiIiP74GO1fDpdOCUoY4t8DHhGkpFBQVM2l8Mo3r1QlCJ6tOQS8itVthHsx+GDr0D0oZYuccf377R1ZszeLFKwfQrVXDIHQyMBq6EZHayTn4eTZMGemVIT79nqCUIX5+zjr++8MW/nTW0Qw7pnUQOho4ndGLSO2zaR58cT9s/AYaJ8L5zwWlDPHsVen83yerOOeEdtw0NLilEwKhoBeR2mPLEvjyAW+ufIPWcPYj0P8qiAv8AqZ1O7L53Wvfc2zbxjxy8Qkhqy1fFQp6EYl+6atg9oOw8gPvqtcz7vXmysc3CMrus/K82vJ1YmOYOL4/9eMjK1ojqzciIsGUsR6++icsfQPq1Iehd8CJN0O94F2dWlzsuP21JWzclcMr1w0isVn9oO07WBT0IhJ9MtNgzqPw/csQE+eF+8m3Q4MWQT/U47N+4otV6dw/qheDuwZ//8GgoBeR6JG9A755wruhtyuG/lfDKX8IahXKkj78cQtPz17LmIEduWJwp5AcIxgU9CJS8+Xu8e7vOu9ZKMqF3mNh6J+hWejCd8WWLP705o/079SMf5x3XER9+Xo4Bb2I1Fz52TD/Ofju/0Fepnd/12F/hZbdQ3rYjH0FXD8thSYJdXj2in7Ex0X2JUkKehGpeQrzIGUyfP045OyEHiNg2N+g3QmhP7SvmN+8uoid2fm8eeOJtG5UL+THDJSCXkRqDl+h9wXr/x6FvVugy1A47W7oOKDauvDAhyuYty6DJy7rzQmJwatZH0oKehGJfMU+7/Z+Xz0MuzdA4kC48Hno8qtq7cbrCzcxde5Grj+lCxf0TazWYwdCQS8ikcs57yKn2Q/BjlXQ9ngY+wZ0Hx6UujSVsWjjbu56bxmndG/JHSOOqdZjB0pBLyKRxzlY+zl8eT9s/QFa9oBLpsCxoyCm+r/43JaZx42vLKJ90wSeGtOXuGq8sXcwKOhFJLJs+MYrOLZ5HjRNgvOfheMvDUqd+KrIK/Rxw8sp5OQX8ep1g2haPz4s/QiEgl5EIkPqIu8Mft1saNQOfv0Y9B0PceELVuccf313KT+kZjJxXH96tGkUtr4EQkEvIuG1bZk3Br/6I6jfAoY/AAOuC9p9WwPx4rcbeGdxGref0YPhvdqGuztVVm7Qm9mLwDlAunPuOP+yR4FzgQLgZ+Bq59yeUrbdAOwFfECRcy45eF0XkbBzDoryvatRC/P8zyVfl/O8ay2s/BDqNvLmwQ++yXsdAb5Zs5MHP1rBiF5t+e1pR4W7OwGpyBn9FOBpYFqJZbOAO51zRWb2f8CdwB1lbD/MObczoF6KSPAU+2DtF7B3qxfK5QZynj+8S2ubB7iq9SM2Huo2hiG3wUm/g/rNg/pjBmLjrn3cPH0x3Vs34rFLexMTxht7B0O5Qe+cm2NmnQ9b9lmJt/OAi4PbLREJuuJib6riVw97UxUPF1sX6tSDuIRfPtdrAo3aQly9stvUqeeVAo6r5w27lPccE1v9n0EFZOcXcf20FMxg0vhkGtSt+SPcwfgJrgFeL2OdAz4zMwc875ybGITjiUhlOAc/fQJfPgjbl3pTFS9+EToOPhi8cfXCMm0x0hQXO/7wxhLWpmcz7ZpBJLWIvNryVRFQ0JvZ34Ai4NUymgxxzqWZWWtglpmtcs7NKWNfE4AJAElJSYF0S0TAf/PrL7wvOtMWQbMucMFEOP7iiD2bDrenvlzLp8u3c/c5PRnSvWW4uxM0VQ56M7sK70va051zpQ7SOefS/M/pZvYuMBAoNej9Z/sTAZKTk6s46CciAKz/2rt13qa50CQJznsKeo+B2Drh7lnE+nT5Np74/Ccu6pfINSd3Dnd3gqpKQW9mI4A/A0OdczlltGkAxDjn9vpfDwfuq3JPRaR8m+bD7Adg/ZyImYteE/y0fS+/f30JvTs25cELIru2fFVUZHrlDOBUoKWZpQL34M2yqYs3HAMwzzl3o5m1B15wzo0E2gDv+tfHAdOdc5+E5KcQqe3SFntDNGtnQYNWcNbDkHx1RMxFj2S+YsfyLZn8dsb31K8bx/NX9Kdenegb1qrIrJsxpSyeXEbbLcBI/+t1QO+AeiciR1byYqOEZnDGP2Dg9RDfINw9i0hFvmJWbM1i3rpdzF+XwYL1GezNL6JuXAzTrx9M2yaRX1u+Kmr+vCGR2mjHam+a5PJ3oW4T72KjQTdCvcbh7llEKRns89ZlsNAf7ABdWzXgnN7tGdy1OSd1a0mrRnXD3NvQUdCL1CS7fob/PQJL3/DmrJ/yRzjpFu9sXijyFbN8i/+Mff0vg/3cPu0Z3LUFg7s0p3Xj6Dx7L42CXqQm2LPJC/gl070rSk+8BU6+DRq0CHfPwqpksM9bt4uFG3aT7Q/2brU42A+noBeJZFlb4OvHYNFU70YbA6+HIb+HRm3C3bOwKPIVs2z/GXspwT7KH+yDujavEfdyrS4KepFIlJ0O3zwBCyeD80HfcfCrP0KTmnP7umAoGezz1u0ipUSwH9W6oYK9ghT0IpEkJwO+/TcsmOgVDOs9Fob+CZp1DnfPqsXhwb5wfQb7CnyAF+zn9/WCfWAXBXtlKOhFIkHuHpj7DMx7FgqyvTIFQ/8CLWt2edzyFPqKWZaWybx1Gf4z9oPB3r11Qy7o18E7Y+/SIqpnxYSagl4knPL3wvzn4LunIC8Teo6CU++E1seGu2eAd4elAl8xuQU+cgt9B57zCn3kFJTyvtBH3v62/mV5JbbLLbEut6CYrLxCCoqKAS/YL+yXeOCMXcEePAp6kXAoyIGFL8C3T0LOLuhxNgz7K7Q7IeSH3r2vgMnfrGdLZu4hgZ1XWDKED4ZycSUrT5lBQp1YEurEUq9OLAnx3uuE+Fia1o+n3f5l8bE0rBtH78SmCvYQU9CLVKfCPFg0Bb55HLK3Q7fTYNhdkNg/5IcuLna8tSiVhz9eSVZeEW0b1yMhPpb68V4g7w/h+vGx1NsfzocF9f7wrh9fynt/m7pxMVFXK6amU9CLVJd9u2DKSO+mH52GwCVTodOJ1XLolVuzuPu9ZaRs3E1yp2Y8cMFxHNNWV9HWFgp6kepQsA+mXwoZ62HM69DjLG+MI8Sy84t4ctZPvPTdBpok1OHRi0/gon6JNf7WeFI5CnqRUPMVwZtXw5bFcOk0OHpEyA/pnGPm0m3c9+Fy0vfmM3pAEneMOJqm9VWuuDZS0IuEknPw4a2w5lP49eNw7LkhP+T6nfv4+/vL+HrNTnq1b8yzV/SnX5Jq4dRmCnqRUJr9IHz/Cgy9AwZcG9JD5RX6+M9XP/PcVz9TNy6Ge8/tyRWDOxEXq3vB1nYKepFQWTAJ5jwK/cZ7c+NDaPbqdO79YDkbd+Uwqk97/jby2FpdxEsOpaAXCYUV78PMP3nz43/9RMi+eN2yJ5f7P1zBx8u20bVVA6ZfN4iTjoqem1pLcCjoRYJtwzfw9nWQOAAufhFig//PrNBXzEvfrufJz9fgK3b86ayjue6ULtSNi77b4EngFPQiwbR9OcwY6xUhG/s6xNcP+iEWrM/grveW8tP2bE4/pjX3nteLjs2DfxyJHgp6kWDZsxleucgL9yvegfrNg7r7ndn5PDxzFW8vTqVD0wQmjU/mzJ61sy69VI6CXiQYcjLglQu9GjbXfAxNOwZt18XFjhkLN/HIJ6vZl1/ETad247enHUX9eP3zlYqp0G+Kmb0InAOkO+eO8y9rDrwOdAY2AJc653aXsu2VwF3+tw8456YG3m2RCFKQA9Mvg90bYdw70KZX0Ha9LC2Tv723jB8272Fw1+bcP+o4urdpFLT9S+1Q0Qm2U4DDL+f7C/CFc6478IX//SH8fwzuAQYBA4F7zExXbkj08BXB29dC6kK4aBJ0HhKU3WbmFnLP+8s47+lvSNudy5OX9WHG9YMV8lIlFTqjd87NMbPOhy0eBZzqfz0V+Aq447A2ZwGznHMZAGY2C+8Pxowq9VYkkjgHH90Oq2fCyH95teQD3qXj/SVbeOCjlezal8+4wZ34w/CjaZJQJwgdltoqkEG+Ns65rf7X24DSvhXqAGwu8T7Vv0yk5vvqYVg8DU75o3fT7gCtTd/L3e8tZ+66XfRObMJLVw3g+MQmQeio1HZB+TbHOefMrJK3JziUmU0AJgAkJSUFo1siobNwMvzv/6DvFXDaXeW3P4LcAh9PfbmGSV+vI6FOLA+cfxxjBiYRqwqTEiSBBP12M2vnnNtqZu2A9FLapHFweAcgEW+I5xeccxOBiQDJyckB/dEQCamV/4WZf4QeI+Ccfwd01eusFdu594PlpO3J5aJ+idw58hhaNtSdliS4Agn6D4ArgX/6n98vpc2nwEMlvoAdDoS26IdIKG2cC29dCx36w8UvVfmq180ZOfzjvyv4fOV2erRpyOsTBjOoa4sgd1bEU9HplTPwzsxbmlkq3kyafwJvmNm1wEbgUn/bZOBG59x1zrkMM7sfWOjf1X37v5gVqXHSV8KMy6BpknfzkCpc9VpQVMykr9fx1JdrMIw7zz6Ga4Z0oY4qTEoImXORN0qSnJzsUlJSwt0NkYMyU2HycCj2wbWfQbNOld7Frux8xk1ewIqtWYzo1Za/n9uT9k0TQtBZqY3MbJFzLrm0dbq0TqQ8ORnw8oWQvxeunlnlkL/8hfms37mP58f156xebUPQUZHSKehFjqQwF2aMgd3r4Yq3oe3xld7Fzux8Lp80nw279jH5ygEM6a4ywlK9FPQiZfEVeV+8bp4Pl7wEXX5V6V3szM5n7KR5bMrI4cWrBnCyasVLGCjoRUrjHMz8A6z+CM5+BHpdUOldHBLyVw7QDUEkbBT0IqX53yOwaAoM+T0MuqHSm+/Y64X85t0KeQk/Bb3I4RZNga8egt5j4fS/V3rzQ0L+qgGc1E0hL+GloBcpadVH8OHtcNSZcN7/q/RVr+l78xg7aT5pu3N56aqBnNhNF0FJ+CnoRfbbNA/eugba9YFLp0Js5SpGpu/NY8zEeWzZk8dLVw9gsK50lQihoBcBSF/l3TykcQe4/E2Ib1C5zbPyGDNJIS+RSUEvkpnm3es1rq53h6gGlRtTT8/KY/SkeWzLzGPK1QNUs0YijoJearfc3fDqxZCX6b/qtXOlNj805AcysEtwbwguEgwKeqm9CvNgxljYuca76rXdCZXafHuWNya/LUshL5FNQS+1U7EP3rkONn0HF78IXYdWavP9Ib89K4+p1wxkQGeFvEQuBb3UPs7BzD95NxAZ8U847qJKbb4t0/viNd0f8skKeYlwKoIttc+cf0HKZDj5Vhh8U6U2VchLTaQzeqldFk2F2Q/ACaPh9Hsrten+kN+xN59p1w6kfyeFvNQMCnqpPVZ/DB/eBt1Oh1FPQ0zF/4d2a2YuYybOY2d2AVOvGUj/Ts3K30gkQijopXbYvADevBra9YZLp1XqqtetmbmMnjiPXQp5qaEU9BL90hbDKxdD43Yw9k2o27DCm27Zk8uYSfPIyC5g2rUD6ZekkJeaR0Ev0W3LEnj5fEhoAuM/gIatKr7pHu9Mfvc+L+T7KuSlhtKsG4le25Z6IV+3MVz5ITTtWOFN0xTyEkWqHPRmdrSZLSnxyDKz2w5rc6qZZZZoU/ni3iJVsX05TD0P6jSAK/9bqRt6eyE/l905Bbx83SCFvNR4VR66cc6tBvoAmFkskAa8W0rTr51z51T1OCKVlr7SC/m4enDlB9C8S4U3Td2dw5hJ89iTU8gr1w6id8emIeyoSPUI1hj96cDPzrmNQdqfSNXsWA1Tz4WYOO9MvkW3Cm+aujuH0RPnkZmrkJfoEqwx+tHAjDLWnWhmP5jZx2bWK0jHE/mlnWu8kMe8kG95VIU33ZzhhXxWbiGvXqeQl+gScNCbWTxwHvBmKasXA52cc72Bp4D3jrCfCWaWYmYpO3bsCLRbUtvs+hmmnAOu2Av5Vj0qvOmhIT+YExIV8hJdgnFGfzaw2Dm3/fAVzrks51y2//VMoI6ZlXpXB+fcROdcsnMuuVWrik+BEyFjnRfyxYXeFMrWx1R40/0hn51fxKvXDeb4xCYh7KhIeAQj6MdQxrCNmbU18+6ubGYD/cfbFYRjinh2b4Ap50JRrhfybXpWeNNDQ36QQl6iVkBfxppZA+BM4IYSy24EcM49B1wM3GRmRUAuMNo55wI5psgBezZ5Y/IF2d7smrbHVXjTTbu82TX7Q/64Dgp5iV4BBb1zbh/Q4rBlz5V4/TTwdCDHEClVZqo3XJOXCePf92rYVNCmXTmMnjiXnEKfQl5qBZVAkJona4sX8rm7Yfx70L5vhTfduGsfYybOOxDyvdor5CX6KeilZtm7zQv5fTth3LvQoX+FN924ax+jJ84jVyEvtYyCXmqOvdu9Mfns7XDFO9BxQIU33R/yeYU+pl83mJ7tG4ewoyKRRUEvNUP2Dph2HmSmwRVvQdKgCm/63c87uWX69zjneFUhL7WQqldK5Nu30wv53Rvh8jeg00kV2sw5x+Rv1jNu8gKaN4jn7ZtOUshLraQzeolsORkwbZR3UdTYN6DzkAptllfo4853lvLu92kM79mGxy/rQ8O6+nWX2km/+RK5cjK8M/mda2Dsa9B1aIU2S9uTyw0vp7B8Sxa/P7MHtww7ipgYC3FnRSKXgl4iU+5uePkCrxrl6BnQ7bQKbTb3513cPH0xhUXFvDA+mdOPbRPijopEPgW9RJ68THj5Qu/mIaOnQ/czyt3EOcdL327gwZkr6dKyARPH9adrq4rfG1YkminoJbLkZcErF3m3AbzsZegxvPxNCn389Z2lvOMfj3/s0t40qlenGjorUjMo6CVy5O+FVy+GLd/DJVPg6LPL3UTj8SLlU9BLZCjYB69eCqkpcMlLcOy55W6i8XiRilHQS/gV5MD0y2DzPLhoMvQcdcTmzjmmfLeBBz5aSecW9Zk4PpluGo8XKZOCXsKrIAdmXAYbv4ULJsJxFx6xeV6hj7++u5R3FqdxZs82PK7xeJFyKeglfArz4LWxsP5ruOA5OOGSIzZP25PLjS8vYmlaJref0YPfnqbxeJGKUNBLeBTlw+uXw7qvYNQz0Hv0EZvP/XkXt0xfTEFRMZOv1Hi8SGUo6KX6FeXD6+Ng7edw3lPQ9/Iym2o8XiRwCnqpXkUF8OZVsOZTOOdJ6De+zKYajxcJDgW9VB9fIbx1NayeCSP/BclXl9l0y55cbtB4vEhQKOileviK4O1rYdWHcPYjMPD6MpvOW7eLm1/1xuNfGJ/MGT01Hi8SCAW9hJ6vCN65Hla8D2c9BINuKLWZxuNFQiPgoDezDcBewAcUOeeSD1tvwL+BkUAOcJVzbnGgx5UaotgH790Iy9+BM++HE28utVnJ8fgzjm3DE5dpPF4kWIJ1Rj/MObezjHVnA939j0HAs/5niXapi+DTO2HzfDj9Hjj5d6U203i8SGhVx9DNKGCac84B88ysqZm1c85trYZjSzhkpsLn/4Clb0CD1nD+s9BnbKlN94/H5xcVM2l8MmdqPF4k6IIR9A74zMwc8LxzbuJh6zsAm0u8T/UvU9BHm/xs+PZJ+O4pcA5O+QMMuR3qNvpFU+ccU/3j8Ukt6jNxXDJHtdZ4vEgoBCPohzjn0sysNTDLzFY55+ZUdidmNgGYAJCUlBSEbkm1KfbBkunw5f2QvR2OvwRO/zs0Lf2/Y16hj7+9u4y3F6dyxrFtePyy3jTWeLxIyAQc9M65NInX+vgAAAmcSURBVP9zupm9CwwESgZ9GtCxxPtE/7LD9zMRmAiQnJzsAu2XVJP1c+DTv3o3CkkcCJe9Ch0HlNl8y55cbnxlET+mZnLbGd353WndNR4vEmIBBb2ZNQBinHN7/a+HA/cd1uwD4BYzew3vS9hMjc9HgZ1rYdbfYfVH0CQJLn4Rel0IVnZoazxeJDwCPaNvA7zrzaAkDpjunPvEzG4EcM49B8zEm1q5Fm96ZdmXQ0rky8mAOY/CgokQl+DNphn8G6hTr8xNnHNMm7uR+z9cofF4kTAIKOidc+uA3qUsf67EaweUPnlaag5fISycDP/7p3fz7n7jYdjfoGHrI26WW+Dj7veX8daiVM44tjWPX9ZH4/Ei1UxXxsqROQc/fQKf3QW71kLXU72rW9v0OuJmmzNyeGX+Rl5fuJk9OYXcenp3bj1d4/Ei4aCgl7JtWwqf/g3W/w9a9oCxb0D34WWOwxcXO75Zu5Npczfwxap0YswY3rMN1wzpwoDOzau37yJygIJefmnvdpj9ACx+GRKawtmPepUmY0sfcsnKK+StlFRembeRdTv30bJhPDefehRjByXRvmlCNXdeRA6noJeDCnNh7jPwzRPezUFOvBl+9UdIaFZq89Xb9jJt7gbe/T6NnAIffZOa8uRlfTj7+LbUjYut3r6LSJkU9OKNwy97Gz6/FzI3wzHnwJn3QYtuv2ha6Ctm1ortTP1uA/PXZxAfF8Oo3u0Zf2Jnjk9sUv19F5FyKehru80LvAueUhdC2xO8ujRdTvlFs/S9eby2YDOvzt/I9qx8Epsl8Jezj+Gy5I40axAfho6LSEUp6GurPZtg1j1e+eCGbWHUf6D3GIiJOdDEOcfiTbuZNncjM5dupdDnOKV7Sx48/3iGHdOaWM2gEakRFPS1TV6WNwY/9xmwGBh6B5x8K8Q3ONik0McHS7Ywde4Glm/JolHdOK4Y3IlxgzvRVTcCEalxFPS1RbEPvn8ZvnwA9u2AE0Z7hceadDjQZNMub+77Gyne3Pej2zTiwQuO4/w+HWhQV78qIjWV/vXWBj/P9ubDpy+HpBO9+fAd+gHe3Pc5a3bw8tyNfLnam/s+oldbxp3YiUFdmmNHqF0jIjWDgj6a7fgJZt3tXdnarDNcOg2OPQ/MyMwt5K1Fqbw8dwMbduXQsmFdfjvsKMYO6kTbJmXXrRGRmkdBHy3ys2H7Mtj6w8FH+gqIb+jdq3XQDRBXl5Vbs5g2dyPvfZ9GbqGP/p2acfuZPTj7uHbEx8WUfxwRqXEU9DVR7m7Y+uOhob5rLd7NvoAGraBdH28+/MDrKazXnE+Xb2Pa3I0sWJ9B3bgYzu/TgXEnduK4Dpr7LhLtFPSRLnuHP8yXHAz1PRsPrm+cCO16e3d1atfbezRqC2ak781jxtzNTF+whO1Z+XRsnsBfRx7DpckdaVpfc99FagsFfaRwDrK2HHqWvvUH2LvlYJvmXb0vUftfBe16U9DqeLYWNSBtdy6pe3JJ3ZRL2g/ppO3ZQOruXLbsyaXYwdAerXj4wk4M7aG57yK1kYI+HJyD3ev9YV5iCCZnp7feYqBlD3ydhpDR5FjS6vVgjXVhfXYsaXtySV2WS9rXuWzfuxBX4qaLZtC2cT06NE2gf6dmXNQvkfP7dqBLywal90NEagUFfagV+7zx80PO1H+E/EwAXEwc2Y27s63ZENa1Ooqlvs4syG3Hut2OnZsLSuxoA3ExRrumXpAP6d6SDk0T6NAsgcRmCSQ2rU/bJvX0haqI/ELtDfriYiguhOIi/8PnPftKWVZc5G9b8r1/fWntC3Nw6SvxpS0hJn0ZMUW5ABRaPKnxXVnhTiKFJBbmd+Qn15GCHK/8b3xcDIn+8D6zQ0KJIK9Ph6YJtGlcT0MvIlJpURX0a+/vR12XR6zzEUuJxyHvi4nFRwyu/B0GIMfVY5nrzPLioSwr7swy14VtdTrStl7DA8F9jv9sfH+gt2xQV3dgEpGgi6qg39OgCxT78Jk/zi2W4gOv4yg2f9QfeI7DZzEUl1i2v92B7UvZdv+6/esP2a/FUWx1iG/WnsTmDenQNIGB/kBvklBHV5qKSLWLqqBP/v3b4e6CiEjEqfI3d2bW0cxmm9kKM1tuZreW0uZUM8s0syX+x98D666IiFRWIGf0RcAfnHOLzawRsMjMZjnnVhzW7mvn3DkBHEdERAJQ5TN659xW59xi/+u9wEqgw5G3EhGR6haUSddm1hnoC8wvZfWJZvaDmX1sZr2CcTwREam4gL+MNbOGwNvAbc65rMNWLwY6OeeyzWwk8B7QvYz9TAAmACQlJQXaLRER8QvojN7M6uCF/KvOuXcOX++cy3LOZftfzwTqmFnL0vblnJvonEt2ziW3atUqkG6JiEgJgcy6MWAysNI593gZbdr622FmA/3H21XVY4qISOUFMnRzMjAOWGpmS/zL/gokATjnngMuBm4ysyIgFxjtnAvtJakiInIIi8TcNbMdwMZyG5auJbAziN2pyfRZHEqfx6H0eRwUDZ9FJ+dcqePeERn0gTCzFOdccrj7EQn0WRxKn8eh9HkcFO2fhWraiohEOQW9iEiUi8agnxjuDkQQfRaH0udxKH0eB0X1ZxF1Y/QiInKoaDyjFxGREqIm6M1shJmtNrO1ZvaXcPcnnMzsRTNLN7Nl4e5LJKhISe3awszqmdkCf/2p5Wb2j3D3KRKYWayZfW9mH4a7L6EQFUFvZrHAM8DZQE9gjJn1DG+vwmoKMCLcnYgg+0tq9wQGAzfX4t+PfOA051xvoA8wwswGh7lPkeBWvAq8USkqgh4YCKx1zq1zzhUArwGjwtynsHHOzQEywt2PSKGS2gc5T7b/bR3/o1Z/UWdmicCvgRfC3ZdQiZag7wBsLvE+lVr6D1mOrJyS2rWCf5hiCZAOzHLO1drPwu9J4M9Acbg7EirREvQi5SqnpHat4ZzzOef6AInAQDM7Ltx9ChczOwdId84tCndfQilagj4N6FjifaJ/mQhQfknt2sg5tweYTe3+Pudk4Dwz24A35Huamb0S3i4FX7QE/UKgu5l1MbN4YDTwQZj7JBGiIiW1awsza2VmTf2vE4AzgVXh7VX4OOfudM4lOuc64+XGl865K8LcraCLiqB3zhUBtwCf4n3R9oZzbnl4exU+ZjYDmAscbWapZnZtuPsUZvtLap9mZkv8j5Hh7lSYtANmm9mPeCdIs5xzUTmlUA7SlbEiIlEuKs7oRUSkbAp6EZEop6AXEYlyCnoRkSinoBcRiXIKehGRKKegFxGJcgp6EZEo9/8BYE09G4vZpQIAAAAASUVORK5CYII=\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"\n", | |
"\n", | |
"plt.plot(t,np.vstack((x,xh[::2])).T)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"### Runge Kutta 4: \n", | |
"\n", | |
"\n", | |
"\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 10, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"(10,)\n" | |
] | |
} | |
], | |
"source": [ | |
"Δt = 0.5\n", | |
"tr = np.arange(0,5,Δt)\n", | |
"xr = np.zeros(tr.shape)\n", | |
"print(xr.shape)\n", | |
"\n", | |
"#ci\n", | |
"xr[0] = 3\n", | |
"\n", | |
"for i in range(1, xr.shape[0]): \n", | |
" \"Apply Runge Kutta Formulas to find next value of y\"\n", | |
" k1 = Δt*dxdt(t[i-1],xr[i-1])\n", | |
" k2 = Δt*dxdt(t[i-1] + 0.5 * Δt, xr[i-1] + 0.5 * k1) \n", | |
" k3 = Δt*dxdt(t[i-1] + 0.5 * Δt, xr[i-1] + 0.5 * k2) \n", | |
" k4 = Δt*dxdt(t[i-1] + Δt, xr[i-1] + k3) \n", | |
"\n", | |
" # Update next value of x\n", | |
" xr[i] = xr[i-1] + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4) " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 12, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"[<matplotlib.lines.Line2D at 0x7fca61c5f748>,\n", | |
" <matplotlib.lines.Line2D at 0x7fca61c5f898>,\n", | |
" <matplotlib.lines.Line2D at 0x7fca61c5f9e8>]" | |
] | |
}, | |
"execution_count": 12, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3SUxf7H8fekdwiEFEIJnYQEQlcBFSkCerGg/sCGoKKC137Fa1fs5WJXWgAromDDQhOlSm8JCZBAIAnpCenJtvn9sRFBKSm72Wz4vs7Zs7vP7jMzWZIPc2bnmVFaa4QQQjgfF0c3QAghRN1IgAshhJOSABdCCCclAS6EEE5KAlwIIZyUW0NWFhQUpCMiIhqySiGEcHrbt2/P01q3+vvxBg3wiIgItm3b1pBVCiGE01NKHTndcRlCEUIIJyUBLoQQTkoCXAghnJQEuBBCOCkJcCGEcFIS4EII4aQkwIUQwklJgAshhB1VmCp4ZcsrpJek27xsCXAhhLCjb5O/5bPEz8guz7Z52RLgQghhJyaLiYUJC+nVqhd9gvvYvHwJcCGEsJOVR1aSUZrBpOhJKKVsXv45A1wp1VYptUYptU8plaCUur/6+LNKqQyl1K7q2xibt04IIZyU1pq4+Dg6NOvA0LZD7VJHTRazMgEPa613KKX8ge1KqZXVr83UWr9hl5YJIYQT23RsE0kFSTx/0fO4KPsMdpwzwLXWmUBm9eMSpVQiEG6X1gghRBMRlxBHsHcwV3S8wm511Oq/BaVUBNAb2Fx96F6l1B6lVJxSKvAM50xRSm1TSm3Lzc2tV2OFEMIZJOQnsDlzMzdH3YyHq4fd6qlxgCul/IAlwANa62LgQ6ATEIu1h/7m6c7TWs/WWvfTWvdr1eof65ELIUSTE7c3Dn93f67ver1d66lRgCul3LGG92da66UAWutsrbVZa20B5gAD7NdMIYRwDkeLj7Lq6Cpu6HYDfh5+dq2rJrNQFDAPSNRa/++k42Enve0aIN72zRNCCOeyMGEhrsqVmyJvsntdNZmFMgi4BdirlNpVfexxYIJSKhbQQCpwl11aKIQQTiKvIo9vk79lbKextPKx/5BxTWahrAdONwP9J9s3RwghnNfniZ9jtBi5rcdtDVKfXIkphBA2UGYsY9H+RQxrN4yIZhENUqcEuBBC2MDXB76mxFDC5OjJDVanBLgQQtST0Wzk430f0z+0PzGtYhqsXglwIYSop58O/0ROeQ6Tekxq0HolwIUQoh4s2sL8+Pl0DezK4PDBDVq3BLgQQtTD2vS1pBSl2G3J2LORABdCiHqYHz+f1r6tuTzi8gavWwJcCCHqaFfOLnbk7ODWHrfi7uLe4PVLgAshRB3Ni59Hc8/mXNP5GofULwEuhBB1kHI8hd/SfmNC9wn4uPs4pA0S4EIIUQcLEhbg5erFhO4THNYGCXAhhKilrLIslh1axjVdriHQ67R72TQICXAhhKilT/d9itaaiT0mOrQdEuBCCFELxYZivjrwFZdHXE64n2O3B5YAF0KIWli8fzHlpnImRTfsZfOnIwEuhBA1VGWu4tN9nzKo9SC6t+ju6OZIgAshRE19n/I9+ZX5Dbpk7NlIgAshRA2YLWYWxC+gR8se9A/t7+jmABLgQghRI6uPruZoyVEmR09u8EWrzkQCXAghzkFrTVx8HO0D2jOs3TBHN+cECXAhhDiHrVlbSchPYGKPibi6uDq6OSdIgAshxDnExcfR0qslYzuNdXRTTiEBLoQQZ5FUkMSGYxu4OepmPF09Hd2cU0iACyHEWcyPn4+Pmw83dLvB0U35BwlwIYQ4g4zSDJanLuf6rtcT4BHg6Ob8gwS4EEKcwcKEhSiluCXqFkc35bQkwIUQ4jQKKwv55uA3XNnxSkJ8QxzdnNOSABdCiNP4IukLKs2VTOrh+EWrzkQCXAgh/qbcWM7nSZ9zadtL6di8o6Obc0YS4EII8TffJH9DUVURt0ff7uimnJUEuBBCnMRoMbIwYSF9gvsQGxzr6OaclQS4EEKcZHnqcjLLMhvFhg3nIgEuhBDVtNbMj59Pp2aduLjNxY5uzjlJgAshRLUNxzZwoPAAk6In4aIafzyes4VKqbZKqTVKqX1KqQSl1P3Vx1sopVYqpQ5W3wfav7lCCGE/cfFxhPiEMKbDGEc3pUZq8l+MCXhYax0FXABMU0pFAY8Bq7XWXYDV1c+FEMIp7c3dy9asrdwSdQvuru6Obk6NnDPAtdaZWusd1Y9LgEQgHLgKWFj9toXA1fZqpBBC2FtcfBz+Hv5c1/U6Rzelxmo1yKOUigB6A5uBEK11ZvVLWUDjvNZUCCHOIbUoldVHVzO+23h83X0d3Zwaq3GAK6X8gCXAA1rr4pNf01prQJ/hvClKqW1KqW25ubn1aqwQQtjDgoQFuLu4c2PkjY5uSq3UKMCVUu5Yw/szrfXS6sPZSqmw6tfDgJzTnau1nq217qe17teqVStbtFkIIWwmtzyX71O+5+rOVxPkHeTo5tRKTWahKGAekKi1/t9JL30PTKx+PBH4zvbNE0II+/os8TPM2szEHhPP/eZGpiY98EHALcBlSqld1bcxwCvACKXUQWB49XMhhHAapYZSFu9fzPB2w2kX0M4+lVjMsOl9MFbavGi3c71Ba70eUGd4eZhtmyOEEA3nqwNfUWIsYXLMZPtVsvZ1+O1laNYGoq6yadGN/1IjIYSwA4PZwKf7PmVg2EB6tOxhn0pS18Pvr0LP/4NI2+9oLwEuhDgv/XjoR3Iqcpgcbafed1k+LLkTAiPgijdBnWkgo+7OOYQihBBNjUVbiIuPI7JFJBeGXWj7CrSG76ZCeR7cvhI8/W1fB9IDF0Kch9akrSG1OJVJ0ZNQdugZs/kjOPALjJgBre23prgEuBDivKK1Ji4+jnC/cEa0H2H7Co7thBVPQdfRMPAu25d/EglwIcR5ZUfODvbk7mFij4m4udh4FLmqBL6eDL6t4OoP7DLufTIZAxdCnFfi4uMI9Azk6s42Xn9Pa1j2EBSmwsRl4NPCtuWfhvTAhRDnjYOFB1mbvpYbI2/E283btoXv/gL2LoZLHoOIQbYt+wwkwIUQ540FCQvwdvNmfLfxti047yD8+DC0HwwXP2Lbss9CAlwIcV7ILM3kp0M/Ma7LOJp7NbddwcZK+Oo2cPOCcXPAxdV2ZZ+DjIELIc4LH+/7GIBbo261bcErnoTseLhxMQS0tm3Z5yA9cCFEk1dUVcSSg0sY3WE0YX5htis48QfYOgcuvBe6Xm67cmtIAlwI0eQtSlpEhamC26Jvs12hx4/Cd9MgLBaGPWO7cmtBAlwI0aRVmir5POlzhoQPoWtgV9sUajbBkjvAYoHr4sDNwzbl1pKMgQshmrTvkr+joLLAtotW/fYypG2GcfOgZSfblVtL0gMXQjRZJouJBQkL6BnUk74hfW1T6KHfYN2b0PtmiHHsDvYS4EKIJmvVkVWkl6YzOXqybRatKs2FpVMgqAuMfq3+5dWTDKEIIZqkcmM5c/fOJSIggqHthta/QIsFvr0bKo7DzUvBw7f+ZdaT9MCFEE3OzpydjPt+HAcKDzCt9zRclA2ibtN7kLwKRr0EodH1L88GpAcuhGgyDGYD7+16jwXxC2jt15r5o+bbZuw7fTusfs66LVq/2+tfno1IgAshmoSkgiT+u+6/JB9P5rqu1/FIv0fwdbfBMEdlEXw9Cfxbw9h37L5EbG1IgAshnJrJYmJ+/Hw+2P0BgZ6BvD/sfS5uc7FtCtcafrgfitJh8i/gHWibcm1EAlwI4bRSi1J5YsMT7Mndw6iIUTwx8AnbLlS1YyEkfGO90rLtANuVayMS4EIIp2PRFhYlLWLm9pl4uHrw2sWvMbrDaNtWkpMIP0+HjpfCoAdsW7aNSIALIZxKVlkWT214ij8y/2BQ+CCev+h5gn2CbVuJoRy+mmTdTf6a2eDSOCfsSYALIZyC1pplh5bx8uaXMWkTT13wFNd3vd4+u8ov/y/kJlrne/uH2L58G5EAF0I0egWVBczYNINVR1fRJ7gPLwx6gbYBbe1TWfxS2L7AOmzSeZh96rARCXAhRKP269FfeW7Tc5QYSnio70PcGnUrrvba9aYw1TrrpE1/uOxJ+9RhQxLgQohGqcRQwqtbXuW7lO/o3qI7c0bOsd1ysKdjNsLXkwFlXWXQ1d1+ddmIBLgQotHZkrmFJzc8SXZ5NlN6TuHunnfjbu9A/XUGZGyH6xdCYHv71mUjEuBCiEaj0lTJ2zve5tPET2kf0J6PR39Mr1a97F9x8irY8Db0mww9rrZ/fTYiAS6EaBTi8+J5fP3jHC46zITuE3iw74N4u3nbv+KSLFh6FwRHweUv2b8+G2qckxuFEOcNo8XIezvf4+afbqbcWM7sEbN5fODjDRPeFot1fW9DGVw3H9xtX+ex4xXcMm8z6YXlNi9beuBCCIdJLkzm8fWPk1iQyNhOY5k+YDoBHgEN14ANM+Hw7zD2XQjubvPi92eVMDFuC2VVJrKKKmkT6GPT8iXAhRANzmwx82nip7yz4x183X1569K3GNa+gedcH/0Dfn0RosdB71tsXvzmQ/nc8fE2vN1dWXz3hUSG2f4/pnMOoSil4pRSOUqp+JOOPauUylBK7aq+jbF5y4QQTVJ6STq3r7idN7a9waDwQSy9amnDh3d5gXVX+eZt4cq3bL5E7M97M7klbgvB/p4snXqRXcIbatYDXwC8B3z8t+MztdZv2LxFQogmSWvN0oNLeW3ra7goF14Y9AJjO421z6XwZ28IfP9v65eXt68AL9uG68KNqTz7QwJ92gUy99Z+BPp62LT8k50zwLXWa5VSEXZrgRCiycstz+XZTc+yNn0tA0MHMmPQDML8whzTmK1zIWkZjHwRwvvYrFitNa8v388Hv6UwIiqEdyf0xsvdTleMVqvPGPi9SqlbgW3Aw1rrwtO9SSk1BZgC0K5du3pUJ4RwRstTlzPjjxlUmip5bMBjTOg+wTZ7VNZF1l5Y/gR0GQkXTLVZsUazhceW7GXJjnRuHNiO58f2wM3V/j9jXWv4EOgExAKZwJtneqPWerbWup/Wul+rVq3qWJ0QwtkUVRXx6NpHeeT3R2jr15bF/1rMTZE3OS68DWXWJWK9A+HqD222RGxZlYk7Fm5jyY50HhrRlRevjm6Q8IY69sC11tl/PlZKzQGW2axFQginYbKYKDOWUWIoodRYSomhhBJDCXkVeczaPYuCygKmxU7jjpg7cHNx8KS3nx6F/GSY+D34BtmkyLzSKiYv2ErCsWJeuTaG8QMadpShTp+oUipMa51Z/fQaIP5s7xdCND5mi5lSY6n1Zig9JYRP3BtKKTH+7f6k4xWmijOW36lZJ94d9i5RLaMa8Kc6gz2LYdencPGj0ME2+2UeyS/j1rgtZBdXMvuWvgyLbPh1w88Z4EqpL4BLgSClVDrwDHCpUioW0EAqcJcd2yiEqAGj2ciu3F2kHE85NYANJSfC9+SALjOWnbNMDxcP/D388ffwx8/dDz8PP4J9gk889vfwx9/d3/q4+v7Px639Wju+1w2QnwLLHoR2F8El021S5J7040yavxWL1nx+5wX0aeeYzY5rMgtlwmkOz7NDW4QQtZRVlsX6jPWsS1/H5qzNp4Syu4v7KcHr7+5PkHeQNWTd/U685u/h/48w/vO4h6v9psA1CFMVfD3JujTsuDngWv//UH4/kMs9n24n0MeDj28fQKdWfjZoaN00gv8ehRA1ZTQb2Zmz0xraGetIPp4MQKhvKKM7jGZI+BBigmII8AzA09XTwa11MJMBlj8Ombth/BfQrE29i1yyPZ3pS/bQJcSfhZP6ExzgZYOG1p0EuBCN3J+97PUZ6/kj8w/KjGW4ubjRN7gvV/W9isHhg+nUvFPDXxDTGGkNWXtg52ew9yuoKICB90D3+l0srrXmw99TeO2X/Qzq3JKPbu6Lv5fjN3yQABeikTFajOzK2cW6jHWsSz99L3tg2EB83X0d3NJGpCyv+ovKzyA7Hlw9ofsVEHsTdLqsXkWbLZrnf0hg4aYjjO3Vmjeu74WHW+NYyFUCXIhGQHrZdWA2wsEVsOtzOPALWEzQug9c8aZ1gSrv+n+xWGk089DiXfy0N4s7h3Tgv6MjcXFpPP8GEuBCOMDJvez1Ges5WHgQ+KuXPTh8MBeEXSC97NPJTrCG9p4voSwXfIPhgnusve3gSJtVU1Rh5M6Pt7HlcAFPXhHJHUM62qxsW5EAF6KBZJVlsSFjA+sz1rMpc9OJXnaf4D481PchhoQPkV72mZQXwN6vrUMkmbvAxR26jbaGdufhNpldcrLMogpui9vKobxS3h4fy1Wx4TYt31YkwIWwE+ll15PZBCm/Wi/A2f8zmA0Q2hNGvQox14NvS7tUezDbuglDcaWJBZMGMKizba7atAcJcCFsKLss+5Sx7FJj6Sm97MHhg+ncvLP0ss8md7+1p737SyjNAp+W0P8OiL0RQmPsWvXW1ALuWLgNDzcXvrzrAnq0bmbX+upLAlwIG/g97Xfe2fkOBwoPABDiE8KoDqOkl11TFcchfol1bDtjGyhX6Hq5dYiky0hws/8FRb/EZ3H/op2EB3qzcNIA2raw7fZn9iABLkQ9Ld6/mBc3v0jHZh2ll10bFjMc+s3a205cBuYq687wI1+EnjeAX3CDNeWTP47wzHfx9GrbnHkT+9PCjpsw2JIEuBB1pLXmvV3vMXvPbC5pcwmvX/J6w+yk7uzykmH357B7ERRngFdz6HMr9L4JwmJtvr3Z2WiteXPFAd5bk8yw7sG8d2MfvD3suwmDLUmAC1EHRouR5zY+x3cp3zGuyzievODJxrFwU2NVWQwJ31iHSNL+AOVinT1y+YvQbQy4Nfxl/0azhSe+2cvibemM79+WFxpwHW9bkd84IWqp3FjOQ78/xIaMDUyNncrdPe+W4ZLTsVggdZ11iGTf92CqgKCuMPxZ6DkeAhy0pRpQbjAx7bMdrNmfy33DuvDg8C5O+W8oAS5ELeRV5DFt9TT2F+znuYue49ou1zq6SXWntfVqRnOV9d5U9bfHhr9uJkP1a38+rn5uMvztfdXvMVZAyhooOgqezaDXeOsXkm36NegQyenkl1YxeeE29qYf58VrorlpYHuHtqc+JMCFqKEjxUe4e+Xd5Ffm885l73BxG9tsDGATWltX3du9yLoWyFmD2fhXGNuUsg6FuHpal28NjYHhz1jXJHFvHN8NHM0vZ+L8LRw7XsFHN/dlZI9QRzepXiTAhaiBvbl7mbZ6GgDzRs4jppV95yPXWEk27F0Mu76AnARw9bB+EejuDZ7+1YHqbg1VNw/r638GrJtn9XOPU9/n6lH93jO9z+Ok5x5/nWPjqyFtLT6jiNvmb8VksfD5nQPp276Fo5tUb437ExeiEVibvpZHfn+Ell4tmTViFu0CGnbfw38wVsKBn62hnbwKtBnC+1kXcepxLfg4fzDZ2rqDudz9yXaa+3iwaPJAOgf7O7pJNiEBLsRZLDmwhBl/zKB7i+68N+w9grwddFm11pCxw/qFYPwSqDwO/q1h0H3Q60Zo1dUx7XIC3+xM5z9f7aFzsB8LJw8gxMGbMNiSBLgQp6G15qPdH/HB7g8YHD6YNy95Ex93B1yZV3zMOq69+wvIOwBuXhD5L+g1ATpeCi7OM2e5IRlMFlJyS/k5Pot3Vh/kgo4tmH1rPwIawSYMtiQBLsTfmCwmXvjjBZYcXMLVna/m6Qufxt2lAf/wDeWQ9KP1YpeUNYCGdhfCv96BHleDV+Nen6Oh5ZRUkpRZQmJmMUlZ1vvknFJMFg3AlT3DePOGXni6Nb3/7CTAhThJubGc/6z9D2vT13JXz7uYFjutYeYHaw1H/7CGdsK3UFUMzdrBxf+xTsFr2cn+bWjkqkxmknNKTwnrpKxi8kr/mk0T1syL7qH+DO0eTGRYAJGh/nQO9nPKOd41IQEuRLWCygLuXX0vCfkJPHXBU9zQ7Qb7V3r86F9DJAWHwN0Xoq6C2AnQfjC4ONeVgbagtSanpOqUHnVSZgkpuX/1qj3dXOgW6s9l1UHdPTSA7qH+BDrJGia2IgEuBJBWnMbdq+4mpzyHty59i6HthtqvsqpSSPzeell56jrrsYgh1t525Fjw9LNf3Y1MpdHaqz4lrLNKKCj7q1fdupkXkWEBDI/6K6wjWvo43WXv9iABLs57CXkJTF09FYu2MGfkHGKDY21ficUCR9Zbp/7t+w6MZRDYAYY+AT3/DwKd92rAmtBak11s7VUnZhWfGAY5lFeGubpX7eXuQrcQf0ZGhdA91P9EWDfzaVpfPNqSBLg4r63PWM9Dvz1EC68WfDj8Qzo062DbCvJTqodIFlVfVh4AMeOsU//aXeDwy8rrymS2UGWyYDBZMJgtVBktGMxmKo1/PU8rLD8x/JGYVczxcuOJ88ObexMZ5s+o6FC6hwYQGeZP+5a+uDaiDYOdgQS4OG99m/wtz258lq6BXflg+Ae2m+NdWWT9InL3F3B0E6Cg01AY9rT1snIP+01HzCqqZF9mUXWgWkP2z6CtMpmtgXu6YydC2HpfZT719RNhXX2sutN8Tt7urnQL9Wd0dOiJHnW3UH+aeUuv2hYkwMV5R2vNnL1zeHfnu1wYdiEzh86s/445FgscWmMN7cQfwFRpXXlv2DPWIZJm9t0Ud3faceatP8yPezNPDEmciburwsPVBQ83FzzdXKvvrc//fNzMw/3EMc8/X3d1wdPd1Xp/0msefyvjz/uwZt60b+GDi/Sq7UYCXJxXzBYzL21+icUHFvOvjv/iuYuew921nr3B8gJYcrt1A16vZtZV92JvhPC+dh0iMVs0KxKymLf+MNuOFOLv6cbkQRGMig7Fx8PNGsiuLni6u+Dp6oqnu/W5BGrTIQEuzhsVpgqmr53OmrQ13BFzB/f1vq/+84Oz98GiG6EoHUa/bt1Zxt2+l2qXVBpZvC2dBRsPk1ZQQdsW3jx9ZRQ39G+Ln6f8SZ9P5F9bnBeOVx7n3l/vZU/uHh4f+DgTuk+of6GJP8DSu6zT/m77EdoNrH+ZZ5FWUM7Cjal8uTWNkioT/SMCeWJMJCOiQuXLv/OUBLho8jJKM7h75d0cKz3GzEtnMqz9sPoVaLHA76/A769ah0n+71MIaG2bxp7G9iOFxK0/zM/xmSiluCImjNsHd6BX2+Z2q1M4Bwlw0aQl5icydfVUDGYDc0bOoU9In/oVWFkM39wF+3+yjnVf8T+7DJmYzBZ+Schi7rrD7Eo7ToCXG3de3JGJF0bQunnj2BxBOJ4EuGiyNh7byINrHqSZZzPmjZxHx+Yd61dgXrJ1vDs/GUa9CgPvsvmXlEUVRr7cepSFG4+QcbyCiJY+PH9VD8b1aYOvjG+Lv5HfCNEk/ZDyA09veJqOzTvy4fAPCfYJrl+BB1fC17dbl2+99VvoYNvt1I7mlxO34TBfbUujzGBmYIcWPDu2B5d1D5bxbXFGEuCiSdFaExcfx1s73mJg6EBmDp2Jv0c9dl/RGtbPhNXPQ0g0jP/MZpe9a63ZmlrIvPWHWLEvG1elGNurNZMHdyA6XJaMFed2zgBXSsUBVwI5Wuvo6mMtgC+BCCAVuEFrXWi/ZgpxbmaLmVe3vsoXSV8wusNoXhz0Yv3meBvK4Lt7IWGpdauyq94Dj3pe8AMYzRZ+2pvJvPWH2ZNeRHMfd6Ze2olbL4xoUrvFCPurSQ98AfAe8PFJxx4DVmutX1FKPVb9fLrtmydEzVSaKvnvuv+y6ugqbutxGw/2fRAXVY/V6gqPwKKbrDu8D38WBj1Q7/HuonIjn285ysKNqWQVV9KxlS8vXB3NuD5t8PZoepsNCPs7Z4BrrdcqpSL+dvgq4NLqxwuB35AAFw5SVFXEfb/ex86cnUzvP52bo26uX4GH18LiiWAxw01fQZcR9Ssur4z5Gw7z1bZ0KoxmBnVuyUvXRnNp12C5KlLUS13HwEO01pnVj7OAkDO9USk1BZgC0K6dg3fzFk1OVlkWd628i7SSNF675DVGRYyqe2Faw+ZZsPxxaNkZxn8OQZ3rWJTmj0MFzFt/iNVJObi7uDA2tjWTB3UgqnVA3dsoxEnq/SWm1lorpc64eo7WejYwG6Bfv341XMNMiHNLK0njzhV3UlRVxKwRs+gf2r/uhRkr4ceHrLu+dxsD18wCr9oHrcFkYdmeY8xdd5h9mcW08PXg30M7c/OF7Qn2l/FtYVt1DfBspVSY1jpTKRUG5NiyUUKcy6GiQ9y54k6qzFXMvXwuPVr2qHthxcfgy5shYztcMh0ueazWW5mVVBr5eNMRFm5MJaekii7BfrxybQxX9w7Hy13Gt4V91DXAvwcmAq9U339nsxYJcQ77C/YzZeUUFIq4y+PoGti17oWlbbGGd1Up3PAJRI2tdRHbjxRy/6KdpBdWMKRLEK9f34uLuwQ12Y10ReNRk2mEX2D9wjJIKZUOPIM1uBcrpW4HjgANsPurENbtz6asnIKXmxdzR86t3w46Oz6GHx+2rmNyy7cQElWr080WzQdrknlr9UHCmnmx5J4L6du+Rd3bI0Qt1WQWypmWbavnikBC1M7OnJ1MXTWVZp7NmDtyLm3829StILMRfnkMts6FTpfBuHngU7vgzSyq4IFFu9h8uICxvVrzwjXRBHjJLjOiYcmVmMIpbM7czL9//TchPiHMGTmHUN/QuhVUmgtfTYQjG+Cif8OwZ8G1dn8GyxOymL5kDwaThTeu78W4PuEyXCIcQgJcNHpr09fy4JoHaRfQjjkj59R978pju6wX55TnwbVzoGftRv4qjWZmLNvHZ5uPEhPejLfHx9KxlV/d2iKEDUiAi0Zt5ZGVPLr2UboGdmXW8Fk096rjGth7voLv7wWfIJi8HFrH1ur0pKxi7vtiJweyS5lycUceGdkND7d6XOkphA1IgItGa9mhZTy5/kligmL4YPgHdVuUymKGVc/Axneh/SC4fiH4tarx6VprPvnjCC/8mEiAlzsfTx7AxV1rfr4Q9iQBLhqlJQeW8Nym5+gf2p93L3sXH3ef2hdy8mbD/e+EUS9DLRa3Kigz8OjXe1iVmM3Qbq14/fpeBPl51r4dQtiJBLhodD5L/OQkticAABL9SURBVIxXtrzC4PDBzLx0Jl5udbiC8eTNhv/1DvSdWKvTNybn8eDiXRSWGXn6yigmDYqQLypFoyMBLhqVuXvn8vaOtxnWbhivXfwaHq4etS/k5M2GJ/0EbQfU+FSj2cL/Vh7go99T6BDkS9xt/enRWtbmFo2TBLhoFLTWvL/rfWbtmcWYDmN4cfCLuLnU8teznpsNH8kv475Fu9iddpzx/dvy9L+i8PGQPxHReMlvp3A4rTVvbnuThfsWcm2Xa3n6gqdxdanl+iGnbDZ8M1zxZq02G/52ZwZPfhuPUvD+jX24omdYLX8KIRqeBLhwKIu28NLml/hy/5fc2P1Gpg+YXvuNGE7ebHj0azBgSo03XyitMvH0t/Es3ZlBv/aBvDU+ljaBdfjCVAgHkAAXDmO2mHlm4zN8l/Idk6Mn80CfB2r/RWHSj/DNPXXabHh32nHuW7STtIJy7h/WhX9f1hk3V5nbLZyHBLhwCKPFyOPrHueX1F+YGjuVu3veXbvwNlXBiqdgyywIi4UbPq7xZsMWi2b2ukO8sXw/wf6eLJpyIQM6yCJUwvlIgIsGZzAbeOT3R1iTtoaH+z7MbdG31a6AvGT4ehJk7YELpln3rHSr2WyVnOJKHlq8m/XJeYyJCeXla3rSzEcWoRLOSQJcNKgKUwUPrHmAjcc28sTAJxjffXztCtj9JSx70BrYE76EbjXfQm11Yjb/+XoP5QYTr1wbw//1bytzu4VTkwAXDabMWMa01dPYkb2D5y96nmu6XFPzk6tK4af/wO7Pod1FMG4uNAuv0amVRjOv/JzEgo2pRIYF8O6EWDoH1+GyfCEaGQlw0SCKDcXcs/IeEvITeGXIK4zpOKbmJ2ftha8mWWeZXDIdLn60xkvAJueUcO/nO0nKKmHSoAimj+ouW5yJJkMCXNhdYWUhd628i+Tjybx56ZsMa1fDvUC0tm66sPwJ8A6Eid/XeJaJ1povtqTx/LIEfD3cmH9bf4Z2D67HTyFE4yMBLuwqtzyXO1fcSXppOu9e9i6DwgfV7MSK4/D9vyHxe+g8Aq7+sMarCB4vN/DYkr38kpDFkC5BvHl9L4IDZEd40fRIgAu7ySzN5I4Vd5BbkcuHwz+kf2j/mp2YthW+ngwlx2DEDLjw3hrvEr/5UD4PfLmL3JIqHh/TnTsGd8TFRb6oFE2TBLiwi7TiNO5YcQclhhJmj5hNbHANNlCwWGDj27B6BjRrA5NXQJu+NarPZLbwzuqDvLcmmXYtfFg69SJ6tqnj5g9COAkJcGFzh4oOcefyOzFYDMy9fC5RLWuw23tprnUtk5TVEHU1jH0HvGq2CmBaQTkPfLmL7UcKGdenDc9d1QM/T/nVFk2f/JYLm9pfsJ8pK6egUMRdHkeXwC7nPunQb7B0ClQWwZVvQd/barSWidmi+XzLUV77JQmt4e3xsVwVW7OphUI0BRLgwmbi8+K5a+VdeLt5M3fkXCKaRZz9BLMJfnsZ1r0JQV3hlm8gpEeN6tqUks9zPySQlFXCwA4teP26XrRrKYtQifOLBLiwiR3ZO5i6eirNPZsz7/J5hPudoydclA5f3w5pf0DvW2D0q+Dhe8560gvLeemnRH7am0V4c2/ev7EPY2JC5YpKcV6SABf1tunYJu5fcz8hPiHMHTmXEN+Qs5+Q9CN8OxUsJhg3D2KuO2cdFQYzH/6ewqzfU1AKHhzelbsu6SgX5YjzmgS4qJe16Wt5cM2DtG/WntkjZhPkHXTmN5+ygmAvuG4+tOx01vK11izbk8nLPyVyrKiSK3uG8d8xkYQ397bxTyKE85EAF3W2InUF09dOp1uLbnw0/COae51l2l5+Cnx1W/UKglOrVxA8+w7v8RlFPP/DPrakFhAVFsBb43vLsq9CnEQCXNRabnku7+x8h++SvyM2OJb3h72Pv8dZFofas9i6gqCrO0xYBN1Gn7X8/NIq3lhxgEVbjxLo48FL11hXDnSVC3KEOIUEuKixClMFCxMWEhcfh9FiZGKPidzT6x583M8w+8NQZl1BcNdnNVpB0Gi28MmmI7y16gBlBjO3XRTBA8O6ynrdQpyBBLg4J4u28NPhn3hr+1tkl2czov0IHuzzIG0D2p75pKx466YLeQetqwdeMv2sKwiuPZDL88v2kZxTypAuQTx9ZRRdQmTJVyHORgJcnNXOnJ28tuU14vPjiWoZxasXv0rfkLNc3q41bJsHvzxeoxUEj+SXMWNZIqsSs2nf0oc5t/ZjeGSwTAsUogYkwMVppZekM3P7TFYcWUGwdzAvDn6RKzteefYd409ZQXA4XP3RGVcQLK0y8f6aZOatO4ybq+LRUd24fXAHPN1kWqAQNSUBLk5Raihlzt45fLLvE9xc3JjaayoTe0w88zj3n2q4gqDFovl2Vwav/JxETkkV1/YJZ/qo7oTIcq9C1JoEuADAZDGx9OBS3t/1PgWVBYztNJb7et937otyLBbY+A78OgMCWsPk5dCm32nfuivtOM/9kMDOo8fp1aYZH93Slz7tAu3w0whxfqhXgCulUoESwAyYtNan/8sVjdrGjI28vu11ko8n0zekLx8M/4AeLWuwJsnfVxD819vg/c+54Dkllbz2y36+3p5OkJ8nr1/Xk3F92sg63ULUky164EO11nk2KEc0sEPHD/H6ttdZn7GeNn5tmHnpTIa1G3buLxCLM2HHQut2Z1UlcOVM6DvpHysIGkwW5m84zLu/JlNlMnPXJR25d2hn/L1kWqAQtiBDKOehwspCPtj1AV8d+AofNx8e6fcIE7pPwMPV48wnWSxw+HfrDJOkn0CbodMwGPE8hEaf8latNb8m5fDCj4kczitjWPdgnrwyig5B516sSghRc/UNcA2sUEppYJbWevbf36CUmgJMAWjXrl09qxP1YTAb+CLpC2btnkW5qZzru17PPbH30MLrLJenlxfArs9hWxwUpIB3C7hwGvSbBC06/uPtyTmlzFi2j98P5NKxlS8LJvXn0m6ymbAQ9lDfAB+stc5QSgUDK5VSSVrrtSe/oTrUZwP069dP17M+UQdaa1YfXc3/tv+PtJI0hoQP4eF+D9Op+RkWktIaMrbD1nmQsBRMldB2oPVinKirwP2fM0aKK428s+ogCzam4u3uypNXRDLxogjcXWu2l6UQovbqFeBa64zq+xyl1DfAAGDt2c8SDSkhP4HXt77O9uztdG7emVnDZ3FR+EWnf7OhDPZ+ZQ3urD3g4QexN0K/yRAac9pTzBbN19vTeH35fvLLDPxfv7Y8cnk3gvzOvlCVEKL+6hzgSilfwEVrXVL9eCTwvM1aJuoluyybd3a+ww8pPxDoFchTFzzFtV2uxc3lNP/kOYnWIZLdi6CqGIJ7wBVvQs//A8/TX86utWZraiHPL0sgPqOYfu0DmX/bAGLa1GwfSyFE/dWnBx4CfFM9Y8EN+Fxr/YtNWiXqrNxYzsKEhcxPmI/JYmJS9CTuiLnjn6sFmqog8QdrcB/ZAK4e1qmA/W+3DpecZiaK0Wxh6+ECViZmsyoxm7SCCkIDvHh7fCxje7WWy9+FaGB1DnCt9SGglw3bIurBoi0sO7SMt3e8TU55DpdHXM4DfR6gjX+bU99YeAS2z4cdn0B5HgRGwPDnoPfN4PvPzRiKK438vj+XlfuyWbM/h5JKEx5uLgzuHMQ9l3Tm6t6t8fGQyUxCOIL85TUB27O38/rW10nITyC6ZTRvXPIGvYN7//UGixkOrrROATy40tq77jraOrbd6bJ/XPKeXljOqn3ZrErM4Y9D+Zgsmpa+HozqEcrwqBCGdAmS0BaiEZC/QieWVpLGzO0zWXlkJSE+Ibw85GXGdBjz14JTpTmw42PYvhCKjoJfCFz8H+g7EZr91TO3WDTxx4pYtS+blYk5JGYWA9CplS+3D+nAiMgQercLlA0VhGhkJMCdUImhhNl7ZvNZ4me4ubhxb+y93NrjVrzdvK1TAFPXW2eSJP4AFqN1OdeRM6D7FdZdcYBKo5lNh/JZuS+b1YnZZBdX4aKgX/sWPDEmkmGRwXRs5efgn1QIcTYS4I2c1pq0kjT25u09cUvKT8JoMXJV56v4d+9/E+wTDJVFsP1j65eSuUng1QwG3GkdJgnqAkBBmYFfk9JZtS+btQdzKTeY8fFw5ZKurRgeGcLQ7sG08D3L1ZhCiEZFAryRyavIIyEvgb15e4nPi2dv3l6KDdYhDW83b6JaRjGh+wSu6HgFkS0j4dguWDUD9n4NxnJo3Qeueh96XAsePhzKLWXV2hRW7cth25ECLBpCAjy5pnc4w6NCuLBjS7zcZQ1uIZyRBLgDlRvL2Ze/70RQx+fFc6zsGACuypXOzTszov0IYoJiiA6KplPzTtZ53MYKiF8KS6dZr5h084aY66D/7ZhDY9l5tJCVq4+wal82KbllAESGBXDvZV0YERlCdHiATPkTogmQAG8gJouJlOMppwyFpBxPwaItAIT7hdOzVU9ujLyRmKAYujfvjE9pLhQcst5S4qz3+Slw/AhYTBDUFUa9SnnU9axLM7JyQza/Jq2ioMyAm4viwk4tufXCCIZFBtMm8BwbMgghnI4EuB1orckozTilZ70vfx+V5koAmnk2IzoommFtLiXGK5ho7U6LkhxrQB/9DgpmwvGj1hX//uThBy06WFf+ixpLYdgQfi7pzKqkHNYv24LBZCHAy42h3YMZHhnCJd1aESDLtgrRpCmtG259qX79+ult27Y1WH0NpbCykPi8+FMCu7CqEABPV08iAyKI9mxFDJ7EVFbS5ngmqvAwFKVBdQ8cAA9/LC06UhUQQYlPO/I9wsl0a81Rwkir8iWvzEB+qYHs4koO5pQC0LaFNyMiQxkeFUz/iBayeJQQTZBSavvpNsyRHngtVZoqSSpIsg6D5FqHQtJL0wFQKDp5NOdSvIimJTFF+XQuSMadgyfON7oHUODTjjyP7mSGDOWIDuWAKZh9lUEcKvOiKNV0mlpz8XYvIMjfgyA/TyKCfLkqtjUjokLpGuIn49lCnKecOsC11hgsBqrMVRjM1vuTHxvMBqqMlVSZKzD8eW+qpMpcSZWpCoO5qvq+svpcg/Ucy59lGKrLN1BlMVJlNpBTdRwz1l5ziHalR5WR68uLiKkyEFVlwFcfodQlgGMuYaTqDvyqB3DAGMwRHUqqDuF4pR+UWAM3wMuNIH9Pgnw9CWvhQYyfJ0HVt5Z+1rBu5edJkL+HXPkohPgHp0iFJ+dfy3bzAQwKjAoM1TejDXqenhYLHho8tcZTazz+du+vNS1POhZiNhNTZaB1pQelpmBSdTipllA+0aEcIYQir7Z4+rc80VsO8vOkm58Hg/4M4+pwbunngaebTN8TQtSdUwR4oHcI4UXHcMMFd+2CG9abq3bFDVfctAtuuOFafcwVt+rnrrhgfe6KO67aDRflhot2w1W5oXBHK1csuGDB5a/H6s/nLlhOOqZxweQZyM4WHUhpHkSQnyfBfh5EVQdzC18PudxcCNFgnCLAHx7/oaObIIQQjY5MWRBCCCclAS6EEE5KAlwIIZyUBLgQQjgpCXAhhHBSEuBCCOGkJMCFEMJJSYALIYSTatDVCJVSucCROp4eBOTZsDnOTj6Pv8hncSr5PE7VFD6P9lrrVn8/2KABXh9KqW2nW07xfCWfx1/ksziVfB6nasqfhwyhCCGEk5IAF0IIJ+VMAT7b0Q1oZOTz+It8FqeSz+NUTfbzcJoxcCGEEKdyph64EEKIk0iACyGEk3KKAFdKjVJK7VdKJSulHnN0exxJKRWnlMpRSsU7ui2OppRqq5Rao5Tap5RKUErd7+g2OZJSyksptUUptbv683jO0W1yNKWUq1Jqp1JqmaPbYg+NPsCVUq7A+8BoIAqYoJSKcmyrHGoBMMrRjWgkTMDDWuso4AJg2nn+u1EFXKa17gXEAqOUUhc4uE2Odj+Q6OhG2EujD3BgAJCstT6ktTYAi4CrHNwmh9FarwUKHN2OxkBrnam13lH9uATrH2q4Y1vlONqqtPqpe/XtvJ2loJRqA1wBzHV0W+zFGQI8HEg76Xk65/EfqTg9pVQE0BvY7NiWOFb1kMEuIAdYqbU+nz+Pt4BHAYujG2IvzhDgQpyVUsoPWAI8oLUudnR7HElrbdZaxwJtgAFKqWhHt8kRlFJXAjla6+2Obos9OUOAZwBtT3repvqYECil3LGG92da66WObk9jobU+Dqzh/P2+ZBAwVimVinXY9TKl1KeObZLtOUOAbwW6KKU6KKU8gPHA9w5uk2gElFIKmAckaq3/5+j2OJpSqpVSqnn1Y29gBJDk2FY5htb6v1rrNlrrCKyZ8avW+mYHN8vmGn2Aa61NwL3AcqxfUi3WWic4tlWOo5T6AtgEdFNKpSulbnd0mxxoEHAL1t7VrurbGEc3yoHCgDVKqT1YOz4rtdZNcvqcsJJL6YUQwkk1+h64EEKI05MAF0IIJyUBLoQQTkoCXAghnJQEuBBCOCkJcCGEcFIS4EII4aT+Hzz7TDCovI+UAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"%matplotlib inline\n", | |
"\n", | |
"from matplotlib import pyplot as plt\n", | |
"\n", | |
"plt.plot(t,np.vstack((x,xh[::2],xr)).T)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Python 3", | |
"language": "python", | |
"name": "python3" | |
}, | |
"language_info": { | |
"codemirror_mode": { | |
"name": "ipython", | |
"version": 3 | |
}, | |
"file_extension": ".py", | |
"mimetype": "text/x-python", | |
"name": "python", | |
"nbconvert_exporter": "python", | |
"pygments_lexer": "ipython3", | |
"version": "3.7.3" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 4 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment