{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "wound-approach", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from numpy import exp, sin\n", "from scipy.integrate import solve_ivp" ] }, { "cell_type": "code", "execution_count": 2, "id": "perfect-helping", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEICAYAAAB25L6yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0JUlEQVR4nO3deVyVZdrA8d8lIoqau2gqYJlWbim2aFqgleWU47QnmZZFm00108xb0Zg1QzPTNFPZYtmeEurYa2Zlm0K9mmZopmmZGyBuuaKIK9zvH/cDHPAA58DZgOv7+ZzPOedZr/Oc48Xt/dyLGGNQSikVuhoEOwCllFKV00StlFIhThO1UkqFOE3USikV4jRRK6VUiNNErZRSIU4TtfI5EUkUkc+DHUcxEWkiIvNEJE9E/hvseAJFRN4Wkb/56FiTRGS6L46lvKeJOoSJyGgRyRSRfBHZLiLzRWRwsOOqijEm1RhzWbDjcHEtEAW0McZc524DETlbRD50kvlBEUkXkUEu62NFxDjfRb6I7BSRj0Tk0nLHyRKRwy7b5YvIi/79eCAi40Rkkb/Po4JDE3WIEpE/AM8BT2GTTDTwMvDbIIZVJRFpGOwY3IgBfjHGnHC3UkROBxYDq4GuwKnAHOBzERlYbvOWxphmQF/gC2COiIwrt81VxphmLo8JPvwsqj4yxugjxB5ACyAfuK6SbSKwiXyb83gOiHDWxQO5wJ+BX4HtwChgBPALsBd41OVYk4DZwEzgILAC6Ouy/mFgo7NuLfA7l3XjsEnuWWAP8Ddn2SJnvTjrfgUOYJNhL5fP+S6wC8gGHgMauBx3EfAMsA/YDFxRyfU4C8gA9gNrgJHO8ieAY8Bx55qOd7PvNOATN8unAF87r2MBAzQst81DwE6XuLOASzz8nicB/wWmO9d2NdAdeMS5XluAy8r9Lt5wvs+tzrUOcz77EaDQ+Yz7ne3fBl4CPnaO/y1wusvxBgHfAXnO8yCXdV2Br5z9vgBeBKY76xo7Me9xrvd3QFSw/93U5UfQA9CHmy8FLgdOlE8K5bZ5ElgKtAfaAd8Af3XWxTv7TwTCgTucZPge0BzoCRwGujrbT3IS2bXO9g85iTHcWX8dtpTZALgBOAR0dNaNc851H9AQaELZRD0cWA60xCbts1z2fReY68QUi/0jMt7luMed2MOAu7F/kMTNtQgHNgCPAo2AoU6C6eHy+aZXci13ALe6WZ7gJL8mVJyoT3OWn+W8z8K7RH3EuUYNneuxGUh2+d42u2w/B3gVaOp878uAO12u16Jyx38bm0zPc46fCsxw1rXG/gEc46y7yXnfxlm/BPgPtkBwkXM9ixP1ncA8INL5buKAU4L976YuP4IegD7cfCmQCOyoYpuNwAiX98OBLOd1PDYRhznvmzvJ5HyX7ZcDo5zXk4ClLusaYEttQyo490rgt87rcUBOufUlScNJmr8AF+CUOp3lYdiS7tkuy+4EMlyOscFlXaTzGTq4iWcINtm6Hj8NmOTy+SpL1CeAy90sP9M5ZycqTtSNneUXOu+zcEq1Lo87KjjvJOALl/dXOfuW/95aYqu/jgJNXLa/CUgvf81d1r8NvO7yfgTws/N6DLCs3PZLnONEO9ekqcu69yhN1LdhCwZ9gv1vpb48QrE+UdlSUFsRaWgqqFfFlnCzXd5nO8tKjmGMKXReH3aed7qsPww0c3m/pfiFMaZIRHKLjycitwB/wCYrnP3autu3PGPMQudm2ktAjIj8L7bE3gRbaiz/GTq5vN/hcpwCESk+d3mnAluMMUWVHKsyu4GObpZ3BIqwJc32FexbfI69LstGGWO+9PDc5b+T3W6+t2bYzxgObHeuA9g/qBVee8cOl9cFlF6/8r8fKL1mpwL7jDGHyq3r4rye5ryeISItsdUgycaY41XEoqpJbyaGpiXY0tOoSrbZhr1JVizaWVZdxf8IEZEGQGdgm4jEAK8BE7D/LW4J/IitxihW6RCMxpjJxpg44GxsHeyfsMnxuJvPsLUasW8DujhxV+dYX2Krd8q7HlhijCmoZN/fYeuT13l4ruragv1NtDXGtHQepxhjejrrvR0Gs/zvB0qv2XaglYg0LbfOnsiY48aYJ4wxZ2Prua8EbvHy/MoLmqhDkDEmD1u//JKIjBKRSBEJF5ErRORpZ7M04DERaScibZ3ta9LONU5ErnZabTyATQpLsfWhBlvHjYjcCvTy9KAicq6InC8i4di67SNAkVNqnAWkiEhz5w/CH6r5Gb7Flhb/7FyneGw1wgwP938CGCQiKSLS2onnPmzy+Z8KPleUiEwAHgceKVea9zljzHbgc+DfInKKiDQQkdNF5GJnk51AZxFp5OEhPwG6O01AG4rIDdg/pB8ZY7KBTOAJEWnkNAm9qnhHEUkQkd4iEoa9QXwc+z8P5SeaqEOUMebf2MT1GDZJbsGWaj9wNvkb9h/TKmxrgRXOsuqai71RWHyD6Wqn5LQW+De2lL8T6I1t5eGpU7Al8n3Y/z7vAf7lrLsPm7w3YVt4vAe86W3gxphj2ERyBbak/jJwizHmZw/3Xw8Mxja5y8KWKK8Bhhtjyn/W/SJyCHvNR2Bb5pSPeV65dtRzvP1MFbgFe7N0LfZ6zqa0ymYhtrXLDhHZXdWBjDF7sCXhP2K/kz8DVxpjivcdDZyPrdJ5HHujs1gH59wHgJ+wrUOm1eSDqcqJc3NA1WMiMgnoZoy5OdixKKVOpiVqpZQKcZqolVIqxGnVh1JKhTgtUSulVIjzS4eXtm3bmtjY2Grte+jQIZo2bVr1hgGmcXlH4/KOxuWduhjX8uXLdxtj2rld6Y/ujnFxcaa60tPTq72vP2lc3tG4vKNxeacuxgVkmgpyqlZ9KKVUiNNErZRSIU4TtVJKhbiAjZ53/PhxcnNzOXLkSKXbtWjRgp9++ilAUXnOV3E1btyYzp07Ex4e7oOolFL1QcASdW5uLs2bNyc2NhaXYRpPcvDgQZo3bx6osDzmi7iMMezZs4fc3Fy6du3qo8iUUnVdwKo+jhw5Qps2bSpN0nWdiNCmTZsq/1ehlFKuAlpHXZ+TdDG9BkrVTUv+7/9InTaNJUuW+PzYOsOLUkrV0JIlSxg2dCjHTpwg9b33WLBwIQMHlp/AvvrqXauPDz74ABHh558rH6r4ueeeo6Cgsok9Kvf2228zYcKEau+vlKo9MqZO5diJExQCx44dIyMjw6fHr3eJOi0tjcGDB5OWllbpdjVN1EqpemL7duI/+IBGQJgIjSIiiI+P9+kp6lWizs/PZ9GiRbzxxhvMmGFnaSosLOShhx6iV69e9OnThxdeeIHJkyezbds2EhISSEhIAKBjx9K5T2fPns24ceMAmDdvHueffz79+vXjkksuYefOnSedVylVRxUWws03s6nLflqMaUThUEOLpBZsarbJp6cJTh11JTfUatQAroohW+fOncvll19O9+7dadOmDcuXL2fZsmVkZWWxcuVKGjZsyN69e2ndujX/+c9/SE9Pp23btpUec/DgwSxduhQR4fXXX+fpp5/m3//+d00+hVKqtvjnP0ndtZCkkVAQfgxOhx3sIGleEgCJvRN9cpp6dTMxLS2N+++/H4Abb7yRtLQ0Nm/ezF133UXDhvZStG7d2qtj5ubmcsMNN7B9+3aOHTum7aOVqi8WL4aJE0m+DwrK9V8rOF5A8oLkWp6oKyn5+qvDy969e1m4cCGrV69GRCgsLEREOPfccz3a37VZnWs76Pvuu48//OEPjBw5koyMDCZNmuTr0JVSoWbvXhg9GgoLyWnhfpOcvByfna7e1FHPnj2bMWPGkJ2dTVZWFlu2bKFr16707duXV199lRMnTgA2oQM0b96cgwcPluzfrl07fvrpJ4qKipgzp3RS6by8PDp16gTAO++8E8BPpJQKCmPg9tshJwfOP5/oFtFuN6toeXXUm0SdlpbG7373uzLLrrnmGrZv3050dDR9+vShb9++vPfeewAkJSVx+eWXl9xMfOKJJ7jyyisZNGhQmRuLkyZN4rrrriMuLq7K+mylVB0wZQrMmQOnnAJpaaRc8hSR4ZFlNokMjyRlWIrvzlnRQNU1ebibOGDt2rUeDZ594MCB6oy57Xe+jMvTa+GJujiAuj9pXN7RuMpZudKYiAhjwJhZs0oWT1813cQ8G2NkkpiYZ2PM9FXTvT40lUwcUGUdtYj0AGa6LDoNmGiMec53fy6UUirEHToEN9wAR49CUhJcd13JqsTeiST2TiQjI8PnbajBg5uJxph1wDkAIhIGbAXmVLaPUkrVOffdB+vWQc+e8OyzAT21t3XUw4CNxphsfwSjlFIhKTUV3noLmjSBmTMhMrLqfXxITBWdRMpsLPImsMIY86KbdUlAEkBUVFRccc+/Yi1atKBbt25VnqOwsJCwsDCPYwoUX8a1YcMG8vLyfHKs/Px8mjVr5pNj+ZLG5R2NyzuBjKvJ1q3E3XEHDQ8fZt0f/8j2K6/0S1wJCQnLjTED3K6sqPK6/ANoBOwGoqraVm8mVk5vJgaPxuWdeh/XkSPG9O9vbx5ef70xRUV+iwsfzUJ+BbY0rYNZKKXqh0cegRUrIDYWpk6tdPgLf/ImUd8EVD7kXIgLCwvjnHPOKXn84x//qHR7HapUqXrso4/sTcOGDWHGDGhRQRfEAPCoC7mINAUuBe70bzilUlMhOdl2/omOhpQUSKxht/kmTZqwcuVKn8TnzokTJ0rGDFFK1WJbt4IzQiZPPQXnnx/UcDwqURtjDhlj2hhjfHMHrAqpqbaZYna27a2ZnW3fp6b653yxsbHs3r0bgMzMTLftIHfv3s0111zDueeey7nnnsvixYsB2zNxzJgxXHjhhYwZM8Y/ASqlAqew0JYK9+yB4cPhj38MdkShOXpecjKUH7O/oMAur0mp+vDhw5xzzjkl7x955BFuuOEGj/b985//zIMPPsjgwYPJyclh+PDh/PTTTwCsXbuWRYsW0aRJk+oHp5QKDX/7G3z1FXToAO++Cw2CP9JGSCbqnAoGnapouadqUvWRkZHB+vXrS94fOHCA/Px8AEaOHKlJWqm64Kuv4Mkn7U3DadOgfftgRwSEaKKOjrbVHe6W+0PDhg0pKioCyg5h6qqoqIilS5fSuHHjk9Y1bdrUP4EppQJn9277X/aiInj0UbjkkmBHVCL4ZXo3UlJO7vgTGWmX+0NsbCzLly8H4P3333e7zdChQ3nhhRdK3vvzpqRSKsCMgVtvtTcRBw2CEBtXPiQTdWKibbIYE2P/BxITY9/XtNVHcR118ePhhx8G4PHHH+f+++9nwIABFfY+/Ne//kVmZiZ9+vTh7LPP5pVXXqlZMEqp0DF5sm2O17IlvPcehIdXuUsghWTVB9ikXNPEXF5hYaHb5UOGDOGXX345afm4ceNKJrFt06YNM2fOPGkbndFFqVpu+XL405/s6zfesCXDEBOSJWqllAqIgwfhxhvh+HG45x64+upgR+SWJmqlVP1kDNx9N2zYAH36wL//HeyIKqSJWilVP737ru1FFxlpu4i7adEVKjRRK6Xqn3XrbFUHwIsvwllnBTeeKmiiVkrVL0eO2Cm1Cgpg9OjSMT1CmCZqpVT98qc/wQ8/wOmn2xnFgzR0qTfqVaIuP/OCDmOqVD3zwQe2qiM83NZLn3JKsCPySMgm6tTVqcQ+F0uDJxoQ+1wsqav9NHSeUqp+yMmB226zr//5TxjgftarUBSSiTp1dSpJ85LIzsvGYMjOyyZpXpJfk/WuXbsqHMb0mWeeKdmuV69eZGVlkZWVxVlnncUdd9xBz549ueyyyzh8+LDf4lNK1cCJE7Y+et8++M1v4IEHgh2RV0IyUScvSKbgeNlxTguOF5C8ILlGxy3fhXzixIkl6+6//34efPBBvvvuO95//31uv/32Ko+3fv167r33XtasWUPLli0rHCdEKRVkTzwBixfDqafa2cRrQb20q5DsQp6T534804qWe6r8MKdvv/02mZmZAHz55ZesXbu2ZJ3rMKYV6dq1a8n41nFxcWRlZdUoPqWUHyxcaEd0E7Htptu18/kpSmekuthnM1K5CslEHd0imuy8k8c5jW7hp3FOqXgYU9chUKHsMKgRERElr8PCwrTqQ6lQ8+uvNmMaAxMngpvZm2qqeEYqO9mJlMxIBb5L1iFZ9ZEyLIXI8LLjnEaGR5IyzE/jnAKXXXaZ22FMY2NjWbFiBQArVqxg8+bNfotBKeVDRUW2jfSOHTBkCPzlL345TWUzUvlKSCbqxN6JTL1qKjEtYhCEmBYxTL1qKom9fTycnovJkye7Hcb0mmuuYe/evZx33nm8+OKLdO/e3W8xKKV86NlnYf58aN3aFnv9NPG0v2akchWSVR9gk7WvE3P5OmfXYUzbtm3rdhjTJk2a8Pnnn3Pw4EGaN29eZt2PP/5Y8vqhhx7yaaxKqRr47jtwxpvnrbegSxe/nSoQM1KFZIlaKaWqLS/PDl164gT8/vcwcqRfTxeIGak8StQi0lJEZovIzyLyk4gM9F0ISinlG6mrUun4+ygkdhMdE8NIHdvf7+csOyOV8dmMVK48rfp4HvjUGHOtiDQCIqvaQSmlAil1dSrjJ4/l6HuFUAg7wgoZ3+xOuLehX+9vQemMVBkZXxHvh5YlVZaoRaQFcBHwBoAx5pgxZr/PI1FKqRpI/vB+jmbZJI0BCuHohqM17igXCsQYU/kGIucAU4G1QF9gOXC/MeZQue2SgCSAqKiouBkzZpQ5TosWLejWrVuVARUWFlY4wWww+TKuDRs2kJeX55Nj5efnnzTYVCjQuLyjcXmnfFytli2jX8H/YHKBd7DJOgwYC9JFWHjxwqDE5Y2EhITlxhj3A5AYYyp9AAOAE8D5zvvngb9Wtk9cXJwpb+3atSctc+fAgQMebRdovozL02vhifT0dJ8dy5c0Lu9oXN4pE9d33xnTtKmJeQDDJAzjMQxznidhYp6NCU5cXgIyTQU51ZObiblArjHmW+f9bMD/NfR+EBYWxjnnnEOvXr246qqr2L9/PwBZWVn06tWrZLvXXnuNuLg49u3bV7LshRdeQETYvXs3ABkZGbRo0aJk3JAnn3wyoJ9FKYWd73DECDh0iJQjF9qOcl2AIUAX/3eUC5QqE7UxZgewRUR6OIuGYatBap3isT5+/PFHWrduzUsvvXTSNtOmTeOFF17gs88+o1WrVgBs2bKFBQsWEF2uYeSQIUNYuXIlK1euLDPAk1IqAHbuhMsvh1274LLLSHx+YcA7ygWKp60+7gNSnRYfm4Bb/RdSqSVLlpCRkUF8fDwDB/q2ReDAgQNZtWpVmWWzZs3iH//4BwsWLKBt27Ylyx988EH++te/Mnr0aJ/GoJSqnrDDh+1wpRs3QlwczJ4NjRr5paNcKPAoURtjVmLrqgNmyZIlDBs2jGPHjtGoUSMWLFjgs2RdWFjIggULGD9+fMmy7OxsJkyYwPfff0+HDh1Kls+dO5dOnTrRu3dvtzH27duXU089lWeeeYaePXv6JD6lVCWOHaPnxImwfLmdTuvjj6Fcr+G6JmR7JmZkZHDs2DEKCws5duwYGRkZNT5m8XjUHTp0YOfOnVx66aUl69q1a0d0dDSzZs0qWVZQUMBTTz3ltv65f//+ZGdn88MPP3DfffcxatSoGsenlKpCURGMH0/rzEw7XOmnn0JUVLCj8ruQTdTx8fE0atSIsLAwGjVq5JNG5MV11NnZ2RhjytRRR0ZG8sknn/DKK6+Qmmpnktm4cSObN2+mb9++9OrVi9zcXPr378+OHTs45ZRTSprhjBgxguPHj5fcaFRK+ckjj8D06RQ2bgyffAIeNPmtC0J2UKaBAweyYMECv9RRR0ZGMnnyZEaNGsU999xTsrx9+/Z8+umnxMfH07ZtW4YPH86vv/4KwMGDB+nduzeZmZm0bduWHTt2EBUVhYiwbNkyioqKaNOmjc9iVEqV8/zz8PTT0LAhPz7xBH1r0ZyHNRWyiRpssvb1TcRi/fr1o0+fPqSlpTFkyJCS5V27duXDDz9kxIgRzJkzh/POO8/t/rNnz2bKlCk0bNiQJk2aMGPGDKSWTe+jVK0xcyY8+KB9/eab7PPjaHihKKQTta+VH+Z03rx5Ja9dhyzt27cvW7duPWl/16m2JkyYwIQJE3wfpFKqrPR0uOUWO0vLP/8JY8aAD+5Z1SYhW0etlFL88AOMGgXHjtkhS//0p2BHFBSaqJVSoSkrC664Ag4cgOuvtzO21NPqxYAmalPFAFD1gV4DpTywZ4/tdbh9u52Q9t13oUH9LVcG7JM3btyYPXv21OtEZYxhz549J810rpRyUVAAV14J69ZBnz7wwQcQERHsqIIqYDcTO3fuTG5uLrt27ap0uyNHjoRkIvNVXI0bN6Zz584+iEipOujECbjhBli61E46OH8+tGgR7KiCLmCJOjw8nK5du1a5XUZGBv369QtARN4J1biUqjOMgbvugo8+sjOHf/YZnHpqsKMKCfW30kcpFVomTYI33oAmTWyyPvPMYEcUMjRRK6WC75VX4Mkn7Q3DmTPBTx3daitN1Eqp4PrgA7j3Xvv61VfhqquCGk4o0kStlAqeRYvgppvsqHhPPAG33x7siEKSJmqlVHCsWWNLz0eOwJ13wl/+EuyIQpYmaqVU4OXm2g4t+/fDb38LL71Ub3sdekITtVIqsPbts13Dc3Nh0CBIS4OwMJ8cOjUVYmPtPcnYWPu+LqhXo+cppYLsyBE7yNKPP8JZZ8G8ebY5ng+kpkJSku3YCJCdbd8DJNbyaRS1RK2UCozCQrj5Zvj6a+jUyU6j1bq1zw6fnFyapIsVFNjltZ0maqWU/xkD998P779vu4TPn2+7iPtQTo53y2sTTdRKKf/7+9/tDcOICJg7F3r39vkpKsr7Pv57EBQeJWoRyRKR1SKyUkQy/R2UUqoOefttW/8gYiuSL77YL6dJSYHIyLLLIiPt8trOm5uJCcYYnWZbKVWp1NWpJC9IJicvh+hGbUmZsZtEgMmT4Zpr/Hbe4huGycm2uiM62ibp2n4jEbTVh1LKh1JXp5I0L4mC4/auXvaxXSRdCYwcSWIA5hhNTKwbibk88WQgfxHZDOwDDPCqMWaqm22SgCSAqKiouBkzZlQroPz8fJo1a1atff1J4/KOxuWduhLXjUtvZOfRnbAFyAJigS4QFdGeGRfMDFpcgVKTuBISEpYbYwa4XWmMqfIBdHKe2wM/ABdVtn1cXJyprvT09Grv608al3c0Lu/UlbhkkhjGY2iIQZzn8RiZJEGNK1BqEheQaSrIqR7dTDTGbHWefwXmAOdV60+GUqpOi45ob0vShdj/fxcCWRDdog40vQiiKhO1iDQVkebFr4HLgB/9HZhSqpaZO5eUWbuJ6ASEAWKfI7pFkDKsDjS9CCJPbiZGAXPEDpjSEHjPGPOpX6NSStUub74Jd9xBYlERXDSUh+5dw44fd9KhVweeGf8Mib3r4B2+AKoyURtjNgF9AxCLUqq2MQaefhoefti+nziRxEmTSNSR8HxKm+cppaqnqAgeegiefdZ2Zpk8GQLQBK8+0kStlPLe8eNw220wfTqEh8O0aXDDDcGOqs7SRK2U8s6hQ3DddXZgpaZNYc4cuPTSYEdVp2miVkp5bu9e+M1vYOlSaNsWPvkEzj032FHVeZqolVKeyc2F4cNh7Vo7kMbnn0OPHsGOql7QYU6VUlX7+Wc7bdbatdCzJ3zzjSbpANJErZSq3LffwuDBsGWLTdbFM7SogNFErZSqUKtly2DoUNizx9ZNf/GFT6fPUp7RRK2Uci8tjd6PPmonHrzlFtu6o/zI/CogNFErpU42eTKMHk2DwkLbqeWtt2x7aRUUmqiVUqWMgccesxPRAhvvvBP+9S9ooKkimLR5nlLKKiyEu++G116DsDB4/XW2xMZyerDjUlqiVkoBR47Y3oavvQaNG9v66HHjgh2VcmiiVqq+y8uDyy+3ybllS9uy46qrqn241FSIjbW1JbGx9r2qGa36UKo+27EDrrgCVq6Ejh3hs8+gd+9qHy41FZKSbEMRgOxs+x7q5qSzgaIlaqXqq40b4cILbZI+4wzb27AGSRogObk0SRcrKLDLVfVpolaqPlq50ibpTZsgLg4WLbL1FDWUk+PdcuUZTdRK1TdffQUXXww7d8KwYZCeDu3b++TQ0RXMYVvRcuUZTdRK1Sdz5tgR8A4csK08Pv4Ymjf32eFTUk7uvBgZaZer6tNErVQdlbo6ldjnYmnwRANin4sl9fnb4dpr4ehR2146LQ0iInx6zsREmDoVYmLs7FwxMfa93kisGW31oVQdlLo6laR5SRQct3f2svOySTr2BvSExGsmwcSJNpP6QWKiJmZf0xK1UnVQ8oJkm6S3AP8HbIGCRpB8XSt4/HG/JWnlHx6XqEUkDMgEthpjrvRfSEqpmsrJy7FJ+h2gEAgDxkJOl/1BjUtVjzcl6vuBn/wViFLKd6IjO0IWNkkb5zkLolto84vayKNELSKdgd8Ar/s3HKVUjU2bRsrM3UR0wpakxT5HdIsgZZg2v6iNxBhT9UYis4G/A82Bh9xVfYhIEpAEEBUVFTdjxoxqBZSfn0+zZs2qta8/aVze0bi844u4Ghw5whmTJ9Nx/nwAXr7xbB5ruY196/fT6oxW3BN/D5dEXRLwuPyhLsaVkJCw3BgzwO1KY0ylD+BK4GXndTzwUVX7xMXFmepKT0+v9r7+pHF5R+PyTo3jWrPGmJ49jQFjGjc25rXXjCkqCn5cflIX4wIyTQU51ZObiRcCI0VkBNAYOEVEphtjbq7Wnw2llG+98w7cc48dVOPMM2HWrBqP2aFCS5V11MaYR4wxnY0xscCNwEJN0kqFgEOH7JjR48bZJD1mDHz3nSbpOkg7vChVG/34I1x/Pfz0EzRpAi+9ZBO2to+uk7xK1MaYDCDDL5EopapmjJ1odsIEOHwYzjoL/vtf6Nkz2JEpP9KeiUrVFvn5cMstMH68TdLjxtmqDk3SdZ5WfShVG6xebUe7W7fODkf38sswdmywo1IBoiVqpUKZMXbC2fPOs0m6Z09bitYkXa9oolYqVB08CDffbCcdPHLEVnksWwZnnx3syFSAaaJWKhT98AMMGADvvQdNm8K0afD66yePyl8JnQ287tA6aqVCiTHw6qvwwAN2gP/evW0HljPP9OowOht43aIlaqVCxYEDcNNNdvaVo0fhjjvg22+9TtKgs4HXNVqiVioENFu/3ibmDRugWTM7f9VNN1X7eDobeN2iiVqpYDIGpkyh/wMPwPHj0Levrero3r1Gh42OttUd7par2kerPpQKlrw82w383ntpcPw43HUXLF1a4yQNOht4XaOJWik/Omkm8NVO04vly6F/f5g9G5o3Z83EiTBlCjRu7JPz6mzgdYtWfSjlJ25nAp+XBJ99RmLyTDh2DPr1g1mz2JWb6/Pz62zgdYeWqJXyk5KZwF0UHC8gees0m6TvvRe++Qa6dQtShKq20BK1Un6Sk+c0sdiCnWg2FugCOS2wI95de22wQlO1jJaolfKT6BbRNkm/Ayx0nrdAdLNTNUkrr2iiVspPUqJGE74JKASMfQ7PCSdl+NNBjkzVNpqolfK17dth9GgSR/+dR7NBwgABaSg8esujJPbWO3zKO1pHrZSvFBbacaIfe8x2B2/cmEnj/8LwwYPJWLyY+Ph4Bg4cGOwoVS2kiVopX1i2zHZY+f57+/7KK2HyZOjalYHAwIsuCmp4qnbTqg+lamLfPjuI0gUX2CTdpQvMmQMffghduwY7OlVHaKJWqjqMgXffhR494JVXICwM/vxnOyv4qFE6G7jyKa36UMpba9bAPffA11/b9xddZOumdZJZ5SdVlqhFpLGILBORH0RkjYg8EYjAlAo5hw7Bww/DOefYJN2uHbzzDmRkVJikdZYV5QuelKiPAkONMfkiEg4sEpH5xpilfo5NqdAxdy78/vd2QGcRuPNOeOopaN26wl10lhXlK1WWqI2V77wNdx7Gr1EpFSqysmDkSFvvnJNjB1FassTWS1eSpEFnWVG+I8ZUnXNFJAxYDnQDXjLG/I+bbZKAJICoqKi4GTNmVCug/Px8mjVrVq19/Unj8k5tj0uOH6fLrFnETJtG2NGjnIiMZPNtt7Ft1ChMWJhH5xo69GKMOfmmoohh4cKvqhVXoGlc3qlJXAkJCcuNMQPcrjTGePwAWgLpQK/KtouLizPVlZ6eXu19/Unj8k6tjmvhQmPOPNMY27bDmBtvNGbbNq/PFRNTegjXR0xMNeMKAo3LOzWJC8g0FeRUr5rnGWP2O4n68mr9yVAqlO3cCWPGwNCh8PPPdqaVL76AtDTo2NHrw+ksK8pXPGn10U5EWjqvmwCXAj/7OS6lAqe463ePHjB9up1l5a9/hVWr4JJLqn1YnWVF+YonrT46Au849dQNgFnGmI/8G5ZSvpW6OpXkBcnk5OUQvTKalGEpdnCkzEzbszAz0254xRXw4otw2mk+Oa/OsqJ8ocpEbYxZBfQLQCxK+YXbKbE+vAPefofEZ7+0VcedO8Pzz8Pvfqe9ClXI0S7kqs5zOyXWicMkF31he6L88Y+26/fVV2uSViFJu5CrOq/SKbG+/x569w5WaEp5REvUqs6LjuzofkqsltGapFWtoIla1V0bNsC4caRM3074RspOibUlnJRhTwU5QKU8o4la1T1OgubMM+Gdd0hc04BH256DNKR0SqwxOiWWqj00Uau6Y8MGuPXWkgQNwPjx8MsvTPrv9yz+6htuH387i79azKSbJwU1VKW8oTcTVe23cSP87W8wbZrtvBIWBrfdZkc/cmkPPXDgQI4eParzFqpaR0vUqvbauNEm5B494O237bLbboNffoE33iiTpIvHhR469GIdF1rVOlqiVrXPxo12wIx33y0tQd96qy1Bn376SZuXHRdadFxoVetoiVrVHps22TrnHj3grbfssltvhXXr4M033SZp0HGhVe2nJWoV+jZtsiXod97xqARdXk6Od8uVCjVaolahq7gE3b27LTGDbXb388+VlqDLi472brlSoUYTtQo9mzfD7bfbKo7yCfqtt6BbN68Op+NCq9pOE7UKqNTVqcQ+F0uDJxoQ+1wsqatdml8UJ+ju3W2rjaIiGDu22gm6WNlxoY2OC61qHa2jVgHjdrjReUnw6y4SZ6yxTexOnLAj2o0da+ugzzjDJ+cuHhc6I+Mr4uPjfXJMpQJFE7UKGLfDjR4vIPmjB0l8HZugb7kFHnvMZwlaqbpAE7UKmEqHG9UErVSFNFGrwDCG6Ij2ZG/YaYcZLQTCgLEQ3ePU0rE5lFIn0ZuJyr/274fJk6FnT1Le2+l+uNHLnw5ujEqFOC1RK79otm6dndE7La2kW2Bix46sH3AGTy7+P8wJo8ONKuUhTdTKdwoKYOZMmDKFAd99V7p82DA70/fIkZwxK5yo+UvYsSODqDbxnCE6kp1SVdFErWpu3Tp45RXbvG7/fgCON29O+O23w5132o4ruA6ONBAYyI4dOjiSUp6oMlGLSBfgXSAKW7M41RjzvL8DUyHu+HGYOxemTIGFC0uXn3ce3H03Szp25KLhw8vsUtngSJqolaqYJyXqE8AfjTErRKQ5sFxEvjDGrPVzbCoUbdkCr70Gr78O27fbZZGRMHq0rd7o3x+AooyMk3bVwZGUqp4qE7UxZjuw3Xl9UER+AjoBmqjri6Ii+OILW3qeN8++BzjrLJucx4yBli2rPEx0NGRnu1+ulKqYGGM831gkFvga6GWMOVBuXRKQBBAVFRU3Y8aMagWUn59Ps2bNqrWvP9XHuMLz8ugwfz6nzptHk23bAChq2JDdQ4awdeRI8vr2BRGP4/ryy/Y880wPjh4NK1kWEVHIQw+t45JLfvXLZ/AkrlCgcXmnLsaVkJCw3BgzwO1KY4xHD6AZsBy4uqpt4+LiTHWlp6dXe19/qktxTV813cQ8G2NkkpiYZ2PM9FXTS1cWFRmzeLExN99sTESEMWAf0dHGpKQYs2NHjeKaPt2YmBhjROzz9OluN/ObuvQ9BoLG5Z2axAVkmgpyqketPkQkHHgfSDXG/G+1/lyokFDhwEiHj5C44rit3li1ym4sAiNG2OqNK66wA/bXUPHgSEopz3nS6kOAN4CfjDH/8X9Iyp9KBkZyGW+joEsByTPuIPFZpxqsXTs7YH9SEnTtGrxglVKAZ60+LgTGAKtFZKWz7FFjzCd+i0r5TU5ejk3S5cbbyOlsYMgQW3q++mqIiAhuoEqpEp60+lgEuL9jpGqPvDyYO5foo43JzjpcZrwNsiC6e0f4+uvgxqiUcksHZarLDhyw422MHAnt28PYsYz46Bzo1MiWpAX73DmCEc3/FdxYlVIV0i7kdc3Bg7at86xZ8OmncPSoXd6gASQk8MkPc4EP4eqHYPcOaNsBVjzDJxmJcHdQI1dKVUATdV1w8CB89JFNzvPnlyZnEbj4Yrj+elvv3KEDOQ2AvYmwumzTixyt3FIqZGmirq3y8+Hjj+n50kvw3Xdw5IhdLmJvCl5/PVxzDXTsWGY37R2oVO2jibo2OXQIPvnElpw//hgOH6Zd8boLLyxNzp06VXiIlJTiEexKl0VG2uVKqdCkiTrUFRTY6oxZs2z1hmuGHTiQDf370+3hh6FzZ48OV9zZJDnZDoYUHW2TtHZCUSp0aaIOotTVqSQvSCYnL4foFtGkDEuxs50cPmxvBM6aZW8MHjpUutMFF9iS87XXQpcu5GZk0M3DJF1MewcqVbtoog4St12554yHl14iMXW1rYMudt55pck5JiZIESulgkUTdZC478p9lOQmS0jMBwYMKE3O2o1bqXpNE3UgFRXBypXw2WdkH82GXE7qyp3dWWDjBjjttKCGqpQKHZqo/W3bNjvo/uef2+dduwAIe6AThVlbT+rKHdYqWpO0UqoMTdS+dvgwLFoEn31mk/Pq1WXXd+kCw4dTuCAe+o2HsKOlJepOERR+lgLam1sp5UITdU0ZA2vW2KT82Wd2YKPizidgGyknJMBll9lHjx4gQkwsZH9P2a7c3z9DzAFtjqGUKksTdXXs3l1anfH557Z6w1W/fjB8uE3Mgwa5HTLUdjxJpMClK3dkJKRM9XfwSqnapt4n6grbMrs6dowWK1eWJuYVK2xJuliHDjYpDx8Ol1xiR6qrgnY8UUp5ql4n6gqnpTKQ2Pjc0sScnk4/13bNERFw0UWl1Rm9e1c4yWtltOOJUsoT9TpRVzgt1bSxJD5TWGbbQ7GxNL36apuYhwyx9RRKKRUA9TNRGwM5OWTvz6mgLXMhtGkDl15qqzMuvZTv1q8nPj4+qGErpeqn+pGojx6F77+Hb76xjyVLYNu2itsyt4yGXzfbwfaLrV8fnNiVUvVe3UzU27fbZFyclDMz4dixstu0bk3hgn9Av9tPbsv8+VPwjM5SppQKDbU/UZ84AatWlSblb76BrKyTt+vZEwYOtM3lBg2C7t2J6Spkfy/allkpFdJqX6Les8cm5OKkvGxZ2TGaAZo3h/PPtwl54ED7ulWrkw6lbZmVUrVBlYlaRN4ErgR+Ncb08lcg90xJZeqmZAqb5hD2cTRJp6Xw8p03wdq1pUn5m2/gl19O3rlbt9KkPGiQLT2HhVV5Tm3LrJSqDTwpUb8NvAi8668g7pmSypStSbCvAL6HwthspuQmwbmv8fKKr8pu3LgxnHtuaRXGBRd41MGkItqWWSkV6qpM1MaYr0Uk1p9BTN2UbJN0mWZyBUy9aAMv7+pSmpQHDoS+faFRI3+Go5RSIUWMa1foijayifqjyqo+RCQJSAKIioqKmzFjhsdBJGQMhUUGFmKbyQkwFBgspMcv9Pg4/pSfn0+zZs2CHcZJNC7vaFze0bi8U5O4EhISlhtjBrhdaYyp8gHEAj96sq0xhri4OOONsIdiDOMxNMQgzvN4TNhDMV4dx5/S09ODHYJbGpd3NC7vaFzeqUlcQKapIKeGRGPhpNNSoEMkjMWWpMcCHSLtcqWUqudCIlG/fHcid3eaSlirGBgshLWK4e5OU3n5br3Lp5RSnjTPSwPigbYikgs8box5w9eBvHx3Ii+TSEZGho6poZRSLjxp9XFTIAJRSinlXkhUfSillKqYJmqllApxmqiVUirEaaJWSqkQ51HPRK8PKrILyK7m7m2B3T4Mx1c0Lu9oXN7RuLxTF+OKMca0c7fCL4m6JkQk01TUjTKINC7vaFze0bi8U9/i0qoPpZQKcZqolVIqxIViog7V+VU0Lu9oXN7RuLxTr+IKuTpqpZRSZYViiVoppZQLTdRKKRXiApaoReRyEVknIhtE5GE36yNEZKaz/lvX6b9E5BFn+ToRGR7guP4gImtFZJWILBCRGJd1hSKy0nl8GOC4xonILpfz3+6ybqyIrHceYwMc17MuMf0iIvtd1vnzer0pIr+KyI8VrBcRmezEvUpE+rus8+f1qiquRCee1SLyjYj0dVmX5SxfKSKZAY4rXkTyXL6viS7rKv0N+DmuP7nE9KPzm2rtrPPn9eoiIulOLlgjIve72cZ/v7GKZhTw5QM7C+JG4DSgEfADcHa5be4BXnFe3wjMdF6f7WwfAXR1jhMWwLgSgEjn9d3FcTnv84N4vcYBL7rZtzWwyXlu5bxuFai4ym1/H/Cmv6+Xc+yLgP5UMBMRMAKYj53o7QLgW39fLw/jGlR8PuCK4ric91lA2yBdr3js9Hs1+g34Oq5y214FLAzQ9eoI9HdeNwd+cfNv0m+/sUCVqM8DNhhjNhljjgEzgN+W2+a32OltAWYDw0REnOUzjDFHjTGbgQ3O8QISlzEm3RhT4LxdCnT20blrFFclhgNfGGP2GmP2AV8AlwcprpuANB+du1LGmK+BvZVs8lvgXWMtBVqKSEf8e72qjMsY841zXgjc78uT61WRmvw2fR1XIH9f240xK5zXB4GfgE7lNvPbbyxQiboTsMXlfS4nf8iSbYwxJ4A8oI2H+/ozLlfjsX8xizUWkUwRWSoio3wUkzdxXeP8F2u2iHTxcl9/xoVTRdQVO2VxMX9dL09UFLs/r5e3yv++DPC5iCwXO3l0oA0UkR9EZL6I9HSWhcT1EpFIbLJ732VxQK6X2GrZfsC35Vb57TdW5cQByhKRm4EBwMUui2OMMVtF5DRgoYisNsZsDFBI84A0Y8xREbkT+7+RoQE6tyduBGYbYwpdlgXzeoU0EUnAJurBLosHO9erPfCFiPzslDgDYQX2+8oXkRHAB8AZATq3J64CFhtjXEvffr9eItIM+8fhAWPMAV8euzKBKlFvBbq4vO/sLHO7jYg0BFoAezzc159xISKXAMnASGPM0eLlxpitzvMmIAP7VzYgcRlj9rjE8joQ5+m+/ozLxY2U+2+pH6+XJyqK3Z/XyyMi0gf7Hf7WGLOneLnL9foVmIPvqvyqZIw5YIzJd15/AoSLSFtC4Ho5Kvt9+eV6iUg4NkmnGmP+180m/vuN+aPi3U1FfENsBXpXSm9A9Cy3zb2UvZk4y3ndk7I3Ezfhu5uJnsTVD3vz5Ixyy1sBEc7rtsB6fHRTxcO4Orq8/h2w1JTeuNjsxNfKed06UHE5252JvbEjgbheLueIpeKbY7+h7I2eZf6+Xh7GFY297zKo3PKmQHOX198Alwcwrg7F3x824eU4186j34C/4nLWt8DWYzcN1PVyPvu7wHOVbOO335jPLq4HH3QE9k7pRiDZWfYktpQK0Bj4r/OjXQac5rJvsrPfOuCKAMf1JbATWOk8PnSWDwJWOz/U1cD4AMf1d2CNc/504EyXfW9zruMG4NZAxuW8nwT8o9x+/r5eacB24Di2DnA8cBdwl7NegJecuFcDAwJ0vaqK63Vgn8vvK9NZfppzrX5wvufkAMc1weX3tRSXPyTufgOBisvZZhy2gYHrfv6+XoOxdeCrXL6rEYH6jWkXcqWUCnHaM1EppUKcJmqllApxmqiVUirEaaJWSqkQp4laKaVCnCZqpZQKcZqolVIqxP0/uU+kmCWgjLQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "h = .2\n", "x = np.linspace(0,2, num=int(2/.2)+1)\n", "y=exp(x)\n", "\n", "def m_y(t, x):\n", " return x\n", "\n", "yv3 = solve_ivp(m_y, (0,2), (1,), method='RK45', t_eval = x, max_step=h).y[0]\n", "xv = []\n", "yv1 = []\n", "yv2 = []\n", "\n", "y0 = 1\n", "xv.append(0)\n", "yv1.append(y0)\n", "yv2.append(y0)\n", "ya0 = yb0 = y0\n", "for ii in np.arange(h,2+h,h):\n", " ya1 = ya0 + h * ya0\n", " yb1 = yb0 + h * (yb0 + h)\n", " xv.append(ii)\n", " yv1.append(ya1)\n", " yv2.append(yb1)\n", " ya0, yb0 = ya1, yb1\n", "\n", "\n", "# plot results\n", "plt.plot(x,y,'r-',linewidth=2, label='Actual')\n", "plt.plot(xv, yv1,'bo',linewidth=2, label='Euler')\n", "plt.plot(xv, yv2,'go',linewidth=2, label='Heun')\n", "plt.plot(x, yv3,'k.',linewidth=2, label='RK45')\n", "\n", "plt.grid()\n", "plt.title(\"Comparison of ODE methods\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d2119885", "metadata": {}, "source": [ "Comparison of SciPy ODE methods" ] }, { "cell_type": "code", "execution_count": 8, "id": "df8e3774", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEICAYAAACqMQjAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0MUlEQVR4nO3de3wU1fn48c9Dwv0mBEEEIaB+f9zkXiReE6mitICKRYQiKBWlBVv94gWoiihKUeT7UhFLW4oVBFFrBaUghKRKCVelylVuCYRyv4UIISE5vz9mdpksu8lu9pZsnvfrta+dOXNmzrOTzT47c2bOijEGpZRSypsq0Q5AKaVU+aVJQimllE+aJJRSSvmkSUIppZRPmiSUUkr5pElCKaWUT5okVMwQkSEi8mW043ARkZoislhETovIR9GOJ1JEZI6IvByibU0Ukbmh2JYqG00S6hIiMlhENohIrogcFJF/ishN0Y6rNMaYecaYO6Idh8N9QBMgwRjzC28VRKSdiCyyE8kZEUkTkRscyxNFxNh/i1wROSwin4vI7R7byRSRc456uSLydnhfHojIcBFZFe52VPRoklDFiMiTwP8Br2B9wLUA3gH6RzGsUolIfLRj8KIl8IMx5oK3hSJyNfBv4HugFXAl8CnwpYgkeVS/zBhTB+gELAc+FZHhHnX6GmPqOB6jQ/haVGVljNGHPjDGANQHcoFflFCnOlYS+a/9+D+gur0sGcgGngaOAAeBu4E+wA/ACWC8Y1sTgY+BD4EzwDdAJ8fyZ4Hd9rKtwD2OZcOxPmCnA8eBl+2yVfZysZcdAXKwPog7OF7n34CjQBbwe6CKY7urgNeBk8Be4K4S9kdbIB04BWwB+tnlLwL5QIG9T0d4Wfd9YImX8pnAV/Z0ImCAeI86Y4HDjrgzgZ/6+XeeCHwEzLX37ffA/wDj7P21H7jD433xF/vvecDe13H2a88DCu3XeMquPweYAXxhb38tcLVjezcA64HT9vMNjmWtgH/Z6y0H3gbm2stq2DEft/f3eqBJtP9vYv0R9QD0UX4ewJ3ABc8PJI86k4A1QGPgcmA18JK9LNle/3mgKvCI/UH8AVAXaA+cA1rZ9SfaH6L32fXH2h/KVe3lv8D6dl0FuB/4EWhqLxtutzUGiAdqUjxJ9AY2ApdhJYy2jnX/Bnxmx5SIlcBGOLZbYMceB4zCSobiZV9UBXYB44FqwG32h9v/c7y+uSXsy0PAQ17KU+wP3pr4ThKt7fK29nwmgSWJPHsfxdv7Yy8wwfF32+uo/ynwR6C2/XdfBzzq2F+rPLY/B+uDvIe9/XnAAntZQ6zkO9Re9oA9n2AvzwDewPoycou9P11J4lFgMVDL/tt0A+pF+/8m1h9RD0Af5ecBDAEOlVJnN9DHMd8byLSnk7GSQJw9X9f+ILveUX8jcLc9PRFY41hWBevb6s0+2t4E9LenhwP7PJa7P7DsD+wfgJ7Y37bt8jisb/jtHGWPAumObexyLKtlv4YrvMRzM9YHvXP784GJjtdXUpK4ANzppbyN3WYzfCeJGnb5jfZ8Jva3ecfjER/tTgSWO+b72ut6/t0uwzrleB6o6aj/AJDmuc8dy+cAf3bM9wG229NDgXUe9TPs7bSw90ltx7IPuJgkHsb6UtIx2v8rlelRHs/jqug5DjQSkXjj4zw61jf7LMd8ll3m3oYxptCePmc/H3YsPwfUcczvd00YY4pEJNu1PRF5EHgS64MSe71G3tb1ZIxZaXfczgBaisjfsY5UamJ9W/Z8Dc0c84cc2zkrIq62PV0J7DfGFJWwrZIcA5p6KW8KFGF9w27sY11XGyccZXcbY1b42bbn3+SYl79bHazXWBU4aO8HsJK5z31vO+SYPsvF/ef5/oGL++xK4KQx5kePZVfZ0+/b0wtE5DKsU08TjDEFpcSigqAd18opA+tb490l1PkvVoesSwu7rKxcHwCISBWgOfBfEWkJ/AkYjXUq4jJgM9apI5cShzA2xrxpjOkGtMM65/4U1gdzgZfXcKAMsf8XuMqOuyzbWoF1Ss3TQCDDGHO2hHXvweo/2OFnW2W1H+s90cgYc5n9qGeMaW8vD3QYac/3D1zcZweBBiJS22OZ1ZAxBcaYF40x7bD6NX4OPBhg+ypAmiSUmzHmNFZ/wgwRuVtEaolIVRG5S0Sm2tXmA78XkctFpJFdP5jr2LuJyL321Um/w/pAWoN1/ttg9WkgIg8BHfzdqIj8RESuF5GqWH0ZeUCR/W15ITBZROrayejJMr6GtVjfkp+291My1qmbBX6u/yJwg4hMFpGGdjxjsD74nvHxupqIyGjgBWCcx1FMyBljDgJfAtNEpJ6IVBGRq0XkVrvKYaC5iFTzc5NLgP+xL7OOF5H7sZL458aYLGAD8KKIVLMvu+7rWlFEUkTkOhGJw7oYoQDriEuFkSYJVYwxZhrWh+bvsT6g92N9m/+HXeVlrH/k77CuivnGLiurz7A6pV2dmffa3xi3AtOwjm4OA9dhXc3kr3pYRyInsU5ZHAdes5eNwUoce7CuZPoAmB1o4MaYfKwPsbuwjlDeAR40xmz3c/2dwE1Yl7VmYn2THgD0NsZ4vtZTIvIj1j7vg3UFmmfMiz3uk/g00Nfkw4NYHfNbsfbnx1w8TbYS66quQyJyrLQNGWOOYx0B/C/W3+Rp4OfGGNe6g4HrsU6jvYDVqe5yhd12DrAN6yqo94N5Yap0Yoz+6JCKDhGZCFxjjPlltGNRSnmnRxJKKaV80iShlFLKp5AkCRG5U0R2iMguEXnWy/LqIvKhvXytiCQ6lo2zy3eISG9HeaaIfC8im0RkQyjiVOWLMWainmpSqnwL+j4J+0qDGcDtWEMyrBeRRXbHo8sIrOufrxGRQcAfgPtFpB0wCOtO3CuBFSLyP47rtVMcHVpKKaUiLBQ30/XAukN1D4CILMAaDM6ZJPpj3eUJ1tUJb4t1Z05/rNv1zwN7RWSXvb2MsgTSqFEjk5iYWJZV+fHHH6ldu3bpFSNM4wqMxhW48hqbxhWYYOLauHHjMWPM5d6WhSJJNKP43ZfZWJewea1jjLkgIqeBBLt8jce6rjtJDdZomAb4ozFmlrfGRWQkMBKgSZMmvP7662V6Ebm5udSp4+2m2ujSuAKjcQWuvMamcQUmmLhSUlI874J3K8/DctxkjDkgIo2B5SKy3RjzlWclO3nMAujevbtJTk4uU2Pp6emUdd1w0rgCo3EFrrzGpnEFJlxxhaLj+gCOoRWwhlXwHJbAXce+s7Y+1o00Ptc1xriej2CNQtkjBLEqpZQKQCiSxHrgWhFpZd+aPwhY5FFnETDMnr4PWGmsu/gWAYPsq59aAdcC60SktojUBbDHcbkDa9wepZRSERT06Sa7j2E0sAxrGObZxpgtIjIJ2GCMWYT1gyXv2x3TJ7ASCXa9hVid3BeA3xhjCkWkCdYvb7li/MAYszTYWJVSSgUmJH0SxpglWAN3Ocued0zn4X20S4wxk4HJHmV7sMazUUopFUV6x7VSSlVk8+ZBYiK33nYbJCZa8yFUnq9uUkopVZJ582DkSDh71vqhlawsax5gyJCQNKFHEkopVVFNmABnPX6b6uxZqzxENEkopVRFtW9fYOVloElCKaUqqhYtAisvA00SSilVUU2eDLVqFS+rVcsqDxFNEkopVVENGQKzZkHLlhgRaNnSmg9RpzVoklBKqYptyBDIzORfK1dCZmZIEwRoklBKKVUCTRJKKRUs+4Y2qlQJyw1t0aQ30ymlVDAcN7QBYbmhLZr0SEIppYIRgRvaokmThFJKBSMCN7RFkyYJpZQKRgRuaIsmTRJKKRWMCNzQFk2aJJRSKhiOG9oI0w1t0aRXNymlVLCGDImZpOBJjySUUrEjhu9XiBY9klBKxYYYv18hWvRIQikVG2L8foVo0SShlIoNMX6/QrRoklBKxYYYv18hWjRJKKViQ4zfrxAtmiSUUqEVrSuMYvx+hWjRq5uUUqET7SuMYvh+hWjRIwmlVOjoFUYxR5OEUip09AqjmKNJQikVOnqFUczRJKFUrIpGB7JeYRRzNEkoFYtcHchZWWDMxQ7kcCcKvcIo5miSUCoWRbMDecgQyMyEoiLrWRNEhaZJQqlwsk/53HrbbZG9Z0A7kFWIaJJQKlwcp3wkkqd8QDuQVchoklCVQzQ6caN5ykc7kFWIaJJQsS9anbjRPOWjHcgqRDRJqMiqTN/oo33KRzuQVQhoklCRU9m+0espHxUDNElEU7RGy4xWu5XtG73jlI/RUz6qggpJkhCRO0Vkh4jsEpFnvSyvLiIf2svXikiiY9k4u3yHiPT2d5shMXUq3/efQHZ8Irek3EZ2fCLf958AU6eGpbliHniAwuEjin2rLhw+Ah54IDbbBavNQMpDJSmJwvjqxYoK46tDUlJ42506lRVbrySRTOIoJJFMVmy9MjLvr6lTWTEhrdh3gRUT0sLfdrTajWbb5eQ133bbrWFpO+gkISJxwAzgLqAd8ICItPOoNgI4aYy5BpgO/MFetx0wCGgP3Am8IyJxfm4zaN//+zQdFr3C/sKm/IFn2F/YlA6LXuH7f58OdVOXOLckjbgL54uVxV04z7klaTHZLsC5ek0CKg+VFa1HcvZCNQ7RhCKEQzTh7IVqrGg9Mrztnv4JnV4ZSKusNIwRWmWl0emVgaw4/ZOwtntp20Ss7Wi1G822y89rDs97LBS/J9ED2GWM2QMgIguA/sBWR53+wER7+mPgbRERu3yBMeY8sFdEdtnbw49tBq3BF/NYQ096kUo+1ahGPuPpRdHimaS96f2Dyxjj17ZLq2dyDiOAZy2TcxiZPt09v2vXLr799lu/2vQrLrtdPNo2OYeRN97wezu7du3im2++CajtTy/04ud8SDyF7rIC4vj8Qi/umTYtoG35snv3bjZu3FisbPJMaMgQhvI3VtOLJFbzPkM58c5GJjTY6GNLwZv8DjRkAEP5OVVJIokMevJLTryzgfGXbfB7O/6+55xefQcaci+D+Tnx9KQna+jJEE68s45n660rVnfPnj2sXbs2JG3/4R1owN0M4efEcT09Wcv1DObkO2t4us4av7djjGHPnj2sXr3a73Veewca0J9B/Jwq9OB61nE9gzj5zmrG1vJ/O6W95r1797Jq1Sr3/OszoAF93e32YB09uJ9T76zmqdoZAFgfd96fS1pWWt2JM4UEBjGCvvRkFE8yh4EsZO+8FDJD1PUViiTRDNjvmM8GrvdVxxhzQUROAwl2+RqPdZvZ06VtEwARGQmMBGjSpAnp6el+B35L4T7e5wHyqUYh8eRj+IFk6ps1bCphO64/UGlKqlePutTjjFXPUZ5DXXIc/xj5+fkcOXLEr/b84WzX2XYOdTmdkeH3dgoKCjh69GhAba8625yj3E4vVlKDfM5TjRX0YsfZJlzp40MqUAUFBRw7dqxY2cmTzTlJLZbRkSRWsJye7KEOnDrI+vXZIWnXm1OnmnOKOiynI0mksoIk9lAPTh1m48bwtWu1fRWnqE8anbiRlaSRxB4awKnjbNq0v1jdgoICcnJyfG7L3/f7xXYTSKczN5LGv7iBvTSCU6fZvHl/6RvwiGvLli0BtN2CUzRmFV24mXS+5kb2cgWc+pFt2wK7SKGk15yfn8+OHTvc86dPt+A0V7Karna7N5FJUziVy5YtF9t1JR/P59KWuea9LTt58kpOAp/Sls94nZk8RzopyD5Devq/AnrNPrkaL+sDuA/4s2N+KPC2R53NQHPH/G6gEfA28EtH+V/s7ZW6TW+Pbt26mUDsj2tpVtPT1ORHE0e+qcmPZjU9zf64lgFtpyzGJMw1Z6lujNUzYAyYs1Q3YxLmFquXlpYWlXZLU5a4WrY0JpmV5giNzIs8Z47QyCSz0rRsGfCmAoorEu16E612A207lO+xUL7mQOOK1P72jKui/J1LAmwwvj7jfS3w9wEkAcsc8+OAcR51lgFJ9nQ8cAzrS2yxuq56/mzT2yPQJPFdv/GmCMxqeppXeNaspqcpAvNdv/GB7eEyWD5+pcmhrjlIE1OImIM0MTnUNcvHryxWL9RJwt92S1OWuJaPX+l+E8PFN3egbQcaVyTa9SZa7QbadijfY6F8zYHGFan97RlXRfk7l6SkJBGKq5vWA9eKSCsRqYbVEb3Io84iYJg9fR+w0g5sETDIvvqpFXAtsM7PbQbtuhvrs7nfeK6KO8gz/IGr4g6yud94rruxfqibusRP669n7fjP6NnyEPFSRM+Wh1g7/jN+Wn99TLbravs/4xeyt2UKIrC3ZQr/Gb8wIq85+u2aiLV7aduxv6+j2Xb5ec1heo/5yh6BPIA+wA9Yp5Em2GWTgH72dA3gI2AXVhJo7Vh3gr3eDuCukrZZ2iPQIwmnUH9jDxWNKzAaV+DKa2waV2CCiYsSjiRC0XGNMWYJsMSj7HnHdB7wCx/rTgYu6Yf3tk2llFKRpXdcK6WU8kmThFJKKZ80SSillPJJk4RSSimfNEkopZTySZNElGVkwKuvWs9KKVXehOQSWFU2GRnQqxfk50O1apCaGv7Rq5VSKhB6JBFF6emQf95QWGg9BzA2YdD0CEYp5Q89koii5ITvqVZ0NflUpVpRAckJu4Hrwt6uHsEopfylRxJRlHT8c1Kr3MFLPE9qlTtIOv55RNpNT7cSRGGh9RzJIxilVMWiRxLRlJxMUvWXSMpfY32lT34tUs1SrdrFI4nk5Ig0q5SqgDRJRFNSknWuJz3d+qSO0DmfKDWrlKqANElEW1JSVD6lo9SsUqqC0T6Jykovb1JK+UGPJCojvbxJKeUnPZKojPTyJqWUnzRJVEauy5vi4vTyJqVUifR0U2WklzcppfykSaKy0sublFJ+0NNNSimlfNIkoZRSyidNEiry9B4NpSoMTRIqsjIyyEgex6sTcslIHqeJQqlyTjuuVURl/G0nvfKXkE81quXnk/q3j0nSDnSlyi09klARlc6t5FONQuLJpyrp3BrtkJRSJdAkoSIq+cGWVKsuxEkh1apXIfnBltEOSSlVAj3dpCIqKQlS0+L0Pj6lKghNEirionUfX0aG3mSuVKA0SahKQQe+VapstE9CVQo68K1SZaNJQlUKOvCtUmWjp5tUpaAD3ypVNpokVKWRRAZJpAPJgGYJpfyhSUJVDtpzrVSZaJ+Eqhy051qpMtEkoSoH7blWqkz0dJOqHKLYc52RAfPmtaB6dT3DpSoeTRKq8ojCrd6urpDz51sxb552haiKJ6jTTSLSUESWi8hO+7mBj3rD7Do7RWSYo7ybiHwvIrtE5E0REbt8oogcEJFN9qNPMHEqFS2urpCiItGuEFUhBdsn8SyQaoy5Fki154sRkYbAC8D1QA/gBUcymQk8AlxrP+50rDrdGNPZfiwJMk6losLVFVKlSpF2hagKKdgk0R94z55+D7jbS53ewHJjzAljzElgOXCniDQF6hlj1hhjDPA3H+srVWG5ukIefjhTTzWpCkmsz+cyrixyyhhzmT0twEnXvKPOWKCGMeZle/454ByQDkwxxvzULr8ZeMYY83MRmQgMB3KADcD/2gnGWwwjgZEATZo06bZgwYIyvZbc3Fzq1KlTpnXDSeMKjMYVuPIam8YVmGDiSklJ2WiM6e51oTGmxAewAtjs5dEfOOVR96SX9ccCv3fMP2eXdQdWOMpvBj63p5sAcVhHOpOB2aXFaYyhW7dupqzS0tLKvG44aVyB0bgCV15j07gCE0xcwAbj43O11KubjP1N3xsROSwiTY0xB+3TR0e8VDuANQ6CS3Oso4gD9rSz/IDd5mFHG38CPi8tTqWUUqEXbJ/EIsB1tdIw4DMvdZYBd4hIA7vD+g5gmTHmIJAjIj3tU1UPuta3E47LPVhHLkpVTBkZtJg3z7oeVqkKJtj7JKYAC0VkBJAFDAQQke7AY8aYXxljTojIS8B6e51JxpgT9vSvgTlATeCf9gNgqoh0BgyQCTwaZJxKRYd9o0Sr8+fRGyVURRRUkjDGHAd6eSnfAPzKMT8bmO2jXgcv5UODiUupcsO+UUKKii6OGaVJQlUgOnaTUuFk3yhRVKVKVMaMysiAV1/VM12q7HRYDqXCyb5RInP2bFo//HDEx4zS0dFVsDRJKBVuSUnsO3+e1hH+hPY2OromCRUoPd2kVIzS0dFVKOiRhFIxSn/XW4WCJgmlYlgURkdXMUZPNymllPJJk4RSSimfNEkopZTySZOEUkopnzRJKKWU8kmThFJKKZ80SSillPJJk4RSsUxH+FNB0pvplIpVURzhLyND7/SOFZoklIpVURrhT0efjS16ukmpWBWlEf685SZVcemRhFKxKkoj/Llyk+tIQkefrdg0SSgVy6Iwwp+OPhtbYj5JFBQUkJ2dTV5eXon16tevz7Zt2yIUlf8iGVeNGjVo3rw5VatWjUh7Knbp6LOxI+aTRHZ2NnXr1iUxMRER8VnvzJkz1K1bN4KR+SdScRljOH78ONnZ2bRq1Srs7SmlKoaY77jOy8sjISGhxAShQERISEgo9YhLKVW5xHySADRB+En3k1LKU6VIEkoppcpGk4SHefMgMRGqVLGe580LfptxcXF07tyZDh060LdvX06dOgVAZmYmHTp0cNf705/+RLdu3Th58qS77K233kJEOHbsGADp6enUr1+fzp0707lzZyZNmhR8gEop5UPMd1wHYt48GDkSzp615rOyrHmAIUPKvt2aNWuyadMmAIYNG8aMGTOYMGFCsTrvv/8+b731FitXrqRBgwYA7N+/n9TUVFq0aFGs7s0338znn39e9oCUUspPeiThMGHCxQThcvasVR4qSUlJHDhwoFjZwoULmTJlCl9++SWNGjVylz/xxBO89NJL2leglIoaTRIO+/YFVh6owsJCUlNT6devn7ssKyuL0aNH8+WXX3LFFVe4yz/77DOaNWvGddddd8l2MjIy6NSpE3fddRdbtmwJTXBKKeWFJgkHj7M6pZb769y5c3Tu3JkrrriCw4cPc/vtt7uXXX755bRo0YKFCxe6y86ePcsrr7zitb+ha9euZGVl8Z///IcxY8Zw9913BxecUjFGR0cPLU0SDpMnQ61axctq1bLKg+Hqk8jKysIYw4wZMxzbr8WSJUt49913mWf3ku/evZu9e/fSqVMnOnToQHZ2Nl27duXQoUPUq1ePOnXqANCnTx8KCgrcndpKVXauEWife8561kQRPE0SDkOGwKxZ0LIliFjPs2YF12ntVKtWLd58802mTZvGhQsX3OWNGzdm6dKljB8/nmXLlnHddddx5MgRMjMz2bx5M82bN+ebb77hiiuu4NChQxhjAFi3bh1FRUUkJCSEJkClQiVKX+d1BNrQ06ubPAwZErqk4E2XLl3o2LEj8+fP5+abb3aXt2rVikWLFtGnTx8+/fRTevTo4XX9jz/+mJkzZxIfH0/NmjVZsGCBdmyr8iWKPyihI9CGniaJCMjNzS02v3jxYvf05s2b3dOdOnW65MonsO6ncBk9ejSjR48OfZBKhUqUfuwIdATacNAkoZQKrSh/ndcRaENLk4RSKrT063xM0SShlAo9/TofM/TqJqWUUj4FlSREpKGILBeRnfZzAx/1htl1dorIMEf5ZBHZLyK5HvWri8iHIrJLRNaKSGIwcSqllCqbYI8kngVSjTHXAqn2fDEi0hB4Abge6AG84Egmi+0yTyOAk8aYa4DpwB+CjFMppVQZBJsk+gPv2dPvAXd7qdMbWG6MOWGMOQksB+4EMMasMcYcLGW7HwO9JFI3A4RhrPCyDBX+1FNP0aZNG5KSkrjnnnvc66xbt849THinTp349NNPg45PKaV8Edfdu2VaWeSUMeYye1qwvv1f5lFnLFDDGPOyPf8ccM4Y87qjTq4xpo5jfjNwpzEm257fDVxvjLlk/AkRGQmMBGjSpEm3BQsWFFtev359rrnmmlJfS2FhIdU/+YQaY8Yg5865y03NmuS99RYXBg4sdRu+NG3alIMHrVz46KOPcs011/DUU0+RlZXFwIEDWbt2LfPnz+fNN9/k888/JyEhgdTUVG699VZEhBdffBGASZMmcfbsWapVq0Z8fDyHDh3ihhtu4IcffiA+PjTXIOzatYvTp0+XWi83N9c9PEh5onEFrrzGpnEFJpi4UlJSNhpjuntdaIwp8QGsADZ7efQHTnnUPell/bHA7x3zzwFjPerkesxvBpo75ncDjUqLtVu3bsbT1q1bLynzJicnx5iWLY2BSx8tW/q1DV9q167tnp45c6YZNWqUMcaYvXv3mvbt25sPP/zQtGvXzhw8eNBrXH//+9/N4MGDL1m2Z88e07hxY1NQUBBUfE7+7q+0tLSQtRlKGlfgymtsGldggokL2GB8fK6W+vXTGPNTX8tE5LCINDXGHBSRpsARL9UOAMmO+eZAeinNHgCuArJFJB6oDxwvLdaghXmscNdQ4SNGjHCXuYYK//bbb4sNFe40e/Zs7r//fvf82rVrefjhh8nKyuL9998P2VGEUkp5CrZPYhHgulppGPCZlzrLgDtEpIHdYX2HXebvdu8DVtrZLrzCNFZ4oEOFO7322mvEx8czxDGg1PXXX8+WLVtYv349r776Knl5eUHFp5RSvgSbJKYAt4vITuCn9jwi0l1E/gxgjDkBvASstx+T7DJEZKqIZAO1RCRbRCba2/0LkCAiu4An8XLVVFiEaazwQIcKd5kzZw5Lly5l3rx5Xgfxa9u2LXXq1Ck2/pNSSoVSUOcpjDHHgV5eyjcAv3LMzwZme6n3NPC0l/I84BfBxFYmrm/rEyZYp5hatLASRIiGhXUNFX733Xfz61//2l3uGio8OTmZRo0a0bt3b5YuXcrUqVP54osvqOVIXHv37uWqq64iPj6erKwstm/fTmJiYkjiU0opT3oy21OYxwr3d6jw0aNHc/78efr370+VKlXo2bMn7777LqtWrWLKlClUrVqVKlWq8M477xT7XWyllAolTRIRUJahwnft2gXAmTNnqFu3rrvO0KFDGTp0aDjDVUqVQUZGbI5pqElCKaWCFMXfWQo7HeBPKaWCFMs/m6pJQimlguT6naW4uNj72VQ93aSUii1R6ByI5d9Z0iShlIodUewciNXfWdLTTUqp2BHLnQNRoknCaepUSEsrXpaWZpUHwTVUePv27enUqRPTpk2jqKjIvXzVqlX06NGDNm3a0KZNG2bNmuVe9sorr9CsWTP3UOOLFi0CYN++faSkpLjvu1iyZAlgDT9es2ZN93Dijz32mHtbd955J506daJ9+/Y89thjFBYWBvW6lCp3YrlzIEr0dJPTT34CAwfCwoWQkmIlCNd8EFzDcgAcOXKEwYMHk5OTw4svvsihQ4cYPHgw//jHP+jatSvHjh2jd+/eNGvWjJ/97GcAPPHEE4wdO5Zt27Zx8803c+TIEV5++WUGDhzIqFGj2Lp1K3369CEzMxOAq6++2t2e08KFC6lXrx7GGO677z4++ugjBg0aFNRrU6pcieXOgSjRIwmnlBQrIQwcCM8/XzxhhEjjxo2ZNWsWb7/9tnscp+HDh9O1a1cAGjVqxNSpU5kyZcol67Zt25b4+HiOHTuGiJCTkwPA6dOnufLKK0ttu169egBcuHCB/Px8r+NBKVXhJSXBuHGaIEJEk4SnlBQYNQpeesl6DmGCcGndujWFhYUcOXKELVu20K1bt2LLu3fvzpYtWy5Zb+3atVSpUoXLL7+ciRMnMnfuXJo3b06fPn1466233PX27t1Lly5duPXWW/n666+LbaN37940btyYunXrct9994X8tSmlYosmCU9paTBzJjz3nPXs2UcRBdOnT6dz586MHTuWDz/8EBFh/vz5DB8+nOzsbJYsWcLQoUMpKiqiadOm7Nu3j2+//ZY33njDfWrLZdmyZRw8eJDz58+zcuXKKL4qpVRFoEnCydkHMWnSxVNPIU4Ue/bsIS4ujsaNG9OuXTs2btxYbPnGjRtp3769e/6JJ55g06ZNfP311+5BAf/yl78w0P5J1aSkJPLy8jh27BjVq1cnISEBgG7dunH11Vfzww8/FNt+jRo16N+/P5995u3nP5RS6iJNEk7r1xfvg3D1UaxfH7Imjh49ymOPPcbo0aMREX7zm98wZ84cd0fz8ePHeeaZZ3j66UtGUC+mRYsWpKamArBt2zby8vK4/PLLOXr0qPuqpT179rBz505at25Nbm6u+3e2L1y4wBdffEGbNm1C9rqUUrFJr25y8vbBnJISdL+E65fpCgoKiI+PZ+jQoTz55JMANG3alLlz5/LII49w5swZjDH87ne/o2/fviVuc9q0aTzyyCNMnz4dEWHOnDmICF999RXPP/+8eyjxd999l4YNG3L48GH69evH+fPnKSoqIiUlpdjlsUop5Y0miQgo7X6EW265hfU+jlbGjx9fbKhwl3bt2vHvf//7kvIBAwYwYMCAS8qbNGnisw2llPJFTzcppZTySZOEUkopnzRJKKWU8kmThFJKKZ80SSillPJJk4RSSimfNElEgGuo8A4dOtC3b19OnToV0PrJycls2LAhPMEppVQJNElEgGuo8M2bN9OwYUNmzJgR7ZCUUjEkIwPmzWtBRkbot61JIsKSkpI4cOAAAOvWrSMpKYkuXbpwww03sGPHDsC6Q3vQoEG0bduWwYMHc+7cOff6o0aNonv37rRv354XXnjBXZ6YmMixY8cA2LBhA8n6YytKVQquX2ydPbsVvXoR8kRR6e64DsdvKBhj/KpXWFhIamoqI0aMAKBNmzZ8/fXXxMfHs2LFCsaPH88nn3zCzJkzqVWrFtu2bSMjI8M9qB/A5MmTadiwIYWFhfTq1YvvvvuOjh07hvw1KaUqBtcvthYVifsXW0P5UxqVLkn4+kA/c+aM1+EvQsE1dtOBAwdo27Ytt99+O2D9WNCwYcPYuXMnIkJBQQEAX331FY8//jgAHTp0KJYEFi5cyKxZs7hw4QIHDx5k69atmiSUqsRcv9h6/nwR1apVCfkvturppghw9UlkZWW5f40O4LnnniMlJYXNmzezePFi8vLyStzO3r17ef3110lNTeW7777jZz/7mXud+Ph49+9ml7YdpVTscP1i68MPZ5KaGvof5NMkEUG1atXizTffZNq0aVy4cIHTp0/TrFkzAObMmeOud8stt/DBBx8AsHXrVr777jsAcnJyqF27NvXr1+fw4cP885//dK+TmJjo/l2KTz75JEKvSCnllpEBr74a+k4BPyQlwZAh+8Lyi62aJCKsS5cudOzYkfnz5/P0008zbtw4unTpwoULF9x1Ro0aRW5uLm3btmXy5Mnunzft1KkTXbp0oU2bNgwePJgbb7zRvc4LL7zAb3/7W7p3705cXFzEX5dSlZqr9/i55whL73EUVbo+iWjIzc0tNr948WL3tPNX415++WXAOj21YMEC4NK+EucRh9PNN998yS/QKaUixNV7XFhIWHqPo0iPJJRSKliu3uO4OOs5hi5B1yMJpZQKlqv3OD3dShAxchQBmiSUUio0kpJiKjm46OkmpZRSPmmSUEop5ZMmCaWUUj4FlSREpKGILBeRnfZzAx/1htl1dorIMEf5ZBHZLyK5HvWHi8hREdlkP34VTJzR5hoqvFOnTnTt2pXVq1cDkJmZSc2aNenSpQtt27alR48exS5xnTNnDq1ataJz58507tyZBx98MEqvQClVWQXbcf0skGqMmSIiz9rzzzgriEhD4AWgO2CAjSKyyBhzElgMvA3s9LLtD40xo4OMr2wyMkJ6lYJrWA6AZcuWMW7cOP71r38BcPXVV/Ptt98CsGfPHu69916MMTz00EMA3HvvvcyaNSvoGJRSqiyCPd3UH3jPnn4PuNtLnd7AcmPMCTsxLAfuBDDGrDHGHAwyhtAK852TOTk5NGjg9YCL1q1b88Ybb/Dmm2+GtE2llCqrYI8kmjg+5A8BTbzUaQbsd8xn22WlGSAitwA/AE8YY/Z7qyQiI4GRAE2aNCE9Pb3Y8vr163PmzJlSGyssLOTMmTNUW7aMavn5SGEhJj+f/GXLyO/QwY9wfTt37hwdO3YkLy+Pw4cPs3jxYs6cOUNubi5FRUXF4rv22mvZvn07Z86cIS8vj7///e+sWbMGsIbr+OUvfxlULKXJy8u7ZB96k5ub61e9SNO4AldeY9O4AhO2uIwxJT6AFcBmL4/+wCmPuie9rD8W+L1j/jlgrEedXI/5BKC6Pf0osLK0OI0xdOvWzXjaunXrJWXe5OTkWBOrVxtTs6YxcXHW8+rVfq1fktq1a7unV69ebdq1a2eKiorM3r17Tfv27YvVPXHihKlRo4Yxxpi//vWv5pFHHgm6/UD4u7/S0tLCG0gZaVyBK6+xaVyBCSYuYIPx8bla6pGEMeanvpaJyGERaWqMOSgiTYEjXqodAJId882B9FLaPO6Y/TMwtbQ4QybMd04mJSVx7Ngxjh496nX5t99+S9u2bUPaplJKlVWwfRKLANfVSsOAz7zUWQbcISIN7Kuf7rDLfLITjks/YFuQcQYmKQnGjQvL3ZPbt2+nsLCQhISES5ZlZmYyduxYxowZE/J2lVKqLILtk5gCLBSREUAWMBBARLoDjxljfmWMOSEiLwHr7XUmGWNO2PWmAoOBWiKSDfzZGDMReFxE+gEXgBPA8CDjjCrXL9OBdXrvvffecw/nvXv3brp06UJeXh5169bl8ccfZ/jw4dELVimlHIJKEvZpoV5eyjcAv3LMzwZme6n3NPC0l/JxwLhgYitPCgsLvZYnJiZy7tw5n+sNHz6cAQMGhCsspZQqld5xrZRSyidNEkoppXzSJKGUUsonTRJKKaV80iShlFLKJ00SSimlfNIkEQF16tS5pGzHjh0kJyfTuXNn2rZty8iRI93LVq1aRY8ePWjTpg3dunUrNgrsxIkTadasGZ07d+baa6/l3nvvZevWrcW2vWnTJkSEpUuXhu9FKaUqBU0SXmRkwKuvhnwA2GIef/xxnnjiCTZt2sS2bdvcd1kfOnSIwYMH8+6777J9+3a+/PJL/vjHP/LFF1+413Wtt3PnTu6//35uu+22YsN8zJ8/n5tuuon58+eH7wUopSoFTRIewjxSuNvBgwdp3ry5e/66664DYMaMGQwfPpyuXbsCkJCQwNSpU5kyZYrX7dx///3ccccdfPDBB4B1R/dHH33EnDlzWL58OXl5eeF5AUqpSkGThIf0dMjPh8JC6zlcIwI/8cQT3Hbbbdx1111Mnz6dU6dOAbBlyxa6detWrG737t3ZsmWLz2117dqV7du3A7B69WpatWrF1VdfTXJycrEjEKWUCpQmCQ/JyVCtGsTFWc/JyeFp56GHHmLbtm384he/ID09nZ49e3L+/Pkybcsa6dcyf/58Bg0aBMCgQYP0lJNSKiiaJDy4Rgp/6SXrOQwDwbpdeeWVPPzww3z22WfEx8ezefNm2rVrx8aNG4vV27hxI+3bt/e5Hdfw4oWFhXzyySdMmjSJxMRExowZw9KlS/360SWllPJGk4QXYRwp3G3p0qUUFBQAVmf18ePHadasGb/5zW+YM2eO+zexjx8/zjPPPMPTT18yDiIAn3zyCV9++SUPPPAAqampdOzYkf3795OZmUlWVhYDBgzg008/Dd8LUUpFX0YGLebNC0snarBDhSs/nD17tlgn9ZNPPkl2dja//e1vqVGjBgCvvfYaV1xxBQBz587lkUce4cyZMxQWFvLkk0/St29f9/rTp09n7ty5/Pjjj3To0IGVK1dy+eWXM3/+fO65555ibQ8YMICZM2fy4IMPRuCVKqUizr7aptX58zBvXshPgWiSiICioiKv5W+88YbX8ltuuYX1662f3zhz5gx169Z1L5s4cSITJ070ut5f//rXS8r69etHv379AoxYKVVh2FfbSFHRxattQpgk9HSTUkpVZPbVNkVVqoTlahtNEkopVZHZV9tkPvxwWK62qRSnm4wxiEi0wyj3nJfSKqUqkKQk9p0/T+swXG0T80cSNWrU4Pjx4/oBWApjDMePH3d3pCulFFSCI4nmzZuTnZ1dbGwjb/Ly8srlB2Qk46pRo0axq7CUUirmk0TVqlVp1apVqfXS09Pp0qVLBCIKTHmNSylVOcT86SallFJlp0lCKaWUT5oklFJK+SSxdNWPiBwFssq4eiPgWAjDCRWNKzAaV+DKa2waV2CCiaulMeZybwtiKkkEQ0Q2GGO6RzsOTxpXYDSuwJXX2DSuwIQrLj3dpJRSyidNEkoppXzSJHHRrGgH4IPGFRiNK3DlNTaNKzBhiUv7JJRSSvmkRxJKKaV80iShlFLKp0qRJETkThHZISK7RORZL8uri8iH9vK1IpLoWDbOLt8hIr0jHNeTIrJVRL4TkVQRaelYVigim+zHogjHNVxEjjra/5Vj2TAR2Wk/hkU4rumOmH4QkVOOZeHcX7NF5IiIbPaxXETkTTvu70Skq2NZWPaXHzENsWP5XkRWi0gnx7JMu3yTiGwIVUwBxJYsIqcdf6/nHctKfA+EOa6nHDFttt9TDe1lYdlnInKViKTZnwNbROS3XuqE9/1ljInpBxAH7AZaA9WA/wDtPOr8GnjXnh4EfGhPt7PrVwda2duJi2BcKUAte3qUKy57PjeK+2s48LaXdRsCe+znBvZ0g0jF5VF/DDA73PvL3vYtQFdgs4/lfYB/AgL0BNZGYH+VFtMNrraAu1wx2fOZQKMo7q9k4PNg3wOhjsujbl9gZbj3GdAU6GpP1wV+8PL/GNb3V2U4kugB7DLG7DHG5AMLgP4edfoD79nTHwO9RETs8gXGmPPGmL3ALnt7EYnLGJNmjDlrz64BIjGOtz/7y5fewHJjzAljzElgOXBnlOJ6AJgforZLZIz5CjhRQpX+wN+MZQ1wmYg0JYz7q7SYjDGr7TYhcu8tV9ul7S9fgnlvhjquiLy/jDEHjTHf2NNngG1AM49qYX1/VYYk0QzY75jP5tKd7K5jjLkAnAYS/Fw3nHE5jcD6tuBSQ0Q2iMgaEbk7RDEFEtcA+9D2YxG5KsB1wxkX9mm5VsBKR3G49pc/fMUezv0VCM/3lgG+FJGNIjIyCvEAJInIf0TknyLS3i4rF/tLRGphfdh+4igO+z4T6zR4F2Ctx6Kwvr9i/vckYoGI/BLoDtzqKG5pjDkgIq2BlSLyvTFmd4RCWgzMN8acF5FHsY7CbotQ2/4YBHxsjCl0lEVzf5VbIpKClSRuchTfZO+rxsByEdluf8uOlG+w/l65ItIH+AdwbQTbL01f4N/GGOdRR1j3mYjUwUpKvzPG5IRqu/6oDEcSB4CrHPPN7TKvdUQkHqgPHPdz3XDGhYj8FJgA9DPGnHeVG2MO2M97gHSsbxgRicsYc9wRy5+Bbv6uG864HAbhcSogjPvLH75iD+f+KpWIdMT6+/U3xhx3lTv21RHgU0J3itUvxpgcY0yuPb0EqCoijYjy/nIo6f0V8n0mIlWxEsQ8Y8zfvVQJ7/sr1B0t5e2BdbS0B+v0g6uzq71Hnd9QvON6oT3dnuId13sIXce1P3F1weqou9ajvAFQ3Z5uBOwkRB14fsbV1DF9D7DGXOwo22vH18CebhipuOx6bbA6ESUS+8vRRiK+O2J/RvGOxXXh3l9+xNQCq4/tBo/y2kBdx/Rq4M5Q7is/YrvC9ffD+rDdZ+87v94D4YrLXl4fq9+idiT2mf26/wb8Xwl1wvr+Cukfvrw+sHr/f8D6wJ1gl03C+nYOUAP4yP6nWQe0dqw7wV5vB3BXhONaARwGNtmPRXb5DcD39j/J98CICMf1KrDFbj8NaONY92F7P+4CHopkXPb8RGCKx3rh3l/zgYNAAdZ53xHAY8Bj9nIBZthxfw90D/f+8iOmPwMnHe+tDXZ5a3s//cf+G08I5b7yM7bRjvfXGhyJzNt7IFJx2XWGY13M4lwvbPsM6zSgAb5z/K36RPL9pcNyKKWU8qky9EkopZQqI00SSimlfNIkoZRSyidNEkoppXzSJKGUUsonTRJKKaV80iShlFLKp/8Pd2Sui3YrlRIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "h = .2\n", "x = np.linspace(0,2, num=int(2/.2)+1)\n", "y=exp(x)\n", "\n", "methods = ('RK45','RK23', 'DOP853','Radau','BDF', 'LSODA')\n", "colours = ('bo','ro','rx','k-','r.','b.')\n", "def m_y(t, x):\n", " return x\n", "\n", "results = []\n", "for c, meth in zip(colours, methods):\n", " y = solve_ivp(m_y, (0,2), (1,), method=meth, t_eval = x, max_step=h).y[0]\n", " y2 = exp(x)-y\n", " plt.plot(x,y2,c, linewidth=1, label=meth)\n", "\n", "\n", "plt.title(\"Comparison of ODE methods\")\n", "plt.legend()\n", "ax = plt.gca()\n", "ax.grid()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 9, "id": "17efc931", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAo4ElEQVR4nO3de5wU1Zn/8c/DJeI4iAI6IgRm1E1AUC4zJlED2mZdFeMlEVmjS0L0J3F+IRqzxqy43sVN2Bg3GpcsufzcCDoZMUaNGDXQxku8MChRhHhjuMpFiIKjolye3x9VM/YMzUz3dPWl4Pt+veo13aeqznm6uubp06eqq8zdERGR+OpS7ABERCQ3SuQiIjGnRC4iEnNK5CIiMadELiISc0rkIiIxp0QuBWdm55nZo8WOo5mZ7W1mD5rZJjO7p9jxFIqZ3WFmN0ZU17VmNjOKuiR7SuQxZmbnmlmDmTWZ2Roze9jMvljsuDri7rPc/Z+KHUeKcUAF0Mfdz063gJkdbmYPhMn+PTNLmtkxKfMrzczD96LJzNaZ2R/M7MQ29Swzsw9Tlmsys5/l9+WBmU00s6fy3Y4UhxJ5TJnZ94D/Am4iSEIDgf8GzihiWB0ys27FjiGNQcBr7r4t3UwzOxR4GngZqAIOBu4DHjWzo9ssvp+7lwPDgceA+8xsYptlTnP38pRpcoSvRfZE7q4pZhPQC2gCzm5nmb0IEv1b4fRfwF7hvOOBVcDlwHpgDXAmMBZ4Dfg7MCWlrmuB2cBvgfeAF4DhKfP/DXgznLcY+ErKvIkESfAWYCNwY1j2VDjfwnnrgc0EyXJYyuv8DfA2sBz4d6BLSr1PAT8G3gEagVPa2R5DgMeBd4FXgNPD8uuAj4Gt4Ta9IM26dwJz0pRPB54IH1cCDnRrs8xlwLqUuJcB/5jh+3wtcA8wM9y2LwOfAa4It9dK4J/a7Be/Ct/P1eG27hq+9i3A9vA1vhsufwdwO/BQWP9zwKEp9R0DzAc2hX+PSZlXBfw5XO8x4GfAzHBejzDmjeH2ng9UFPv/Zneeih6Apk68aXAysK1t0mizzPXAs8CBwAHAX4AbwnnHh+tfDXQHLgyT5V1AT2Ao8CFQFS5/bZjoxoXLXxYmzu7h/LMJeqldgH8G3gf6hfMmhm19B+gG7E3rRH4SsADYjyCpD0lZ9zfA/WFMlQQfMhek1Ls1jL0rUEvwgWVptkV34A1gCvAp4IQwAX025fXNbGdbrgW+maY8ESbHvdl1Ij8kLB8SPl9Gdol8S7iNuoXboxG4MuV9a0xZ/j7gf4B9wvf9eeBbKdvrqTb130GQbD8X1j8LqAvn9Sb4gJwQzvta+LxPOP8Z4CcEHYYx4fZsTuTfAh4EysL3phrYt9j/N7vzVLyG4dcEvYpFEdU3EHgUWELQK6ws9sbN47Y7D1jbwTJvAmNTnp8ELAsfH0+QqLuGz3uGyebzKcsvAM4MH18LPJsyrwtBr2/0LtpeCJwRPp4IrGgzvyWpECTV14AvEPZaw/KuBD3lw1PKvgU8nlLHGynzysLXcFCaeEYTJOPU+u8Grk15fe0l8m3AyWnKB4dt9mfXibxHWH5s+HwZYa84ZbpwF+1eCzyW8vy0cN2279t+BMNrHwF7pyz/NSDZdpunzL8D+GXK87HA38LHE4Dn2yz/TFjPwHCb7JMy7y4+SeTnE3Qcjiz2/8qeMhVzjPwOgp5lVH4D/Ke7DyHoYayPsO5SsxHo28F488EEwxHNlodlLXW4+/bw8Yfh33Up8z8EylOer2x+4O47CIZmDgYws6+b2UIze9fM3gWGAX3TrduWu88j+Fp+O7DezGaY2b7h+t3TvIb+Kc/XptTzQfgwNeZmBwMrw7h3VVd7NgD90pT3A3YQ9FR3pbmNv6eUnenu+6VMv2hn/bbvyYY071s5wTh/d2BNyvvwPwQ98/asTXn8AZ9sv7b7D3yyzQ4G3nH399vMa3Yn8AhQZ2Zvmdk0M+veQRySg6Ilcnd/gtY7N2Z2qJn90cwWmNmTZjY4k7rM7HCCntBjYd1NKf/Yu6NnCHpfZ7azzFsE/9zNBoZlnfXp5gdm1gUYALxlZoOAXwCTCb527wcsIhgmadbuJTbd/VZ3rwYOJxgD/j5B8tya5jWs7kTsbwGfDuPuTF1/Ihg+ams88EwH+9pXCDoVr2bYVmetJNgn+qZ8QOzr7kPD+dle5rTt/gOfbLM1wP5mtk+beUFD7lvd/Tp3P5xgnP3LwNezbF+yUGpnrcwAvhP+U19GcBZGJj4DvGtmvzOzF83sP82sa96iLDJ330Qwvn27mZ1pZmVm1t3MTjGzaeFidwP/bmYHmFnfcPlczvOtNrOvht8CvkuQNJ4lGI91gjF2zOybBD3yjJjZUWb2+bDH9j7BmPCOsNdZD0w1s57hB8b3OvkaniPobV4ebqfjCYYp6jJc/zrgGDObama9w3i+Q5CcfrCL11VhZpOBa4Ar2nwbiJy7ryEYWrzZzPY1sy5hx+i4cJF1wAAz+1SGVc4BPhOe4trNzP6Z4IP2D+6+HGgArjOzT4WnvJ7WvKKZJczsiPB/cDPBB3JeX/+ermQSuZmVE3x632NmCwm+FvYL533VzBalmR4JV+9GMA56GXAUwQGmiYV+DYXk7jcTJLZ/J0iiKwl6xb8PF7mR4J/tJYKzHV4IyzrrfoIDmc0HwL4a9rwWAzcTfEtYBxxBcJZKpvYl6NG/Q/D1fCPwn+G87xAk96UEZ6jcRXBsJSvu/jFBojmFoKf/38DX3f1vGa7/OvBFglMKlxH0SM8CTnL3tq/1XTN7n2CbjyU4s6htzA+2OY/8vmxf0y58neBg7mKC7TmbT4aE5hGcrbPWzDZ0VJG7byToSf8rwXtyOfBld29e91zg8wTfqq8hGNpsdlDY9maCY1Z/JhhukTwx9+LdWMLMKgk+4YeF46Kvunu6sciO6vkC8CN3Py58PgH4grt/O9KA91Bmdi1wmLv/S7FjEZGdlUyP3N03A41mdjaABYZnuPp8YD8zOyB8fgJBr0REZLdXtERuZncTfB3/rJmtMrMLCE6ru8DM/krwNTCjXymG46mXAXPN7GWCA23tnQkgIrLbKOrQioiI5K5khlZERKRzinIBo759+3plZWWn1n3//ffZZ599Ol6wwBRXdhRXdhRXdko1LsgttgULFmxw9wN2mlGMn5NWV1d7ZyWTyU6vm0+KKzuKKzuKKzulGpd7brEBDV5iP9EXEZEIKJGLiMScErmISMyVzN1atm7dyqpVq9iyZUu7y/Xq1YslS5YUKKrMFTKuHj16MGDAALp31wXlRKSEEvmqVavo2bMnlZWVmNkul3vvvffo2bNnASPLTKHicnc2btzIqlWrqKqqynt7IlL6SmZoZcuWLfTp06fdJC5gZvTp06fDby4i0tq0p6eRbEy2Kks2Jpn29LRdrBEfJZPIASXxDGk7iWTvqIOPYvzs8S3JPNmYZPzs8Rx18FFFjix3JTO0IiKST4mqBPXj6hk/ezy1NbVMb5hO/bh6ElWJYoeWs5LqkRdb165dGTFiBMOGDeO0007j3XffBWDZsmUMG/bJvRJ+8YtfUF1dzTvvfHKHr9tuuw0zY8OG4HLNjz/+OL169WLEiBGMGDGC66+/vqCvRUR2lqhKUFtTyw1P3EBtTe1ukcQhxol81iyorIQuXYK/s2blXufee+/NwoULWbRoEb179+b222/faZk777yT2267jUceeYT9998fgJUrVzJ37lwGDhzYatnRo0ezcOFCFi5cyNVXX517gCKSk2RjkukN07lqzFVMb5i+05h5XMUykc+aBZMmwfLl4B78nTQpmmTe7Oijj2b16ta3dKyvr+eHP/whjz76KH37fnJv4UsvvZQbbrhBY9ciJax5TLx+XD3XJ65vGWbZHZJ5LBP5lVfCB21ud/vBB0F5FLZv387cuXM5/fTTW8qWL1/O5MmTefTRRznooINayu+//3769+/PEUccsVM9zzzzDMOHD+eUU07hlVdeiSY4EemU+W/NbzUm3jxmPv+t+UWOLHexPNi5YkV25Zn68MMPGTFiBKtXr2bIkCGceOKJLfMOOOAAevfuTX19PZdeeikAH3zwATfddBOPPvroTnWNGjWK5cuXU15ezpw5czjzzDN5/fXXcwtQRDrt8mMv36ksUZXYLcbJY9kjbzMU3WF5pprHyJcvX467txojLysrY86cOfz85z9nVjiG8+abb9LY2Mjw4cMZNmwYq1atYtSoUaxdu5Z9992X8vJyAMaOHcvWrVtbDoSKiEQplol86lQoK2tdVlYWlEehrKyMW2+9lZtvvplt27a1lB944IH88Y9/ZMqUKTzyyCMcccQRrF+/nmXLlrFo0SIGDBjACy+8wEEHHcTatWvx8O5Lzz//PDt27KBPnz7RBCgikiKWify882DGDBg0CMyCvzNmBOVRGTlyJEceeSR33313q/KqqioeeOABzj//fJ5//vldrj979myGDRvG8OHDufjii6mrq9PBUBHJi1iOkUOQtKNM3ABNTU2tnj/44IMtjxctWtTyePjw4Tud0QLB+ebNJk+ezOTJk6MNUEQkjUh65Ga2n5nNNrO/mdkSMzs6inpFRKRjUfXIfwr80d3HmdmngLKOVhARkWjknMjNrBcwBpgI4O4fAx/nWq+IiGQmiqGVKuBt4P+Z2Ytm9kszK83bV4uI7Ias+RS5TldgVgM8Cxzr7s+Z2U+Bze5+VZvlJgGTACoqKqrr6upa1dOrVy8OO+ywDtvbvn07Xbt2zSnmfCh0XG+88QabNm3qcLmmpqaW89lLieLKjuLKTqnGBbnFlkgkFrh7zU4z3D2nCTgIWJbyfDTwUHvrVFdXe1uLFy/eqSydzZs3Z7RcoRU6rky3VzKZzG8gnaS4sqO4slOqcbnnFhvQ4Glyas5DK+6+FlhpZp8Ni74ELM613mJovozt0KFDGT58ODfffDM7duxomf/UU0/xuc99jsGDBzN48GBmzJjRMu+mm26if//+LZfBfeCBBwBYsWIFiUSi5bz0OXPmAMGpinvvvXfLZW4vuuiilrpOPvlkhg8fztChQ7nooovYvn17gbaAiMRRVGetfAeYFZ6xshT4ZkT17tqsWcFVslasCH6bP3VqzieWN/9EH2D9+vWce+65bN68meuuu461a9dy7rnn8vvf/55Ro0axYcMGTjrpJPr378+pp54KBFdBvOyyy1iyZAmjR49m/fr13HjjjYwfP57a2loWL17M2LFjW843P/TQQ1vaS1VfX8++++6LuzNu3DjuuecezjnnnJxem4jsviI5j9zdF7p7jbsf6e5nuvs7Ha+VgwJcx/bAAw9kxowZ/OxnP2u57srEiRMZNWoUAH379mXatGn88Ic/3GndIUOG0K1bNzZs2ICZsXnzZgA2bdrEwQcf3GHb++67LwDbtm3j448/1i9CRaRdsfyJft6vYxs65JBD2L59O+vXr+eVV16hurq61fyampq0l6d97rnn6NKlCwcccADXXnstM2fOZMCAAYwdO5bbbrutZbnGxkZGjhzJcccdx5NPPtmqjpNOOokDDzyQnj17Mm7cuEhfl4jsXuKZyPN1Hdsc3XLLLYwYMYLLLruM3/72t5gZd999NxMnTmTVqlXMmTOHCRMmsGPHDvr168eKFSt48cUX+clPftIyjNPskUceYc2aNXz00UfMmzeviK9KREpdPBN5vq5j28bSpUvp2rUrBx54IIcffjgLFixoNX/BggUMHTq05fmll17KwoULefLJJxk9ejQAv/rVrxg/fjwQ3HVoy5YtbNiwgb322qvlaojV1dUceuihvPbaa63q79GjB2eccQb3339/pK9LRHYv8Uzk+b6OLfD2229z0UUXMXnyZMyMb3/729xxxx0tByc3btzID37wAy6/fOeL1acaOHAgc+fOBWDJkiVs2bKFAw44gLfffrvlbJSlS5fy+uuvc8ghh9DU1MSaNWuAYIz8oYceYvDgwZG9LhHZ/cTz6ofNZ6dEfNZK8x2Ctm7dSrdu3ZgwYQLf+973AOjXrx8zZ87kwgsv5L333sPd+e53v8tpp53Wbp0333wzF154Ibfccgtmxh133IGZ8cQTT3D11VfTvXt3unTpws9//nN69+7NunXrOP300/noo4/YsWMHiUSi1amJIiJtxTORQ16uY9vR+dpjxoxh/vz09/ebMmUKPXv23Kn88MMP5+mnn96p/KyzzuKss87aqbyiomKXbYiIpBPPoRUREWmhRC4iEnNK5CIiMadELiIFNe3paSQbk63Kko1Jpj09rUgRxZ8SuYgU1FEHH8X42eNbknmyMcn42eM56uCjihxZfMX3rBURiaVEVYL6cfWMnz2e2ppapjdMp35cPYmqRLFDiy31yFM0X8Z22LBhnHbaabz77rtZrX/88cfT0NCQn+BEdiOJqgS1NbXc8MQN1NbUKonnKJaJPF9jbM2XsV20aBG9e/fm9ttvz6k+EUkv2ZhkesN0rhpzFdMbpu/0/yzZiWUiL8QY29FHH83q1asBeP755zn66KMZOXIkxxxzDK+++ioQ/BL0nHPOYciQIZx77rl8+OGHLevX1tZSU1PD0KFDueaaa1rKKysr2bBhAwANDQ0cf/zxkcUsEgfN/6/14+q5PnF9yzCLknnnxXKMPN9jbNu3b2fu3LlccMEFAAwePJgnn3ySbt268ac//YkpU6Zw7733Mn36dMrKyliyZAnPPPNMy4WyAKZOnUrv3r3Zvn07X/rSl3jppZc48sgjI4lPJM7mvzW/1f9r8//z/Lfma4ilk2KZyKH1GNtVY66KZAdovtbK6tWrGTJkCCeeeCIQ3BDiG9/4Bq+//jpmxtatWwF44oknuPjiiwEYNmxYq0RdX1/PjBkz2LZtG2vWrGHx4sVK5CLA5cfufKG5RFVCSTwHsRxagfyMsTWPkS9fvrzlrkAAV111FYlEgkWLFvHggw+yZcuWdutpbGzkxz/+MXPnzuWll17i1FNPbVmnW7duLfcB7ageEZFMRJLIzWyZmb1sZgvNLO+nbeR7jK2srIxbb72Vm2++mW3btrFp0yb69+8PwB133NGy3JgxY7jrrrsAWLx4MS+99BIAmzdvZp999qFXr16sW7eOhx9+uGWdysrKluua33vvvZHEKyJ7tih75Al3H+HuNRHWmVZ7Y2xRab7r/d13383ll1/OFVdcwciRI9m2bVvLMrW1tTQ1NTFkyBCmTp3aciu44cOHM3LkSAYPHsy5557Lscce27LONddcwyWXXEJNTQ1du3aNLF4R2XPFcow8X2NsTU1NrZ4/+OCDLY9T795z4403AsFQTF1dHQDvvfdeq8vYpvbcU40ePXqnOwGJiOTC3D33SswagXcAB/7H3WekWWYSMAmgoqKiujkBNuvVqxeHHXZYh21t3769JHuyhY7rjTfeYNOmTR0u19TURHl5eQEiyo7iyo7iyk6pxgW5xZZIJBakHfVw95wnoH/490Dgr8CY9pavrq72thYvXrxTWTqbN2/OaLlCK3RcmW6vZDKZ30A6SXFlR3Flp1Tjcs8tNqDB0+TUSMbI3X11+Hc9cB/wuU7WE0U4uz1tJxFJlXMiN7N9zKxn82Pgn4BF2dbTo0cPNm7cqCTVAXdn48aN9OjRo9ihiEiJiOJgZwVwn5k113eXu/8x20oGDBjAqlWrePvtt9tdbsuWLSWZxAoZV48ePRgwYEBB2hKR0pdzInf3pcDwXOvp3r07VVVVHS73+OOPM3LkyFybi1ypxiUiu7/Y/rJTREQCSuQiIjGnRC4iEnNK5CIiMadELiISc0rkIiIxp0QuIhJzSuQie6B83cBcikOJXGQPVIgbmEvhxPJ65CKSm3zfwFwKSz1ykT1U6g3Ma2tqlcRjTIlcZA+VjxuYS3EokYvsgfJ9A3MpLCVykT1QIW5gLoWjg50ie6B83cBcikM9chGRmFMiFxGJOSVyEZGYiyyRm1lXM3vRzP4QVZ0iItKxKHvklwBLIqxPREQyEEkiN7MBwKnAL6OoT0REMmfunnslZrOB/wB6Ape5+5fTLDMJmARQUVFRXVdX16m2mpqaKC8vzyHa/FBc2VFc2VFc2SnVuCC32BKJxAJ3r9lphrvnNAFfBv47fHw88IeO1qmurvbOSiaTnV43nxRXdhRXdhRXdko1LvfcYgMaPE1OjWJo5VjgdDNbBtQBJ5jZzAjqFRGRDOScyN39Cncf4O6VwDnAPHf/l5wjExGRjOg8chGRmIv0Wivu/jjweJR1iohI+9QjFyki3TtToqBELlJEunemREGXsRUpIt07U6KgHrlIkenemZIrJXKRItO9MyVXSuQiRaR7Z0oUlMhFikj3zpQo6GCnSBHp3pkSBfXIRURiTolcRCTmlMhFRGJOiVxEJOaUyEVEYk6JXEQk5pTIRURiTolcRCTmlMhFRGIu50RuZj3M7Hkz+6uZvWJm10URmEih6OYOEndR9Mg/Ak5w9+HACOBkM/tCBPWKFIRu7iBxl/O1Vtzdgabwafdw8lzrFSmU1Js7nNL3FB6e/7Bu7iCxYkEezrESs67AAuAw4HZ3/0GaZSYBkwAqKiqq6+rqOtVWU1MT5eXlOUSbH4orO6UY168bf82dK+5kwsAJnF91frHDaaUUtxcors7IJbZEIrHA3Wt2muHukU3AfkASGNbectXV1d5ZyWSy0+vmk+LKTqnFNW/pPO87ra9P+PUE7zutr89bOq/YIbVSaturmeLKXi6xAQ2eJqdGetaKu78bJvKTo6xXJJ9Sb+5wftX5urmDxE4UZ60cYGb7hY/3Bk4E/pZrvSKFops7SNxFcWOJfsD/huPkXYB6d/9DBPWKFIRu7iBxF8VZKy8BIyOIRUREOkG/7BQRiTklchGRmFMiFxGJOSVyEZGYUyIXEYk5JXIpGboKoUjnKJFLydBVCEU6J4ofBIlEIvUqhLU1tUxvmK6rEIpkQD1yKSmJqgS1NbXc8MQN1NbUKomLZECJXEpKsjHJ9IbpXDXmKqY3TNeFq0QyoEQuJSP1KoTXJ67XVQhFMqRELiVDVyEU6Rwd7JSSoasQinSOeuQiIjGnRC4iEnNK5CIiMadELiISc0rk0oqudyISP1HcfPnTZpY0s8Vm9oqZXRJFYFIcut6JSPxEcfrhNuBf3f0FM+sJLDCzx9x9cQR1S4Hpeici8ZNzj9zd17j7C+Hj94AlQP9c65Xi0fVOROLF3D26yswqgSeAYe6+uc28ScAkgIqKiuq6urpOtdHU1ER5eXmOkUZvd4rrxXde5Lol13F6v9N5YM0DXDPkGkbuP7LocRWC4sqO4speLrElEokF7l6z0wx3j2QCyoEFwFc7Wra6uto7K5lMdnrdfNpd4pq3dJ73ndbX5y2dl/Z5seIqFMWVHcWVvVxiAxo8TU6N5KwVM+sO3AvMcvffRVGnFIeudyISPzkf7DQzA34FLHH3n+QekhSTrnciEj9R9MiPBSYAJ5jZwnAaG0G9IiKSgSjOWnnK3c3dj3T3EeE0J4rg9mT6YY6IZEq/7CxR+mGOiGRK1yMvUfphjohkSj3yEqYf5ohIJpTIS5huRCwimVAiL1G6EbGIZEqJvETphzkikikd7CxR+mGOiGRKPfIO6HxuESl1SuQd0PncIlLqNLTSAZ3PLSKlTj3yDOh8bhEpZUrkGdD53CJSypTIO6DzuUWk1CmRd0Dnc4tIqdPBzg7ofG4RKXXqkYuIxJwSuYhIzCmRi4jEXCSJ3Mx+bWbrzWxRFPWJiEjmouqR3wGcHFFdIiKShUgSubs/Afw9irpERCQ75u7RVGRWCfzB3YftYv4kYBJARUVFdV1dXafaaWpqory8vLNh5o3iyo7iyo7iyk6pxgW5xZZIJBa4e81OM9w9kgmoBBZlsmx1dbV3VjKZ7PS6+aS4sqO4sqO4slOqcbnnFhvQ4Glyqs5aERGJOSVyEZGYi+r0w7uBZ4DPmtkqM7sginpFRKRjkVxrxd2/FkU9IiKSPQ2tiIjEnBK5iEjMKZGLiMScErmISMwpkYuIxJwSuYhIzCmRi4jEnBK5iEjMKZGLiMScErmISMwpkYuIxJwSuYhIzCmRi4jEnBK5iEjMKZGLiMScErmISMwpkYuIxFxUt3o72cxeNbM3zOzfoqgz1bSnp3HrJVeyqlslYxInsKpbJbdeciXTnp4WdVPpzZoFlZXQpUvwd9Ystbu7tR22e9wJJ+xxr3lPfJ+L+Zrzso+5e04T0BV4EzgE+BTwV+Dw9taprq72bPz04ine5/v4vErcCf72+T7+04unZFVPp8yc6V5W5g6fTGVlQXmKZDJZlHY7knVcEbXbqbgK1HbJtJtF27vN/hVh21nHFYP3uSNAg6fLw+kKs5mAo4FHUp5fAVzR3jrZJvKVXQf5vEq87/fxqxLB33mV+Mqug7Kqp1MGDWq98ZunQa3bjvwfLcN2O5J1XBG126m4CtR2ybSbRdu7zf4VYdtZxxWD97kju0rkUdx8uT+wMuX5KuDzbRcys0nAJICKigoef/zxjBsYs30FA5ZBbQPccBxc9WdILIMdrMiqns44bvkKLE25L1/Bn1PabmpqijSW41bsot0VrdvtSLZxRdVuR9LFVai22ypWu9m0vbvsX1G23Z5S2r8K0na67J7NBIwDfpnyfALws/bWiVOPfGXXQWk/Sdu2HXWP6b0+6dt9r8+gjlbNLS71yEuyp6YeeXZKav+KsG120SOP4mDnauDTKc8HhGWR+d23z+Pss6H+Hrg+Gfw9++ygPN9+sH0q71PWqux9yvjB9ql5bXcK6dudQn7bfWps+nafGpvfdgGYOhXKWrdNWVlQvju2W8y29ZoL124h2k6X3bOZgG7AUqCKTw52Dm1vnWx75D966kf+04un+Mqug3w75iu7DvKfXjzFf/TUj7KqpzMGDXL/GjO9kaDtRgb515i50wdp1D0ms/TtmmVXT7ZxZfp6c7WruJ6sndnqfX6ytgAHotyDg06DBvkOs2AjFOIAWJu2vZ22I++RZ9huRzodVwRtt2eXceW53XZFsI+Rr4OdQd2MBV4jOHvlyo6WzzaRp8rLDt2OTA82Rx1XVN8Cs43LLH272X6AdCauYp5U0F5cpUBxZadU43LPLbZdJfJIziN39znu/hl3P9TdC/A9pXDOOw9mzIBBg8As+DtjRlCeT8X6FjhwYHblUbrySvjgg9ZlH3wQlOdb8+nFJ5xwXMFPLxbJlX7ZmYHzzoNly2DHjuBvvpN4c5t70gcIwIoV2ZVHZdYsmDQJli8Hd2P58uC5krnEhRJ5CduTPkCgeN8GivlNQCQKSuSyk2J8gEDxvg0U65tAs2L+alx2D0rkUjKK9W2gmMcFWg/roGEd6RQlcikpxfg2UMzjAhrWkSgokcser/U3AS/ocYFiD+vI7kGJXIRPvgnMm/fngh4XKPawjsbmdw9K5CJFVKxhHY3N716UyEWKqFgHeDU2X1jTnp5GsjHZqizZmIzs5jhK5CJFVowDvBqbL6yjDj6K8bPHtyTzZGOS8bPHc9TBR0VSfxTXIxeRmBk4MBhOSVcu0UtUJagfV8/42eM5pe8pPDz/YerH1ZOoSkRSv3rkInugYp5yCXvmgdZEVYLamlruXHEntTW1kSVxUCIX2SMV81IMe+qB1mRjkukN05kwcALTG6bvNGaeCyVykT1UsS7FsCceaE02Jjlj5ni63lvPzAv+l6731nPGzPGRJXMlchEpqD3xQOuMh+bz8V31rHsugbux7rkEH99Vz4yH5kdSvxK5iBRUMX8EVSzP/PhyPvpb6zHxj/6W4JkfXx5J/UrkIlJQxT7QWgz5/haiRC4iBVXsA63FOFsm399CckrkZna2mb1iZjvMrCaakERkd1eMA63FPFsm399Ccu2RLwK+CjwRQSwiInlTzLNl8n2FzZwSubsvcfdXowlFRCR/in22TD6vsGnunnslZo8Dl7l7QzvLTAImAVRUVFTX1dV1qq2mpibKy8s7tW4+Ka7sKK7sKK7spIvrnHO+wLp1PXZatqJiC3V1zxYqtJy2WSKRWODuOw9ju3u7E/AngiGUttMZKcs8DtR0VFfzVF1d7Z2VTCY7vW4+Ka7sKK7sKK7spItr5kz3sjL3YIQ8mMrKgvJix5YpoMHT5NQOL5rl7v/YqY8OEZES0jyUceWVwXDKwIHBwcZC/aI1n3T1QxHZY5x33u6RuNvK9fTDr5jZKuBo4CEzeySasEREJFM59cjd/T7gvohiERGRTtAvO0VEYk6JXEQk5pTIRURiLpIfBGXdqNnbQJo7BmakL7AhwnCioriyo7iyo7iyU6pxQW6xDXL3A9oWFiWR58LMGjzdL5uKTHFlR3FlR3Flp1TjgvzEpqEVEZGYUyIXEYm5OCbyGcUOYBcUV3YUV3YUV3ZKNS7IQ2yxGyMXEZHW4tgjFxGRFErkIiIxV1KJ3MxONrNXzewNM/u3NPP3MrPfhvOfM7PKlHlXhOWvmtlJBY7re2a22MxeMrO5ZjYoZd52M1sYTg8UOK6JZvZ2Svv/J2XeN8zs9XD6RoHjuiUlptfM7N2UeXnZXmb2azNbb2aLdjHfzOzWMOaXzGxUyrx8bquO4jovjOdlM/uLmQ1PmbcsLF9oZru8qUue4jrezDalvFdXp8xr9/3Pc1zfT4lpUbg/9Q7n5XN7fdrMkmEeeMXMLkmzTP72sXQXKS/GBHQF3gQOAT4F/BU4vM0y/xf4efj4HOC34ePDw+X3AqrCeroWMK4EUBY+rm2OK3zeVMTtNRH4WZp1ewNLw7/7h4/3L1RcbZb/DvDrAmyvMcAoYNEu5o8FHgYM+ALwXL63VYZxHdPcHnBKc1zh82VA3yJtr+OBP+T6/kcdV5tlTwPmFWh79QNGhY97Aq+l+X/M2z5WSj3yzwFvuPtSd/8YqAPOaLPMGcD/ho9nA18yMwvL69z9I3dvBN4I6ytIXO6edPfm27o+CwyIqO2c4mrHScBj7v53d38HeAw4uUhxfQ24O6K2d8ndnwD+3s4iZwC/8cCzwH5m1o/8bqsO43L3v4TtQuH2rUy2167ksl9GHVdB9i0Ad1/j7i+Ej98DlgD92yyWt32slBJ5f2BlyvNV7LwhWpZx923AJqBPhuvmM65UFxB86jbrYWYNZvasmZ0ZUUzZxHVW+DVutpl9Ost18xkX4RBUFTAvpThf26sju4o7n9sqW233LQceNbMFFtwTt9CONrO/mtnDZjY0LCuJ7WVmZQTJ8N6U4oJsLwuGfEcCz7WZlbd9THcIipCZ/QtQAxyXUjzI3Veb2SHAPDN72d3fLFBIDwJ3u/tHZvYtgm8zJxSo7UycA8x29+0pZcXcXiXLzBIEifyLKcVfDLfVgcBjZva3sMdaCC8QvFdNZjYW+D3wDwVqOxOnAU+7e2rvPe/by8zKCT48vuvum6Osuz2l1CNfDXw65fmAsCztMmbWDegFbMxw3XzGhZn9I3AlcLq7f9Rc7u6rw79LCW5SPbJQcbn7xpRYfglUZ7puPuNKcQ5tvvrmcXt1ZFdx53NbZcTMjiR4/85w943N5Snbaj3BDV6iGk7skLtvdvem8PEcoLuZ9aUEtleovX0rL9vLzLoTJPFZ7v67NIvkbx/Lx8B/Jw8WdCMY5K/ik4MkQ9ss821aH+ysDx8PpfXBzqVEd7Azk7hGEhzg+Yc25fsDe4WP+wKvE9GBnwzj6pfy+CvAs/7JwZXGML79w8e9CxVXuNxggoNPVojtFdZZya4P3p1K6wNRz+d7W2UY10CCYz7HtCnfB+iZ8vgvwMkFjOug5veOICGuCLddRu9/vuIK5/ciGEffp1DbK3ztvwH+q51l8raPRbZxI9oYYwmO9r4JXBmWXU/QywXoAdwT7tjPA4ekrHtluN6rwCkFjutPwDpgYTg9EJYfA7wc7swvAxcUOK7/AF4J208Cg1PWPT/cjm8A3yxkXOHza4Eftlkvb9uLoHe2BthKMAZ5AXARcFE434Dbw5hfBmoKtK06iuuXwDsp+1ZDWH5IuJ3+Gr7HVxY4rskp+9azpHzQpHv/CxVXuMxEgpMfUtfL9/b6IsEY/Esp79XYQu1j+om+iEjMldIYuYiIdIISuYhIzCmRi4jEnBK5iEjMKZGLiMScErmISMwpkYuIxNz/BxysuzxG/JkKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "h = .2\n", "x = np.linspace(0,2, num=int(2/.2)+1)\n", "y=exp(x)\n", "\n", "methods = ('RK45', 'DOP853','Radau')\n", "colours = ('bo','ro','gx')\n", "def m_y(t, x):\n", " return x\n", "\n", "results = []\n", "for c, meth in zip(colours, methods):\n", " y = solve_ivp(m_y, (0,2), (1,), method=meth, t_eval = x, max_step=h).y[0]\n", " y2 = exp(x)-y\n", " plt.plot(x,y2,c, linewidth=1, label=meth)\n", "\n", "\n", "plt.title(\"Comparison of ODE methods\")\n", "ax = plt.gca()\n", "plt.grid()\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "85c088dc", "metadata": {}, "source": [ "Now look at y'' = -y for which we expect solutions of sin(x) or cos(x).\n", "We set this up as a coupled system of first order equations as:\n", "y' = -u and u' = y" ] }, { "cell_type": "code", "execution_count": 5, "id": "268b884a", "metadata": {}, "outputs": [], "source": [ "def m_y(t,x):\n", " # x is a vector for (u, y)\n", " y = (x[1], -x[0])\n", " #print(x,y)\n", " return y\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "cbee4d08", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABld0lEQVR4nO2deXwV1fXAv+clJOwh7EtIgIgQkkBYBBEXKKgIAq07UpdWi9raVvuzqWJd6m5ata3aWle0UtRqFSyoFQRUkFVDNmTNyr4kgQCB5L37+2NmksnjZX3LvJc3389nPu/NnTv3nndn3py59557jiilsLGxsbEJXxxWC2BjY2NjYy22IrCxsbEJc2xFYGNjYxPm2IrAxsbGJsyxFYGNjY1NmGMrAhsbG5swx1YENl4jInNE5H9Wy2EgIu1E5GMRKReRf1stT6AQkfki8piPynpYRN72RVm+RkQ+EZGbmpH/NhH5cxPyRYvI9yLSwysBQxBbEQQRInK9iGwUkQoR2avf8OdbLVdjKKUWKKUusVoOE1cBvYBuSqmrPWUQkWEislhXFsdEZIWInGc6PkBElH4tKkRkv4j8V0QudiunQEROmvJViMgL/v15ICI3i8jX/q4nGFFKXaaUerMpeUUkCvg98McmlHsKeB241zsJQw9bEQQJIvIb4M/AE2gPsXjgb8AsC8VqFBGJtFoGDyQA25RS1Z4OikgisBrIBgYCfYEPgf+JyHi37F2UUh2BEcDnwIcicrNbnhlKqY6m7U4f/hYb75gFfK+U2t3E/P8CbhKRaD/KFHwopezN4g2IASqAqxvIE42mKPbo25+BaP3YRKAESAcOAHuBHwLTgG3AEWCeqayHgfeBd4FjwLfACNPxe4Gd+rE84EemYzejPUSfAw4Dj+lpX+vHRT92ADiK9rBNMf3Ot4CDQCHam5rDVO7XwJ+AUiAfuKyB9kgCVgJlQC4wU0//A3AaqNLb9BYP5/4TWOoh/e/Al/r3AYACIt3y3APsN8ldAExp4nV+GPg38LbettnA2cB9ensVA5e43Rev6ddzt97WEfpvrwSc+m8s0/PPB14ElujlrwMSTeWdB2wAyvXP80zHBgKr9PM+B14A3taPtdVlPqy39wagV3OuS1Pkcyun3jr18m9tyn2D9ob/e9P+tXqezvr+ZcA+oIcpz3bgIqufC4HcLBfA3hTAVKDa/aHjlucRYC3QE+gBrAEe1Y9N1M9/EGgD/AztYfsvoBOQDJwEBur5H0Z7UF6l579H/3O00Y9fjfaW7ND/OMeBPvqxm/W6fglEAu2oqwguBTYBXdCUQpLp3LeARbpMA9CU1C2mcqt02SOAO9AUnnhoizbADmAeEAX8QH+wDDH9vrcbaMt9wE88pE9Ce7i2o35FMEhPT9L3C2ieIqjU2yhSb4984H7Tdcs35f8Q+AfQQb/u64HbTO31tVv589EenGP18hcA7+jHuqI9KG/Qj83W97vpx78BnkV74bhQb09DEdwGfAy016/NaPQHaTOvS73yeSir3jo5UxHUe9+gKZCr3cpeoMvSTc97udvxxcCvrH4uBHKzXAB7UwBzgH2N5NkJTDPtXwoU6N8noj3oI/T9TvrDapwp/ybgh/r3h4G1pmMOtLfOC+qpOxOYpX+/GShyO17zUNL//NuAc9HfmvX0CLQ39WGmtNuAlaYydpiOtdd/Q28P8lyA9jA3l78QeNj0+xpSBNXAVA/pQ/U6+1G/Imirp0/Q9wvQ38pN28/qqfdh4HPT/gz9XPfr1gVtePAU0M6Ufzawwr3NTcfnA6+a9qehDYuApgDWu+X/Ri8nXm+TDqZj/6JWEfwU7cVjeCP3aGPXpV75PJRVb52cqQjqvW/Q3u6nup3fBShC65H9w0P5C4AHffHfDpXNniMIDg4D3RsZb++LNpxiUKin1ZShlHLq30/qn/tNx08CHU37xcYXpZQLbWipL4CI3CgimSJSJiJlQArQ3dO57iilvkAbVngROCAiL4tIZ/38Nh5+Qz/T/j5TOSf0r2aZDfoCxbrc9ZXVEIeAPh7S+wAutDfl+jDqOGJK+6FSqotpe6WB892vySEP160j2jxHG2Cv6Tr8A61n0BD7TN9PUNt+7vcP1LZZX6BUKXXc7ZjBP4HPgHdEZI+IZIhIGw91N+W61CefO02ts06ZHu6bUjQFiylPGdoQXQrwjIfyOqEp9LDBVgTBwTdob38/bCDPHrSHg0G8ntZS+htfRMQBxAF7RCQBeAW4E23YoAuQgzbMY6AaKlgp9Vel1GhgGNoY+G/RHr5VHn5DUyfxzOwB+utyt6SsZWjDX+5cA3xjeph44kdo4/lbm1hXSylGuye6mxRMZ6VUsn68wWvgAff7B2rbbC8QKyId3I5pFSlVpZT6g1JqGNo8w+XAjfXU4c11qaEZdTZGFto9WIOIpKH1OBYCf/VwThKwuQV1hSy2IggClFLlaOP7L4rID0WkvYi0EZHLRCRDz7YQ+L2I9BCR7np+b+y8R4vIFXov5C60h85atPFohTbHgIj8BO3NqUmIyDkiMk5/ezuONibu0t963wMeF5FOusL5TQt/wzq0t8l0vZ0mog2zvNPE8/8AnCcij4tIV12eX6I9aH5Xz+/qJSJ3Ag8B97m99focpdRe4H/AMyLSWUQcIpIoIhfpWfYDcbp5ZFNYCpytmyhHisi1aIr6v0qpQmAj8AcRidJNlmcYJ4rIJBFJFZEINAOAKrSekzveXpcamlFnYywFjDZDRIxJ6HnAT4B+IvJz0/F+aPMpa1tQV8hiK4IgQSn1DNqD8fdoD+FitLfyj/Qsj6H9WbPQxja/1dNayiK0iWBjAvEK/S0sD627/A3awyYVzUqoqXRG61GUog0LHKbWhvuXaMphF5qlx7/QrDqahVLqNNoD5jK0nsbfgBuVUt838fztwPloJqEFaG/EVwKXKqXcf2uZiBxHa/NpaBOP7jJ/7LaO4MPm/qZ6uBFt0jUPrT3fp3ZI6ws0q5x9InKosYKUUofR3qr/D+2apKNNkhrnXg+MQxvyeghtItugt173UWALmnXRPz3U4dV1caNJdTaBj4GhImIMoz6JNnz1d6WtG/gx8JiIDNaPXw+8qR8LG4yZdZswQkQeBs5SSv3YallsbPyNiMxFM1K4q5F80WhDQhcqpQ4EQrZgIRgXA9nY2Nj4DKXUy03MdwrNcizssIeGbGxsbMIce2jIxsbGJsyxewQ2NjY2YU5IzhF0795dDRgwoEXnHj9+nA4dOjSe0SKCXT4IfhmDXT4Ifhlt+bwnGGXctGnTIaXUmW62rV7a3JJt9OjRqqWsWLGixecGgmCXT6nglzHY5VMq+GW05fOeYJQR2KhsFxM2NjY2Nu7YisDGxsYmzLEVgY2NjU2YE5KTxTY2NuFFVVUVJSUlVFZWAhATE8OWLVsslqphrJSxbdu2xMXF0aZNfQ5b62IrAhsbm6CnpKSETp06MWDAAESEY8eO0alTp8ZPtBCrZFRKcfjwYUpKShg4cGCTzvHJ0JCIvC4iB0Qkp57jIiJ/FZEdIpIlIqNMx24Ske36dpMv5LGMjAyyZ91PSeQAXOKgJHIA2bPuh4yMxs+18T0ZGSy7fwUDBoDDAQMGwLL7V9jXIwSprKykW7duiEjjmcMcEaFbt241vaem4Ks5gvlo4Rbr4zJgsL7NRYsNi4h0RfN0OA4tfN1DIhLrI5kCTvbqclIWP0GcsxAHijhnISmLnyB7dbnVooUly8rPYdwTs9hU2JVqJWwq7Mq4J2axrPwcq0ULX+q8LAnlEkO5xOASafTFyVYCTae5beUTRaCU+pK6EZvcmQW8pZuyrgW6iEgftHCLnyuljiilStGCZjekUIIP042dsvgJBM2Zv+G4Q4DYJQusky+MWfyPPURRSVdKEaArpbSjghFPXG332Cyi7ssSxHCUGI7iAPvFyUICNUfQj7rhDUv0tPrSz0B3JTsXoFevXqxcubJFglRUVLT4XE+oRTuZuEZzbliMFiD1IrSGVWiKoJ+zkOKIeHacexnxKVA8e3bA5PMHwSxj/4ULqd5cTfGGG/mLq1iLioMWMLka2IQihcM1D55+i59gxaG5yNiVAZUzmNsQ/Cdf4n//STW118SI0dkG7a1UgC7//ScrV15c57yYmBiOHTtWs+90OuvsB4IuXbqQnJxMdXU1CQkJvPzyy3Tp0oXCwkKuueYa1q1bB8D8+fN5/fXX+fDD2rAUzz//PPfffz/5+fl069aNr776itmzZ5OQoAWNmzFjBvfee69P5a2srGz6NfS0yqwlG1qw75x6jv0XON+0vxwYA9wD/N6U/gBwT2N1BdPK4uKIBFUN6iSoKn1zgXKCOgEqG5TSNxeorJnzAiqfPwhmGbNmzlMuvb33gZoH6pTe9tWgXgR1XN83rktxRELA5QzmNlTKf/JVg/pAvxbGf8K4Fs6aTznjvLy8vDr7R48ebbCet99WKiFBKRHt8+23vZe9Q4cONd9vvPFG9dhjjymllMrPz1fJyclKKaXeeustlZqaqg4ePFgjY1FRkbrkkktUfHy8OnjwoFJKa9/p06d7L1QDuLeZUtavLN6NKUYuWnzc3Q2khwx9nUVsROsBGJugvd0UA4nUxtezh4n8T+ySBQjam2Y3tJiUUWhtHwH8HGir7xvDd32dRRZIGia4GVA8j/AfasekRd9caLFS9wKC8mrYbsECmDsXCgs1zVJYqO0v8OFfb/z48ezeXfdR9d577/HUU0/xv//9j+7du9ek33333WRkZAT1HEegFMFi4EbdeuhcoFxpMVk/Ay4RkVh9kvgSPS1k2Ozox8Oc2ZAKLWJ2W/2Y/dAJDH2dRTVzNIZidkfQhomMYYk9EfEectn4AvOcwC4Un6N4BO0amHGgPSTSgdN4N19w//1w4kTdtBMntHRf4HQ6Wb58OTNnzqxJKyws5M477+R///sfvXv3rklftGgR/fr1Y8SIEWeU88033zBixAguu+wycnNzfSNcC/GV+ehCtBi3Q0SkRERuEZHbReR2PctStDi1O9Di2f4cQCl1BHgU2KBvj+hpwYvbG85Trv2kULchFXCUzgA1b6dV+jH7oeNf9kTE8w1ae5sjbSjgdNvONfM2u4ErgDKgr7PQnjj2E0YPTaGNHS8GBgHVROACyulMOdp1uQZ4lVrl3dIedFE971r1pTeVkydPkpaWRu/evdm/fz8XX1w7j9GjRw/i4+N57733atJOnDjBE088wSOPPHJGWaNGjaKwsJDNmzfzy1/+kh/+8IfeCeclvrIamq2U6qOUaqOUilNKvaaUekkp9ZJ+XCmlfqGUSlRKpSqlNprOfV0pdZa+veELefyJ+Q0nH0UVVTyEdkO7EEoiEsiZOY8TSaNrHkTH0CJ6HwRKp8+xTPZw4Mi067kbWEPtW6cCcmbOI+oP95Mzcx4lEQn0B3oBr4FtseJHjB6w0UOL0NMduHAoRYwqJ0aVoxAEbRgvAu960PH1vGvVl95U2rVrR2ZmJoWFhSilePHFF2uOtW/fnqVLl/LSSy+xQB+Dys/PJz8/nxEjRjBgwABKSkoYNWoU+/bto3PnznTs2BGAadOmUVVVxaFDh7wT0AtsX0PNxPyGkwD8G+gIHIuIxaFcxFUXkLrocfrcPLXmodMZoad05L6zxpM6IcZS+Vs7u7vv40C7WAY74uso5tQJMZCeTuqix4mrLmBPRALPoHVNbVNf/7EnIr5m/F+5pbvng9oedHU9+ZrC449D+/Z109q319J9Qfv27fnrX//KM888Q3V1dU16z549+fTTT5k3bx6fffYZycnJHDhwgIKCAgoKCoiLi+Pbb7+ld+/e7Nu3zzCQYf369bhcLrp16+YbAVuArQiaSX1vOGe8uZgeOo6nn+KGyTfxwY51HPjd7+xhCF+jD9cVRyTw8BtvkF7poPzyH+N4+qkaxUx6ep1T+jqL6AxEU3e82p7D8S2l0+fwKXCHKU1xZs+4dPqcGkWxH7gabXivJT3oOXPg5ZchIQFEtM+XX9bSfcXIkSMZPnw4CxcurJM+cOBAFi9ezE9/+lM2btxYz9nw/vvvk5KSwogRI/jVr37FO++8Y+lksu1rqJnsiYgnzllIFbW2zzXp9ZyTvbqci5e9yH3At8BU3X49m3mkBkLoVo4xXLcOmA3cpg4jjbSvcR1Be+AY1l4NXUeb5pM6IYY7153FjAOlKHWE3RHxlE6fc0bPOHVCDNnMI3bJAvo6i9hNFC+eczl3tbAHPWeObx/8oK2tMPPxxx/XfM/JqfWuM2LECHbv3n3GOoeCgoKa73feeSd33nmnbwX0ArtH0ExKp8/hKHA5cFJP8/SGY8YYTroHmELtQjN7GMI3GO2bBvyK2oVJDbWv8QbqAF4Hvqbx62jTREw+nuR3V/H1wTKG3vPPBnto7j3oi2fdw33ZZTjuTWfAAKgsq4R9+6z4NWGBrQiaSeqEGP6adjnQjnbuY9D1YAw3OKhdcWxOt/GOvs4iTqK1rcMtvT5SJ8TUzOEo4GnaN3odbZqG2cfTEhJ50nWMi/54bZN9PC0rP4ffLPoLkyq/YK/uIyqq/ABHne0bP9mmRdiKoLmkp7OqeyU/e3d+ncnhM95wTJgnvFzYpqS+Zk9EPCuoXRdgTq8X0xvotUeO8FXnSPq/+dsGr6NN0zD7eLoEuIdTRFHJ4n/safL5nTjFRyi6A90oRVCUHaxq9FyblmErgmaye/duNm3axIwZM5p8jnki7CiaB77T2MMQvqJ0+hxeBr4wpTVnmCc2NpYpU6bw/vvv+0O8sOM3h+8nmqoagwoHEE0VvznctBVdxvltqDUlFaBXdUg5HQgpbEXQFEyLyBbExXFpWRU7rnusyVY/5mGIzgj7iOLN8bPtYQgfMXBMW5ZHRtPfEYcLodjRv3nDPBkZXFgsvPqzO22vpD4gHm1Iropa9yrm9Kaeb6wAN8qI4rSvRLRxw7YaagKGVYqgmRv+VlVoi4+aavWTnq7n0wyZr3vqKb4rKrKHIXzEp0lJTJg8kZRPPwVg58qVTJw4scnnZ68u5/YNH/AJUE6tV1LbqqtlnOrck+ij+/kR8DbQ1ZTeronntzu6H4AVwDbgB4CKaHOGawob32D3CJqAYZXiAn4BjMQ7q58rrriCDz/8EJfL1Xhmm0b54IMPuOKKK1p8fuySBUSjuT8w+hC2VVfLaTdtEusj2lBArRJwRkbTbtqkJp/vjIwGoBPwEqAQHDHWhqaMiIggLS2NlJQUZsyYQVlZGaCZhaakpNTke+WVVxg9ejSlpaX89re/ZejQoQwfPpwf/ehHNeesX7+etLQ00tLSGDFiRB2X1VZgK4ImYLY+MezN3dObw9lnn0337t1Zu3at98KFOZWVlXzyySfMmjWrxWUY17E5Vkc2DbBwIZ/NmsG0zp1rVnRFzH8N3BZfNXR+xPzXICGBc4C9DgfO2C4waFDTZViwgDoxSn3getRwMZGTk0PXrl3ruJgw+Oc//8nzzz/PZ599RmxsLBdffDE5OTlkZWVx9tln8+STTwKQkpLCxo0byczM5NNPP+W2226rs0o50NiKoAkY1ienqTvm2WKrn4wMzj/dhfnnT7XHpFuKPm/zbod4hpeXU9VvXIvb0OzewOyszrbqajlLS0qY9uGH4HJBQUHzV3fNmQMFBUQoxdTZsznZnFW3AfBD3VQ31JdccgmRkdoI/LnnnktJSQmguakw0isrKy13UW0rgiZQOn0OJ4HpaD5TwLvFR9mry7l929ccUscQO7ZxizDmbXa7DvIk3jmNM6y6BG1x2TrsxWXecODAAb7//nvOP/98n5Q3bdo0Tp482XhGAz/7oW6OG2ozr7/+OpdddlnN/rp160hOTiY1NZWXXnqpRjFYga0ImkDqhBj+dd71lBJFdBMXkTVE7JIFDAfM7yf2mHTzMNrqbuA8Pa2lbWi26toHzJfO9uIyL/jss8+YPHkyUVFRPinv0ksvpbKysulzan7yQ91cN9RmHn/8cSIjI5lj6hmNGzeO3NxcNmzYwJNPPkllZaVX8nmDrQiaQno6+ZMGctm8e5q8iKwh+jqLalzu2g7PWkZfZxEH0Pw9ed2GpsVlU77+mnUjBnl1fcOdpUuXMm3aNJ+V161bN6Kios7w9VMvfvJD3Vw31Abz58/nv//9LwsWLPA4BJSUlETHjh3r+CsKNLYiaCLLli1jypQpPinLGHtW1F0Na49JN509EfGsppmriZvA2LFj2bVrl6W+4UMSfc6mICKB/73zDmm3Pey7ea+MDKJccGTbDtTGjZzemIXafaR+30N+9kPdVDfUAJ9++ikZGRksXryY9iaZ8vPza84tLCzk+++/Z8CAAT6RryX4KkLZVBHZKiI7ROReD8efE5FMfdsmImWmY07TscW+kMfXlJWVkZuby/jx431SnjEmfQh4CmpCK9pj0k2ndPoc/g18bkrzRRu2adOGCy64gC+++KLxzDY1GHM2R11FPAyMce322bxX9upyOladwolL70mfplPFQU5UuL8G6ATAD3VT3VDfeeedHDt2jIsvvpi0tDRuv10L2vj1118zYsQI0tLS+NGPfsTf/va3OnGOA43XsxMiEgG8CFwMlAAbRGSxUirPyKOUutuU/5dopvgGJ5VSad7K4U9WrVrF+PHjadu2rU/KM7vcfcVZyARHX7pdfjOpE2JY6ZMaWj8p53Vm+Wcd+UVVZ1yuvewxuTde6WXZU6ZMYfny5VxzzTW+EDUsMNbaDAWS9bTaORvv3sRjlyzg6B0XMJBadxMCRJYdBvp5PskPfqhb4oZ6x44dHsu64YYbuOGGG3wqnzf4Ypp6LLBDKbULQETeQXOnk1dP/tnAQz6oN2D4clgIqLPS+KKbbmLbhAnMnTtXO7Zype/qacXsvPJKop5/ngnFxYgIcVAbR8DLNpw8eTLPP/+8lxKGF8bcTAS1wZrM6d6WfYxaBWDQxnY54TN8oQj6AcWm/RJgnKeMIpIADKSuf7C2IrIRza3IU0qpj+o5dy4wF6BXr16sbOGfvaKioknn9l+4kKIcOGvtJyx3FfNH6cWKRTuJT4Hi2bNbVLcnevbsyXvvvcfZZ5/dLPmsJBhkXLp0KUOGDGHVqlVnHPNGvv4LF3IwW3F0VwE7RYhy9GfHuZf5/LoHQxs2RHPlS3TE0cdVzGmgHbUP7N2OOHZ6+TsTHZqKN4cUBagiilNuwV+CCafTeUZwmkBSWVnZ9GuolPJqA64CXjXt3wC8UE/e3wHPu6X10z8HAQVAYmN1jh49WrWUFStWNClf1sx5ygXqAKhHQFWDcoHKmjmvxXV7Yvv27apfv37K5XI1Sz4rCQYZb7zxRvXSSy95POaNfMZ1fxDUem05kl+uezC0YUM0V76smfPUBlDX623my3bLmjlP5X3yiTqxYYPat2GDUhs2KNeGDer49hKvy/YnR48etbT+vLy8M9KAjcrDM9UXk8W7gf6m/Tg9zRPXAXVmV5RSu/XPXcBK6s4fWIYx5hkLzEPr7vrD1j8xMRGXy0V+fr5Py23tfPnll1x44YU+L9e47g8Ao/Q0e41H46ROiOGDlIsR6YjLB2tt3Ms+3S4GB1HsBk4SxbGOPWjfMaLRc22ahi+GhjYAg0VkIJoCuA643j2TiAxFe65+Y0qLBU4opU6JSHdgAhAUfhaMsU3B92OeZkSEiy66iC+//JJBzfGlEsYUFRVx/Phxhg4d6vOyzWPd9hqPZpCeTv6333LZPS/guOmmunM2Pig7assWopOSaP/991T17avZ43ey1glda8LrHoFSqhq4E/gM2AK8p5TKFZFHRGSmKet1wDt698QgCdgoIpvRPM4+pUzWRlZi9i+kPKT7kgsvvNDjWLeNZ4zegD/8s5ivr9mfvr3Go3HWrFnDeeed13hGL+jYsWPTF5bZNBmfrCNQSi1VSp2tlEpUSj2upz2olFpsyvOwUupet/PWKKVSlVIj9M/XfCGPLyidPodjwDRqFy35xdY/I4Pe733HF/PfxiUOEiffaDug84QpONCqG25g+H++8Es7mf0OLQC+wl7j0RSKi4uprKzkrLPO8ms9ViqC+txQN5WJEyeyceNG/wjnJfbK4npInRDD++fN4ShROHw85mkme3U5s1a+wjCqOYaiv6vYdkDnAWPBUpyzkKHAbFXql3Yy+x0qAN6z/Q41CaM34G8vmh06dOD48ePUHVhwIyMDVqyom7ZihdcvDU1xQx2q2IqgPtLTOThrOBf8+g6f+Beqj9glC3AAH6EF4QB7ctITxiSuAn4NnIWf2snkd2jc0qXkXTTK9jvUBAIxLATayu82bdpw6tSp+jOdcw5cc02tMlixQts/5xyfyWF2Q71+/XrGjx/PyJEjOe+889i6dSugOam77rrrSEpK4kc/+lEdD6p33HEHY8aMITk5mYceql1WNWDAgBr3Jhs3bmxWpD1vsENVNsC6deu46qqr/FqHMQnZBjsoSkMY7WEERHdP9wdjx45l06ZNOJ1OIiJsC5WGWLNmDX/+858DUlfHjh0b9tQ5aRK895728L/jDvj737X9SU2LkNYYhhvqW265BYChQ4fy1VdfERkZybJly5g3bx4ffPABr732Gu3bt2fLli1kZWUxatSomjIef/xxunbtitPpZPLkyWRlZTF8+HCfyNcS7B5BA6xbt45x4zyujfMZZgd01R7SbTSM9nASOEd93bp1o1evXmzZssVvdbQGjh8/Tl5eHqNHjw5IfR06dGg8PsGkSZoSePRR7dMHSqA+N9Tl5eVcffXVpKSkcPfdd5ObmwvA6tWr+fGPfwzA8OHD6zzo33vvPUaNGsXIkSPJzc0lL89aGxlbEdTD7t27OXXqFAMHDvRrPcbk5AngUTRlYE9OnonRTvcB3+tpgWincePGsW7dOr/WEbLoE/hLOycy/MQJDnUc6n9Dh337aHPkBJVHj9Z4Ij2xY/eZnkhXrNB6Ag88oH26zxm0gPrcUD/wwANMmjSJnJwcPv7440bjCuTn5/OnP/2J5cuXk5WVxfTp02vOiYyMrIm7EMj4BLYiqId169YxduxYv09+GZOT5REJLARWSi97ctIDqRNiyJ5xH2/hoCP4bfLeHVsR1I8xgX/CtZ/f412UuKZyosJJzLGDtEUz7Y3iNO3K9tb1RGrMCbz3HjzySO0wkQ+UAZzphrq8vJx+/TTnd/Pnz6/JN2HCBP71r38B1MQtBjh69CgdOnQgJiaG/fv388knn9ScM2DAADZt2gTABx984BN5m4KtCOohEMNCQJ3JyTGzZ7Ms/WZ7ctIT6el0/dvPkZ7diXf5b/LenXHjxrF+/Xq/1hGqGBP41wFGAEZ/GzpElh1G0IwFjIdXrSdSnQ0b6s4JGHMGGzb4TA6zG+r09HTuu+8+Ro4cWSc+wS233EJFRQVJSUk8+OCDNUNnI0aMYOTIkQwdOpTrr7+eCRMm1Jzz0EMP8etf/5oxY8YEdF7Kniyuh3Xr1nHfffcFtM4xY8bw1VdfBbTOUGLjxo2MGTMmoIG+R4wYwfbt2zl+/DgdOnQIWL2hgHkVdqAMHQyPow16IvX0cjBpktfzBA25od62bVvN98ceewzQhpLeeecdj2WZew5mLrjggjplBQq7R+ABp9PJpk2bGDt2bEDrHTNmTI3pmc2ZGIogkERHR5OSklLTXbepZU9EPIcIzOp7gypq4yCretJtmo+tCMzok19fRPWnb0UFx3uMDOgq35EjR5Kfn09VVVVA6gs1rFAEZGRw1kEX/5s4C5c4KIkcYK/81imdPoe1aBP4xtu5vyfwq7t0q7Gw20dtdL/qLt38Vmc4YCsCE8bk10nXXh4lMJNfZjp16kSvXr1qzM9salFKsXHjxoCZKBpkry7nsvyNKFWGAxXweyKYSZ0Qw6dDLqRaOvvc42h9tO8YwckufXARxX7gBG042aWP7YnUS+w5AhPG5NdUahvGV+H2msrZZ5/Nxo0bSUtLC0h9oUJhYSHR0dH07ds3oPXGLlnAROBKU1qg74mgJT2dPWvXct0jr+C45hrfehytj969aQ8cO9aZ9nv3UtWzJ126dPF3ra0eu0dgwpjkcmDdKt+hQ4cGrWMqK9mwYUPgh4XQrn0/tJXf7uk28O233zJypDUhRNq3b8+JEycsqbu1YSsCE3si4qkisJNf7gwZMsRWBB6wZH4A7doL2pi0yy093Dl8+DBHjhwhMTHRkvptReA7bEVgonT6HPKAWwjc5Jc7iYmJ5OXlNexUKwyxShEYK5pfAr7W0+yV3xqZmZmMHDkSh8Oax0igFYHhhnrEiBGMGjWKNWvWAFBQUEC7du0YOXIkSUlJjB07to556Pz58+nRowdpaWmkpaVx4403BkzmpmLPEZhInRDD48UzOJm5HJc6yZ6IeEqnzwncKt+MDKIX7WTAyWo2t21L34iE2vrDcYFZRgbZq8uJ+e/bbHIV0WflNrIvvzGg7ZE6IYZs5uH6+B+8pA4zyHxNwhwrh4VAM+11Op1UVVXRpo374J3vMVxMAHz22Wfcd999NQGlEhMT+e677wDYtWsXV1xxBSdPnuSOO+4A4Nprr+WFF17wu4wtxSeqXESmishWEdkhIvd6OH6ziBwUkUx9u9V07CYR2a5vN/lCnhaTns7hiWcx/skH/ep6uj6yV5czcc3L/JgqFIG3Wgo2DCsucRVxH5DmKgl8e+grv3/w7TI2DxsW8HsimPn222/reNQMNCJC+/btG3dA5weOHj1KbGysx2ODBg3i2Wef5aWXXgqwVC3Ha0UgIhHAi2irzIcBs0VkmIes7yql0vTtVf3crsBDwDhgLPCQHsfYMjIzMy2z2DGsln4LGIMg4RybwGiPPsA9eppV7TFs2DDy8/MteegEK1YrAgjs8JDhfXTo0KHceuutPPDAA/XmHTVqVJ0Vwu+++27N0NAbb7wRCHGbhS+GhsYCO5RSuwBE5B1gFtAUv6qXAp8rpY7o536OZr250AdyNRullKWKwLBEicQOnA61v1vQ3Bi4pweSqKgoBg8eTG5uriVzFcHGsWPHKCkpYejQoZbU37lzZ5+X2WDUM+oODX3zzTfceOON5OTkNKmsYB8a8oUi6AcUm/ZL0N7w3blSRC4EtgF3K6WK6zm3n6dKRGQuMBegV69erFy5skXCVlRU1HvugQMHEBG2bNliiQ/6REcc/V1ac1RRqxB2O+LY2cLf6w8aakNfYrRHFVo7GN3XxtrDX/L17t2bd9991ycxcwPVhi2lMfmys7OJj4/n66+/rjePL4mJieHYsWM1+6WlpURERHDq1Cn27NnjE3fx5vIby5OSksLBgwfJz8/nxIkTuFyuOuevXr2as88+m2PHjlFZWcnp06ebVL4vqaysbPI9FqjJ4o+BhUqpUyJyG/Am8IPmFKCUehl4GWDMmDGqpSHcVq5cWW/4t48//pixY8cGLDycO9mX30Dc4icQ4M/ANUA8UHb5DZbJ5ImG2tCXGO1xDTAf6IpmsdNYe/hLvu+++46dO3f6pOxAtWFLaUy+rKwsJk6cGLDfsGXLFjp16lSzf+zYMTp16kTHjh0pKiqiffv2AfHWacjw/fff43K5SEhIoLi4GIfDUXOsoKCABx98kNtvv51OnTrRtm1boqKi6sgfCNq2bdvkyXxfKILdQH/TfpyeVoNSyuQjllcBw1HLbmCi27krfSBTi8jMzGTEiBFWVU/qhBhWHJrL4HWfsdpZSEfpzvkz5oathUrqhBi+On0XKz79C51RlFhssZOWlsZ//vMfS+oOGnRLrq8//gs/UMcpeeUT6yzb9u3jZIWTti7Fie++I5ooqrt009xN9O7t8+qMOQLQhn7efPPNGuWzc+dORo4cSWVlJZ06deJXv/oVV155ZQOlBRe+UAQbgMEiMhDtwX4dcL05g4j0UUrt1XdnAsa4y2fAE6YJ4kvQfFhZQmZmJldffbVV1UN6OjJ2JXET/8GIhx5it9NJqu7SNixJT8c5diXDj64ncvXqwLgwaIARI0awefNmXC6XZbbzVmNYckUDl6NZtvVb/ATZzCM1wLKcqHDSrmwvXdGGDqM4TZuyvZygD+39UJ/T6fSYPmDAAI9GBMZQ0M0338zNN9/sB4l8h9d3s1KqGrgT7aG+BXhPKZUrIo+IyEw9269EJFdENgO/Am7Wzz2CFqFxg749YkwcW4GVE8XuDB8+vCaiUThjdVBvM127dqVLly7k5+dbLYplxC5ZQBXwD2on86yy5DKC1PQCjEgRZwSpsWkSPpkjUEotBZa6pT1o+n4f9bzpK6VeB173hRzecPToUfbv38/gwYOtFgWwFYFBVlZWUFnppKWlkZmZaZlbBavp6yxiP9AD6y3bzMFo6g1SY9MkwrN/64GsrCxSUlICGh6uIRITEzl06BDl5eG5mMwgmHoEoCmCzZs3Wy2GZeyJiGcz4D5IYoXvJXMwGhe1/sHsIDXNx1YEOlZPFLvjcDhISUkhOzvbalEsw+l0kpubS0pKitWi1GD0CMKV0ulz+AL4lynNKt9LRpAaAQ5Q6yzSDlLTfGxFoEclW/2r+xj+8stBFYEq3IeHduzYQa9evfyyeKhFZGTQ9m/L2PTx0rCNVpY6IYZveibicvQIWDCa+jCC1JwmimPAUSLtIDUtJOydzhlWEO2w3grCneHDh4f1MESwDQtlry7n0uV/Jwk4SXDdKwEjPZ2C559n0va1OAYNstaSSw9SA/1ov3s3p4Ee/TyuR7VphLDvEcQuWYATeAFt8RYEj3+fcO8RZGVlBdVwXeySBUQA/wXa6WnBcq8EitLSUsrKyhgwYIDVotShXbt2fvcDZbihTk5OZsSIETzzzDO4XLVRKr7++mvGjh3L0KFDGTp0aB2fQg8//DD9+vUjLS2NlJQUFi9eDEBRURGTJk1i5MiRDB8+nKVLNZsbw7W14Z/o9ttvrylr6tSpjBgxguTkZG6//fZ6zVqbQ9j3CPo6izgMxGK9FYQ7w4cPJycnJ2zt1rOysoLKd7txT7TBugh2VpOTk0NycnLQ3Y+BUARmX0MHDhzg+uuv5+jRo/zhD39g3759XH/99Xz00UeMGjWKQ4cOcfHFF5OYmMj06dMBuPvuu7nnnnvYsmULF1xwAQcOHOCxxx7jmmuu4Y477iAvL49p06ZRUFAAaAYjnuaj3nvvPTp37oxSiquuuop///vfXHfddV79tuC6mhawJyKeHILDCsKdLl260LVr17C1Ww+2oSHjnnBR934JhnslUGRnZ5OaGnwDYdHR0VRVVdV9O/7mG3jySe3Tx/Ts2ZOXX36ZF154AaUUL774IjfffHONN9bu3bvzyCOP8NRTT51xblJSEpGRkRw6dAgR4ejRowCUl5c3KSa3MWdWXV3N6dOnEZFGzmicsFcEpdPnsAowO4YNpghU4To8VF5ezsGDBxk0aJDVotRgRCvLAV7R04LpXgkEOTk5QakIHA4H0dHRVFZWagnffAOTJ8MDD2ifflAGgwYNwul0cuDAAXJzcxk9enSd4yNHjiQ3N/eM89atW4fD4aBHjx48/PDDvP3228TFxTFt2jSef/75mnz5+fmMHDmSiy66iK+++qpOGZdeeik9e/akU6dOXHXVVV7/lrBXBKkTYljTZyjV0t1yK4g66NZMA5Z8yeYrrggfCxX9d3/RdQjJx4+zNzoxaH536oQYcmbOI8IRx71AsSM+OO6VABKsPQJwGx5auRJOnwanU/sMAk+vzz33HGlpadxzzz28++67iAgLFy7k5ptvpqSkhKVLl3LDDTfgcrno06cPRUVFfPfddzz77LM1w1AGn332GXv37uXUqVN88cUXXssW9nMEpKdT/MYbXPTZFzhSUy33Z2NgWDPNAvYRPhYqxu/eBvyOIPvd6em6DI/TtndvZONqUuOC4W4JDEopsrOzg2pdh5k6imDiRIiK0pRAVJS272N27dpFREQEPXv2ZNiwYWzatIlZs2bVHM/MzCQ5Oblm35gjMPPaa6/x6aefAjB+/HgqKys5dOgQPXv2JDo6GoDRo0eTmJjItm3b6qyyb9u2LbNmzWLRokVcfPHFXv2WsO8RnDp1ioKCAoYMGWK1KHUwonP9AM2LH4SHhYrxu2cAP9TTgvF3p6Sk1BuUpLVSUlJC27Zt6dGjh9WieKRO2Mrx42H5cnj0Ue1z/Hif1nXw4EFuv/127rzzTkSEX/ziF8yfP79mcvfw4cM8+OCDpDfikTU+Pp7ly5cDmqvtyspKevTowcGDB2vmO3bt2sX27dsZNGgQFRUV7N2r+e+srq5myZIlPgkOFPY9gq1btzJo0CCiooJrWbo5Olck1KygbO0WKsbvcxDcljmGIpg6darVogSMYB4WAg+WQ+PH+1QBGG6oq6qqiIyM5IYbbuA3v/kNAH369OHtt9/mZz/7GceOHUMpxe23386MGTMaLPOZZ57hZz/7Gc899xwiwvz58xERvvzySx588EHatGmDw+HgpZdeomvXruzfv5+ZM2dy6tQpXC4XkyZNqmNa2lLCXhHk5OQEZVd3T0Q8cc5CBKhGeyiKkW6taH5lT0Q8/ZyFnEYLTymm9GD63SkpKQGLzhUsBLsiaNOmDS6Xi6qqKtq0aePz8huz17/wwgvZsGFDzb45ItnDDz/s8Zxhw4axevXqM9KvvPJKj/EMevXqVacOXxH2Q0PBqggMCxWA99F8fIeDhUrp9DnsBa6kVgkE4+9OTk4Ou6GhYFcEIhKQ9QStEbtHkJPDT37yE6vFOIPUCTFkM4/YJQvY7ixkn3Qmbsadrd5CJXVCDK8cnE3ZN//BxWn2RMRbGpWsPpKTk9myZUvrX+ynRySLXbKAbGchsxcuJ/ujHdZEJGuIffs4UeEkuuIEJ7dto62fo5W1NnyiCERkKvAXtN78q0qpp9yO/wa4FW2U4yDwU6VUoX7MCRguNouUUjMJIMHaIzBbqKR8+CHrX3+d1EWPWyxUAEhP50T0Xxg1MhbHiy8GjRWXO507d6Z79+7k5+e36tgEhhWXEy204GTXPqItsuJSStW7eMqIVtZF3/d3tLJgRynVeCYTXr/KiEgE8CJwGTAMmC0iw9yyfQeMUUoNRxvpMBuFn1RKpelbQJVARUUF+/btC6pFS55ITk72uDCltZKbm1vH7C5YCQfLIcOKywE8BLTFGiuutm3bcvjw4XofcEa0si76BuEbrUwpxeHDh2nbtm2Tz/FFj2AssEMptQtARN4BZgF5JsFWmPKvBX7sg3q9Ji8vj6FDhwZNMJr6SExMZN++fRw/fpwOHTo0fkKIk5OTw5w5wTUn4AlDEZhtx1sbZmutyHrSA0FcXBwlJSUcPHgQgMrKyjoPOnVoD576CgqQLb6fOG4K7jIGkrZt2xLXjDUuvlAE/YBi034JMK6B/LcAn5j224rIRrRho6eUUh95OklE5gJzQZs5X9nClYIVFRU15y5dupTu3bu3uCx/YJbPTN++fXn77beDYr1DfTL6AqUUWVlZlJWV+eQa+xOHw8EXX3zBhAkTmn1uoGRsKYZ8iY44+ruKcaI9VI1Xpt2OOHZaKH9FRQUdO3as2U+cPpf+Lu0xVIUmpwModvRn5/K3gkLGQFNYWNj0zEoprzbgKrR5AWP/BuCFevL+GK1HEG1K66d/DgIKgMTG6hw9erRqKStWrKj5fvfdd6unn366xWX5A7N8Zq6//no1f/78wApTD/XJ6AuKi4tVz549vSrDn/KZ2bRpk0pJSWnRuYGSsaUY8mXNnKdcoOaBygSlQLlAZc2cFxTyGRhyKlCvgvomCOQMxmsMbFQenqm+6BHsBvqb9uP0tDqIyBTgfuAipdQpkyLarX/uEpGVwEhgpw/k8kj/hQvJfu7zGiuINEdPsleXB58VhBvhMk8QKvMDZGSw99BwcnN3IFJFQkIbXp2zgikxG4L6PmouhvXa+4v/yFVUURKREJRWXGYru83OQgoklg4z7gg6OYMVXyiCDcBgERmIpgCuA643ZxCRkcA/gKlKqQOm9FjghFLqlIh0ByZQdyLZ5xTlwMQ1T9S4b7jEdYBeweLLpgGSk5N5+eWXrRbD74SKIlhWfg4X/nEWV3CaF4iiTWEsUU9Us2zeIqZYLZwvSU9naFUVRZ2fJan0qDb2bLVMnjBZ2SX/4x+sX7+e1NfCwMrOR3htNaSUqgbuRFvztAV4TymVKyKPiIhhBfRHoCPwbxHJFJHFenoSsFFENgMr0OYI8vAjZ639BEEb7/wt0Ivg9GXjTjj1CILSnNeNxf/YQxSVLMBFT6AbpURRyeJ/7LFaNJ+zfft2+vfvb9nEZ3MJl/+KL/HJOgKl1FJgqVvag6bvHl+SlFJrILAv4v1cJTXfrbSCaC4DBw7k4MGDlk9A+Zvc3NygXODnzm8O3080VTipfZuKporfHL4fCH6Lp+aQl5cXEr00g2HDhpGXl9fgugOburTiJZGe2e3QOrbVhFaUqYiICIYMGUJenl87TJailAqZh0482ouDQruX3NNbE7m5uQwb5r40KHjp2rUrHTp0oKSkpPHMNkAYKoId516GQgtWv0ZPC0ZfNp5o7V3e4uJiOnbsSGxsrNWiNMqpzj0BzczteQ/prYlQUc5mjF6BTdMIO0UQnwI5M+fxH9pyDIInIlkTaO2KIFQmigHaTZuEMzKaOOA+4DTgjIym3bRJFkvme0KtRwCt/7/ia8JOERTPnk3qosfZ2SeWlMJC4qoLNB8+wW7yl5FBhw+/ZeMzL+ISR+sKXamHp1w9/XqGLVsWGr9t4UIi5r9G24QE4oEdffoQMf81WLjQasl8SlVVFTt37gyKhYzNwe4RNI+wUwQApaWlVFRU0L9//8YzBwnZq8uZtv7fxFCJA0Wcs5CUxU+QvbrcatG8xnBs1laVMQdC57fNmQMFBQybNYu8v/5V229l7Nixg7i4ONq1a2e1KM3CVgTNIywVwZYtW0hKSgopi4LYJQsYCCyEmjgFoWD22hQMx2b3AOfoaaH021rzQycU5wegdmhINdMLZ7gSloogLy8v5MY8+zqLcABRUMe5VrCbvTaFvs4iFJo5byj+ttauCELtvwLQrVs32rZty549rW9dhz+wFUGIYJi3OgGXh/RQZk9EPPvRnIW5p4cCrVkRhNIEvjv2hHHTCUtFEIpWEEboyneBT/W0UDF7bYzS6XPIBMzT9aH024YOHcr27duprq5uPHOIEYovTQatWUH7mrBUBKF4c6dOiCFn5jyqHD14gdAye22M1AkxfJE8hePSCRcScr+tffv29OnTh127dlktim/IyEDdv5CCiAS2Z2fTccys4Lficicjg4gdkTzwQC4OBwwYAMvuXxFavyGAhF3M4uPHj3PkyBESEhKsFqV56E612m6/mccuuYS4/PzgdP7VEtLTKd2+nbF3Xonj9tuDNjxlQxhvn2effbbVonhN9upyJq55mSLgUWCwqxgVAo4ZzSwrP4epn02ngtO8wquUtlangD4i7HoEhYWFDB06NGQDjg8cOLAmWllrIhR7aWaSk5NbzTCEYcXVD7hbTwslKy7QnAKOpIrncSK0bqeAviA0n4ZeUFhYGLKTXwCRkZGcffbZfP/991aL4jMMH0OhrAiGDRvWaiYmDWstB7URyczpocBvDt9Pb6rrWNnVOgW0cSfsFEFBQUFIP3BAe+hs2bLFajF8xr59+2jTpg3du3e3WpQW05omJg1rrSpC10LNcP5XTd3f0BqdAvqCsFMEhYWFIa8IkpKSWs1DB0J/WAg0y6GtW7fidDobzxzkGBZqPwP26mmhZMUFtc7/nkeLjeueblOXsFMEraVHYCuC4KJTp0706NGDgoICq0XxmtQJMSwffyv/RuhMaFqoGU4B2wFG6PrW6hTQF/hEEYjIVBHZKiI7ROReD8ejReRd/fg6ERlgOnafnr5VRC71hTxnoDs12xoRz9H9+2kzZHLomcOZsBVBkKHfX4OKDpBz1lmh4TSvIdLT2f+LifQekEAnpULHMaMZ3SngsJ49yQNISAhdp4D6/VUSOcBvDie9VgQiEgG8CFwGDANmi4j7v/oWoFQpdRbwHPC0fu4wtBjHycBU4G96eT7FcGrmchXzFJDgKgoNp2b1cNZZZ1FUVERlZaXVoviELVu2hLQiMO6vH1FJF0LIaV4DtIYhVObMYdjmzeR16wYFBSHrFNC4v+KchX5zOOmLHsFYYIdSapdS6jTwDjDLLc8s4E39+/vAZNE8vs0C3lFKnVJK5QM79PJ8imEONxj4hZ4WauZwZqKiohg0aBDbtm2zWhSfkJeXR1JSktVitBjj/vo5cL6eFsr3F7QSRQD06tULp9PJwYMHrRalxRj3l6J24tvX95cvFpT1A4pN+yXAuPryKKWqRaQc6Kanr3U7t5+nSkRkLjAXtIu7cuXKJgt4ockczqz5+jqLmlVOIKioqGiSTD169OD999/nyJEj/hfKjabK2BTKyso4efIk33//PVu3bvVJmb6UrykY91cEZzrNq0+OQMvYXHbs2MHo0aODVsbmtF+/fv1YsGABaWlpfpXJHV9d4wtNZrt+e34ppbzagKuAV037NwAvuOXJAeJM+zuB7mgRI39sSn8NuKqxOkePHq2aQ3FEglJwxlYckdCscgLBihUrmpTvgQceUA888IB/hamHpsrYFFatWqXOO+88n5WnlG/lawrG/eUCdVr/bOz+CrSMzWXIkCHqm2++sVqMemlO+916663qb3/7m/+EqQdfXWPj/qrSN2+eX8BG5eGZ6ouhod2AOcJLnJ7mMY+IRAIxwOEmnus1hjmcmVAzh3OntawlCPmJYmrvLwH+BOwjtO8vl8tFYWFhSA/XmQl14wrj/noDWKGn+fr+8oUi2AAMFpGBIhKFNvm72C3PYuAm/ftVwBe6dloMXKdbFQ1EG8Zf7wOZ6mA4bCuJSAhJp2aeCPWb26A1KALz/fUFsMzRM6Tvr+LiYjp27EhMTGjK706o/1eM+2sR7TmIf8x5vZ4jUNqY/53AZ2jDpK8rpXJF5BG0bshitCGff4rIDuAImrJAz/cekIe2CPAXSinfr8jRHbbB46xcuZKJEyeGnFOzOmRkcPrLw+zMy+OUCAcjEiidPke7MULBxC8jg+zV5cQuWUCes5Cxjp5krzwUOvK7Y7q/ku+6iwP9+5P6f/9nsVAtJy8vjwEDBlgths8IdUVg3F/5yR+RsnAhccOH+/z55RPvo0qppcBSt7QHTd8rgavrOfdx4HFfyBEuZK8uZ/SSDGYA5Wjmiv1CyDukYQ4nwA+Ai10H6BlC8jfEsGHDWLt2beMZg5i8vLzQ887bAHFxcTVeh7t27Wq1OC3i9OnT7Nq1iyFDhvil/LBbWdwaMMzJ3kabcYfQMlc0m8OlAz0JLfkbIuTfPml9ikBESEpKCuk5tR07dhAfH090dLRfyrcVQQhieIGM5ExzslDALKc5TnGoyN8Qhh8oFcJB01vb0BCEvoL291yarQhCEMMLpEKLYeyeHuyY4y+HovwN0a1bN9q3b8/u3T43fgsISncJHh8f+tfCjK0IGsZWBCGIYU62E20hBoSWuaIh/wvAaj0tlORvjFB+6Ozdu5d27dq1Goshg1C+JmArAhsPGOZkDkd/7gUKHPEhZa5oyP8f2nKU0PRu2RAh+dDRHZut6j+SpP37SZx8Y2g7zjOTkYHjuSXk/G+535y2+Rt/K4Kwi1ncKjCZK/ZNTOTU0qWk+smawC/o8u+Ke5PU1auJS0gIbXNeN4YNG0ZmZqbVYjQLw5JrC3AX0N9VTFwrseTKXl3Oxcv/zmDgNKFnZVddXc327dv9ZjEEdo8g5ElOTg7JEIllZWWUl5fTv3//xjOHGKHYIzAsuX4IzNDTWoslV+ySBUQCnwKGzU0o/bZdu3bRt29f2rdv77c6bEUQ4oTiQwc019NJSUk4HK3vFjSuSShZDvVtwDFjqGP8hjac6RQwFAjE6vvW9y8MM0I1aHprcC1RHz169MDhcLB//36rRWkyeyLiUWhDJ8otPdQxfoOL0LRSy83NJTk52a912IogxElOTg7JHkFubm6rVQQiEnLXpXT6HPYCV1L71txaLLkMK7W11IatDKXfZvcIbBpl6NChbNu2jerqaqtFaRZ5eXl+f8uxklDrqaVOiGHJ+NmUEY0LodjRv9VYchlWauLoy0OEnpVaIBSBbTUU4nTo0IE+ffqwa9cuzj77bKvFaTKttkegO9Tr+/G75KpSSu5+JjQcAqancyL6L4waGYvjxRfZqTtnbBXoVmou16Mc7tSJznuziOvc2WqpGka/jzr/9222uoroNO4Ksi+/wW/3kd0jaAWE2tvn0aNHKS0tbVX+bAwMM8ypqpRRhFb84kCMRVuJw+Fg6NChITFkZ9xHEa4iHgaGuIr9eh/ZiqAVEGrj0Xl5eQwdOrRVWgwZZpijgZ/qaaFiqtjaFQGEjrm1cR/1AQyH5v68j1rfPzEMCbUeQWueHzBMEgVt3FW5pQcrSilbEQQR5vsowkO6r7EVQSsg1HoErXZ+gFqTREGLtORySw9W9uzZQ3R0NN27d288cwgTKorAuF+qqL2HzOm+xitFICJdReRzEdmuf8Z6yJMmIt+ISK6IZInItaZj80UkX0Qy9S3NG3nCFcNyyOn0fXA3f9CaewTm+NhLgA8IDVPFcOgNgKYIcnJyrBajUYz76GfAXj3Nn/eRtz2Ce4HlSqnBwHJ9350TwI1KqWRgKvBnEeliOv5bpVSavmV6KU/4kZFB/pwn6XbyNDsiI4PXoZbu1KwkcgC5n35K7KzbglNOLzHHLy4GFknHkDBVzMnJCQtFkJCQQHl5OWVlZVaL0iCpE2LIvPxe/o0Qg/9NXr01H50FTNS/vwmsBH5nzqCU2mb6vkdEDgA9gDIv67ah1rrgGrRVk8HqUMuQ8wTwE2C8aw8ShHJ6jckhYOqqVfzrvvtIXRT8kVhzc3M555xzrBbD7zgcjpo5tQkTJlgtTv2kp9Nu5vfEff8+HbdvpyP41TGjeOMPRUTKlFJd9O8ClBr79eQfi6YwkpVSLhGZD4wHTqH3KJRSp+o5dy4wF6BXr16j33nnnRbJXFFRQceOHVt0biBornyJk2+kv6uY09SNWFbs6M/O5W81cGbLaUkbGnI60bq4xhuIP+QMlmtcXl7OnDlz+Pjjj9H+HrUEi4wGP//5z7n99tsZPnw4EHzyueONfE8//TTDhg1jxowZjWf2Am/bcNWqVSxbtoxHH33UZzJNmjRpk1JqzBkHlFINbsAyIMfDNgsoc8tb2kA5fYCtwLluaYLmFPBN4MHG5FFKMXr0aNVSVqxY0eJzA0Fz5XMiSoFyglKmzYn4R0DVsjYMpJzBdI179+6tioqKzkgPJhldLpfq1KmTOnz4cE1aMMnnCW/k++Mf/6h+9atf+U6YevC2DR966CF1//33+0YYHWCj8vBMbXSOQCk1RSmV4mFbBOwXkT4A+ucBT2WISGe0ubP7lVJrTWXv1eU7BbwBjG1MHpu6mK0IqutJDwYMeaoJTcdfLSUlJSXoJyeLi4vp0KEDXbt2tVqUgBAqlkM5OTmkpKQEpC5vJ4sXAzfp328CFrlnEJEo4EPgLaXU+27HDCViuEIP7n9MEGJYF1QBf0B70AajlYoh5/8BO/S0YJTT14SCIsjNzQ3YAycYsBXBmXg7WfwU8J6I3AIUAtcAiMgY4Hal1K162oVANxG5WT/vZqVZCC0QkR5ow0OZwO1eyhN2pE6IIZt5xC5ZwHvOQi5y9KHX5T8JOisVQ853Fj/Fb3FREpFQ64OnFZOcnMxXX31ltRie0f3ZrP74bwxTZZREDqi9JmNbaec8I4Oyr8uo2LefQyJUmu/DIPIFVVlZSWFhYcD8h3mlCJRSh4HJHtI3Arfq398G3q7n/B94U78NdaxUUq68koNXXcWU2bMtFsoD6en0OnCA6i//Rv8jRxCRVhWesj5SUlJ46aWXrBbDI4Yl1yJgJnUtzlrrIG326nJSP36Sa4FjwMAgtbLbunUriYmJREVFBaQ+e2VxKyI1NTWohyGys7NJTU09w4KmNTNs2DC2bNkSlIv9DH829wCG4Wio+EVqKcZvfhEwXB4G428O5LAQ2IqgVZGamkp2drbVYtRLTk4OqanB9N7lfzp37kyPHj3Iz8+3WpQz6OsswoU2LBCKIRxbgjlsZTCH5Az0Aj9bEbQigl0RGD2CcCNYJ4z3RMSzB83QwD29tWL8NkVwW9nZPQKbFpOYmMiBAwc4duyY1aJ4xFYEwUXp9DmsBx4wpbV2Sy7Deq0ceALt9wbjb7YVgU2LiYiIICkpKSgfOi6XK2wcm7kTrOaKqRNiWDbkAk5JZ1xIyIVwbAmGL6gTEQn8Hdjg6Bdcv3nBAo7Fx7O/oIBBU6bAgsDMXdihKlsZxtvn+PHjrRalDvn5+XTt2pUuXbpYLUpgycggeulWvlv1Aa53HOyJiKd0+hz69z4EVoeCTE9n/7p1XPOHl3Fcey1xmPzZrFxpnVz+xGRlN/zSSznwy19y+eWXWyyUzuzZON//kLzqUyQBEUVFOG++hYj//hcWLvRr1XaPoJURrPME4ToslL26nMtXvU5vqnChakJXFgVJp23z5s2MGDHCajEsYfjw4WRlZVktRg0nl64govoUB6h14xxRfYqTS1f4vW5bEbQybEUQXMQuWUB74BNqI00JcNbaT6wTSqeiooK9e/dy1llnWS2KJQSbIog+qnnouRS40kO6P7EVQSvDUATKC6+y/iBcFYFhlhhFXRPNfq4SS+Qxk5OTQ1JSEpGR4TlCHGyKoAjNcimCug9mI92f2IqgldG7d28A9u3bZ7EkdQnHNQRQa5booq6zvd0O69dVZ2Vl1bidDkeSkpLIz8+nsrLSalEAeLbb41QSySlqY12fog3PdvN/PAtbEbQyRCTohodOnTpFfn4+Q4YMsVqUgGOYK65H87MO2p98x7mXWSeUTrgrgqioKAYPHhw08b5n3taXfKK5GkEBh4nlNG2ZeVtfv9dtK4LWhB4OcsCqTWRfeqn1YSt1eVa2H8igykoOdhjSKsNTNoRhrhjl6Me9QLEjnpyZ84gPAmef4a4IQBse2rx5s9ViADAlZgPvX30fX7S9hEhRjE44wrp5i5gSs8HvdYfn4GArxXAiNgtt1aTVYSsNeUqA3weBPJagmysq9RiqZ08is9aS2qcPK1euJNFCsZRSZGVlheVwnZmgmidIT4dHH+XXg9J46ikjcZK++Re7R9CKMBxqzUAL7gDWOtQy5JkMXBsE8liJiDBixAgyMzOtFgXQgtG0b9+eHj16WC2KpYwYMSJ4FAGQmZlJWlpawOu1FUErwrBQcaB19ZRbulXyuFtBBJuDr0CRlpYWNIognNcPmDGGhoLFym7z5s2hpwhEpKuIfC4i2/XP2HryOUUkU98Wm9IHisg6EdkhIu/q0cxsWohhoSLURiozp1shj4I6VhBWymM1waQI7PkBjd69eyMiQWFld/ToUfbu3cvgwYMDXre3PYJ7geVKqcHAcmoXxLlzUimVpm8zTelPA88ppc4CSoFbvJQnrDEsVAA+QFvEZKVDrdLpcyhCC1Fn2NAHo4OvQJGWlmbtxKQ+eV8SOYDNv/89cX98Lewm7+uQkUHOD3/PkEMVZPbta7lxRVZWFikpKURERDSe2cd4qwhmUWsV9ya1Q9ONoscp/gFgxDFu1vk2Z2JYqJREJFAMfCCdLHWolTohho/OuZITtA0bp2YNMWTIEIqKijh+/Lgl9RuT93HOQs4FpqnDpCx+guzV5ZbIYzVGe1xFJR2gxv2HVe1h1fwAgHgzNiYiZUqpLvp3AUqNfbd81WgxiauBp5RSH4lId2Ct3htARPoDnyilPBrWichcYC5Ar169Rr/zzjstkrmiooKOHTu26NxA4Cv5vv32W9544w2ef/55H0hVl+bI+MYbb+ByubjllsB19oL5Gs+dO5e77rqL+Pj4gMuYOPlG+ruKUWiL2yLQemrFjv7sXP5WnbzB3IbgG/mM9qhGawfjPdxTe7SE5sr4pz/9icGDBzNr1iyv666PSZMmbVJKjTnjgFKqwQ1YBuR42GYBZW55S+spo5/+OQgoABKB7sAOU57+QE5j8iilGD16tGopK1asaPG5gcBX8h05ckR17NhROZ1On5RnpjkyXn755eqDDz7wuQwNEczX+Kc//an6+9//bomMTkQpUC5QyrQ5kTPyBnMbKuUb+ZrTHi2huTKOGTNGrV692id11wewUXl4pjY6NKSUmqKUSvGwLQL2i0gfAP3To3ckpdRu/XMXsBIYCRwGuoiIsZYhDtjdmDw2TSM2NpYePXqwfft2S+X49ttvGTlypKUyBBNWzhMYk/TV1HV3Ea6T9+bfXYW1xhXV1dXk5uZatq7D2zmCxcBN+vebgEXuGUQkVkSi9e/dgQlAnq6dVgBXNXS+TcsZOXIk3377rWX179+/nxMnTjBgwADLZAg2rLQcMowJ/ggYa1XDefLeaA8BXgaysa49tm3bRlxcHJ06dQp43eC9IngKuFhEtgNT9H1EZIyIvKrnSQI2ishmtAf/U0opw7nH74DfiMgOoBvwmpfy2JgYNWqUpYrgu+++Y9SoUWjTRzZkZBDx5EdkrV3HhEmTAm6lYhgTLKANAmE/eW82rsgGPpLYwLeHbsm1POUihm/fbpnlklcuJpRSh9EWjrqnbwRu1b+vAc8eBfShorHeyGBTP6NGjeKZZ56xrH57WKgu2avLGf/Js0wFyrDA5UZ6OmedPEl+t+dIKz1GdHQ01vtAtRBTtLJzXnuNVatWkfqW/z19mjEsl1YBv8A6Nyz2yuJWjDE0pCxaNWn0CGw0DJcb/0Lr/kLgXW5s3ryZpKQkoqOjA1ZnKDBmzBg2btwY8HqNe+JW4CI9zQo3LLYiaMX07t2btm3bUlhYaEn9do+gLoZrjTZY53Jj48aNjBlzpvVguDNs2DAKCws5duxYQOvt6yxCoQ3NWOmGxVYErZxRo0bx3XffBbzesrIy9u/fz9lnnx3wuoMVwxpFoVnuuKcHAlsReKZNmzYMHz484P+VPRHxlKBZLbmnBxJbEbRmMjLot72MlVf+BJc4AjMRpU9+fd5tKCOOH2dvdGJ4uzEwYVipnAAeo9YfVCCtVGxFUD9WDA+VTp/DSrT7wcAKyyVbEbRisleXc+m21USqchyogCyhNya/ql37+R3WL9sPJgwrlfKIBN4BvnD0CaiVSkVFBfn5+SQnJwekvlDDCkWQOiGGTweOQaSLpW5Y7MA0rZjYJQsYD5iDItZORPnHOsKY/LoCbSw8EHWGDCYrlYRLL6Xgyiu5ZO7cgFWfmZlJSkoKUVG2k19PjBkzhscfD/A9mp7Ojg8+4Lb5i3BceCFxYIkll90jaMX0dRbRm9oHsjndn3WCHYOgMZKSkli3bl1A67SHhRpm6NCh7N27l7KysoDVeerUKXJychg9enTA6vSErQhaMXsi4hE0dwIut3R/1lkOnMaOQdAQtiIIPiIiIkhLS2PTpk0BqzMzM5PBgwfToUOHgNXpCVsRtGKMycl30WITgP8nokqnz+FLtCXjdgyC+klMTCQ/P9//5oqmGAQbFywg7tb77cn7+sjIYMCeUyy/+KqAGVesW7eOcePG+a38pmIrglaMMTnZztGLRwmMS4HUCTEsGnwekRJjxyBogMjISNLS0vw+OWlM3sc6C7kWmOTaa0/e10P26nIu3bUBUWUBM64IFkVgTxa3ZvTJycQT93Nzjx50O7SFuHbt/F7njqVLue/5d3Fceqllk1+hwNixY1m3bh2TJk3yWx3G5H1b4AFq//D25P2ZxC5ZwAXAlaY0fxs6rF+/nnnz5vml7OZg9wjCgPbt25OcnBwQ07iqqio2btwYFG85wc64ceP8Pk9gTNILdd/67Mn7M+nrLKI/gTOuOHz4MPv372fo0KF+Kb852IogTDjvvPNYs2aN3+vJzMxk4MCBdOnSxe91hTqGIvCnLyhjkv40gTMYCFX2RMTjQFvoF4i2Wr9+PWPGjLEkRrE7tiIIEwKlCNasWcOECRP8Xk+o03/hQo79+hVO791PkcN/E5Ol0+dQDcxAW9EM9uR9fRjGFe8DS/Q0n7eVafJ+7bRpDFu5KSgm721FECYYisDfnkjXrFnDeeed59c6WgNFOZD68ZPcjIuj+G8FduqEGD648KcU0ob29uR9gxjGFZ0dvXkA/xhXGJP3cc5CBgBXq6NBMXlvTxaHCXFxcbRv354dO3YwePBgv9WzZs2awK/ODEHOWvsJAjxJ7duYXyYm09PZG/VnJg+7Bcff/25P3jeEblwxtOph5nTtSoei74iLjfVpFcbkvQJuQLv2wbDy3qsegYh0FZHPRWS7/nlGq4nIJBHJNG2VIvJD/dh8Eck3HUvzRh6bhjnvvPNYvXq138ovKiri9OnTJCYm+q2O1kI/VwmgrcA2jxD7Y2Jy1apVXHjhhT4vt7XSpk0bzj33XL7++mufl22+vmbX01ZP3ns7NHQvsFwpNRhYru/XQSm1QimVppRKA36ANlT5P1OW3xrHlVKZXspjUx8ZGQzIOcznP/2lbxfLmMY8v05IYOSBY+T88PeWj3kGO7sd2nu54N/A6S6Xi6+++spWBM3koosu4ssvv/R5ucb1rUZb8e+ebhXeKoJZwJv69zeBHzaS/yrgE6XUiUby2fiY7NXlXJnzOagKny6WMY95dgJu4GRQjHkGOzvOvawmcPpbwNf4ZxJ3y5YtdOnShX79+vm03NbOhRdeyKpVq3xerjEhPQ/I0dOCYfJevJk8FJEypVQX/bsApcZ+Pfm/AJ5VSv1X358PjAdOofcolFKn6jl3LjAXoFevXqPfeeedFslcUVFBx44dW3RuIPCXfImTb6SPq5jTQDtq3T8UO/qzc/lbzSrLLGPi5Bvp7yoGtLccQRvqaEm5viLYrzFArzffZN+2KM5a+wlvu4rZTkduOO964lOgePZsn9Xz0UcfsXXrVn73u98167xgb0N/y3f69GlmzZrFBx98QPv27VtUhicZ+y9cSH62iyu/eY2tKKoc/dlx7mU+v+71MWnSpE1KqTMdTimlGtyAZWjKy32bBZS55S1toJw+wEGgjVuaANFoPYoHG5NHKcXo0aNVS1mxYkWLzw0E/pLPiSgFqhqUMm1OpNllmWU0ynXp5bm8KNdXBPs1VqqujOvXr1fDhg3zSz3XXnuteuONN5p9XrC3YSDku/DCC9Wnn37a4vPrk3HNmjVqxIgRLS7XG4CNysMztdGhIaXUFKVUiodtEbBfRPoA6J8HGijqGuBDpVRNVDal1F5dvlPAG8DYxuSxaRnGGKQL345NGuc7qe0R+KLccGLUqFHs3buXvXv3+q7QBQtQCQmsevddLrr/flgQ2GDorYGLunbly2uuAYcDBgzwWRsuX76cyZMn+6QsX+HtHMFi4Cb9+03AogbyzgYWmhNMSkTQ5hdyzjzNxhcYY5OFwLNo45K+GJs0yv0TYDhLCIYxz1AiIiKCiRMnsnz5ct8UOHs2zptvYUdREZHAgD17cN58CwRg6KHVMHs2ExYvYdXRo1oft7DQZ224bNkypkyZ4gMhfYe3iuAp4GIR2Q5M0fcRkTEi8qqRSUQGAP0B99mXBSKSDWQD3akbutPGhxiLZaId8fwF+MrR1yeLZVInxJA94z7+QgTdCIyH09bIlClTfKYITi5dQUT1KbagKWgBIqpPcXLpCp+UHw6cXLqC811VRFEbWN4XbXj8+HE2btzIBRdc4LWMvsSrBWVKqcPAGX0cpdRG4FbTfgFwhtmCUuoH3tRv0wxMYRKnzZ3Lt8OGcdddd/mkXKZm0T73XYbs2IGI2AuWWsCUKVN48sknUUqhdZBbTvRRbYT2Uuo6UDPSbRon+ugBHMBS6q7z8LYNv/76a0aNGhV0E/G2i4kwZNq0aSxdutRn5S1dupRp06Z5/QALWzIyqLxnPs6SvWz1gd+hIuJxcma40CLseZumYrRVFD5oQ9Nam8+nTuWcr7OCwr+QGVsRhCGTJ0/mm2++oaKiwiflGYrApmVkry4n9eMn+TlOyvDe79Cz3R7nGyKoNqWdog3PdrNdfzSVZ7s9zim9P1WNNu/V0jY0r7VJAG5S5UG31sZWBGFIp06dGDduHF988YXXZZWVlZGZmcnEiRO9FyxMMfzP3Auco6fV+p9pPjNv68t8HLxFO1zAYWI5TVtm3tbXNwKHATNv68tp2lJKLH8F1tKpxW1o9i90B5CKd9fXH9iKIEzx1fDQ559/zgUXXEA7f0c+a8UYfmYMv0PKLb25TOq0jv906Mgf+m4mUhSjE46wbt4ipsRs8Im84cCUmA2sm7eI0QlH+D/mMbXzz1vchsZ1VGiTsuKWHgzYiiAcycjgrP9u5eN/vIpTpNlj0v0XLqwZ81xyzTWM/2Rt0I15hhLGmguhrg+alq7F+Ob884kbFMfu3YNxuaCgAKY8Pkmb2LdpGunpTHl8EgUFsGnTlfTs+QGTH5vYojYMheBAtiIIQ7JXl3P5ipc5FycHaf6YdFEOpCx+gr7OQgYDP1FHgm7MM5Qw1mIAZAL3491ajP/85z9cccUVvhHOhpEjR1JVVUVOTsuWOZVOn8MeYCbUzNsE21obWxGEIbFLFuBAW93XU09rzpil4UtfAb9DswsOtjHPUMJY41ESkcBZwEsIX0+9q/lrMfTVxP957jmufPllezWxjxARrkhK4j8XXNCiVcapE2J4LnkyXaQDkUEaHMhWBGGIMTbZBu0GaO6YtOFLH+yA6D4hPZ3URY8TV11AF6WYdt21ZM8Y3LxhCH018XdFRbQBUvbutVcT+4rZs5n1+XI+KC9v2Srj9HSWtTnM7csW41Au4qoLSF30eFAN1dmKIAwxj02aA3U3dcxytyOOY2guY82+a4NpzDOUueGGG/jnP//ZrHOM1cQbgL9iryb2JSeXruACZxXxwCE9rTltm5OTw6FDh4Lass5WBGGIMSYtwHdowzvNGbPcce5lvIrmT8SwgAi2Mc9Q5pLNm9n+XTZfRfRrchCh6KMHcAE3AlPd0m28w1hl/B+gq1t6vWRkoO5fSEnkAP6ZmsqM3UfJ/dEDQWtQYSuCMMQ8Jp0EvI6DTyf/vMljlnHJir906MpIRy9cQTrmGcpsWVvBnFPH+d61p8lBhIqIJw9tuE/c0m28w2hD96HUhto2e3U5E9e8TF9nIbHA/wVJkPr6sBVBOGIak+749NPMGDSGj794G9fv7q3/7dO0TH7nK6/Q4fgxEqf/FMfTTwXlmGcoE7tkAXcAP6b2odPYZPwzXR+t00MDezWxrzBWGRvmvdU03rbmRWT3AIkEt0GFrQjCnOzV5fx+13ry1VGcDbx9mpfJHwCeoIrUj58M2jecUKavs4ghnPl239BkfLdJ+eThoJRYezWxjzFWGR8mlgPADxHKiG6wbfs6i2qUeCgYVNiKIMyJXbKAwcCH1HpZ9PTmYrzhuIDr0Wyig/kNJ5TZExGPULvAzGNge1MPrVqEDz74A7P7nM8rXe61VxP7GPMq47/wNMWOs/kH8IMnflBvD3pPRDxfoC0iCwWDCq/cUNuEPsYbiuFl0ZhEdn9zMd5wjGXy7ufb+I7S6XPot/gJIoB/A7uBXwN9nYWURA6gdPocum//hpQtKxAgFxiNi//b+yW5M8/nvprwUJP0zcYr0tOZAhQ8Dtmzypm9eCu/R1PScc5C+i1+gmzmkZqRQfbqcmKXLKCDs5AbgU+A4XoxhkFFMLpptxVBmLMnIp44ZyFCbbjJNoALB4hwjM4AdEbxHlpPoC11Q1IG440dyqROiCGbecQuWcBFzkJSgcuAIdQ+eICaMeghwGtoilzrodnzAv4idskC4oAPqH14CpC0+GnUYicp+v5i4JdoDuZcwJ6IBEqnzwlagwqvhoZE5GoRyRURl4iMaSDfVBHZKiI7ROReU/pAEVmnp78rIlHeyGPTfMympArtzfNBIBInDiCGo8RwlB3AW0AZtsmo3zFN5p+OSOBVIIG6E8dCre+aSGr/yHYPzb+496CNNTiROGuuixOYhmaWLWhKINgNKrydI8gBrgC+rC+DiEQAL6K91AwDZovIMP3w08BzSqmzgFLgFi/lsWkmZlNSB8KDOMiDOr7sXcBAtLecPmgPJNtkNDD0dRYxnVqvlS59K0Vbx2HMIZh7aDb+w+wgUItRoMUArwYqgG8Jbi+j9eFtqMotQGORqcYCO5RSu/S87wCzRGQL8AO0uUeAN4GHgb97I5NNMzGFsAToLQ5eNB02PGGabxSFEFddYA8JBQDz0J2L2vi5nYD7qPvACeYx6NaCMX9jnswvRrsubdGGgtzfrkNh+DQQcwT90NrKoAQYB3QDypRS1ab0M+IaG4jIXGAuQK9evVi5cmWLhKmoqGjxuYHAavkSHXH0dxXXTAx76jLudsSx025Dr2iqjGrcpfRb83LNgyeauj0A9P0SR392nHsZ8b0P+eS3B3sbWiVf/96HWHneXM5a+wn9XCW0w8Ez+uuSp/+KAraPu5QdQdyWACilGtyAZWhDQO7bLFOelcCYes6/CnjVtH8D8ALQHa2nYKT3B3Iak0cpxejRo1VLWbFiRYvPDQRWy5c1c55yaa61PG4uUF+cN9dSGRvD6jZsCk2W8emnVdbMeao4IkE5QTndrkclbdQvu71tnXwWESzy/bLb26qSNnWuiXGdihz9VdbMeUo9/bTVYtYAbFQenqmN9giUUlO81DW79Ye8QZyedhjoIiKRSusVGOk2FmK2WOnrLKyxGurE0RrLh/jehxopxcZnmIbult2/gnFPzOI0kcRSSimxRFFtLxqzkJm39eX0E22poGOda7Ju3iIiL5agdjRnJhALyjYAg3ULoSjgOmCxrp1WoPUYAG4CFtVThk2gMFmsOJQiRpUTo8pxKFVj+VBsuza2BPPCJnvRWHDQWq6JV3MEIvIj4HmgB7BERDKVUpeKSF+04aBpSqlqEbkT+Axt8errSqlcvYjfAe+IyGNojjBf80YeG5tWjWlhUy32ojFLaeiaBPu8gAlvrYY+RPNO4J6+B82U1thfCpwRKV1plkRjvZHBxsbGxsY7bF9DNjY2NmGOrQhsbGxswhxbEdjY2NiEObYisLGxsQlzRLPiDC1E5CBQ2MLTu1MbgzoYCXb5IPhlDHb5IPhltOXznmCUMUEp1cM9MSQVgTeIyEalVL2eUq0m2OWD4Jcx2OWD4JfRls97QkFGA3toyMbGxibMsRWBjY2NTZgTjorgZasFaIRglw+CX8Zglw+CX0ZbPu8JBRmBMJwjsLGxsbGpSzj2CGxsbGxsTNiKwMbGxibMabWKQESmishWEdkhIvd6OB4tIu/qx9eJyIAAytZfRFaISJ6I5IrIrz3kmSgi5SKSqW8PBko+kwwFIpKt17/Rw3ERkb/qbZglIqMCKNsQU9tkishREbnLLU/A21BEXheRAyKSY0rrKiKfi8h2/TO2nnNv0vNsF5GbAijfH0Xke/0afigiXeo5t8H7wY/yPSwiu03XcVo95zb4n/ezjO+a5CsQkcx6zvV7G7YIT9FqQn1Dc3e9ExgERAGbgWFueX4OvKR/vw54N4Dy9QFG6d87Ads8yDcR+K/F7VgAdG/g+DTgE7TIiecC6yy83vvQFstY2obAhcAoTNH2gAzgXv37vcDTHs7rCuzSP2P177EBku8SIFL//rQn+ZpyP/hRvoeBe5pwDzT4n/enjG7HnwEetKoNW7K11h7BWLQwmLuUUqeBd4BZbnlmAW/q398HJouIEACUUnuVUt/q348BW2ggXnMQMwt4S2msRYs418cCOSYDO5VSLV1t7jOUUl8CR9ySzffam8APPZx6KfC5UuqIUqoU+ByYGgj5lFL/U7Wxw9eCdbHW62m/ptCU/7xPaEhG/RlyDbDQH3X7i9aqCPoBxab9Es580Nbk0f8E5UC3gEhnQh+SGgms83B4vIhsFpFPRCQ5sJIBWuzt/4nIJhGZ6+F4U9o5EFxH/X88q9sQoJdSaq/+fR/Qy0OeYGnLn6L18jzR2P3gT+7Uh65er2doLVja7wJgv1Jqez3HrWzDemmtiiAkEJGOwAfAXUqpo26Hv0Ub6hiBFgXuowCLB3C+UmoUcBnwCxG50AIZGkS08KczgX97OBwMbVgHpY0PBKXNtojcD1QDC+rJYtX98HcgEUgD9qINvQQrs2m4NxCU/6nWqgh2A/1N+3F6msc8IhIJxACHAyKdVmcbNCWwQCn1H/fjSqmjSqkK/ftSoI2IdA+UfHq9u/XPA2iR6NyjyTWlnf3NZcC3Sqn97geCoQ119htDZvrnAQ95LG1LEbkZuByYoyurM2jC/eAXlFL7lVJOpZQLeKWeei2/F/XnyBXAu/XlsaoNG6O1KoINwGARGai/MV4HLHbLsxgwLDOuAr6o7w/ga/RxxNeALUqpZ+vJ09uYsxCRsWjXKpCKqoOIdDK+o00o5rhlWwzcqFsPnQuUm4ZAAkW9b2BWt6EJ8712E7DIQ57PgEtEJFYf+rhET/M7IjIVSAdmKqVO1JOnKfeDv+Qzzzv9qJ56m/Kf9zdTgO+VUiWeDlrZho1i9Wy1vzY0i5ZtaJYE9+tpj6Dd7ABt0YYTdgDrgUEBlO18tOGBLCBT36YBtwO363nuBHLRrB/WAucFuP0G6XVv1uUw2tAsowAv6m2cDYwJsIwd0B7sMaY0S9sQTSntBarQxqlvQZt7Wg5sB5YBXfW8Y4BXTef+VL8fdwA/CaB8O9DG14170bCm6wssbeh+CJB8/9Tvryy0h3sfd/n0/TP+84GSUU+fb9x7prwBb8OWbLaLCRsbG5swp7UODdnY2NjYNBFbEdjY2NiEObYisLGxsQlzbEVgY2NjE+bYisDGxsYmzLEVgY2NjU2YYysCGxsbmzDn/wEV8Zqat2NGyQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "h = .2\n", "xdist = 6*np.pi\n", "x = np.linspace(0,xdist, num=int(xdist/h)+1)\n", "y=exp(x)\n", "\n", "methods = ('RK45','RK23', 'Radau','BDF', 'DOP853')\n", "colours = ('bo','ro','rx','k-','r.','b.')\n", "\n", "\n", "for c, meth in zip(colours, methods):\n", " y = solve_ivp(m_y, (0,6*np.pi), (-1,0), method=meth, t_eval = x, max_step=h).y[1]\n", " #y2 = exp(x)-y\n", " plt.plot(x,y,c, linewidth=1, label=meth)\n", " #print(y)\n", "\n", "\n", "plt.title(\"Comparison of ODE methods on sin(x)\")\n", "plt.legend()\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 7, "id": "7b130b87", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABEm0lEQVR4nO2deZgU1dW438MMqzAgoKOCBjWJCkQGBjVjYjLzERXwk8UYPjWiRiNCJEYNISH5AQIRFB3ilrgk4fNTCUiiICa44kxMAlFABsOiERUVN5QEYYzDMnN/f9yqmZqe7pleqquru8/7PPV0Lbfqnr7VferWueeeI8YYFEVRlNynXaYFUBRFUYJBFb6iKEqeoApfURQlT1CFryiKkieowlcURckTVOEriqLkCarwFd8QkW+LyNOZlsNFRDqLyOMi8omI/D7T8ihNiEitiByXQPnFIjImjnIni8jqlITLYVThhxARuUhE1jl/ivdF5AkR+Wqm5WoLY8wiY8xZmZbDw/lAMdDLGPOtaAVEpL+IrHAeCntFpEpETvcc7ycixrkXtSLyoYj8UUTOjLjOdhH5zFOuVkTuSu/Xy16MMV2NMW/EU1ZETgYGAY/Fcd2Xgd0icm6KIuYkqvBDhohcD9wGzMUqq2OAXwGjMyhWm4hIYaZliMLngH8aYw5GOygixwN/A/4BHAscBSwDnhaRsojiPYwxXbGK5xlgmYhcFlHmXEeRucvkVL9AtHZNtK1Dem8S4SpgkYl/lugi5xwlEmOMLiFZgO5ALfCtVsp0xD4Q3nOW24COzrFyYAcwFdgJvA+MAUYC/wT+BfzUc60bgD8ADwN7gZeAQZ7jPwFed45tAcZ6jl2GVZa/AHYBP3f2/dU5Ls6xncAerFId6PmeDwAfAW8B/w9o57nuX4FbgX8DbwIjWmmPk4BqYDewGRjl7J8F7AcOOG16RZRzHwRWRtl/N/C8s94PMEBhRJkpwIceubcD34jzPrfztO0uYCnQM6K+K4C3gedjtHVbbRhZ/vPAn4FPgI+Bh1uRb5TTlrudtj3Jc2y7891fdq71MNApxnVi1ul8x8876/cDvwT+hP2tvQAc7yn7BvDViPvziGf7ZmAVIM52H+AznP+FLp57kmkBdPHcDBgOHIxULhFlZgN/Bw4HDgNWA3OcY+XO+TOA9sCVjkL4HdANGOD8EY51yt+AVYjnO+WnYBVse+f4t7C93nbA/wCfAkc6xy5z6vo+UAh0prnCPxtYD/TAKv+TPOc+gH0974ZVcP/EUcjONQ44shcAk7APNonSFu2BbcBPgQ7AfzkK4wTP93uolbb8APhOlP0VQL3znfoRXeEf5+w/ydneTvwK/wfOPeyLfYDfCyx2jrn1PQAc4mnXyLZuqw0jyy8Gfubcy054FGiEbF907vOZTvtOddq4g+d7vuj8LnoCW4GJMa4Vs05aKvxdwKmOvIuAJc6xQ5yyh3nO7eJ838uAM7APk74Rde8BTs70fzpsS8YFaFNAWIjtJW7y6XrHAE87P9QtQL9Mf0ePbN8GPmijzOvASM/22cB2Z70cq9ALnO1uzp/lNE/59cAYZ/0G4O+eY+2wbwVnxKi7BhjtrF8GvB1x/DKaFP5/OX/KL+P0PJ39Bdied3/PvquAas81tnmOdXG+wxFR5DkDq7S9118M3OD5fq0p/IPA8Cj7T3Tq7ENshd/J2f8VZ3s79k1it2e5Mka9W4Fhnu0jsQ+5Qk99x0W069ue7XjaMPLePADcR4RijCLbdGBpxG/iXaDc8z0v9hyfD9wT41ox66Slwv+N59hI4BVnvY9TtlPE+adh31jfAi6Mcv13ga+l+z+bbUs22PDvx/Z8/eIB4BZjzEnYHsVOH6+dKruA3m3YXI/C/shd3nL2NV7DGFPvrH/mfH7oOf4Z0NWz/Y67YoxpwJqEjgIQkUtEpEZEdovIbmAg0DvauZEYY54D7sK+qu8UkftEpMg5v32U79DHs/2B5zr/cVa9MrscBbzjyB3rWq3xMVbZRnIk0IA1KcXCreNfnn1jjDE9PMuvY5z7OewYgNuuW7FvFMWeMpFt692Opw0jz5+KfdN6UUQ2i8jlMWRr9vty2vYdYtwf4D9EvzeJ1NnaNXc7n928hY0xL2BNPYI1iUXSzXOu4hB6hW+MeZ7mfypE5HgReVJE1ovIX0TkxHiuJSL9sT21Z5xr13oUShhYA+zD2t1j8R5WYbgc4+xLlqPdFRFphzUzvCcinwN+DUzGern0ADZh/2AuprULG2PuMMaUAv2xpoIfYZXsgSjf4d0kZH8PONqRO5lrPYs1W0UyDljTxm9jLLaz8GqcdXl5Bzsu4X04dDLGeOWObFvvdjxt2Ox8Y8wHxpgrjTFHYd8GfiUin48iW7Pfl4gI9jeS8P1JoM7WrvEp9q32i979InI11hz2HvbB4j3WB2viS+be5DShV/gxuA/4vqNMpmC9WOLhi1iXrUdFZIOI3CIiBWmTMkGMMZ9g7e+/FJExItJFRNqLyAgRme8UWwz8PxE5TER6O+UfSqHaUhE5z3mruBb7wPk7TbbTjwBE5DvYHn5ciMgpInKaiLTH2oTrgAbn7WMpcKOIdHMeLNcn+R1ewPYGpzrtVA6cCyyJ8/xZwOkicqOI9HTk+T5wCfDjGN+rWEQmAzOBaRFvF/FyD/b7f8655mEiMjrek5NpQxH5loj0dTb/jb230WRfCpwjIsOce/dD7G8iYd/2BOpsi5XA1z3X/SJ2IPpiYDz2/pd4yn8deM4Ysy+JunKarFP4ItIVOB34vYjUYAe8jnSOnScim6IsTzmnF2LtvlOAU7ADb5cF/R1awxhTif3z/j+ssn0H28te7hT5ObAO6yXxD6xnzc9TqPIx7IDsv7F/nvOMMQeMMVuASuxbx4fAl7CeH/FShH1D+DfWRLALuMU59n3sQ+ANrEfO77BjNQlhjNmPVfAjsL3eXwGXGGNeifP814CvYl0tt2PHL74JnG2Mifyuu0XkU2ybj8R6UkXK/HiEH/6yGFXfDqzAun/uxT5gT4tHZg+JtuEpwAsiUuvU/QMTxQ/eGPMqVpHeiW3Tc7HupvsTlC/uOuPgPuDbYinEPthuNsZsdO7hT4EHRaSjU/7b2IeqEoHrxhRqRKQf8EdjzEDHDvyqMSaa7bWt63wZ+0P5urM9HviyMeZqXwXOEkTkBuzA2cWZlkVRWkNEfocdTF7eRrmTgXuNMZHzKBSysIdvjNkDvCki3wJrYxSRQXGevhboISKHOdv/hfXUURQlxBhjLmpL2TvlXlZlH5vQK3wRWYw1K5wgIjtE5ArsK9sVIrIRO0EkLvunY/ucAqwSkX9gByBjeVIoiqLkFFlh0lEURVFSJ/Q9fEVRFMUfQh1UqXfv3qZfv35Jnfvpp59yyCGH+CuQj4RdPgi/jCpf6oRdRpUvcdavX/+xMeawqAczPdW3taW0tNQkS1VVVdLnBkHY5TMm/DKqfKkTdhlVvsQB1pksDq2gKIqi+IAqfEVRlDwhZYUvIkeLzRK0xQmQ9IMoZURE7hCRbSLysogMSbVeRVEUJTH8GLQ9CPzQGPOSiHQD1ovIM8ZOzXcZAXzBWU7DJjBIdCo5AAcOHGDHjh3U1dW1Wq579+5s3bo1mSoCIQj5OnXqRN++fWnfvn1a61EUJTtIWeEbY97HxiDBGLNXRLZiQ6l6Ff5o4AFnQOHvItJDRI50zk2IHTt20K1bN/r164cN5BedvXv30q1bt5jHM0265TPGsGvXLnbs2MGxxx6btnoURckefLXhOzFvBmOjGHrpQ/P43DuIP2Z5M+rq6ujVq1eryl4BEaFXr15tvgkpipI/+OaH70SxfAS41th4N8leZwIwAaC4uJjq6upmx7t3705tbW2b16mvr2fv3r3JipF2gpKvrq6uRRvGS21tbdLnBoHKlzpByli0eTM9amrYXVLCngED4jon7G0YdvlaEMtfM5EFm33nKeD6GMfvxZOGDJuY4Mi2rhvND3/Lli1x+aLu2bMnrnKZIij54m2vaITRx9iLypc6aZVx9Wpj5s61n6tXG9O5szEFBfZz9erMy+cDYZSPdPrhOxlxfgtsNcYsiFFsBXCJ463zZeATk4T9PiwUFBRQUlLCwIEDOffcc9m9ezcA27dvZ+DAphwhv/71ryktLeXf/27KlFdZWYmIsGvXLgCqq6vp3r07JSUllJSUMHv27EC/i6KkhTVrYNgwmD7dfj7wAOzfD/X19jObesU5hB82/K9gE2f8l9j8pzUiMlJEJorIRKfMSmyihm3Y6JTf86HeNpk/H6qqmu+rqrL7U6Fz587U1NSwadMmevbsyS9/+csWZR588EHuvPNOnnrqKQ499FAA3nnnHZ5++mmOOeaYZmXPOOMMampqqKmpYcaMGakJpyhhoLq6uYIH6NABCgrsZ69eMG+efTAogeGHl85faZ7nNFoZAwSeZOSUU2DcOFi6FCoqrLJ3t/2irKyMl19+udm+pUuXctNNN7Fq1Sp6927K+X3dddcxf/58Ro+OO5udomQn5eVWse/fbz8vucQu1dVW2V97bdOxVaugTEPYB0FOz7StqLDKfdw4mDGjufL3g/r6elatWsWoUaMa97311ltMnjyZp59+miOOOKJx/2OPPUafPn0YNKhlrpY1a9YwaNAgRowYwebNm/0RTlEySVmZVeRz5jQp9LIymDYNdu1S806GCHW0TD+oqIBJk+zvbvp0f5T9Z599RklJCe+++y4nnXQSZ555ZuOxww47jJ49e7J06VKuu+46AP7zn/8wd+5cnn766RbXGjJkCG+99RZdu3Zl5cqVjBkzhtdeey11IRUlE6xZYxV4eXmTko8ksvdfXh6sjHlMTvfwwZpx7r7bKvu7725p008G14b/1ltvYYxpZsPv0qULK1eu5J577mHRokUAvP7667z55psMGjSIfv36sWPHDs444ww++OADioqK6Nq1KwAjR47kwIEDfPzxx6kLqShBEzlQG8s+H633rwRCTvfwvTb7igq7+GnW6dKlC3fccQdjxozhe99rGoc+/PDDefLJJykvL6d3796cffbZ7Ny5s/F4v379qK6u5ogjjuCDDz6guLgYEeHFF1+koaGBXr16pS6cogRN5EBtdXVsZR6r96+klZzu4a9d21y5uzb9tWv9q2Pw4MGcfPLJLF68uNn+Y489lhUrVnD55Zfz4osvxjz/D3/4AwMHDmTQoEFcc801LFmyRGcRK9mJa6pxPXHUVBM6crqHP3Vqy31uTz8VImf6Pv74443rmzZtalwfNGgQ7777bovzt2/f3jjLdvLkyUyePDk1gRQlDLimGq8NPx4i7f5K2shpha8oSsAkaqpx7f7qohkIOW3SURQl5ESz+ytpQxW+oiiZQ+3+gaImHUVRUiMVG3yydn8lKVThK4qSPH7Y4NVFMzDUpKMoSvKoDT6rUIWfBLHCI8dLeXk5L730UnqEU5QgURt8VqEKPwniCY+sKHmBhknIKvJD4a9Zk7bY22VlZY2Tq1588UXKysoYPHgwp59+Oq+++ipgg61dcMEFnHTSSYwdO5bPPvus8fxJkyYxdOhQBgwYwMyZMxv39+vXrzGmzrp16yjXnpMSVtwomKrsQ0/uD9qmcWKHGx75iiuuAODEE0/kL3/5C4WFhTz77LP89Kc/5ZFHHuHuu++mS5cubN26lZdffpkhQ4Y0XuPGG2+kZ8+e1NfXM2zYMF5++WVOPvlkX+RTlKxDZ92mFV8UvogsBP4b2GmMGRjleDnwGPCms+tRY0wwufwSCegUJ7HCI3/yySdceumlvPbaa4gIBw4cAOD555/nmmuuAeDkk09uptCXLl3Kfffdx8GDB3n//ffZsmWLKnwlP4nWOVN8xS+Tzv3A8DbK/MUYU+IswSVuTcOgUqzwyNOnT6eiooJNmzbx+OOPU1dX1+p13nzzTW699VZWrVrFyy+/zDnnnNN4TmFhIQ0NDQBtXkdRcgL1+Ek7vih8Y8zzwL/8uJbvpHFQyQ2PXFlZycGDB/nkk0/o06cPAPfff39jua997Wv87ne/A2xwNTcl4p49ezjkkEPo3r07H374IU888UTjOf369WP9+vUAPPLII77JrCihRT1+0k6QNvwyEdkIvAdMMcZEzeUnIhOACQDFxcVURzzlu3fv3hhpsjXq6+ubyg0caBeAOM6NB/fan//85+nfvz8LFy7k6quvZuLEicyePZuzzjoLYwx79+7l4osvZtKkSZxwwgmccMIJlJSU0NDQwHHHHcfAgQP54he/SN++fTnttNOoq6tj7969/OhHP+Lqq6+mqKiIr371q82/TwLU1dW1aMN4qa2tTfrcIGhNvqLNm+lRU8PukhL2DBgQrGAOYW8/CJ+MRbfc0nTf9u0LnXyRhF2+FhhjfFmAfsCmGMeKgK7O+kjgtXiuWVpaaiLZsmVLi33R2LNnT1zlMkVQ8sXbXtGoqqryT5A00EK+1auNmTvXmHvvNaZzZ2MKCuzn6tXhkC+EhF1GlS9xgHUmhk4NpIdvjNnjWV8pIr8Skd7GGM3lp/iDd8BPBBoa7OLTQL0SgXrTZCWBKHwROQL40BhjRORU7NjBriDqVvIE74Bfu3bWDizSZAtWBeUfGsM+a/HLLXMxUA70FpEdwEygPYAx5h7gfGCSiBwEPgMucF49FCU1XEXeq5dVPq4Suu022LWraeBPFZR/pMHVWQkGXxS+MebCNo7fBdzlR12K4lK0eTP86EfRlbxXAc2bpwrKT1xvGrfd1Zsma8j9mbZKztKjpqa5It+1y07xj0QVlL9oDPusRRW+krXsLimJT5GrgvIfjWGfleRH8DSfccMjDxgwgEGDBlFZWdk4Kxbgr3/9K6eeeionnngiJ554Ivfdd1/jsRtuuIE+ffrwla98hYEDB7JixQoA3n77bSoqKhg8eDAnn3wyK1euBGD79u107tyZkpISSkpKmDhxYuO1hg8fzqBBgxgwYAATJ06kvr4+oBYIB3sGDIh/Up03wFcag+kprTN/PlRVNV+vqrLr0Hxd8R/t4SeBG1oBYOfOnVx00UXs2bOHWbNm8cEHH3DRRRexfPlyhgwZwscff8zZZ59Nnz59OOeccwC47rrruOqqq9ixYwdnnHEGO3fu5Oc//znjxo1j0qRJbNmyhZEjR7J9+3YAjj/++Mb6vCxdupSioiKMMZx//vn8/ve/54ILLgioFUJCoj1N9TDJKKecAuPGwdKldn3sWDAGli+3yt49pqQH7eGnyOGHH859993HXXfd1RhX57LLLmuMiNm7d2/mz5/PTTfd1OLck046icLCQj7++GNEhD177HSFTz75hKOOOqrNuouKigA4ePAg+/fvR0R8/GYhxumhF22OOlm7dTReS0apqLAKfdw4q+CNsd6zN94IY8bYYxUVtmxVFbxw20f6NuYj2sP3geOOO476+np27tzJ5s2bufTSS5sdHzp0KJujKKcXXniBdu3acdhhh3HDDTdw1llnceedd/Lpp5/y7LPPNpZ78803GTx4MEVFRfz85z/njDPOaDx29tln8+KLLzJixAjOP//89H3JsODpoQ8qLIQhQxLroesAbsapqIBJk6wlbvp0u2/OHOjSpalMVRXMH7uGFf8ZD388qG9jPpH1PXwRiboUFRXFPNbWkm5+8Ytf8JWvfIUpU6bw8MMPIyIsXryYyy67jB07drBy5UrGjx9PQ0MDRx55JG+//TYbNmxgwYIFjeYjl6eeeor333+fffv28dxzz6Vd9ozj6aHLgQOJ99A1Q1PgeO32YNfvuMM+t2+/3a5Pnw6FhdbEM2OGfQO461vVFNYf0LcxH8l6hR8rZsSePXtSiQuUEG+88QYFBQUcfvjh9O/fvzHKpcv69esZ4Angdd111/G3v/2Nv/zlL4299d/+9reMGzcOsFm06urq+Pjjj+nYsSO9evUCoLS0lOOPP55//vOfza7fqVMnRo8ezWOPPZaw7FmHJ6Kiad8+uR66ZmgKFNdu7w7Qunb7kSOtOccY2+tfvtzq9Tlz7BvA8ZeX09C+vUbP9JGsV/iZ5qOPPmLixIlMnjwZEeHqq6/m/vvvbxxk3bVrFz/+8Y+ZOnVqq9c55phjWOUkfNi6dSt1dXUcdthhfPTRR43eN2+88QavvfYaxx13HLW1tbz//vuAteH/6U9/4sQTT0zfFw0Lnh76xspKVdpZgNduf+ONTYO0Bw/CsmV2fe1aW7ZDB9vzr6yEBWvK7D123saq6srUgydF1IafBG7GqwMHDlBYWMj48eO5/vrrATjyyCN56KGHuPLKK9m7dy/GGK699lrOPffcVq9ZWVnJlVdeyS9+8QtEhPvvvx8R4fnnn2fGjBm0b9+edu3acc8999CzZ08+/PBDRo0axb59+2hoaKCioqKZy2ZO43jm7NFX/Kwh0m5fUdE0OOsybpx9AFRUwIIFMGUKTJx4Nr/61efVg8cnVOEnQVv+7l/72tdY63ZZIrjhhhsAWsS279+/P3/7299alP/mN7/JN7/5zRb7i4uLY9ahJIAGVYuPFNupqgruvtsq+7vvbqnw165t7qHj9J/42c+OpXdve473uJIcqvCV/EV98uMjwXaaP9/a7b3ulWPHwv/8D8yebfe7vXW3TDSL5/XXQ03NO8yZ06/xrUBJDbXhK/mL+uTHR4Lt5B2kBViyxNrt3TmBrk2/rRfUqipYseKoxrcCr6ePkhxZqfA1snJ85Ew7pSsUguZQjY8E28k7SDtjBjz6qB2Y9fbQKyqi9+pdXJv9zJlbmD27+WQtJXmyzqTTqVMndu3aRa9evfJnZmkSGGPYtWsXnTp1yrQoqZFOs4sGVYuPJNop2iBtIrg2fZHdjddz3wrUtJM8Wafw+/bty44dO/joo49aLVdXVxdqZReEfJ06daJv375prSPtpDvZhkZ9jI8E26mtQdq2cHv/1dXNxwS84wJr17b+lqC0xK+MVwuB/wZ2GmMGRjkuwO3YBOb/AS4zxryUTF3t27fn2GOPbbNcdXU1gwcPTqaKQAi7fKFBQyFkHV4XSldJRw7SJoI34FpFhQZZSwW/bPj3A8NbOT4C+IKzTADu9qleJdfRUAhZR6SLZbyDtLGIHBNI5eGR7/iV4vB5EenXSpHRwANOHtu/i0gPETnSGPO+H/UrOY6aXbKKaGaWRE060c5PZUxAsYhfnhyOwv9jDJPOH4GbjDF/dbZXAT82xqyLUnYC9i2A4uLi0iVLliQlT21tLV27dk3q3CAIu3wQfhlVvtTxS8bFi4/mxBP3Mnjw7sZ9Gzb04JVXunHhhe+kLN+GDT2YNas/o0a9x4oVRzFz5pZmdWWKMN7jioqK9caYoVEPJhtgLErAsX7AphjH/gh81bO9Chja1jVLS0tNslRVVSV9bhCEXT5jwi+jypc6fsn43HPG9O5tP6NtJ0tVVVXjtdbdudqYuXPNujtX+3JtPwjjPQbWmRg6NSgvnXeBoz3bfZ19ihIeNMxC0njt7JMm+RsKYe1aeHLmGkqnWvfc0g4deHL+KlatLVPTToIEpfBXAJNFZAlwGvCJUfu9EiY0zELKpMvOPnUqMK+6mXtu6d5qSqfp/UkUX7x0RGQxsAY4QUR2iMgVIjJRRNzwjSuBN4BtwK+B7/lRr5LDBJ1oXMMspEyk772vs2J1VrQv+OWlc2Ebxw1wtR91KXlAJnrb6u+fEn773rfAcc/98+xqup5TTqnn96CTsOInK2PpKDlOJnrb6u+fEn773kelrIyGqdMYPqus8e3BfdCccoqP9eQwWRdaQckDMtXbVn//5iQwiJ0O3/topHNwOB9Qha+EDw1qlnlCPIitk7CSR006SjjRROOZJQ6z2vz5LQdmq6pIe97ZtA4O5ziq8JW8IVMKKiuJwysmMtFJEPZ07+CwxslPHFX4Sk7jVfKuglqwAEaOtJ9eBaXK30Mcg9iZCGoWyOBwDqM2fCWniQytO20aTJkC3/iG/bz1Vg25G5M4BrGDtqcHNTicq2gP3y+CniikRCXSbOMq+XPPtb3QefPg4ovhmWfs57x5VvmPGdO856i9/fhQe3p2oT18P4j0aLjtNti1Sz1MMkC0ZBnz5sF559le6Pjx8MQTTQpqxAh48EHo0qXpGtrbj4+0T7ZSfEcVfiq4fspvv93k0bBvH0yeDA0NqvwzQDQ/7WnTrNIfPx4eesiaca6/Hnr0sGad8ePhscdg7Fi45hp7zpMz11D692roVK73LQat2dNV4YcTVfjJ4u3VFxRAodOUIlbxNzS0VP4h8mUOJUlGq/TmPAX7OWJEU49+3rwmRXTrrXYb7Oett8LBg/Cd78A559hz7vtOU2RGvW+xybQ9PfK+g4ZZaAtV+IngVUheP2WAK6+EY46BXr3g2mvtMa/yT0cC7lwihYk+kWacBQtsT378eHjkEavEvYpo8GC45ZaWNvsOHaBPH3jvd9U0HNhPuwbrg/76wmoe+UuZKpGQobluE0cVfrxEs9N7p/9fckmTgvrSl6xy9yp/DcjVOtEm+sSp8L1mnBEjmpttvvMdu3/w4OZvAJG9wnHjYNkyuz333HLq9nWgU7v9mMIOTP59OVOX+fptFR/QMAuJowo/XiIV0q5dsaf/e93ZXOXvlnHeEoqKivQB4CXF+Dle98Dx462yd/e3ZVeOtEXzeBmjz13FBUdU84ePy5m6TBNtuITGjOL8jyrKy5k0qUzDLMRLrFRYYVhCleJw9WpjOnc2pqDAfq5endI1DnbsmNw1AiTw9G2rbQq7eNvFK5+bBm/6dH9S602fbgzYz2QJY/q7SBKVMV2pDGMRVb5m/6POZnj31b7dd1/kyzC0kuLQrwQow0XkVRHZJiI/iXL8MhH5SERqnOW7ftQbCK5/PaQePtfzliAHDmiSjUgSiJ8zf75Nkg1NJplp06Br19Sn26tveWwyMbu2BZ7/UcO+/dz1rWoNsxAnKZt0RKQA+CVwJrADWCsiK4wxWyKKPmyMmZxqfYESbSBx2rTkr+cxW5jCQjXppMApp8DYsf0pKbHmBNf10lU+yboHRvqWf/ihnZS1fHnzAd589gTJeLRK539UX7efdh06cPzl5Y1yqVto6/jRwz8V2GaMecMYsx9YAoz24bqZx+9EHJ74JBsrK9VjJwUqKmDmzC2MGwe1tc2VvXs8GYUcac+/4ALrbLVkid3WhBsheANy/kcFN86hoKr523ay9z1fEGvySeECIucDw40x33W2xwOneXvzInIZMA/4CPgncJ0x5p0Y15sATAAoLi4uXeL+0xKktraWrl27JnWuS9HmzQz64Q+RAwcw7duzsbKSPQMGpHRNr3xHvfUWPWpq2F1S4tt1/cSPNkwntbW1LF06kAcf7Mf48du5/PLtaalnw4YezJrVn1Gj3mPFiqOYOXMLgwfvjku+MLcfJC6j2xZuG0RuZ1q+oAmjfBUVFeuNMUOjHoxl3I93Ac4HfuPZHg/cFVGmF9DRWb8KeC6ea2ds0NY7eJjgQGK8rL/rrtQHgdNMGAekvCxYsMHXgdrWSGYQN+ztZ0wUGdv4vd98c8t2fu45uz8Q+UJGGOWjlUFbP9wy3wWO9mz3dfZ5Hyq7PJu/AcIblspvu30MetTUJO13nq94XQKrqmDWrP7MmGFnyroDdukYQIw0YeRsdMY4Jr9lenatkhp+2PDXAl8QkWNFpANwAbDCW0BEjvRsjgK2+lBveggogfbukpI2E0wozfEm3Fi7Fi666G3mzWt6CKQjLnpeJdzIRPJ4JVBSVvjGmIPAZOAprCJfaozZLCKzRWSUU+waEdksIhuBa4DLUq03bcSR6ccP9gwYkLqbZy6QQFhpr0tgbS387nfH+DJQ2xp5lXAjoN++kjl8mWlrjFkJrIzYN8OzPg3w3y6SDoJMoB1HgomcJon4Oc1n1L5HRUW/tIrY+ABZs4Y/z66m6znlVEwuy00XzSxPHh+aWcAhRkMrRCPfFXFQJBE/x2tPv+OOo6iqCsB+7DyYzti3n7onO7CeVZROLsvNYF1Z/NvXYGptoxmvXDKdsSrT9WeCBE0IkfZ01w8/7fZ058HUrqGeTu3288SPqzM3yzRgsinxeyhmAYccVfjQZFqYPt1+Bq10M11/pmgjUXaksnFn1Lr288GDdwdjT/c8mNp17EDx/5QzZ441LeW6MvEOlEP4J555TX75cH8SRRU+ZN47IdP1Z5JW4udEKptTTqHRK8clkJmVngfT+vmr+OnjZXkTZyfbes0ZnwUcctSGDymH5s36+kNKqOKdl5VRVVeWlzlcMx47J040x27baA8f2jQt+IHXPOFGeqyqgpEjoaqujPXzV/HnM239VXVlobSRZoIwvaLnlYumh2zpNefr/UmIWFNww7CkPbRCmsImGNNyCvpzzxnTvbsxEybY9UMOOWCKioyprLT7i4rsscrKlvHG0zVtvS3CMG28tTj3mZQvnhADYWi/tmhLxqDj30cS9jYMo3ykOx5+VpLmgdJI+zOAMfDww+4+gwjs3m33i8Bnn8GUKdak7XUrC+sAWboJ8yzXbBvMTBbtNecW+WvDTyGHaiy8Ez/cP8bYsTB0KGzcaGOqV1W5k4bepV+/fo12UWhKz/fM7DV86Y/VzF1dzrS5zdPr5fpEEm8busrG3T91anjinYdqfCGNaOyc3CJ/e/hpmEYerVe/f78dFpg0yW67ttBHH+3DHXfY9dtvp3H9oxVrWF47jIqq6Tx5cBjPzF6T871IL942dJWN9zuHKd55mMYX8pp8nMOSJPnbw/dpGnlrvfq1a+2zZMqUJqW+bJk9b8ECwRjo0cOac4yx55+5vpqClfsppB5Tv5+ffqWa88aVMWiQvV6uZ17Kpp5z3kTRDDNJhOfIZ/Krhx/ZE0ggh2osWuvVHzxoFfzs2TZzkptrZu1amDNnE8uXw7PP2jLLl9usSnP/Vo5p34F6sZN85v6tnBEjmq7nkrW9/Th6Y9nQcw7z+EJekc9zWJIh1mhuGBZfvXQ8me5TTTgSywNn2DDrbeNdjyznenFEG91vvK7He6iy0phDDrFeKu61g0j4EUvGlIjzHrTmmZNW+RIgl710gk5y0hpttqGP/+tkCOM9phUvnYwr9dYWXxX+3Ln2RwH2c+7cpK8dzVWtc2d76S5d4nNhi+eH0lo948e3LOv3H9L3H3Mc9yARN8Aw/tm8hF0+YxwZI9yTM+2K2UK+tojiXh3UQyuM97g1hZ8/Jp0UBmkjY7pUVFhL0Lnn2unmY8bYSw4bBoWFzcul4sIW6RIHtp7Bg+Ghh2DBArsva8w7cdwDdQMMlqLNm1u4J2dbOIVoptl8cZtNmFhPgjAsvk+8SnKiVawez/jxifXq25QvARkqK40RsTKkq/eVlt6Lj72xjPWuIr5DLPknTNiWAeES4/XvfjfmW1cyeXz9JpV7HK9pMBWyrYfvi2IGhgOvAtuAn0Q53hF42Dn+AtAvnusmqvC9fzz3Rqy7c7WpHp64km/NTt+7d9OM2LZs9bFI9IcSTam4D5x0/SGD+jEna0LIyJ8tis04lvwLFmwIXr4EWX/XXVFt4EEoy3hI9R6n+6EVVKcoEdKq8IEC4HXgOKADsBHoH1Hme8A9zvoFwMPxXDtRhe/+SNfdudq8/t3vmq3X32s+pbOpb5f4gE5b9vNUbZwJ/VCi/ADcOpN94PguY4oko2AyovBjjENEkz+Mvb9IcsKGH4Os7OH7MAidboVfBjzl2Z4GTIso8xRQ5qwXAh8D0ta1kzHprLtztfmUzuYg7cw+Ck2DtIv6uhqNeLxvpk+3JpzKyubnJqpk4/6htNGjdGV0lb6ff86gFVaivbGw9PBdIuXPGoUfQVZ56cQgqIdWJhwb2qI1hS/2ePKIyPnAcGPMd53t8cBpxpjJnjKbnDI7nO3XnTIfR7neBGACQHFxcemSJUsSkueYRYv43G8WUkhDsl9JURQl4xzs2JGNlZXsGTAgofMqKirWG2OGRj0Y60kQ7wKcD/zGsz0euCuizCagr2f7daB3W9dOqYcvBeZTOput19/b6iCh++SvrLT7vYOh8frUJ0MqPfxouL3LYcP8652pDT8+stmGH5Y2jEVeypdGG74fbpnvAkd7tvs6+6KWEZFCoDuwy4e6m1FVBcNnlbH1zlWsGHotj31/FWc8MIGrtk+jqq6sMRfnKafY8AdXXdXkYjllCmzaZCeBXnwxPPhg00zZZ5+1M2G9bl6BxXSJI1a/d4r/2rX2u4XdHc3r6hoZJC3bXDFjuZK+8kq3zAqmZCc+RACISawnQbwL1ib/BnAsTYO2AyLKXE3zQdul8Vw7FS+dBQs2NPbcJ0xoaeeOnLXqery4A7LpHAw1xr+eQbTepV8zctPZu/LDxpqXvT+fCbuMKl/ikM4evjHmIDAZOzC71VHmm0VktoiMcor9FuglItuA64GfpFpvNKZObepluQmu582D4uKmmPNuj3f5crjmGttxHjECnnjChiZ+6CH7cM1orz4BovUuly+3wdvCHIsm6yb3JEDkRD2g8e0yrGSjzEri+DLT1hiz0hjzRWPM8caYG519M4wxK5z1OmPMt4wxnzfGnGqMecOPetvCG4TrBz9oUvDeUMVeJT9wINx6q31IVFVlh2nB+5DzsnGjnTh5++3N/8gZ+RPHCJiWDUHSksGd5blhQw8gvGY1LzozNU+I1fUPw5LqTFuvH67XzOGuuyaaTKQVTNeroJ8um77I2Mqgc6p+0mF8nXaxbb8v4xOX2sLbhmGZbOUlzPfYmHDKRz7G0tmwoUczM4E35rw3VPHUqXD99c178mE03TTSRnhhr4mnosIOOovAjTdmyGziCV/bsK8pfG1VlR1cPu+83AwvXFEBo0a9l1VvL7n6xqU0kbMK/5VXujUqt7Vrm2LOr10L997btO4SaiXvEkce3kgTT0WFNWW5WbcC/xM7AdMa2hVQ19CB9d3KARv73xj78HXlDLv5LBGqqmDFiqMak6Nkw4MsMqFLNsgcDzo+4SFW1z8Mi+/B00JEUvIlMQsvlRAMvrWh41e87s7VvpoMwnqPI/3wMxmaoC3cNgxTOAUvftzjdH63MP4GacWkk78pDrMRN7ywm86tjRDP3qxMYE0oY8bYtxtofiytlJVBWRmlwKSdNCZuz1WTgWtWE9kNNH97Cet3bi0sdVhljsmaNc1Sl2ZT2sx0owo/m0gwD2/kn3jZMqv0b7zRevEE/aPPlxywrmnQm20v7N81mjkz7DJHJUaOW+/4RC53NtoiZ234OUsCs/Aybc/32k7dt41p06Br1xwapI0jR28oceQu2rw505L4S4wct7k6PpEoqvDzCPdHH5R/vte3e+1aq+znzbP7c2KQNo5B9FAOGHrkHvTDH2bfw6o1omRV04TzTajCzxO8P/qf/cy6ao4ZY/ena5KN13ZaW2uVfaSdOPSeUa0RozfpJZQTmjxyy4EDUeXOWqLEntK0mU2oDT9PCMqe7wan8/65RozIUdtpHIPooRwwdOSur9tPQ0Eh7Txyu29jWf0gdpwEXHJmfMIHtIefzSRgP27Nnj9oUMvyVVWwePHRLQ+0QWSPdsECG7pi/PgctJ3GEckUQjihyZF7+xVzGNnhaarqrNyhePtQ0kssf80wLOqH3woppkKLFnbCr3juzz1nzPDuq83DJXNNGasbs4P57dudLfc4jCELXNyosmGUzZjsucdhgnwMrZDzxGE/jkXkINby5XbW69ix8I1vWNv+0qU24qhbvrVBxsiByYpOa1heO4zzaqZTVTCM68vsG0g+2k7DPmA4ePDucL19KGlFFX62EsUbIV5aC6m8apVN/OISz2t+pBnnLz+vpqB+P4XU065+P68vrG4sm/UDtQkS9gHDDRt6qLtiHqGDttlKgpOwvMRSuBs32j/+7bfb3v7xx5/Mtm0t4w7dcgv86Ed23R3gmzYNzj3XBkPbtrKcP7fvAA37aVfYgcm/L2fqRfnZewzzgGFVFcya1Z9ly5pkyqW8BFlFxOzgdKEKP5uJ8EZIFq/Zwf3jn3MOvPRST7p0sWXctJDGwMyZTevLl9vz582zyv7BB2H8+DLaT7IPo4LycqbWlWXnFP0cZ+1amDlzCxUVJUCWh1PIZmLMDk4HKZl0RKSniDwjIq85n4fGKFcvIjXOsiKVOhX/iTQ7gP3dDRnyLwoLm3LkGidr2O7dLTOITZtms4ZNn24/q+qaZgTnmxknFmGbhDV1atM4jUuu36uw3QMgpfG4REm1h/8TYJUx5iYR+Ymz/eMo5T4zxpSkWJeSJrx/cFeB2zj6L2NMOeec0+RHDy3Xx49vPqlKTQPRccc63HaJDG6npJ9Q3oMEgyKmQqoKfzRQ7qz/H1BNdIWvZAne3r7b0ejQAU4/3dr2RZrs/O56ZaVV/DkRaTGNhHISVp4RynuQwnhcoohxUz8lc7LIbmNMD2ddgH+72xHlDgI1wEHgJmPM8lauOQGYAFBcXFy6ZMmSpGSrra2la9euSZ0bBH7LV7R5Mz1qathdUsKeAQN8uebq1R2ZP7+UmTO3ADB9+kDAcMklb/HAA/0Aw5w5NvjWrFn9mTlzSwsTQTrJ1nu8cGE/HnywH+PHb+fyy7cHKtPixUdz4ol7G+9TbW0tr73Wl1de6caFF74TqCzxkK577Nc98FO+yHsD1osq0XtTUVGx3hgzNOrBWA767gI8C2yKsowGdkeU/XeMa/RxPo8DtgPHt1Wv0YlX8ZPiJKxYTJiwrXEizs03N+XFHTGiad1NoBJEHuBIsvEeZ3oSVuTkN3fiVdgmXLmk4x77eQ/8lM+vRC20MvEqpZmwwKvAkc76kcCrcZxzP3B+PNdXhR8nSWTCioe8asM0EClfWLJKeRVe9+77QqvsjfH/Hvt9D9IlXyoPo9YUfqoTr1YAlzrrlwKPRRYQkUNFpKOz3hv4CrAlxXoVLylMwvKVbI0N7xcJJJiHzE3C8sb2GTXqvbwaQwjLPYhFuuMupTpoexOwVESuAN4CxgGIyFBgojHmu8BJwL0i0oB1A73JGKMK308CHPSJSYC+xKEk2vePwOsN5Y0q6v6pg4pU6U0GcscdR1FVlcMDxxETmsI8EQ7SnxUuJYVvjNkFDIuyfx3wXWd9NfClVOpR4sCnSVhJE82XOJ8UfoLfP3D3QEfxre9WzrhZZY31HnroFsaNK8m8p0o6yLJOSLQJkH67N2ssHcUfwmJWyhQJfn+ve+CMGWmet+DJcHXy9cN4cuaaxnoGD94dKpOGrwQ4ockPgjA3aWgFxR/CYFbKJNG+fxsKJrDE2h7F1579lO6tBpruT5hMGr4S4IQmPwjC3KQKX/GPTJuVMk2C3z/d9lqwYwXDupVT6lF867uVs2p+bodQALQTEgVV+LlKQNH3lOQIwl4Ldqxg+Lgynpy/itK91oY/3LHh5wX53gmJQG34uYjHZsuwYfnrJhlignIPdK87fFYZMz6b1qjsc9KEkwShDKaWRlTh5yJZNliVj0TmGHYVTGQgOz8UT+hy6oaIyOQ9uZ7XVxV+LhKUx0y+T7TykXQqnsixAs1q1USg3lIhQG34uUgQg1VZ5uMcdtIVxTGosYJsJjBvqRCgPfxcpawpAUlaULOR7/hlevHapd2xAnd/2EIJhIF8egNSha8kR75PtEoDfiker3nIHRPwmodyPatVInjfgGbPbnrLCkTpZ8Akqgo/H0jHD8s1G82Zo+YcH/Aqnq5d7ctZpE0/3gHcfLNLp0LGgqk5JtGGn02nvqK5J106vYTUhp/rpNPWrj7OrVK0ebNt/zjGUSIVj5sn2FU8icbZySe7dCpkLJiaYxJtZ+o5sG8/2xdWc3xZWdpjKqnCz3XyPahZplizhkE//CEcPBjXg9areCIHcJPpoQcxi1dJAU/Yh3aFHZj8+3JOOTL9KRfVpJPr+G1rV1fM+Kiupt2BA0kPaic6gOsdqHV7idOmWfNQoHZpJT48JtGCqlWcck1ZIPMkVOHnOn7a2nUGb/yUl9PQvn3SD1q3hz5smE0Y71XW0Wy83oHatWutsn9m9hq++do8KjqtUc+cMOJ40lXVlQXmJZSSSUdEvgXcgE1ycqoTBz9aueHA7UAB8BtjzE2p1KskiF+2djUPxU9ZGRsrKxmyZ0/CcyEi7bhjx8KYMbB8ud12j0UmUVm61JYdOhTar1vDH+uGUfC/+2FRBypWraJiqt6rsBH0PIlUe/ibgPOA52MVEJEC4JfACKA/cKGI9E+xXiVZUjHJqCtmQuwZMCCpuRDeAdyKCli2DERgwgSr+N1jp5xiFfxVVzWdu3+/fZG7fkg1BQd1nkTYCdpLKNWMV1sBRKS1YqcC24wxbzhllwCj0by2wZOqx46Gmw2ESM+Rigq45hprlevSpfkxY+Dhh6G42Jp+OnSAKVNgwR3l/FdhBwrIjljwgdBGBFnvG5NLutNOBu0lJDbJeYoXEakGpkQz6YjI+cBwJ78tIjIeOM0YMznGtSYAEwCKi4tLlyxZkpRMtbW1dO3aNalzgyAT8h2zaBHHLlyINDTQ0K4d2y+/nLe//e2Y5V0ZizZvpkdNDbtLSmyvNSTkyz3esKEHs2b1Z9So93j00T6AcN55O1ix4ihmztzChg09ePDBfnTsWM+8ef9g8ODdbNjQg5XT9zDz67+ny8gTY963fGnDos2bGfTDH9LuwAEa2rdnY2VlizZx23nmzC2NbejdTqd8flJRUbHeGDM06kFjTKsL8CzWdBO5jPaUqQaGxjj/fKzd3t0eD9zVVr3GGEpLS02yVFVVJX1uEGREvtWrjenc2ZiCAvt5773GzJ1r90ehqqqq5TkxymaCfLjHzz1nTO/e9tPd7tzZGDBm+vSm48OGGVNU1FTOLXvzzemXMZ34Jt/cufY3DPZz7tyoxdz2nD69ebunXT4fAdaZGDq1TZOOMeYbqT1veBc42rPd19mnBI3XJNOrF1x7bdvmHR2ozSiRNl6wt+v0060J5447rI0/MhG6axZQ33uHONMd5vqEtSDcMtcCXxCRY0WkA3ABsCKAepVouEHVdu1qqci9A7pr1nDMokX2waADtRnDGzffVejLlsGzz8IFF1gbvosGRmuFON2Tcz2QWqpumWOBO4HDgD+JSI0x5mwROQprxhlpjDkoIpOBp7BumQuNMZtTllxJjcgeT69eTQO6BQUgwrEHDsCiRXDbbfYBoQO1GSWyt3/vvVbpr13b3Msj13qlvtGGe3JQLpKZGBx2SdVLZxmwLMr+94CRnu2VwMpU6lJ8JtLjxmu6aWgAQIyx+3btsm8FSkbJWNyXPKE1F0m/8wx7HyTpjp/jRWPp5DORPR63x+/08BsOHKCdmnGUPCGoB2q6kt3Egyp8xRLZ4we2L1zIcZdfrmYcP2jDB1zJLzI1OKwKX2kiosf/9r59HKfKKXU0HaQSQaaimWrwNEVJN5oOMuvwRh91STkxieMFt/6uNRnLsqUKX1HSTZAxiDR8tS94o49C08CqmyYyYTyRZk++fhhPzlwTfJYt1KSjKOknqBhEajryDd8HVj1vee3ZT+neaqDp3gRl0lGFryhBEEQ6SJ0V7Su+DqzGOdM33ahJR1FyBQ1f7Su+zrr1MxFRCmgPX1FyBQ1f7Rt+zrptmlnb9JYX1MzaSLSHryi5hBsrSZV967QxuO1nYhLfB4BTQHv4iqLkF3EMbvs56zaTM2sj0R6+oij5RYLzIvzwyfcOAE+alLn4R6rwFSVo1Fc+syQ4uO2HSSYsYZfVpKMoQaK+8pknwcHt1kwy8UyaDirscjxoD19RgkTDLISDBAe3UzHJ+DkAnCqq8BUlSNLhK68morSTqEnGa/d3s5Z57f4VFcG7ZELqGa++BdwAnAScaoxZF6PcdmAvUA8cNLEyqitKruO3r7yaiNJOayYZkejnZDLJSWuk2sPfBJwHPB9H2QpjTIkqeyXv8dNXXk1EacdrknF76F6TTDSPHa/df8aMzNnsI0lJ4RtjthpjXvVLGEVREkTDKaQdbyJ5t+fu7t+woUejx06k+2ZFBXx/6BoOzJnH3HPXZFzZA4jxpr1P9iIi1cCUVkw6bwL/BgxwrzHmvlauNQGYAFBcXFy6ZMmSpGSqra2la9euSZ0bBGGXD8Ivo8pnKdq8mR41NewuKWHPgAEJnattmDgbNvRg1qz+jBr1HkuX9uXyy7czbtyOxv0XXfQ29fXC5z96iauXXUpH9rOPDqy45m6OGPu5tMtXUVGxPqYlxRjT6gI8izXdRC6jPWWqgaGtXKOP83k4sBH4Wlv1GmMoLS01yVJVVZX0uUEQdvmMCb+MKl/qhF3GsMo3fboxYMyZZ75vevc25rnn7P7KSmNEjDnzTGOmMdfUS4ExYOrbFZg5XeY2lksnwDoTQ6e2adIxxnzDGDMwyvJYvE8cY8y7zudOYBlwarznKkpOox42WYfXY+fFF3sybVqTrX7ePLj4YnjmGeg0vJx2nay5rV3HDoy4uTwjrphe0j7xSkQOAdoZY/Y662cBs9Ndr6KEHvWwyToiPXYOPXQLc+eWMGKE9dEfPx6eeMI+DO68u4xz5q+yyU7KyyktK6M0w/KnNGgrImNFZAc2dcufROQpZ/9RIrLSKVYM/FVENgIvAn8yxjyZSr2KkhMk62GjbwX+E2ebRk6iGjx4N9OmwaOPWmX/0EPWAcvNVTt8VhlVXw5P9NKUevjGmGVYE03k/veAkc76G8CgVOpRlJwkmSxI+lbgPwm0aeRkqQ0bejBvHjz+uH0Y3HqrfW4MHtx8Rm0YPHRAY+koSuZIZhKWpjH0nxTa9JVXujWbkAVW2btKPqhctfGiCl9RMkmiuW5Dkhs1p0ihTS+88B3Ky49vti9sSt6LKnxFCRNr1rTe49c0hv6TR22qCl9RwkJrtuTIB0EOK6WMkCdtqgpfUcJCpC35gQfsvl694NprdaBWSRlV+IoSFry25IIC+N//hYMHbUjGhga76ECtkgKq8BUlLHhtyW+/Db/+te3tt2tnHwAiOlCbDbQ1DpNBVOErSphwbclr1sD//V+TGee222DXrlAqEcVDyOdJqMJXlDCSR54jOUXI50mowleUsJInniOhIxWTTMjnSajCVxRFcUnVJBPyNzNV+IqiKC5+mGRC/GaWak5bRVGU3CHHU0ZqD19RFMUl5CaZVFGFryiK4iUZk0yIfe+9pKTwReQW4FxgP/A68B1jzO4o5YYDtwMFwG+MMTelUq+iKEpoCLnvvZdUbfjPAAONMScD/wSmRRYQkQLgl8AIoD9woYj0T7FeRVGUcJBs5rIMkJLCN8Y8bYw56Gz+HegbpdipwDZjzBvGmP3AEmB0KvUqiqIEQjypD7NooNdPG/7lwMNR9vcB3vFs7wBO87FeRVEU/4nXVJNFA71tKnwReRY4IsqhnxljHnPK/Aw4CCxKVSARmQBMACguLqY6ydej2trapM8NgrDLB+GXUeVLnbDLmEn5jlm0iGP37UMaGmjYt4/tCxfy9r59ABRt3kyPmhoKTziBRunKymDfvlCbdDDGpLQAlwFrgC4xjpcBT3m2pwHT4rl2aWmpSZaqqqqkzw2CsMtnTPhlVPlSJ+wyZlS+1auN6dzZmIIC+3nvvcbMnWs/nf0HO3a05UIEsM7E0KmpeukMB6YCXzfG/CdGsbXAF0TkWOBd4ALgolTqVRRFSTteU403CY0nP4EYE7oAaa2RqpfOXUA34BkRqRGRewBE5CgRWQlg7KDuZOApYCuw1BizOcV6FUVR0k9ZGUybZkNTu544DQ12gLagANO+fagHaSNJqYdvjPl8jP3vASM92yuBlanUpSiKkjEio2A6+Qk2FhUxJEt696AzbRVFUdomhifOnjAP0EZBFb6iKEo8hDgKZrxotExFUZQ8QRW+oihKnqAKX1EUJU9Qha8oipInqMJXFEXJE1ThK4qi5AliQy+EExH5CHgrydN7Ax/7KI7fhF0+CL+MKl/qhF1GlS9xPmeMOSzagVAr/FQQkXXGmKGZliMWYZcPwi+jypc6YZdR5fMXNekoiqLkCarwFUVR8oRcVvj3ZVqANgi7fBB+GVW+1Am7jCqfj+SsDV9RFEVpTi738BVFURQPqvAVRVHyhKxX+CIyXEReFZFtIvKTKMc7isjDzvEXRKRfgLIdLSJVIrJFRDaLyA+ilCkXkU+cjGE1IjIjKPmc+reLyD+cutdFOS4icofTfi+LyJCA5TvB0zY1IrJHRK6NKBNoG4rIQhHZKSKbPPt6isgzIvKa83lojHMvdcq8JiKXBizjLSLyinMfl4lIjxjntvqbSKN8N4jIu577ODLGua3+59Mo38Me2baLSE2Mc9PefkkTK9ltNixAAfA6cBzQAdgI9I8o8z3gHmf9AuDhAOU7EhjirHcD/hlFvnLgjxlsw+1A71aOjwSeAAT4MvBChu/3B9iJJRlrQ+BrwBBgk2fffOAnzvpPgJujnNcTeMP5PNRZPzRAGc8CCp31m6PJGM9vIo3y3QBMieM30Op/Pl3yRRyvBGZkqv2SXbK9h38qsM0Y84YxZj+wBBgdUWY08H/O+h+AYSIiQQhnjHnfGPOSs74Xm9O3TxB1+8ho4AFj+TvQQ0SOzJAsw4DXjTHJzr72BWPM88C/InZ7f2f/B4yJcurZwDPGmH8ZY/4NPAMMD0pGY8zTxuaYBvg70DcddcdDjDaMh3j+8ynTmnyO/hgHLPa73nST7Qq/D/COZ3sHLRVqYxnnx/4J0CsQ6Tw4pqTBwAtRDpeJyEYReUJEBgQrGQZ4WkTWi8iEKMfjaeOguIDYf7JMtiFAsTHmfWf9A6A4SpkwteXl2De3aLT1m0gnkx2T08IYZrEwtOEZwIfGmNdiHM9k+7VKtiv8rEBEugKPANcaY/ZEHH4Ja6IYBNwJLA9YvK8aY4YAI4CrReRrAdcfFyLSARgF/D7K4Uy3YTOMfa8Prb+ziPwMOAgsilEkU7+Ju4HjgRLgfazZJIxcSOu9+9D+p7Jd4b8LHO3Z7uvsi1pGRAqB7sCuQKSzdbbHKvtFxphHI48bY/YYY2qd9ZVAexHpHZR8xph3nc+dwDLsK7OXeNo4CEYALxljPow8kOk2dPjQNXU5nzujlMl4W4rIZcB/A992HkwtiOM3kRaMMR8aY+qNMQ3Ar2PUm9E2dHTIecDDscpkqv3iIdsV/lrgCyJyrNMDvABYEVFmBeB6Q5wPPBfrh+43jq3vt8BWY8yCGGWOcMcURORU7D0J5IEkIoeISDd3HTuotymi2ArgEsdb58vAJx7TRZDE7FVlsg09eH9nlwKPRSnzFHCWiBzqmCvOcvYFgogMB6YCo4wx/4lRJp7fRLrk844NjY1Rbzz/+XTyDeAVY8yOaAcz2X5xkelR41QXrBfJP7Ej9z9z9s3G/qgBOmHNANuAF4HjApTtq9hX+5eBGmcZCUwEJjplJgObsd4GfwdOD1C+45x6NzoyuO3nlU+AXzrt+w9gaAbu8SFYBd7dsy9jbYh98LwPHMDakK/AjgutAl4DngV6OmWHAr/xnHu581vcBnwnYBm3Ye3f7m/R9V47CljZ2m8iIPkedH5jL2OV+JGR8jnbLf7zQcjn7L/f/d15ygbefskuGlpBURQlT8h2k46iKIoSJ6rwFUVR8gRV+IqiKHmCKnxFUZQ8QRW+oihKnqAKX1EUJU9Qha8oipIn/H80h76v9VW/6AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "h = .2\n", "xdist = 6*np.pi\n", "x = np.linspace(0,xdist, num=int(xdist/h)+1)\n", "y=exp(x)\n", "\n", "methods = ('RK45', 'Radau', 'DOP853')\n", "colours = ('bx','r.','k-','r.','b.')\n", "\n", "\n", "for c, meth in zip(colours, methods):\n", " y = solve_ivp(m_y, (0,6*np.pi), (-1,0), method=meth, t_eval = x, max_step=h).y[1]\n", " y2 = sin(x)-y\n", " plt.plot(x,y2,c, linewidth=1, label=meth)\n", " #print(y)\n", "\n", "\n", "plt.title(\"Comparison of ODE errors on sin(x)\")\n", "plt.legend()\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "2965340d", "metadata": {}, "source": [ "So in conclusion, with the above methods, the errors are very small but they do eventually get amplified." ] }, { "cell_type": "code", "execution_count": null, "id": "c304ec35", "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 }