{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "understanding-monte", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.integrate import solve_ivp\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "scenic-buddy", "metadata": {}, "outputs": [], "source": [ "def myderiv(x,y):\n", " return y" ] }, { "cell_type": "code", "execution_count": 3, "id": "gross-spank", "metadata": {}, "outputs": [], "source": [ "def euler(initialx, initialy, stepsize, derivFn):\n", " x = initialx\n", " y = initialy\n", " while True:\n", " #print(y)\n", " yield y\n", " x += stepsize\n", " y = y + stepsize * derivFn(x,y)\n", " \n", " " ] }, { "cell_type": "code", "execution_count": 4, "id": "lightweight-playlist", "metadata": {}, "outputs": [], "source": [ "def heun(initialx, initialy, stepsize, derivFn):\n", " x = initialx\n", " y = initialy\n", " while True:\n", " yield y\n", " x += stepsize\n", " yest = y + stepsize * derivFn(x,y)\n", " y = y + stepsize * (derivFn(x,y) + derivFn(x+stepsize, yest))/2\n", " \n", " \n", " " ] }, { "cell_type": "code", "execution_count": 5, "id": "accredited-chapter", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "stepsize = 0.1\n" ] } ], "source": [ "numsteps = 50\n", "limit = 5\n", "stepsize = limit/numsteps\n", "print (\"stepsize = \",stepsize)\n", "x = np.linspace(0,limit,numsteps)\n", "y1 = np.exp(x)\n", "y2 = np.fromiter(euler(0, 1, stepsize, myderiv), float, len(x))\n", "y3 = np.fromiter(heun(0, 1, stepsize, myderiv), float, len(x))" ] }, { "cell_type": "code", "execution_count": 6, "id": "monetary-channels", "metadata": {}, "outputs": [], "source": [ "# Result is in y4.y[0]\n", "y4 = solve_ivp(myderiv, (0, limit), (1,), t_eval = x, vectorized=True)" ] }, { "cell_type": "code", "execution_count": 7, "id": "developing-stereo", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA90ElEQVR4nO3deXiU1fXA8e/JRthRohFlCVqsCrIIqKNoJ0YRqSKK4gparVgsVqUu2IqAC0sEtajFrVRQhFIQFeuCJhkVGRRQFAT9ARL2HUXCNlnO74/3TTIJ2bdJJufzPPPMzH2XuXdCTi7nve+9oqoYY4wJLxGhroAxxpiqZ8HdGGPCkAV3Y4wJQxbcjTEmDFlwN8aYMGTB3RhjwpAFd1MriMhNIrIg1PXIJSINRWS+iOwTkf+Guj41RUReE5Enquhco0Xkjao4lyk/C+5hRkRuFJGlIpIhIttE5AMR6RXqepVGVWeoau9Q1yPINUA80FJVry1qBxE5Q0Tedf8A7BeRNBE5L2h7goio+7PIEJEdIvKeiFxS6DzpInIoaL8MEXm+epsHInKriCys7s8xoWHBPYyIyHDgWWAsTmBqC/wTuDKE1SqViESFug5FaAf8n6pmFbVRRE4BvgBWAO2BE4F5wAIR8RTavYWqNgG6AB8D80Tk1kL7XKGqTYIew6qwLaY+UlV7hMEDaA5kANeWsE8DnOC/1X08CzRwt3mBzcCDwE5gG9Af6Av8H7AX+FvQuUYDc4D/APuBr4EuQdtHAOvcbauAq4K23YoTGJ8B9gBPuGUL3e3ibtsJ/IoTQDsFtXM6sAvYADwCRASddyEwEfgZWA9cVsL3cTrgA34Bvgf6ueVjgACQ6X6ntxdx7OvA+0WUTwE+c18nAApEFdrnfmBHUL3TgYvL+HMeDfwXeMP9blcApwIPu9/XJqB3oX8X/3J/nlvc7zrSbfthINtt4y/u/q8BLwD/c8//JXBK0PnOA5YA+9zn84K2tQc+dY/7GHgeeMPdFuvWeY/7fS8B4kP9exPOj5BXwB5V9IOEPkBW4UBSaJ/HgMXA8cBxwCLgcXeb1z3+USAauMMNoG8CTYGOwCGgvbv/aDf4XePuf78bTKPd7dfi9GYjgOuAA0Ard9ut7mfdDUQBDSkY3C8FlgEtcAL96UHHTgfeceuUgPOH5/ag82a6dY8EhuL8EZMivotoYC3wNyAGuMgNSr8Nat8bJXyX24E/FFGe6AbMhhQf3E92y09336dTvuB+2P2OotzvYz3w96Cf2/qg/ecBLwGN3Z/7V8CdQd/XwkLnfw0nAJ/tnn8GMMvddizOH81B7rYb3Pct3e1+4GmcTsSF7veZG9zvBOYDjdyfTXegWah/b8L5EfIK2KOKfpBwE7C9lH3WAX2D3l8KpLuvvTjBO9J939QNQOcE7b8M6O++Hg0sDtoWgdM7vKCYz14OXOm+vhXYWGh7XqBxA+3/Aefi9m7d8kicHvUZQWV3Ar6gc6wN2tbIbcMJRdTnApwAHXz+mcDooPaVFNyzgD5FlJ/mfuZJFB/cY93y89336bi956DHHcV87mjg46D3V7jHFv65tcBJzR0BGgbtfwOQVvg7D9r+GvBq0Pu+wA/u60HAV4X297vnaet+J42Dtr1JfnC/Dacz0TnUvyv15VEbc52mYvYAcSISpcXkiXF60huC3m9wy/LOoarZ7utD7vOOoO2HgCZB7zflvlDVHBHZnHs+ERkMDMcJcLjHxRV1bGGqmupeUHwBaCcib+H8z6AhTu+0cBtOCnq/Peg8B0Uk97MLOxHYpKo5JZyrJLuBVkWUtwJycHq0xxdzbO5n7A0q66+qn5Txswv/THYX8XNrgtPGaGCb+z2A80e42O/etT3o9UHyv7/C/34g/zs7EfhZVQ8U2tbGff26+3qWiLTASdH8XVUzS6mLqSC7oBo+/Di9tP4l7LMV50JhrrZuWUXl/uIiIhFAa2CriLQDXgGG4fyXvQWwEifFkqvE6UhVdbKqdgfOwMkpP4ATUDOLaMOWCtR9K9DGrXdFzvUJTuqpsIGAX1UPlnDsVTj58R/L+FkVtQnn30ScqrZwH81UtaO7vbxTwhb+9wP539k24BgRaVxom/NBqpmqOkZVz8DJ218ODC7n55tysOAeJlR1H06+/AUR6S8ijUQkWkQuE5Fkd7eZwCMicpyIxLn7V2YccncRudod7XIvTiBZjJPfVZycPSLyB6BTWU8qIj1F5BwRicbJ1R8Gctze6WzgSRFp6v4RGV7BNnyJ0yt90P2evDgpjlllPH4McJ6IPCkix7r1uRsnYD1UTLviRWQYMAp4uND/Gqqcqm4DFgCTRKSZiESIyCki8jt3lx1AaxGJKeMp3wdOdYfbRonIdTh/fN9T1Q3AUmCMiMS4w2+vyD1QRBJF5EwRicS5SJ6J8z8cU00suIcRVZ2EE+wewQmsm3B6z2+7uzyB8wv4Hc4oi6/dsop6B+diae5FtqvdHtoqYBLO/yZ2AGfijI4pq2Y4Pf+fcf5rvwd4yt12N07A/wlnZMybwNTyVlxVAzjB5zKc/xH8Exisqj+U8fg1QC+c4Y3pOD3XAcClqlq4rb+IyAGc77wvzoimwnWeX2ic+7zytqkYg3EuGK/C+T7nkJ9OSsUZJbRdRHaXdiJV3YPT4/4rzs/kQeByVc099kbgHJx00yici725TnA/+1dgNc6omtcr0zBTMnEvdhhTLiIyGviNqt4c6roYY45mPXdjjAlDFtyNMSYMWVrGGGPCkPXcjTEmDNWKm5ji4uI0ISGhQsceOHCAxo0bl75jGLE21w/W5vqhMm1etmzZblU9rqhttSK4JyQksHTp0god6/P58Hq9VVuhWs7aXD9Ym+uHyrRZRArfMZzH0jLGGBOGLLgbY0wYsuBujDFhqFbk3IuSmZnJ5s2bOXz4cIn7NW/enNWrV9dQrUIjNjaW1q1bEx0dHeqqGGPqiFob3Ddv3kzTpk1JSEggaLrSo+zfv5+mTZvWYM1qlqqyZ88eNm/eTPv27UNdHWNMHVFr0zKHDx+mZcuWJQb2+kBEaNmyZan/gzHGmGC1NrgD9T6w57LvwZjw5F+kzJjRFr+/6s9da9MyxhgTzvwfZ5B0aSQBbceMGUpKiuDxVN35a3XPvTZ4++23ERF++KHkab6fffZZDh4safGdkr322msMGzaswscbY+oW36trCWg02UQSCAg+X9We34J7KWbOnEmvXr2YOXNmiftVNrgbY+oX747/EEOASMkmJgaq+sZcC+4lyMjIYOHChfzrX/9i1ixn9bXs7Gzuv/9+OnXqROfOnXnuueeYPHkyW7duJTExkcTERACaNMlfk3nOnDnceuutAMyfP59zzjmHbt26cfHFF7Njx46jPtcYE+YOHsTz1T9IIYk7rl9FSgpVmpKBupJzL+GCYqUGQZYy3fE777xDnz59OPXUU2nZsiXLli3jq6++Ij09neXLlxMVFcXevXs59thjefrpp0lLSyMuLq7Ec/bq1YvFixcjIrz66qskJyczadKkyrTCGFPXfPghHDqE5xzlyJA9VR7Yoa4E9xCZOXMm99xzDwDXX389M2fOZP369fzpT38iKsr56o499thynXPz5s1cd911bNu2jUAgYGPXjamP3nrLeb766mr7iLoR3EvoYVfXTUx79+4lNTWVFStWICJkZ2cjIvTs2bNMxwcPXwweo3733XczfPhw+vXrh8/nY/To0VVddWNMbRYIwPz5zuurr4bNm6vlYyznXow5c+YwaNAgNmzYQHp6Ops2baJ9+/Z06dKFl156iaysLMD5IwDQtGlT9u/fn3d8fHw8q1evJicnh3nz8hey37dvHyeddBIA06ZNq8EWGWNqhZQU+PVX6NwZfvObavsYC+7FmDlzJldddVWBsgEDBrBt2zbatm1L586d6dKlC2+++SYAQ4YMoU+fPnkXVMePH8/ll1/OeeedR6tWrfLOMXr0aK699lq6d+9ean7eGBOGaiAlAzhzl5T0AKYCO4GVRWz7K6BAnPtegMnAWuA74KzSzq+qdO/eXQtbtWrVUWVF+fXXX8u0X10X/H2kpaWFriIhYm2uH8K+zVlZqnFxqqC6YoWqVq7NwFItJq6Wpef+GtCncKGItAF6AxuDii8DOriPIcCUCvy9McaY8PT557B7N3ToAB07VutHlRrcVfUzYG8Rm54BHsTpuee6Epju/lFZDLQQkVZFHGuMMfVPcEqmmueMqlDOXUSuBLao6reFNp0EbAp6v9ktM8aY+i0nJz+4DxhQ7R9X7qGQItII+BtOSqbCRGQITuqG+Ph4fIUmVmjevHmB0SfFyc7OLtN+dd3hw4fzvqOMjIyjvq9wZ22uH8K5zU1XraL7li0cPv54FmdkQDX/PldknPspQHvgW3csd2vgaxE5G9gCtAnat7VbdhRVfRl4GaBHjx5aePXv1atXl2n8ergv1pErNjaWbt26AbZCfH1hbQ4zH3wAQOwNN+B1R9VB9bW53GkZVV2hqseraoKqJuCkXs5S1e3Au8BgcZwL7FPVbVVbZWOMqWNUa24IpKvU4C4iMwE/8FsR2Swit5ew+/vATzhDIV8B7qqSWoZIZGQkXbt2zXuMHz++xP1t2l5jTJFWrIC1a+H44+H882vkI0tNy6jqDaVsTwh6rcCfK1+t8klOhp49Ieh/OqSlwZIl8OCDFT9vw4YNWb58eaXrV5ysrKy8OWqMMWEst9fevz9ERtbIR4bFHao9e8LAgU5AB+d54ECnvDokJCSwe/duAJYuXVpkvmzXrl0MGDCAnj170rNnT7744gvAuUN10KBBnH/++QwaNKh6KmiMqV1qOCUDdWXisFIkJsLs2U5AHzoUpkxx3gf35Cvi0KFDdO3aNe/9ww8/zHXXXVemY++55x7uu+8+evXqxcaNG7n00ktZvXo1AKtWrWLhwoU0bNiwchU0xtR+a9bgX9EYX4PReBteRDXM7luksAju4ATyoUPh8cdh5MjKB3aoXFrmk08+YdWqVXnvf/31VzIyMgDo16+fBXZj6gn/pEUkkULgSANi+kQWWJjDv8nPjI0zaLCpAZ42VRv2wya4p6U5PfaRI53nxMSqCfBFiYqKIicnByg4nW+wnJwcFi9eTGxs7FHbGjduXD0VM8bULqr45u0lQIy7VqozvN3jcQJ70vQkjmQdYcamGaQMTqnSAB8WOffcHPvs2fDYY/kpmtwcfFVLSEhg2bJlAMydO7fIfXr37s1zzz2X9746L8waY2qpb7/Fu3O2s1ZqpBZYK9WX7iOQHSCHHALZAXzpvir96LAI7kuWFMyx5+bglyyp3Hlzc+65jxEjRgAwatQo7rnnHnr06EFkMVe+J0+ezNKlS+ncuTNnnHEGL774YuUqY4ype954Aw+LSRkwhccflwIpGW+Cl5jIGCKIICYyBm+Ct2o/u7jpImvyYVP+ls6m/E0LdRVqnLW5jsvKUj3xRGd630WLitxl0cZF+sdpf9RFG4veXhpKmPI3bHLuxhhTq3z6KWzdCiefDOeeW+QunjYejrQ9UuUXUyFM0jLGGFPrzJjhPN94Y7VP71sUC+7GGFPVDh+GOXOc1zfdFJIqWHA3xpiq9t57ziLY3bvDaaeFpAoW3I0xpqrlpmRC1GsHC+7GGFO1fv4Z3n8fIiLg+utDVg0L7iVo0qRJgfc2pa8xplRz5kAgABddBK1Ct4S0BXdjjKlKb7zhPIcwJQNhFtz9m/yM+3wc/k3+av+skqb0nThxYt5+nTp1Ij09nfT0dE4//XTuuOMOOnbsSO/evTl06FC119MYU4M2boTPPoPY2Bqd3rcoYXMTU+4kPIHsADGRMVUyCU/hKX/37t1Lv379gJKn9C3OmjVrmDlzJq+88goDBw5k7ty53HzzzZWqozGmFpk503nu1w+aNQtpVcImuOdOwpOt2XmT8FQ2uBee8ve1115j6dKlQMlT+hanffv2eX8sunfvTnp6eqXqZ4ypZWrBKJlcpQZ3EZkKXA7sVNVObtlTwBVAAFgH/EFVf3G3PQzcDmQDf1HVj6qn6gXlTsKT23Ov8kl4CiluSt/g6YCh4JTADRo0yHsdGRlpaRljwoj/9bX4Vvweb9NWePr0CXV1ypRzfw0oXNOPgU6q2hn4P+BhABE5A7ge6Oge808RqZEFAz1tPKQMTuHxxMerfF7kohQ3pW9CQgJff/01AF9//TXr16+v1noYY0LP74ek29oyksdJOjgf/7KYgttr8HpgrrIskP2ZiCQUKlsQ9HYxcI37+kpglqoeAdaLyFrgbKBGWuRp46n2oJ5r8uTJ/PnPf6Zz585kZWVx4YUX8uKLLzJgwACmT59Ox44dOeecczj11FNrpD7GmNDxpWQTyIogmygCmpO3IAdUz/XAshBn1shSdnKC+3u5aZlC2+YD/1HVN0TkeWCxqr7hbvsX8IGqziniuCHAEID4+Pjus2bNKrC9efPm/OY3vym1btnZ2cXOqR5O1q5dy759+wDIyMg4agx+uLM21w91tc2bZ/zEkFevJ0AMUQ2ESZO+pWPHXwGYsXEGU9dPJYccIojgtva3cVPb/Jx8ZdqcmJi4TFV7FLWtUhdUReTvQBYwo7zHqurLwMsAPXr0UG/u8iSu1atX07Rp01LPs3///jLtV9fFxsbSrVs3AHw+H4W/r3Bnba4f6mybJ0zgFF7B1/cpvI/0wuM5K29Tg00NmLFpRl7P/bbE2wr03KurzRUO7iJyK86F1iTN7/5vAdoE7dbaLTPGmPC0YQN89BGemGg8006DuIKbc68H+tJ9eBO8NZY6rlBwF5E+wIPA71T1YNCmd4E3ReRp4ESgA/BVpWtpjDG11dSpoOrctBQXV+QuNXk9MFdZhkLOBLxAnIhsBkbhjI5pAHwsziT0i1X1T6r6vYjMBlbhpGv+rKrZ1VV5Y4wJqexsJ7gD3HFHaOtSSFlGy9xQRPG/Stj/SeDJylTKGGPqhA8/hM2b4ZRToJZdKwiruWWMMaZGvfKK8/zHPzpT/NYitas2tUxkZCRdu3alU6dOXHHFFfzyyy8ApKen06lT/qjQV155he7du/Pzzz/nlU2aNAkRYffu3YBzRbx58+Z07dqVrl278thjj9VoW4wxVWzbNmfFpagouPXWUNfmKBbcS5A7t8zKlSs59thjeeGFF47a5/XXX+e5557jo48+4phjjgFg06ZNLFiwgLZt2xbY94ILLmD58uUsX76cRx99tEbaYIypJv/+t5Nzv+IKOOGEUNfmKGEV3P1+GDfOea5qHo+HLVsKjuqcPXs248ePZ8GCBcQFXSW/7777SE5ORkKw4rkxpgbk5MCrrzqva9mF1FxhMyuk3w9JSc4CKDExkJKSf/tvZWVnZ5OSksLtt9+eV7ZhwwaGDRvGN998wwlBf7XfeecdTjrpJLp06VJEHf106dKFE088kYkTJ9KxY8eqqaAxpmalpsL69dC2LfTuHeraFClseu4+nxPYs7OdZ5+v8ufMnc/9hBNOYMeOHVxyySV524477jjatm3L7Nmz88oOHjzI2LFji8ynn3XWWWzYsIFvv/2Wu+++m/79+1e+gsaY0Mi9kHrbbVBLpz8Jm+Du9To99shI57kqRiXl5tw3bNiAqhbIuTdq1Ij333+fF198kRnuHM7r1q1j/fr1dOnShYSEBDZv3sxZZ53F9u3badasWd78EX379iUzMzPvYqsxpg7ZtQvmzXNGx9x2W6hrU6ywSct4PE4qxudzAntVpWTACeSTJ0+mf//+3HXXXXnlxx9/PB9++CFer5e4uDguvfRSdu7cmbc9ISGBpUuXEhcXx/bt24mPj0dE+Oqrr8jJyaFly5ZVV0ljTI3wP/Yxvsy/4vUE8LTJn23Fv8lf41MMlCRsgjs4Ab0qg3qwbt260blzZ2bOnMkFF1yQV96+fXveffdd+vbty7x58zj77LOLPH7OnDlMmTKFqKgoGjZsyKxZs+yCqzF1jP+LHJJeuIoAA4lZBil+J+aEalrfkoRVcK9qhZfNmz9/ft7rlStX5r3u0qXLUSNpgALL6A0bNoxhw4ZVfSWNMTXG98oaAnqKM297tubN217eZT6Tk6FnT0hMzC9LS4MlS+DBB6umrmGTczfGmOrmXT2FGAJESjYxMZJ3bS93mc9IiSzTMp89e8LAgU5AB+d54ECnvKpYz90YY8ri++/xfPUPUhp8i++B/+Ht2ygvDVzeaX0TE2H2bCegX3ZZAh984LwP7slXVq0O7qpqeWmc78EYE2KTJwPgue10PI83Ompzeaf1TUyEoUPh8ccTGDmyagM71OK0TGxsLHv27Kn3gU1V2bNnD7GxsaGuijH11549MH268/ovf6mSU6alwZQpMGhQOlOm5Kdoqkqt7bm3bt2azZs3s2vXrhL3O3z4cNgHvtjYWFq3bh3qahhTf738Mhw+DJddBqedVunT5ebYZ88GkXT+8IeEvPdV1YOvtcE9Ojqa9u3bl7qfz+fLW1vUGGOqXGYm5N7AeM89VXLKJUvyA7nPl5+DX7KkHgR3Y4ypFebOhS1b4PTTq2wemaKGOyYmVm3evdbm3I0xplZ49lnn+Z57oA4N8Cg1uIvIVBHZKSIrg8qOFZGPRWSN+3yMWy4iMllE1orIdyJyVnVW3hhjqtXixfDll3DMMTBoUKhrUy5l6bm/BvQpVDYCSFHVDkCK+x7gMqCD+xgCTKmaahpjTAjk9trvvBMaHT38sTYrNbir6mfA3kLFVwLT3NfTgP5B5dPVsRhoISKtqqiuxhhTczZtgjlznKlm//znUNem3Cp6QTVeVbe5r7cD8e7rk4BNQfttdsu2UYiIDMHp3RMfH4+vghOwZ2RkVPjYusraXD9Ym0Pr57Fp/JD9AJ27Z9B47VpYuxaA7/d9z/J9y+navCsdm1d+wZ1qa7OqlvoAEoCVQe9/KbT9Z/f5PaBXUHkK0KO083fv3l0rKi0trcLH1lXW5vrB2hw6iz45oA05oJFkasMGWbpokVu+cZE2fKKhRo6J1IZPNNRFGxeVeJ4JE1RTUwuWpaY65bkq02ZgqRYTVys6WmZHbrrFfc6dxHwL0CZov9ZumTHG1Bm+ScsIEOPM/pgVmbeyW1GzP5akJiYIK05Fg/u7wC3u61uAd4LKB7ujZs4F9ml++sYYY2q/Q4fwfjnBmf0xIqfAym7lnf0xeIKwRx+lyu9CLUmpOXcRmQl4gTgR2QyMAsYDs0XkdmADMNDd/X2gL7AWOAj8oRrqbIwx1WfqVDx7/0fKqXfhu+XfeBOp8OyPEDxBGNUyQVhxSg3uqnpDMZuSithXgbp3WdkYYwACAZgwAQDPuH54rj76pqXyzv6YO0HYyJHOc1XfiVocu0PVGGNyvf66MwTyjDOgf/9Kny54grDHHstP0VT1DJBFseBujDEAWVkwbpzz+m9/g4jKh8fgCcKg4ARh1c0mDjPGGHCi7rp1cMopcN11VXLKmpggrDjWczfGmJwcePJJ5/XDD0NU3e/3WnA3xpi334ZVq6BNmzo3QVhxLLgbY+o3VXjiCef1Qw9BTExo61NFLLgbY+q3Dz6Ab76B+Hi47ba8Yv8mP+M+H4d/kz+Elau4up9YMsaYilLF/+A8fIzAe00nPA0bAk5gT5qeRCA7QExkDCmDU4od256c7EwnEHyRNC3NGRFT1AXVmmI9d2NMveWfvISk7//BSB4naeqN+N1OennmkAnl/DElseBujKmfcnLwJX+VP0FYQPImCCvPHDKhnD+mJJaWMcbUT//9L96tM4jhdgKRkcTESN4EYeWdQyZU88eUxIK7Mab+ycyERx7Bw1pSHvwIX4v+eL35E4RB+eaQCdX8MSWx4G6MqX9efdVZWenUU/E8eTmeSkTC4PljcoN6bUjNWM7dGFO/ZGTAmDHO6yefrPTdqKGcP6Yk1nM3xtQvzz4LO3bA2WfDgAGVPl0o548pifXcjTH1x+7dzsB0gPHjQY6erz1cWHA3xtQfY8fC/v1w6aWh71pXMwvuxpj6YcMGeOEF53XuvO1hrFLBXUTuE5HvRWSliMwUkVgRaS8iX4rIWhH5j4iExyw8xpg6zT90OuMCw/Ff8ih065ZfXsY5ZJKTj15BKS0tP8tT21T4gqqInAT8BThDVQ+JyGzgepwFsp9R1Vki8iJwOzClSmprjDEV4H9jHUkf/JUAMcR8HkGK3xnTXp45ZHKnGcgdGRM8BLI2qmxaJgpoKCJRQCNgG3ARMMfdPg3oX8nPMMaYilPFN9qXP81AZkTeNAPlmUOmtk4zUJwK99xVdYuITAQ2AoeABcAy4BdVzXJ32wycVNTxIjIEGAIQHx+PL/fbLqeMjIwKH1tXWZvrB2tz1Yj79FO86z4khhs4EhFBVJTSrNm3+Hy/0mxfM6IkClUlSqJotrdZiZ8vApddlsDjjycwaFA6IulUtrrV9nNW1Qo9gGOAVOA4IBp4G7gZWBu0TxtgZWnn6t69u1ZUWlpahY+tq6zN9YO1uQocOKDapo0q6KIH3tKxY1UXLSq4y6KNi3TsZ2N10cZFRZ8jSGqqalyc6siRznNqauWrWJk2A0u1mLhamZuYLgbWq+ouABF5CzgfaCEiUer03lsDWyrxGcYYU3HjxsGmTdCtG55x/fBEHr1LWeeQqa3TDBSnMjn3jcC5ItJIRARIAlYBacA17j63AO9UrorGGFMBa9fmD2V5/nmILCKyl0NtnWagOJXJuX8pInOAr4Es4BvgZeB/wCwRecIt+1dVVNQYY8rlvvsgEIBbboHzzqv06WrrNAPFqdTcMqo6ChhVqPgn4OzKnNcYYyrlvfecR7NmzjQD9ZDdoWqMCS+HD8O99zqvx4yBE04IaXVCxYK7MSa8TJoE69ZBx47w5z/nFYfrnajFsSl/jTFhwz9vO77RR/ByLp7nxkJ0tFMexneiFsd67saYsOBfpCRdcwwjsx4lKdKHPzb/Smc434laHAvuxpiw4HtuBYGcSGeKAWIK3DnqTfASExlDpEQSExmDN8Fb4rmCF7weOrTuBXawtIwxJhzs2oX3g4eIYS6BiAhiYiLwevM3e9p4SBmcgi/dhzfBW+pNS7VxwevysuBujKn7/vIXPPs+JKX7Q/iunow30Zn1MVi43olaHAvuxpi67e23YdYsaNQIz+z78JxcuaXzSroT1YK7McbUhJ9/dpLi4Mwjc/LJlT5lXbsTtTh2QdUYU3fddx9s3w7nnw/DhpXr0HAZz14cC+7GmLrpgw9g2jSIjYWpUyEiosw3KkH+ePbcAJ+ba+/Zs5rrXUMsLWOMqXP8nxzAd/23zs1Kj10Np55arhuVoOB49qFDnVExde2iaUksuBtj6hS/H5L6RBPIvp8YuZeUc6PwUPSNSqWNjgkezz5yZPgEdrC0jDGmjvH983sC2RHOzUoRDfAtdPqo5b1RCY4ez144B1+XWc/dGFN3bNqE953hxDCPgBS8WakiNyqFw3j24lhwN8bUDdnZcPPNePZ/RopnJL7LJx51s1JZb1SC8BnPXhwL7saYuuHJJ+Gzz+CEE/C8/RCe4yt3s1K4jGcvjuXcjTG138KFzsIbIvD663D88WU+NNzHsxenUsFdRFqIyBwR+UFEVouIR0SOFZGPRWSN+3xMVVXWGFMP/fwz3Hgj5OTAQw/BxReX6/BwH89enMr23P8BfKiqpwFdgNXACCBFVTsAKe57Y4wpP1X44x9h0yY4+2x47DGg7KsqQfjMz15eFc65i0hz4ELgVgBVDQABEbkS8Lq7TQN8wEOVqaQxpn7yj3gH31un4m2UhGfmyxAdXe6blSC8x7MXR1S1YgeKdAVeBlbh9NqXAfcAW1S1hbuPAD/nvi90/BBgCEB8fHz3WbNmVageGRkZNGnSpELH1lXW5vqhvrf5pw/3MWzCJQSIIToqh4nPrqRjx1+ZsXEGU9dPJYccIojgtva3cVPbm0o87zfftGDMmDPo128r7757IqNGraJbt19qoEWlq8zPOTExcZmq9ihyo6pW6AH0ALKAc9z3/wAeB34ptN/PpZ2re/fuWlFpaWkVPrausjbXD/W6zTt26Njm4zWSTAXVyEjVsWOdTYs2LtKGTzTUyDGR2vCJhrpo46ISz5maqhoX5zwX9T7UKvNzBpZqMXG1Mjn3zcBmVf3SfT8HOAvYISKtANznnZX4DGNMfZOZCQMH4t33NjERWURGKjExHHWz0uOJj5cpJVPSePZwVuGcu6puF5FNIvJbVf0RSMJJ0awCbgHGu8/vVElNjTH1w/Dh8OmneFq1IuXF/fi+j8XrrfjNSuE+nr04lb2J6W5ghojEAD8Bf8AZgTNbRG4HNgADK/kZxph64oT334fnn4eYGHjrLTznHoenX+nHJSc7QxuDA3ZamtM7Lyq41weVCu6quhwn915YUmXOa4yph778klOffdZ5/cILcO65+Df5yzRXTO5Y9tz0S/C8MfWVTT9gjAm9bdvg6quJyMyEu+6CP/6xXEMew31u9oqw6QeMMaF15Aj+3qMYt3UwC065GZ55Bih6fvaSBI9lHzq0fgd2sJ67MSaUcnLwX/4kSSufJUAMUZsgbVkUHk/+/Oy5PffS5mcvPDd7fbhoWhIL7saY0HngAXyfxBAghmyi0KwcfD5nZEx55mcP97nZK8KCuzEmNJ5+Gp5+Gm9kL2KihEAWREVp3nh2OHrIY3GjYp56KrznZq8IC+7GmJo3axb89a8AeKb9iZSTI/H5oFmzb/F4zir2sJJGxRQO4paWMcaYmpSaCoMHO6+Tk+Gmm2CTH3r5YG8znBvdi2ajYsrOgrsxpkb4/eCbtQ3vq2PxZGbCX/4C999fYMhjlERx1llnlZhfr48zPFaEDYU0xlQ7vx+SLsph5OTjSDr4Lv6L/uYMeRQpMOQxMyez1CGPhUfFFF5lyTgsuBtjqp3v7V8IHM4hmygCxODzjoYIJ/zkDnmMlEiiI6JLHPIYnGN/7LH8FI0F+KNZWsYYU722bME7415imEYAiGkYiffi/MWtg4c8NtvbrMSUTEkzPFp6piAL7saY6rN1KyQm4tmyhpRTG+O7bgrey6IKzPAI+UMefT4fUPyQR7BRMWVlwd0YUz22bYOLLoI1a6BrVzwpT+M5tmGZDrWJwCrPcu7GmKq3Y4cT2H/8ETp3hk8+wX/gR1vUugZZz90YU6X8/9uL77Y38e5sgadTJyewH/w/W9S6hlnP3RhTZfzztpN0RUNG7rybJEnFP/5TOO64cs/wCDbksbIsuBtjqsaqVfhufY2ARjtDHiNi8X13LFBwuGPwDI/JyUcH7bQ0uPNOG/JYWZUO7iISKSLfiMh77vv2IvKliKwVkf+4S/AZY8LZ4sVwwQV4f30naFFrKXVR69wLp7lB+5tvWjDQXZizPi5qXZWqIud+D7AaaOa+nwA8o6qzRORF4HZgShV8jjGmNvrgA7jmGjh4EM/lcTw76HvmpuxnwGUt8XjOzNutqEWtC88VM3nyGcybV3R+3YY8lk+leu4i0hr4PfCq+16Ai4A57i7TgP6V+QxjTO3k98O465bjv2IsHDwIt96K/7kHuXfNBaScdDH3fn9OmUfG5F447ddvqwXwKlLZtMyzwINAjvu+JfCLqma57zcDJ1XyM4wxtYx/kZL0u0xGzu5EUvZH+G98DqZOxbdpYaUunL777omWV68iFU7LiMjlwE5VXSYi3gocPwQYAhAfH593Z1p5ZWRkVPjYusraXD/U1jZLZiYf3r2VQOZNzoVTEaY26suRTz+l2b5mREkUqkqURNFsbzPuvHMdp522n27dfsk7xzfftOCHH5py2mn7GTPmDEaNWkW3br8QG9uAq67qnve+Pqi2n7OqVugBjMPpmacD24GDwAxgNxDl7uMBPirtXN27d9eKSktLq/CxdZW1uX6olW3etUv1wgt1EedqQw5oZES2NmyoumhR/i6LNi7SsZ+N1UUbncLUVNW4OOe58PsJE/LLVZ0255bXF5X5OQNLtZi4WuGeu6o+DDwM4Pbc71fVm0Tkv8A1wCzgFuCdin6GMaYWWbEC+vWD9HQ8rVrx7AMpzP2haakXTktaYMMunFaf6hjn/hAwXETW4uTg/1UNn2GMqUnvvgvnnQfp6dCjB/73XuTeg9eV+cJp8EXToUMteNeEKgnuqupT1cvd1z+p6tmq+htVvVZVj1TFZxhjQkAV/9DpjLtyMf6MTnDDDfDZZ/gOfF+uC6d2t2nNs7lljDFF27MH/5XjSPriMQLEEBOtpAyLwtNQ8u44zZ0rxpvgLXaa3lmz4K23CqZibCKw6mfTDxhjjvbFF9C1K74vogkQ44yKyYnG96mzyEZRd5wWvts0d5pesLtNQ8F67saYfDk5zoQvjzwC2dl4O+4mai3kZGYTFQ1eb2TeruW5cFqYXTStftZzN8YAzlS9406bhv/hdyA7Gx58EOYPRm+5GBIfRQcnQWu7cFpXWM/dGIP/2S9Juq8zAQYRw3WkTFyO56/n4ft8HNknLURP/JRsicSX7uPzmZ4ic+tLljg59+ALp9ZDDx3ruRtTn+3fD0OG4LvvbQK4U/VGNsQXOA8oeqre4nLrUVE2TW9tYj13Y+qrlBS4/XbYsAFv1AVESQ452QVz67kXTn3pPrwJXifH3qbo3PqSJcVfOLXee82z4G5MPeL3g++jI3hXPIfnrQecwu7dYdIQ9M3esO589JQvoPU4nNlDip+qt/ASeHbhtHax4G5MPeH3Q1JiNoEjkcRwFymR7+AZfSk89BC+xRMttx5mLOduTH2wdSu+O2cSOKJOXp0YfHfNdoY8Rkdbbj0MWXA3JpxlZcEzz8Bpp+FdMZkoCSCSRVSs4L2hVd5uRd2UFDxu/dFH8wN6VpbdlFQXWFrGmDDk94NvWjreT0biWfeGU3hDFNqkL/xUMK+eP21Afm49N/3y4IOWW6+rrOduTJjxv72DpAuOMPKl1iStewn/CVfBu+/iG9qX7NYL0QvGkn3SwrzJvopLv/TsaRN+1WUW3I0JF/v2wYgR+K55nkB2pLtKUgN8f5oFV1xRZF4dKDb9ApZbr8ssuBtT1wUC8I9/wCmnwIQJeLM/ISoyC4nIJqoBeHvHkJwMh9cWzKsfXuukZKDoaQNKGrduaj/LuRtTV2Vn438yFd+zy/H+/BYe9sCFF8LfBqNv9SkwZr1nT4/bC/fwcKInL/WS20MvnH5JTHTy7YVZbr3usJ67MXVNdjbMnIn/lJtJGnU+I3++jyRJxZ/8Ofh8+BrtdMasB+XWi0u9JCZSINBb+iV8WHA3pq7IyoIZM6BTJ7jxRnwbEjjizrV+JKIBvqxeJD8lNN1bMLfedK+zkEZxMzZa+iU8VTgtIyJtgOlAPKDAy6r6DxE5FvgPkACkAwNV9efKV9WYeioQwP/EJ/j+uRrvnjl4+AESEmh5dVtyngtAlpITkUnL09fRofmZDBzoIfnVFPYf66PpXi9j/uhh9uyiUy+Wfglflcm5ZwF/VdWvRaQpsExEPgZuBVJUdbyIjABG4CyabYwph6j9+2HCBPwTvyBxzywC2puYiLtIG5GCZ/Sl7Fk8kYhfe5Oz/gIi2n/Onpa/Z8gFZ7ppFQ9Dh3p4esrRI19sqbv6ocJpGVXdpqpfu6/3A6uBk4ArgWnubtOA/pWsozH1ht8P4x7Yi//ap/EMHAgjRjA90JEjxKBEcUSjmJ51JsnPRNN0r5cGCV8TeeFTNEj4usT0i6Ve6h9R1cqfRCQB+AzoBGxU1RZuuQA/574vdMwQYAhAfHx891mzZlXoszMyMmjSpEmFjq2rrM1hSJUN87bzpxeuIjMnihgCpJDEb3tk8aezevHfqU9CdjREZnLFw/8g8ZjLGDPmDG75+0wOxi+k0Y5eTHvyBkaNWgXAmDFn0K/fVt5990RGjVpFt26/hLZ9ZRT2P+ciVKbNiYmJy1S1R5EbVbVSD6AJsAy42n3/S6HtP5d2ju7du2tFpaWlVfjYusraHEZ27VKdOFH11FN1LCM0gkwF1QjJ1Luu/konTFB97u1FGjPkdypJf9OYIb/T595epBMmqKamqsbFqY4c6TynpuaXpaY6py/8vrYL259zCSrTZmCpFhNXKzXOXUSigbnADFV9yy3eISKtVHWbiLQCdlbmM4wJN/5F6sz7kj4Nj2+ccxMS0PLU5eT8FIBsJScqk9izN9CzZ08GDvQw6dVx7gXScXkXSIuaUz052RbMMI7KjJYR4F/AalV9OmjTu8AtwHj3+Z1K1dCYcLF+Pf6xaSROvYFAThti5EHSdAGey1qQfPxE1njnEfFF/gXSH7Qb8UuuKfICae74dLvxyBSnMuPczwcGAReJyHL30RcnqF8iImuAi933xtQ7fj+Me/QQ/hHvOHeOnnwy02ev4YhGOxdHiWb60Jfh/ffpecsZ/HfSRUS3dS6QRrf9mk+nXZu3WEbhC6R245EpTYV77qq6EJBiNidV9LzG1Hn79+N/ZjGJj11AIDuaGLmENB3P51F/Y8v5GfCxk3ohMpMtLb8mObkTDz4I8yZ7uOovKZx1tY+v3/IyakjDYnvotl6pKY3NLWNMJfn94PvwMF75FM/yKfDhh0xvcB9HchKBKI6gTO//GAP/eD5jR3xL9K19yVp/PlHtv+Czt8dx32TnPImJ8JerPDz+mIeRI6FbN1+BHnpp49Mt/WKCWXA3pqI2b8b/3FISJ/UhkB1FjFxAmo7mc+5hy3l74POgHnqnbSz5vpHbOx/n9s7HMW+yJy8gF+6hH3NMCzIzrYduKsaCuzFl5P8iB9+MLXgPfYBn+RSSl1/CwmZtOJITRW4PfVyvcbRq15PPvv2uyB564d55cGAv3EO/6qozmDfPeuimYmziMGNKsnkz/Pvf+HuPIvGCw/x9SisSp92Mf3ksPWNXsuDkWIgMgGRCZCYLMjdy/e2NmTfZQ6PF4/Ce34RGi/N76MWtbFRUDn3UqFV2B6mpMOu5GxPE/9Gv+GZsxnvwAzwrX4EffySZB5weOjHk9dB7T6bX77oxseMS7p7XFzadD22+YOJV4/ICdOEeekn586KGMHbr9gteb0223oQTC+6mXvL7wecDb+e9eA6lkvxMNI02fM/9W+8loKcSIwlM1J84GHMtPXsojx5pDN/m59AX7FvNfZ6eJCZ6+HrZOP69zscf2o5j2JX5C0zbCBcTShbcTf2RnU3y/TtptH0d9/+3J4HsSGIklomaxjo6M63ZgAK98/uPu58PZrQh8ZIoJr7jL7KHnpYG86d4GDnU46RZ3CBtI1xMqFlwN2EpeYLSM2EXa374iLmfHmZA1iI6LNnBusNXMq2ZlyM5keQF8Rb38kF/P9ti/8f8V/+U1zvvfefnJF4yGIBhVx7dQy8uzXL11dZDN6Fnwd3UacnJ0LN7DmvSP2Xu/J0MaPEjHb7PYN3KTjzW9Dcc+HkAZMewIPImGh+znPnNn2PbWd8z/+Px+UF8mJ/ExwcTu8nP+4G+ZK8/n8j2X/DwkHF5n1NUD708aRbroZuaZsHd1AkzZ7ZBFdZs/5K587Yy4MTNdFiXybolJ/OYxnNg7zluEA84Qfzw3/n6xCtYuvts0CjIVk4/7wCJb890gvjoo4P44bUemnyVPwb98FoPtCn5QqgFcVNbWXA3tUZyMvTsCWt++Y6572xjwOn76bAtmiVLlN/v+oAr3vkDB/Z0huzu+UF819/5ul1flmb3zA/inbeT+Oq/ueObX1k6MP8i6B23ngAUHcTT1joBe95sD4mJHtK8+QHcLoSausiCu6k2ucE6OACmpcFTT8ED9ytrtn7J3Pd2MKD9LjpsiWTdojY89kwTDuzuDNlnFAjgifg4vd1vWJrdIz+I/3YjiTMf4Y6tsSy9LSiID+sACQl0WA+N264gPmEPO9Jb0qH5mXm98MJBvKQ8uc20aOoiC+6m0nKDeGysO7zQC4cPw7rVAZ4aH8Hga//Hyo376MQ+pqcO5uG4l7nihl4c2NulYColrxceFMDb/kjiA335Lqs3dxzXnaVDg4L4XztC0pl0SCs+iM+ffWaBMebludhpAdzUZRbcTQHF9bZz75Ts2RPW7FvB3Pk7GHC20iGjJes+aMzYR9txMCuCrJwIoiSbRjmHmMdVNDmuHU9P/WdeEB9+zF0M3zqNme2SC6ZSTvmJxORbuGN3C5Y+HBTAH+kB/c9kr89HB/VaEDemjCy4h4kSUyAPFF++Zt8K5n6whwGXtaRDkzNYt+IIT42PYfD1H7Jy0346NTjE9PcHMtvzLOzZwxVPX+fmvU9nwfTcHvdf2dbsauYHhoJGkYlyYbPXScz8kvHHJcHumLwgvtLTBx75M3f8mM3SPwQF8Ye6QP+ie+EA33zTgrFjLYgbU1YW3EOovAE5t/ccHd2iwG3paWmwbp1z3OD7V7Dypz10an8s0yd24uE79zLw6uZOsN6SQaeGAaa/dw0Pd3mfK25K4MCujk6gnhqg8TFfMn/X353e9isFe9uJqdMAOL3diQXTJvHfkXhxK2ZHr4I384P1SaNj4N4DdEpeyYJH8ss7ndcRep5Jh4yyp1Jmz4YffmhqQdyY8ihucdWafIRygezchYaDpaaqXnZZ+conTCj/uYYMcRYvHj7+O+09JE2Hj/9O4+JUJ01SjYvL0eGjv9Leg9/X4X9N0bjmAU0d69fUv32sxzTcr8MHvKq9vS/p8KRJGtdgn6b2eECHn/mQEnVAkUwl6oAOP+4WVdDhx91SZHmPdslOGapIQHu0S1aNi9PeHZ8oUN77kqmq06erfvyxvvTCh+65AkrUAX1p3neq6izizC2/Uy76m3KLs4hz7uLMRbfv6EWchwwp/ru1hZPrB2tz+VDCAtkhD+xageAeHERzv5jqCbBlKB/7jca1zNbU/+7W1De2aNwxmTr8no+09zVzdfgd8zSu6SGddJ1f4xof1OFXv6q9E1/W4X3+oXGxv2qqd4wOP2dMwcDbbqhqo0bFBuSSgnXvdiMKBuVTHlFNSNDeZ44tWJ70qurzz+tLf/9XwUA99xtVddoWXD58/HcFvrPGJ3+nJ1+Upo1P/k5TU/OD83NvL9Kxn43V595eVGKwLukPZHHsl75+sDaXT0iCO9AH+BFYC4woad/yBve8HuHjy/Ts/rN0+N+/0LhjszT1tQ2a+spajWsR0OF3f+AE2D+9q3HNDuukQd9oXNNDOvzmGdr7smk6/Np/a1zjA5p6y2s6/PJ/FgyWF45XvfrqowPvb+9TbddOh7cbWu7AW+aA3G6EU55QqLzLeNU+fVSvu05/12NSwW1XvKn6/vs6/J6PigzKxQXrkgJ14T9qwdvK0+OuKvZLXz9Ym8unpOBeLTl3EYkEXgAuATYDS0TkXVVdVRXnT0x0cstPP3IaZHfmq/cCDD/mdhJvdfLCg4+7haenBOeM72T469PYctwtPD2rUC552jTGtxsB2UEX/Tb8Ap+9xcp2pxYsP9wANmxgZbvmBcubtYaIeGjUiJWNOhS8gHhyV/i9sHLd2bAwqLznxXDzpXRKa8iCfwflpAdfCQ/+nU7P/8SCkUHlN/SFhx4C4Phh78PyoG3ndyIt9kymz4DhT7g595NbMn3imZwUDdMnnllk+bhxJV2gPDPv+7787Px8v+W9jakbquuC6tnAWlX9CUBEZgFXAlUS3AFW/rQHsk/PD5YtT4ZjfwsxMazktwUDbLvOcOEAVm7pBl8GlXf0QK82dFrZlgVbgoLl77rCFf+l08fRLJgaVH7dpfCnIXT69w4WjAsqv+M6eOgJADpNWFHwAuJVSfDQvU65P6j8wi6knXAm0+cVDsjnclILmD6p81EB+fKznbYveONihj/xY4FtGcUE5aeeKq3cKatooLYgbkztJE7PvopPKnIN0EdV/+i+HwSco6rDgvYZAgwBiI+P7z5r1qxyfcY/5xzkvy95ITsaIjO59k4fd13TqMRtRZV7TjmRMWPOoPfNn7B+x2Hax8ey4I2LufHGjbz5Ztsyl48a5fzdKs+5LrhgNxddtJNu3X7Ja9c337Rg1qw2XH/9pqPKf/ihKQDt2u3kvPOOHLXthhs2les7rEsyMjJo0qRJqKtRo6zN9UNl2pyYmLhMVXsUubG4fE1lHsA1wKtB7wcBzxe3f4Vz7uO/0x7XzC0yL1zWC6TlvehXlaNlKpqTtrxk/WBtrh/qVM4d2AK0CXrf2i2rEvkTOZ2Jz7cHr/fMIvLCZUtPlCcNUd7yih5jjDGVVV3BfQnQQUTa4wT164Ebq+rk5Z3IyQKsMaa+qZbgrqpZIjIM+AiIBKaq6vfV8VnGGGOOVm3TD6jq+8D71XV+Y4wxxYsIdQWMMcZUPQvuxhgThiy4G2NMGKqWm5jKXQmRXcCGCh4eB+yuwurUBdbm+sHaXD9Ups3tVPW4ojbUiuBeGSKyVIu7QytMWZvrB2tz/VBdbba0jDHGhCEL7sYYE4bCIbi/HOoKhIC1uX6wNtcP1dLmOp9zN8YYc7Rw6LkbY4wpxIK7McaEoTod3EWkj4j8KCJrRWREqOtT3URkqojsFJGVoa5LTRGRNiKSJiKrROR7Ebkn1HWqbiISKyJfici3bpvHhLpONUFEIkXkGxF5L9R1qQkiki4iK0RkuYgsrfLz19Wcu7tO6/8RtE4rcINW0TqttZGIXAhkANNVtVOo61MTRKQV0EpVvxaRpsAyoH+Y/5wFaKyqGSISDSwE7lHVxSGuWrUSkeFAD6CZql4e6vpUNxFJB3qoarXctFWXe+5567SqagDIXac1bKnqZ8DeUNejJqnqNlX92n29H1gNnBTaWlUvd5GdDPdttPuom72wMhKR1sDvgVdDXZdwUZeD+0lA8KKhmwnzX/r6TkQSgG7AlyGuSrVzUxTLgZ3Ax6oa7m1+FngQyAlxPWqSAgtEZJm7pnSVqsvB3dQjItIEmAvcq6q/hro+1U1Vs1W1K84SlWeLSNim4UTkcmCnqi4LdV1qWC9VPQu4DPizm3atMnU5uFfrOq2m9nDzznOBGar6VqjrU5NU9RcgDegT4qpUp/OBfm4OehZwkYi8EdoqVT9V3eI+7wTm4aSaq0xdDu5567SKSAzOOq3vhrhOpoq5Fxf/BaxW1adDXZ+aICLHiUgL93VDnEEDP4S0UtVIVR9W1daqmoDze5yqqjeHuFrVSkQauwMEEJHGQG+gSkfB1dngrqpZQO46rauB2eG+TquIzAT8wG9FZLOI3B7qOtWA84FBOL255e6jb6grVc1aAWki8h1OJ+ZjVa0XwwPrkXhgoYh8C3wF/E9VP6zKD6izQyGNMcYUr8723I0xxhTPgrsxxoQhC+7GGBOGLLgbY0wYsuBujDFhyIK7McaEIQvuxhgThv4fr7qz1cn1cWcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(x,y1,'r-',linewidth=2, label='Actual')\n", "plt.plot(x,y2,'bx',linewidth=2, label='Euler')\n", "plt.plot(x,y3,'g.',linewidth=2, label='Heun')\n", "plt.plot(x,y4.y[0],'b.',linewidth=2, label='RK45')\n", "\n", "plt.grid()\n", "plt.title(\"Comparison of ODE methods\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 8, "id": "medieval-concept", "metadata": {}, "outputs": [], "source": [ "y2err = abs(y2-y1)\n", "y3err = abs(y3-y1)\n", "y4err = abs(y4.y[0]-y1)" ] }, { "cell_type": "code", "execution_count": 9, "id": "adjustable-vinyl", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAm2ElEQVR4nO3dfZRU1Znv8e9Dg7yIQgTFF16ajHEUUUDQK6POpX2XKGGNmUTTYfCK6RmiWcZEHQ1XEXMxSEw0aq4uEl0S7RFZOokvYxIVuuNoNAKCimiub92IRlAUFFHk5bl/nFNNdfU51dXdVV11qn6ftWp11z6nztm7lad3P3ufvc3dERGR5OlR7AqIiEjnKICLiCSUAriISEIpgIuIJJQCuIhIQimAi4gklAK4FJSZ/d7Mphe7HlHMzM3s4Dxdq8nMTu7E52aa2Xoz22Jmg/JRF6kcCuBlLAwqn4XBIfW6tTvr4O5nuPvC7rxnFDNrNLMLil2PdGbWC/g5cKq793f3jQW+X6OZTSrkPaR79Sx2BaTgznL3J9o7ycx6uvuOjLIqd9+Z6406er4wBOgDvNzRD5qZAebuu/JeK0kM9cArlJmdZ2ZPm9mNZrYRuMbM7jKz28zsUTP7FKgxs8PCntsmM3vZzKakXaPN+RH3aen5hvd8ysxuMLOPzOwtMzsjSx2bzOwyM3vRzD41szvMbEiYlvnEzJ4wsy+lnX+smf05rOsLqd6mmc0FTgBujfgr5GQzey38zC/DwIiZ9TCz/21mzWa2wcx+Y2YD0u41LTy20cxmZdT7GDNbbmYfh+mRn0e07RDgr+HbTWa2NCz/BzNbZmabw6//kPGznGtmTwNbgS9nXPMyM3sgo+xmM/tF3M84GzPb38y2pqd2zOwoM3s//OtBis3d9SrTF9AEnBxz7DxgB/A9gr/E+gJ3AZuB4wh+ue8FvA78CNgDOBH4BPj78BqZ5/eJuE8jcEHaPbcD3wGqgJnAuwQ9ybj6P0vQUz0I2AA8D4wj6LkuBWaH5x4EbAQmh3U5JXy/b2Y90q7vwCPAQGA48D5wenjs/LDtXwb6A/8J3B0eGwVsAf4R6E2QBtmR+lkDzwDTwu/7A8fGtK86rEPP8P0+wEfAtPC/ybnh+0FpbVgLHB4e75VxvQOAT4GB4fue4c9sfMS9vwVsyvIaHp73KDAz7XM3ArcU+/9tvYKXeuDl73dh7zL1+k7asXfd/RZ33+Hun4VlD7r70x78aT6WIADNc/cv3H0pQcA7N+0aLee7++c51KfZ3X/lQaplIUHQGZLl/Fvcfb27vwP8N/AXd18Z3uu3BMEc4NvAo+7+aFiXx4HlBAE9m3nuvsnd1wINYZsBaoGfu/ub7r4FuBI4x8x6Al8HHnH3J919G3AVkJ7K2A4cbGaD3X2Luz+bw88F4KvAa+5+d/jf5F7gVeCstHPucveXw+Pb0z/s7n8DngT+OSw6HfjA3Vdk3sjd/8PdB2Z5rQ1PXUjws8XMqgj+29+dY3ukwBTAy9/UjH+Yv0o79nbE+ellBwJve+s8azNBbzfbNbJ5L/WNu28Nv+2f5fz1ad9/FvE+9dkRwD+n/7ICjif4BZFTfQjSEqnrHUjQ1pRmgh7tkPBYS7vd/VOC3n7KDOAQ4NUwDXJmO3VIybxn6r4d+Xm3BNzwa1eD7YPAKDMbSfBXzWZ3f66L15Q8UQCvbFFLUaaXvQsMM7P0/0+GA++0c41ieJsgxZH+y2pPd58XHu9oPd8l+KWQMpwgTbIe+BswLHXAzPoBLXlid3/N3c8F9gOuB+43sz07cc/UfTvy8/4dcKSZjQbOBOqjTjKz2ozZSZmv4WFbPgcWE/wymIZ63yVFAVyy+QtBr/RyM+sVDgqeBSwqZqVi3AOcZWanmVmVmfUxs0lmNjQ8vp6MQb923AtcYmYjzaw/cB1wnwczde4HzjSz481sD+Ba0v4tmdm3zWzf8C+XTWFxLrNFHgUOMbNvmVlPM/smQb79kVwrHQbc+4H/AJ5LS4VknlfvwdTFuFf6535DMH4xBQXwkqIAXv4ezuhZ/TbXD7r7FwQB+wzgA+D/Av/i7q8WqK6d5u5vA18jGHB9n6BHfhm7/x//BfD1cPbLzTlc8k6CYPUk8BbwOcGAL+7+MnAhQZD8G8FA47q0z54OvGxmW8L7npM2xpCtDRsJes0/JEjJXA6c6e4f5FDfdAuBI8hTsHX3pwl+AT3v7pkpHikicy+Vv4BFJB/C9MerwP7u/nGerrkU+A93/3U+rif5oQAuUkbC8YqfA3u7+/l5uubRwOPAMHf/JB/XlPzQk5giZSIcKF1PMHPl9DxdcyEwFbhYwbv0qAcuIpJQGsQUEUmobk2hDB482Kurqzv12U8//ZQ998xlKm35UJsrg9pc/rra3hUrVnzg7vtmlndrAK+urmb58uWd+mxjYyOTJk3Kb4VKnNpcGdTm8tfV9ppZ5PRNpVBERBJKAVxEJKEUwEVEEqro88C3b9/OunXr+Pzz7CuRDhgwgFdeeaWbalVYffr0YejQofTqpTXxRaTzih7A161bx1577UV1dTXhZiiRPvnkE/baa69urFlhuDsbN25k3bp1jBw5stjVEZEEK3oK5fPPP2fQoEFZg3c5MTMGDRrU7l8cIpJs8+dDQ0PrsoaGoDxfih7AgYoJ3imV1l6RSnT00fCNb+wO4g0Nwfujj87fPYqeQhERKUc1NbB4cRC0zzijmt//Pnhf02br784riR54sVVVVTF27NiW17x587Kef9ddd3HRRRd1U+1EJKlqamDmTLj77mpmzsxv8IaE9cDnzw/+/Ej/ITQ0wLJlcPnlnb9u3759WbVqVZfrF2fHjh307JmoH7WI5EFDA9x2G0yb1sRtt1VTU1PBPfDuyCmlq66u5oMPgs1Qli9fHvko7Pvvv8/ZZ5/N0UcfzdFHH83TTz8NwDXXXMO0adM47rjjmDZtWmEqKCIlKxWfFi+G889vakmnZA5sdkWiuoXpOaWZM4PfbPnIKX322WeMHTu25f2VV17JN7/5zZw+e/HFF3PJJZdw/PHHs3btWk477bSW+epr1qzhqaeeom/fvl2roIgkzrJlu+NTY+Pu+LVsWf564YkK4LA7p/TjH8NVV+XnB9GVFMoTTzzBmjVrWt5//PHHbNmyBYApU6YoeItUqKi0br5TKDkHcDOrApYD77j7mWY2kmB38kHACmBauAluQaVySlddFXzN9w8kXc+ePdm1K9hMPG7e9q5du3j22Wfp06dPm2OVtFymiHS/juTALwbSn2W/HrjR3Q8m2JV7Rj4rFiU9p3TttRQkp5SuurqaFStWAPDAAw9EnnPqqadyyy23tLwv5GCoiEi6nAK4mQ0Fvgr8OnxvwInA/eEpqX3zCio9pwStc0pdkcqBp15XXHEFALNnz+biiy9mwoQJVFVVRX725ptvZvny5Rx55JGMGjWK22+/vWuVERHJUU57YprZ/cBPgL2AS4HzgGfD3jdmNgz4vbuPjvhsHVAHMGTIkPGLFi1qdXzAgAEcfPDB7dZh586dsUE0iV5//XU2b96c9ZwtW7bQv3//bqpRaVCbK0Oltbmr7a2pqVnh7hMyy9vNgZvZmcAGd19hZpM6emN3XwAsAJgwYYJnTsV75ZVXclqkqlwWs0rp06cP48aNy3pOpe1aAmpzpai0NheqvbkMYh4HTDGzyUAfYG/gF8BAM+vp7juAocA7ea+diIjEajcH7u5XuvtQd68GzgGWunst0AB8PTxtOvBgwWopIiJtdOVJzH8HfmBmrxNMJbwjP1USEZFcdOhBHndvBBrD798Ejsl/lUREJBeJWgtFRER2UwCHNtN7tFysiCRB4gJ4/Uv1VN9UTY85Pai+qZr6l+qLXSURkaJIVACvf6meuofraN7cjOM0b26m7uG6ggbxbMvF3nDDDS3njR49mqamJpqamjjssMP4zne+w+GHH86pp57KZ599VrD6iUjlSlQAn7VkFlu3b21VtnX7VmYtmdWl62Y+Sn/11Ve3HEstF7ts2TIeeOABLrjggnav99prr3HhhRfy8ssvM3DgwNh1VEREuiJRy8mu3by2Q+W5ylxO9q677mL58uVA9uVi44wcObJlffHx48fT1NTUpfqJiERJVAAfPmA4zZubI8sLJW652PSlZqH1crO9e/du+b6qqkopFJEyVqitHnORqBTK3JPm0q9Xv1Zl/Xr1Y+5Jcwt2z7jlYqurq3n++ecBeP7553nrrbcKVgcRKV3dvdVjukQF8Nojallw1gJGDBiBYYwYMIIFZy2g9ojagt0zbrnYs88+mw8//JDDDz+cW2+9lUMOOaRgdRCR0pW+1ePVV+/es6BQG82kS1QKBYIgnu+AnZnTPu+88zjvvPMAGDx4MPfdd1+bz/Tt25fHHnss8nqrV69u+f7SSy/NX0VFpCQVYqvHXCSqBy4iUooyt3os1C5hmRTARUS6oLu3ekynAC4i0gWF2uoxF4nLgYuIlJKoqYI1Nd2TB1cPXEQkoRTARUQSSgGc4GnJsWPHMnr0aM466yw2bdoEQFNTE6NHj24571e/+hXjx4/no48+ain72c9+hpnxwQcfAMHmpQMGDGhZV+Xaa6/t1raISOVQAGf3WiirV69mn3324Ze//GWbc+6++25uueUW/vjHP/KlL30JgLfffpvHHnuM4cNbP8p/wgknsGrVKlatWtVqYSwRkXxKZAB/5hn4yU+Cr/k2ceJE3nnnnVZlixcvZt68eTz22GMMHjy4pfySSy5h/vz5mFn+KyIi0o7EzUJ55hk46ST44gvYYw9YsgQmTszPtXfu3MmSJUuYMWNGS1lzczMXXXQRK1euZP/9928pf/DBBznooIMYM2ZMRB2fYcyYMRx44IHccMMNHH744fmpoIhImsT1wBsbg+C9c2fwtbGx69dMrQe+//77s379ek455ZSWY/vuuy/Dhw9n8eLFLWVbt27luuuui8xvH3XUUTQ3N/PCCy/wve99j6lTp3a9giIiERIXwCdNCnreVVXB10mTun7NVA68ubkZd2+VA+/Xrx+PPvoot99+O/X1wc4/b7zxBm+99RZjxoyhurqadevWcdRRR/Hee++x9957t+yxOXnyZLZv394ywCkikk+JS6FMnBikTRobg+Cdr/QJBMH65ptvZurUqXz3u99tKd9vv/34wx/+wKRJkxg8eDCnnXYaGzZsaDleXV3N8uXLGTx4MO+99x5DhgzBzHjuuefYtWsXgwYNyl8lRURCiQvgEATtfAbudOPGjePII4/k3nvv5YQTTmgpHzlyJA899BCTJ0/mt7/9Lcccc0zk5++//35uu+02evbsSd++fVm0aJEGOUWkIBIZwPMtcznZhx9+uOX79KVhx4wZ02aGCtBqy7SLLrqIiy66KP+VFBHJkLgcuIhIMcyf33aFwYaGoLxYFMBFRHJQzK3T4pRECsXdKypP7O7FroKIdFD61mkzZwYbN3TX1mlxit4D79OnDxs3bqyYoObubNy4sc0u9yJS+tK3Tps5s7jBG0qgBz506FDWrVvH+++/n/W8zz//vGyCXp8+fRg6dGixqyEiHZS5dVp3rfsdp+gBvFevXowcObLd8xobGxk3blw31EhEpK30rdNSgbs7d6CPUvQUiohIEhRz67Q4Re+Bi4gkQTG3ToujHriISEIpgIuIJJQCuIhIQimAi4gklAK4iEhCKYCLiCRUuwHczPqY2XNm9oKZvWxmc8LykWb2FzN73czuM7M9Cl9dERFJyaUHvg040d3HAGOB083sWOB64EZ3Pxj4CJgRfwkREcm3dgO4B1I7HvQKXw6cCNwfli8EphaigiIiEs1yWQXQzKqAFcDBwC+BnwLPhr1vzGwY8Ht3Hx3x2TqgDmDIkCHjFy1a1KmKbtmypWWz4EqhNlcGtbn8dbW9NTU1K9x9QpsD7p7zCxgINADHA6+nlQ8DVrf3+fHjx3tnNTQ0dPqzSaU2Vwa1ubRcf7370qWty5YuDco7q6vtBZZ7REzt0CwUd98UBvCJwEAzS62lMhRou1mkiEjClOLOO3FymYWyr5kNDL/vC5wCvEIQyL8enjYdeLBAdRQR6TbpO+9cfXXxl4zNJpce+AFAg5m9CCwDHnf3R4B/B35gZq8Dg4A7CldNEZHuU2o778RpdzlZd38RaLOTgru/CRxTiEqJiBRTqe28E0dPYoqIpEnfeefaa3enU1I58VKiAC4ikqYUd96Jox15RETSlOLOO3HUAxcRSSgFcBGRhFIAFxFJKAVwEZGEUgAXEUkoBXARkYRSABcRSSgFcBGRhFIAF5GKNH9+28fjGxqC8qRQABeRipSkdb/j6FF6EalI6et+z5wZrDpYqut+x1EPXEQqVlLW/Y6jAC4iFStz3e9SXDI2GwVwEalISVr3O44CuIhUpCSt+x1Hg5giUpGStO53HPXARUQSSgFcRCShFMBFRBJKAVxEJKEUwEVEEkoBXEQkoRTARaRslcOKg9kogItI2SqHFQez0YM8IlK2ymHFwWzUAxeRspb0FQezUQAXkbKW9BUHs1EAF5GyVQ4rDmajAC4iZascVhzMRoOYIlK2ymHFwWzUAxcRSSgFcBGRhFIAFxFJKAVwEZGEUgAXkcQr9zVP4rQbwM1smJk1mNkaM3vZzC4Oy/cxs8fN7LXw65cKX10RkbbKfc2TOLn0wHcAP3T3UcCxwIVmNgq4Alji7l8BloTvRUS6XfqaJ1dfvfvhnXKZLhin3QDu7n9z9+fD7z8BXgEOAr4GLAxPWwhMLVAdRUTaVc5rnsQxd8/9ZLNq4ElgNLDW3QeG5QZ8lHqf8Zk6oA5gyJAh4xctWtSpim7ZsoX+/ft36rNJpTZXBrU5P1auHMicOaOYMuVdHnroQGbPXsO4cZvyeo/O6mp7a2pqVrj7hDYH3D2nF9AfWAH8U/h+U8bxj9q7xvjx472zGhoaOv3ZpFKbK4Pa3HVLl7oPHhx8jXpfbF1tL7DcI2JqTrNQzKwX8ABQ7+7/GRavN7MDwuMHABs6/etFRKQLyn3NkzjtroUSpkfuAF5x95+nHXoImA7MC78+WJAaioi0o9zXPImTy2JWxwHTgJfMbFVY9iOCwL3YzGYAzcA3ClJDERGJ1G4Ad/enAIs5fFJ+qyMiIrnSk5giIgmlAC4iiVGpj8zHUQAXkcSo1Efm42hHHhFJjPRH5mfODDYproRH5uOoBy4iiVKJj8zHUQAXkURpaAh63lddFXwtlx3mO0MBXEQSI5XzXrwYrr12dzqlUoO4AriIJEalPjIfR4OYIpIYlfrIfBz1wEVEEkoBXEQkoRTARaTk6InL3CiAi0jJ0ROXudEgpoiUHD1xmRv1wEWkJJXDE5f1L9VTfVM1J/7pRKpvqqb+pfq8Xl8BXERKUtKfuKx/qZ66h+to3tyM4zRvbqbu4bq8BnEFcBEpOeXwxOWsJbPYun1rq7Kt27cya8msvN1DAVxESk45PHG5dvPaDpV3hgYxRaTklMMTl8MHDKd5c3Nkeb6oBy4iUgBzT5pLv179WpX169WPuSfNzds9FMBFpGjK+YGd2iNqWXDWAkYMGIFhjBgwggVnLaD2iNq83UMBXESKptwf2Kk9opam7zex9H8upen7TXkN3qAALiJFlP7AztVX7555kqRcN+ye791jTo+CzPeOo0FMESmq9Ad2rroqmcG77uG6limDqfneQN573JnUAxeRokr6AzvdMd87jgK4iBRNOTyw0x3zveMogItI0ZTDAztx87rzOd87jgK4iBRU1FTBlSsHMn9+8MBOZs67pib6QZ5S1R3zveMogItIQUVNFZwzZ1RZTRUs9HzvOJqFIiIFFbW29+zZa6ipGVvsquVN7RG13RKwM6kHLiIFl7m297hxm4pdpU4p1nzvOArgIlJwmVMFV64cWOwqdVh3rO/dUQrgIlJQUVMF58wZlaipglDc+d5xFMBFpKCipgrOnr0mUVMFobjzveMogItIXsStLAhtpwqOG7cpUVMFobjzveMogItIXpT7yoLFnO8dRwFcRPKiXFYWhOjZJsWc7x1H88BFJG+SvrIgtL+6YDEDdqZ2e+BmdqeZbTCz1Wll+5jZ42b2Wvj1S4WtpogkQdJXFoTSnG0SJ5cUyl3A6RllVwBL3P0rwJLwvYhUsHJYWRBKc7ZJnHYDuLs/CXyYUfw1YGH4/UJgan6rJSKlKm62yU9/mvyVBaE0Z5vE6ewg5hB3/1v4/XvAkDzVR0RKXNxsk8suS/7KglCas03imLu3f5JZNfCIu48O329y94Fpxz9y98g8uJnVAXUAQ4YMGb9o0aJOVXTLli3079+/U59NKrW5MiSxzStXDmTOnFFMmfIuDz10ILNnr+nQ+ial0uYn1j/Br9/6NRu2bWC/3vtxwcgLOHnIybHlndXV9tbU1Kxw9wltDrh7uy+gGlid9v6vwAHh9wcAf83lOuPHj/fOamho6PRnk0ptrgxJbfNVV7lD8LWjSqHN97x4j/eb28+5hpZXv7n9/J4X78n7vbraXmC5R8TUzqZQHgKmh99PBx7s5HVEJIE026Q05DKN8F7gGeDvzWydmc0A5gGnmNlrwMnhexGpAJptUjrafZDH3c+NOXRSnusiIiVk/vxgwDJ9YLK92SZJenBn+IDhNG9ujixPCj1KLyKRymW2SdwmDEmabRJHj9KLSKSordCStrZJe4/FQ5ALX7t5LcMHDGfuSXNL6lH59iiAi0ispK9tkm2gMrWuSZICdialUEQkVtJnm5TDQGU2CuAiFSzusfhUedJnmyTpsfjOUAAXqWDZNmGI2gqtlNc2iRqsLIeBymwUwEUqWLZNGC6/PDmzTeJ2jAdKbhOGfNIgpkiFS/pAJWQfrGz6flPZBOxM6oGLVID2ct1JHqiE8h+sjKMALlIB4nLdPXsmb6AyKtdd7oOVcRTARSpAXK57x47kDVRG5bonf2VyWQ9WxlEAF6kQ6bnumTOTN1AJ8bnuR197tKwHK+NoEFOkjMQtQLVsWVCenuuuqSntAcv6l+rbPOaeLded9KcqO0M9cJEyUi657rhUyT5994k8v9xz3XEUwEXKSLnkuuNSJUBF5rrjKICLJFC2aYFJy3VHzSqJS5V8+NmHFZnrjqMcuEgCpVIlqV51+rolmfO6SznXHbfc6z5992HjZxvbnD98wPCKzHXHUQ9cJIHiUiWQrFy3UiVdowAuUsI6miop1QWoMtMkT6x/Aoh/UlKpktwohSJSwjqaKonKaRc7hRKVJrnhkxs47KXDsu5LqVRJ+9QDFymye+8dlrWXnfRUSVSaZNuubcxaMqvsl3stNAVwkSI79NBPYtfkhmSnSrLNKEk9fKNUSecphSJSZOPGbcq6eXCSUyXtzSgBlCrpAvXARbpJ3IDkvfcOi+xlp46XYqokqqfdkRklvXv0VpokDxTARbpJ3GPuhx76Seya3KWYKol7zD1qMBKiZ5Rcesil6nXngVIoInmWbUGpqFTJqlWtZ5rU1Ox+X+xUSdSCUnE97SqrYqfvbHONqBkljY2Nha56RVAPXCTPsm0UHJUqefXVvYray45Kh6TKO9LT3uk7NaOkmymAi3RSXE47vaeduVFwVKrk3HPfLto6JXFBOltOu8qqIq+VmkGiGSXdRwFcpB1xgfqNNzrW044bkFy5cmC3tKMjA4+ptEmUbD3t2iNqafp+E7tm7yrrzYRLhQK4SDviUiLnnNOxnnbcgOSrr+6V1/pGBeqOpkNSOe8o6mmXDg1iipB94PHyy6MHH1PnpnraV13VtqedOSgZlSoxexv4uw7XOWqAEYici923Z98ODzzOPWluq2tB6562AnbxqQcuFaUz6RCITomkzsu1p53PQcm4HvXFv784MlBHPUgD7adD1NMubQrgUpY6GqizpUNS52UG6ricdmZPHtoflIybCRJ3LC53HReo47SXDlFOu7QphSKJFpf6eOMN+OlPo1fxO+ec3NMhqetFpUT+6Z/ie9px87Sj0h6vrH+FG/98Y5u0R0pUSiQzeLdnUN9BfLbjM6VDyowCuJSM+fOhV6+BTJq0u6yhIQjEl10Wv9N63HKrnQnUUTvZZEuJRPW0a2o6lp/u6T3ZuiN6Jkjq+8xjcbnruED9izN+AdCmTgrcyaYALt0uW6/5vvtGMXZs62B85ZXxQTp9udWuBGrI/jRkkMZoHfwuv7w2LwOJceKm8cHu3HVHA7UCdnlRAJecZJulAdHH4nrO2dIbhx66hm98Y2ybYDxuXHyQTh9gTA/UsxbVM+/TWTB7LXM/Gc7ORXM5eUgtU2bV0++Hs/g/29ay7w+HM2XWXL51RC11t9Tzv16Yxdong8BXd8tcli2r5d3B0avsPb32aRa+sDDnQN3RtEdqGl/UVL8RA0a0PNauQF25NIhZ5uIG8yZP7lh5tlkacfOkTz654wOGGw+8n+0XVvPjHj3YfmE17w4OBvPeHVwfWQ5BoJ77aTXM7sHcT6uZtaieWYvquW51Hbv2bgZzdu3dzHWr67jw0e+y/Yw6NmwLZm9s2NbM9jPqeHHod7npzdazOm56s46DzogfMFywYkFeBhL3rto7diZItg0PNMAo6oHnUVwvNa4n2l65WW7np3rBm4fXU//e7h5Z7f5z+eCNWn78u3r6TZnF+9vWsm/v4Wx9aC5zvh72RHMsf2huLYNr6jn5kVns+tNaenwynCtumUtNTRA06m5pe+wH59Sy8cDoz9S/VM/2C2fxY9Yy4MLhvDt4LvUvwfxXf8Z2+zxoD83M+F3Q073z+YVss62tygHWvAzXra6DvYNjqUDdu0df6JXR4+21lVd7LcB3tc4db9u1lWW2gJ3bW5fn8kRiR8Tlp7/3d9/jsFGHZc1PK3ctUczdO/9hs9OBXwBVwK/dfV628ydMmODLly/P+frz50cHpZW/qWXcv+RePmBt8D97Pq6V7R5HH01sUJx9f8fLe51+GZv8vXbPf2huLU+sD3qcrYLW9n5MGTadP25YyLZdu8t79+jH+UdND4JiB8rT0wUQBJ8FZy0AiHzgY/qY6M9MHxN9j949+vLxjra91x5UsYu2wXK/3iP44gvY5NFPE+aDYbH7NnZ0IDH1s8oMxgdtPIhJ6SO3FaCxsbGi2tzV9prZCnef0Ka8swHczKqA/wecAqwDlgHnuvuauM90NICn/gzmvSOhaRJUN8L+LzKxz3Se+XxhzuU/Gh38w8nHtbLdY9ThMON3dWxr3n2s94gXdwfFApXfMXUBs5bMCoLM28furtewZ3cHmQKVjxgwAqBD924JyBnlOGDkXo5hgOMRx0J5aF/LE4lv7v659/vyi7t/QWWUpwL1pXc8wHurD2X/0a9yw4yzW3rNzzwDjY0waRJMnLj7H3dmebq4Y4UuL9Q9tm1rHdDKrX2Z5Xfe+Sbnn//lNvfOVVwAx9079QImAn9Me38lcGW2z4wfP947YsSNI5wZxzo9P3Vse/B1xrFeNaeqQ+UjbhyRt2sV8x7Z7m3XWOQxrqGg5XaNdfzeszt2j878zAddP8h7/+ukVuW9/3WSz3xkZofK73nxHnd3n333H9x6bXVsu1uvrT777j9kLf/zn9379nWvqgq+/vnPHlve0NAQe35Hr5XP8kLe49ZbV5R1+zLLe/TY1ebeHQEs94iY2pUc+EHA22nv1wH/I+I3Rx1QBzBkyJAOLeS+dvNaaDoXdu4B3hN2OjRNYuewZ4MeT47la4ddH1wwD9cq5j2y3Xu/3vuxPuJYj2HPsSui3IYuw/NQvucBQX54S8QxDloWWd+4a+1dvYZPm05mZ1p5VfPJnHnMEB757yPalH/7lOCvx3n/ZW2O/dspJ/H0yn+kMa184uYr+caee/DK5i9yLj9o4x40Njbyzp++gu3sg7thO3vwzp++QuPQ+PL6+uFs2zaSXbuMbdt2ceedTWzbtjay/Gtf20J9/ZuR5wMdulY+ywt57+ee68vhhzeWbftyuXdeREX1XF7A1wny3qn304Bbs32maz3wL2J6ZO2Xt+2pdf5a2e6x30+i79HjmujP5Kt8v5+M8B/de48zfVLrY9Mn+WGXzvReF7Qu73XBJB/2bx0rP/ba9B7qF616qPe8eE/ksda92vbL73nxHp9+7e1hj/aLmJ5u6/Jsx5LSU1MPvPza17YHvrMgPfCu5MAnAte4+2nh+yvDXwg/iftMuefA33gTfrezbQ58v3ens+HAtrnrfJVPrVrAUT1r2Ty8njuf3Z13Pf/YsxmwtpaDzqiPzMfWv5SfciBv12psbKR370kll8cs5D2UAy+/9nVXDrwrAbwnwSDmScA7BIOY33L3l+M+U+6zUFJP7EVN+epM+Q//64ds2LYhp/PLRaXNTgC1uRIUahZKp1MoYeCfTBDE3wBmtXd+R1Mo6RoaGjr92aRSmyuD2lz+utpeCjCIibs/CjzalWuIiEjn6FF6EZGEUgAXEUkoBXARkYRSABcRSaguLWbV4ZuZvQ80d/Ljg4EP8lidJFCbK4PaXP662t4R7r5vZmG3BvCuMLPlHjUPsoypzZVBbS5/hWqvUigiIgmlAC4iklBJCuALil2BIlCbK4PaXP4K0t7E5MBFRKS1JPXARUQkjQK4iEhCJSKAm9npZvZXM3vdzK4odn0KzczuNLMNZra62HXpDmY2zMwazGyNmb1sZhcXu06FZmZ9zOw5M3shbPOcYtepu5hZlZmtNLNHil2X7mBmTWb2kpmtMrPc19PO5dqlngPvzObJSWdm/whsAX7j7qOLXZ9CM7MDgAPc/Xkz2wtYAUwt8//GBuzp7lvMrBfwFHCxuz/bzkcTz8x+AEwA9nb3M4tdn0IzsyZggrvn/cGlJPTAjwFed/c33f0LYBHwtSLXqaDc/Ungw2LXo7u4+9/c/fnw+0+AVwj2XC1b4TLPW8K3vcJXafem8sDMhgJfBX5d7LqUgyQE8KjNk8v6H3clM7NqYBzwlyJXpeDCVMIqYAPwuLuXfZuBm4DLgV1Frkd3cuAxM1sRbvKeN0kI4FIhzKw/8ADwfXf/uNj1KTR33+nuY4GhwDFmVtbpMjM7E9jg7iuKXZdudry7HwWcAVwYpkjzIgkB/B1gWNr7oWGZlJEwD/wAUO/u/1ns+nQnd98ENACnF7kqhXYcMCXMCS8CTjSze4pbpcJz93fCrxuA3xKkhfMiCQF8GfAVMxtpZnsA5wAPFblOkkfhgN4dwCvu/vNi16c7mNm+ZjYw/L4vwSD9q0WtVIG5+5XuPtTdqwn+HS91928XuVoFZWZ7hgPzmNmewKlA3maXlXwAd/cdwEXAHwkGtxZ7lp3vy4GZ3Qs8A/y9ma0zsxnFrlOBHQdMI+iRrQpfk4tdqQI7AGgwsxcJOimPu3tFTKurMEOAp8zsBeA54L/c/Q/5unjJTyMUEZFoJd8DFxGRaArgIiIJpQAuIpJQCuAiIgmlAC4iklAK4CIiCaUALiKSUP8fFA0ztfB8yvoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(x,y2err,'bx',linewidth=2, label='Euler')\n", "plt.plot(x,y3err,'go',linewidth=2, label='Heun')\n", "plt.plot(x,y4err,'b.',linewidth=2, label='RK45')\n", "\n", "plt.grid()\n", "plt.title(\"Error in methods for y'=y\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "included-criticism", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.8.10" } }, "nbformat": 4, "nbformat_minor": 5 }