{ "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": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEICAYAAAB25L6yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAy4UlEQVR4nO3deViVZfrA8e8tIoqa5hKaClimlZopbZoWqGU5ZU22Spll0T7VZPOraEybaGaapsWpdGxPSXJqrKwsS6EmsxKXXDM3QHLJFUVUFJ7fH88LHPAA58DZgPtzXefinHe9z8vrzePzPosYY1BKKRW6GgU7AKWUUlXTRK2UUiFOE7VSSoU4TdRKKRXiNFErpVSI00StlFIhThO18jkRSRSRucGOo4SINBOR2SKSJyL/CXY8gSIib4nIUz461gQRme6LYynvaaIOYSIySkQyRSRfRLaKyBwRGRjsuKpjjEk1xlwc7DhcXA1EAW2NMde420BETheRj51kvl9E0kVkgMv6WBExzu8iX0S2i8gnInJRheNkichBl+3yReQl/349EJExIvKtv8+jgkMTdYgSkT8CLwBPY5NMNPAKcEUQw6qWiDQOdgxuxAC/GGOOulspIicDC4AVQFfgRGAWMFdE+lfYvLUxpgXQB/gSmCUiYypsc7kxpoXL614ffhfVEBlj9BViL6AVkA9cU8U2EdhEvsV5vQBEOOvigVzgT8BvwFbgSmA48AuwG3jM5VgTgPeB94D9wBKgj8v6R4ANzrrVwO9d1o3BJrnnneM+5Sz71lkvzrrfgDxgOdDL5Xu+A+wAsoHHgUYux/0WeBbYA2wCLq3iepwGZAB7gVXACGf5RKAQOOJc07Fu9p0GfOZm+WTgG+d9LGCAxhW2GQdsd4k7Cxjq4e95AvAfYLpzbVcA3YFHneu1Gbi4wn3xuvP7/NW51mHOdz8EFDnfca+z/VvAy8CnzvF/AE52Od4AYJHze1kEDHBZ1xX42tnvS+AlYLqzrqkT8y7nei8CooL976Y+v4IegL7c/FLgEuBoxaRQYZsnge+BE4D2wHfAX5x18c7+44Fw4HYnGb4LtAR6Ov+wT3K2n+Aksqud7cc5iTHcWX8NtpTZCLgOOAB0dNaNcc51H9AYaEb5RD0MWAy0xibt01z2fQf4yIkpFvtHZKzLcY84sYcBd2H/IImbaxEOrAceA5oAg50E08Pl+02v4lpuA25xszzBSX6RVJ6oT3KWn+Z8zsK7RH3IuUaNneuxCUh2+b1tctn+Q+DfQHPn9/4jcIfL9fq2wvHfwv7xPMc5fiqQ5qxrg/0DeJOz7gbnc1tn/ULgOWyB4ALnepYk6juA2c51CQPigOOC/e+mPr+CHoC+3PxSIBHYVs02G4DhLp+HAVnO+3jgIBDmfG7pJJNzXbZfDFzpvJ8AfO+yrhG21DaoknMvA65w3o8BciqsL00aTtL8BTgPp9TpLA8DDgOnuyy7A8hwOcZ6l3WRznfo4CaeQdhk63r8GcAEl+9XVaI+ClziZvmpzjk7UXmibuosP9/5nIVTqnV53V7JeScAX7p8vtzZt+LvrTW2+usw0Mxl+xuA9IrX3GX9W8BrLp+HAz87728Cfqyw/ULnONHONWnusu5dyhL1rdiCwRnB/rfSUF6hWJ+o7H8p24lIY1NJvSq2hJvt8jnbWVZ6DGNMkfP+oPNzu8v6g0ALl8+bS94YY4pFJLfkeCIyGvgjNlnh7NfO3b4VGWPmOw/TXgaiRWQWtsTeDFv6rfgdOrl83uZynAIRKTl3RScCm40xxVUcqyo7gY5ulncEirElzRMq2bfkHLtdll1pjPnKw3NX/J3sdPN7a4H9juHAVuc6gP2DWum1d2xzeV9A2fWreP9A2TU7EdhjjDlQYV0X5/00532aiLTGVoMkG2OOVBOLqiF9mBiaFmL/S3xlFdtswT4kKxHtLKupkn+EiEgjoDOwRURigFeBe7H/LW4NrMRWY5SocghGY8wkY0wctsqlO/AwNjkecfMdfq1B7FuALk7cNTnWV9jqnYquBRYaYwqq2Pf32PrktR6eq6Y2Y0vU7YwxrZ3XccaYns56b4fBrHj/QNk12wocLyLNK6yzJzLmiDFmojHmdGw992XAaC/Pr7ygiToEGWPysPXLL4vIlSISKSLhInKpiDzjbDYDeFxE2otIO2f72rRzjRORq5xWGw9gk8L32PpQg63jRkRuAXp5elAROVtEzhWRcGzd9iGgyCk1zgRSRKSl8wfhjzX8Dj84x/6Tc53isdUIaR7uPxEYICIpItLGiec+bPL5v0q+V5SI3As8ATxaoTTvc8aYrcBc4J8icpyINBKRk0XkQmeT7UBnEWni4SE/A7o7TUAbi8h1wOnAJ8aYbCATmCgiTZwmoZeX7CgiCSLSW0TCgH3YP7hFbs6hfEQTdYgyxjyHTVyPY5PkZmyp9kNnk6ew/5iWY1sLLHGW1dRH2AeFJQ+YrnJKTquBf2JL+duB3thWHp46Dlsi34P97/MubEsOsA8gDwAbsS083gXe8DZwY0whMAK4FFtSfwUYbYz52cP91wEDsU3usrAlypHAMGNMxe+6V0QOYK/5cGzLnIoxz67QjnqWt9+pEqOx1UWrsdfzfcqqbOZjW7tsE5Gd1R3IGLMLWxJ+CPs7+RNwmTGmZN9RwLnYKp0nsA86S3Rwzr0PWINtHaKdYfxInIcDqgETkQlAN2PMjcGORSl1LC1RK6VUiNNErZRSIU6rPpRSKsRpiVoppUKcXzq8tGvXzsTGxtZo3wMHDtC8efPqNwwwjcs7Gpd3NC7v1Me4Fi9evNMY097tSn90d4yLizM1lZ6eXuN9/Unj8o7G5R2Nyzv1MS4g01SSU7XqQymlQpwmaqWUCnGaqJVSKsQFbPS8I0eOkJuby6FDh6rcrlWrVqxZsyZAUXnOV3E1bdqUzp07Ex4e7oOolFINQcASdW5uLi1btiQ2NhaXYRqPsX//flq2bBmosDzmi7iMMezatYvc3Fy6du3qo8iUUvVdwKo+Dh06RNu2batM0vWdiNC2bdtq/1ehlFKuAlpH3ZCTdAm9BkrVTwv/9z9Sp01j4cKFPj+2zvCilFK1tHDhQoYMHkzh0aOkvvsu8+bPp3//ihPY11yDa/Uxa9YsRISff656qOIXXniBgoKqJvao2ltvvcW9995b4/2VUnVHxtSpFB49ShFQWFhIRkaGT4/f4BL1jBkzGDhwIGlpVU/+UdtErZRqILZuJf7DD2kChInQJCKC+Ph4n56iQSXq/Px8FixYwOuvv16aqIuKihg3bhy9e/fmjDPO4F//+heTJk1iy5YtJCQkkJCQAEDHjmVzn77//vuMGTMGgNmzZ3PuuefSt29fhg4dyvbt2485r1KqnioqghtvZGOXvbS6qQlFgw2tklqxscVGn54mOHXUVTxQq1UDuGqGbP3www+55JJL6N69O23atGHJkiX88MMPbNq0iaVLl9K4cWN2795NmzZteO6550hPT6ddu3ZVHnPgwIF8//33iAivvfYazzzzDP/85z9r8y2UUnXF3/9O6o75JI2AgvBCOBm2sY2k2UkAJPZO9MlpGtTDxBkzZvDAAw8AcP311zNjxgw2btzInXfeSePG9lK0adPGq2Pm5uZy3XXXsXXrVgoLC7V9tFINxYIFMH48yfdBQYX+awVHCkiel1zHE3UVJV9/dXjZtWsX8+fPZ+XKlYgIRUVFiAhxcXEeNZlz3ca1HfR9993HH//4R0aMGEFGRgYTJkzweexKqRCzezeMGgVFReS0cr9JTl6Oz07XYOqo33//fUaPHk12djZZWVls3ryZrl270q9fP6ZMmcLRo0cB2L17NwAtW7Zk//79pfu3b9+eNWvWUFxczKxZZZNK5+Xl0alTJwDefvvtAH4jpVRQGAO33QY5OXDuuUS3ina7WWXLa6LBJOoZM2bw+9//vtyykSNHsmXLFqKjoznjjDPo06cP7777LgBJSUlceumlpQ8TJ06cyGWXXcbgwYPLPVicMGEC11xzDYMGDaq2PlspVQ9MngyzZsFxx8GMGaQMfZrI8Mhym0SGR5IyJMV356xsoOravNxNHLB69WqPBs/et29fTcbc9jtfxuXptfBEfRxA3Z80Lu9oXBUsW2ZMRIQxYMzMmaWLpy+fbmKejzEyQUzM8zFm+vLpXh+aKiYOqLaOWkR6AO+5LDoJGG+MecF3fy6UUirEHTgA110Hhw9DUhJcc03pqsTeiST2TiQjI8PnbajBg4eJxpi1wJkAIhIG/ArMqmofpZSqd+67D9auhZ494fnnA3pqb+uohwAbjDHZ/ghGKaVCUmoqvPkmNGsG770HkZHV7+NDYqrpJFJuY5E3gCXGmJfcrEsCkgCioqLiKnbRbtWqFd26dav2HEVFRYSFhXkcU6D4Mq7169eTl5fnk2Pl5+fTokULnxzLlzQu72hc3glkXM1+/ZW422+n8cGDrH3oIbZedplf4kpISFhsjDnL7crKKq8rvoAmwE4gqrpt9WFi1fRhYvBoXN5p8HEdOmRMv3724eG11xpTXOy3uPDRLOSXYkvTOpiFUqphePRRWLIEYmNh6tQqh7/wJ28S9Q3ADH8FEghhYWGceeaZpa+//e1vVW6vQ5Uq1YB98ol9aNi4MaSlQatKuiAGgEddyEUkErgIuMO/4ZRJTYXkZNv5JzoaUlIgsZbd5ps1a8ayZct8Ep87R48eLR0zRClVh/36KzgjZPL003DuuUENx6MStTGmwBjT1hjjmydg1UhNtc0Us7Ntb83sbPs5NdU/54uNjWXnzp0AZGZmum0HuXPnTkaOHMnZZ5/N2WefzYIFCwDbMzEpKYmLL76Y0aNH+ydApVTgFBXZUuGuXTBsGDz0ULAjCs3R85KToeKY/QUFdnltStUHDx7kzDPPLP386KOPct1113m075/+9CcefPBBBg4cSE5ODsOGDWPNmjUALF68mG+//ZZmzZrVPDilVGh46in4+mvo0AHeeQcaBX+kjZBM1DmVDDpV2XJP1abqIyMjg3Xr1pV+3rdvX+mgTSNGjNAkrVR98PXX8OST9qHhtGlwwgnBjggI0UQdHW2rO9wt94fGjRtTXFwMlB/C1FVxcTELFy50m5CbN2/un8CUUoGzc6f9L3txMTz2GAwdGuyISgW/TO9GSsqxHX8iI+1yf4iNjWXx4sUAfPDBB263GTx4MC+9VNbPx58PJZVSAWYM3HKLfYg4YACE2LjyIZmoExNtk8WYGPs/kJgY+7m2rT5K6qhLXo888ggATzzxBPfffz+DBg2qtPfhP/7xDzIzMznjjDM4/fTTmTJlSu2CUUqFjkmTbHO81q3h3XchPLzaXQIpJKs+wCbl2ibmioqKitwuHzRoEL/88ssxy8eMGVM6iW3btm157733jtlGZ3RRqo5bvBgefti+f/11WzIMMSFZolZKqYDYvx+uvx6OHIG774arrgp2RG5polZKNUzGwF13wfr1cMYZ8M9/BjuiSmmiVko1TO+8Y3vRRUbaLuJNmwY7okppolZKNTxr19qqDoCXXoLTTgtuPNXQRK2UalgOHbJTahUUwKhRZWN6hDBN1EqphuXhh+Gnn+Dkk+2M4kEautQbDSpRV5x5QYcxVaqB+fBDW9URHm7rpY87LtgReSRkE3XqilRiX4il0cRGxL4QS+oKPw2dp5RqGHJy4NZb7fu//x3Ocj/rVSgKyUSduiKVpNlJZOdlYzBk52WTNDvJr8l6x44dlQ5j+uyzz5Zu16tXL7KyssjKyuK0007j9ttvp2fPnlx88cUcPHjQb/EppWrh6FFbH71nD/zud/DAA8GOyCshmaiT5yVTcKT8OKcFRwpInpdcq+NW7EI+fvz40nX3338/Dz74IIsWLeKDDz7gtttuq/Z469at45577mHVqlW0bt260nFClFJBNnEiLFgAJ55oZxOvA/XSrkKyC3lOnvvxTCtb7qmKw5y+9dZbZGZmAvDVV1+xevXq0nWuw5hWpmvXrqXjW8fFxZGVlVWr+JRSfjB/vh3RTcS2m27f3uenKJuR6kKfzUjlKiQTdXSraLLzjh3nNLqVn8Y5pfJhTF2HQIXyw6BGRESUvg8LC9OqD6VCzW+/2YxpDIwfD25mb6qtkhmp7GQnUjojFfguWYdk1UfKkBQiw8uPcxoZHknKED+NcwpcfPHFbocxjY2NZcmSJQAsWbKETZs2+S0GpZQPFRfbNtLbtsGgQfDnP/vlNFXNSOUrIZmoE3snMvXyqcS0ikEQYlrFMPXyqST29vFwei4mTZrkdhjTkSNHsnv3bs4//3wmT55M9+7d/RaDUsqHnn8e5syBNm1ssddPE0/7a0YqVyFZ9QE2Wfs6Mefn55f77DqMabt27dwOY9qsWTPmzp3L/v37admyZbl1K1euLH0/btw4n8aqlKqFRYvAGW+eN9+ELl38dqpAzEgVkiVqpZSqsbw8O3Tp0aPwhz/AiBF+PV0gZqTyKFGLSGsReV9EfhaRNSLS33chKKWUb6QuT6XjH6KQ2I10TAwj9eZ+fj9n+RmpjM9mpHLladXHi8DnxpirRaQJEFndDkopFUipK1IZO+lmDr9bBEWwLayIsS3ugHsa+/X5FpTNSJWR8TXxfmhZUm2JWkSOAy4AXgcwxhQaY/b6PBKllKqF5I/v53CWTdIYoAgOrz9c645yoUCMMVVvIHImMBVYDfQBFgP3G2MOVNguCUgCiIqKiktLSyt3nFatWtGtW7dqAyoqKqp0gtlg8mVc69evJy8vzyfHys/PP2awqVCgcXlH4/JOxbiO//FH+hb8HyYXeBubrMOAm0G6CPMvnB+UuLyRkJCw2BjjfgASY0yVL+As4ChwrvP5ReAvVe0TFxdnKlq9evUxy9zZt2+fR9sFmi/j8vRaeCI9Pd1nx/Iljcs7Gpd3ysW1aJExzZubmAcwTMAwFsMQ5+cETMzzMcGJy0tApqkkp3ryMDEXyDXG/OB8fh/wfw29H4SFhXHmmWfSq1cvLr/8cvbu3QtAVlYWvXr1Kt3u1VdfpV+/fuzZs6d02aRJkxARdu7cWbpPs2bNSscNufPOOwP6XZRS2PkOhw+HAwdIOXS+7SjXBRgEdPF/R7lAqTZRG2O2AZtFpIezaAi2GqTOKRnrY+XKlbRp04aXX375mG2mTZvGv/71L+bOncvxxx8PwObNm5k/fz7RFRpGnnzyySxbtoxly5aVdpBRSgXI9u1wySWwYwdcfDGJL84PeEe5QPG01cd9QKrT4mMjcIv/QiqzcOFCMjIyiI+Pp39/37YI7N+/P8uXLy+3bObMmfztb39j3rx5tGvXrnT5gw8+yF/+8hdGjRrl0xiUUjUTdvCgHa50wwaIi4P334cmTfzSUS4UeJSojTHLsHXVAbNw4UKGDBlCYWEhTZo0Yd68eT5L1kVFRcybN4+xY8eWLsvOzubee+9l6dKldOjQoXT5xx9/TKdOnejdu/cxx9m0aRN9+/bluOOO46mnnmLQoEE+iU8pVYXCQnqOHw+LF9vptD79FCr0Gq5vQrZnYkZGBoWFhRQVFVFYWEhGRkatj1kyHnXbtm3ZvXs3F110Uem69u3bEx0dzcyZM0uXFRQUkJKSwpNPPnnMsTp27EhOTg5Lly7lueeeY9SoUezbt6/WMSqlqlBcDGPH0iYz0w5X+vnnEBUV7Kj8LmQTdXx8PE2aNCEsLIwmTZr4pBF5SR11dnY2hYWF5eqoIyMjmTNnDlOmTCE11c4ks2HDBjZt2kSfPn3o1asXubm59OvXj23bthEREUHbtm0BOxb1ySefzC+//FLrGJVSVXj0UZg+naKmTeGzz8CDJr/1QcgOytS/f3/mzZvnlzrqVq1aMWnSJK644gruuuuu0uXt27fn888/Jz4+nnbt2jFs2DB+++03APbv30/v3r3JzMykXbt27NixgzZt2hAWFsbGjRtZt24dJ510ks9iVEpV8OKL8Mwz0LgxKydOpE8dmvOwtkI2UYNN1r5+iFiib9++9OnTh7S0tHJ1y127duXjjz9m+PDh/Pe//+Xcc891u/8333zD+PHjady4MWFhYUyZMoU2bdr4JValGrz33oMHH7Tv33iDPX4cDS8UhXSi9rWKw5zOnj279L3rkKV9+vTh119/PWZ/16m2Ro4cyciRI30fpFKqvPR0GD3aztLy97/DTTeBD55Z1SUhW0etlFL89BNceSUUFtohSx9+ONgRBYUmaqVUaMrKgksvhX374Npr7YwtdWz2cF8JaKI21QwA1RDoNVDKA7t22V6HW7faCWnfeQcaNdxyZcC+edOmTdm1a1eDTlTGGHbt2kXTpk2DHYpSoaugAC67DNauhTPOgA8/hIiIYEcVVAF7mNi5c2dyc3PZsWNHldsdOnQoJBOZr+Jq2rQpnTt39kFEStVDR4/CddfB99/bSQfnzIFWrYIdVdAFLFGHh4fTtWvXarfLyMigb9++AYjIO6Eal1L1hjFw553wySd25vAvvoATTwx2VCGh4Vb6KKVCy4QJ8Prr0KyZTdannhrsiEKGJmqlVPBNmQJPPmkfGL73Hvipo1tdpYlaKRVcH34I99xj3//733D55UENJxRpolZKBc+338INN9hR8SZOhNtuC3ZEIUkTtVIqOFatsqXnQ4fgjjvgz38OdkQhSxO1UirwcnNth5a9e+GKK+Dllxtsr0NPaKJWSgXWnj22a3huLgwYADNmQFiYTw6dmgqxsfaZZGys/VwfNKjR85RSQXbokB1kaeVKOO00mD3bNsfzgdRUSEqyHRsBsrPtZ4DEOj6NopaolVKBUVQEN94I33wDnTrZabR8OIZ7cnJZki5RUGCX13WaqJVS/mcM3H8/fPCB7RI+Z47tIu5DOTneLa9LNFErpfzvr3+1DwwjIuCjj6B3b5+forK87+O/B0HhUaIWkSwRWSEiy0Qk099BKaXqkbfesvUPIrYi+cIL/XKalBSIjCy/LDLSLq/rvHmYmGCM2em3SJRS9ULqilSS5yWTk5dDdJN2pKTtJBFg0iTw4/R1JQ8Mk5NtdUd0tE3Sdf1BImirD6WUD6WuSCVpdhIFR+xTvezCHSRdBowYQeK99/r9/ImJ9SMxVySeDOQvIpuAPYAB/m2MmepmmyQgCSAqKiouLS2tRgHl5+fTokWLGu3rTxqXdzQu79SXuK7//nq2H94Om4EsIBboAlERJ5B23ntBiytQahNXQkLCYmPMWW5XGmOqfQEnOj9PAH4CLqhq+7i4OFNT6enpNd7XnzQu72hc3qkvcckEMYzF0BiDOD/HYmSCBDWuQKlNXECmqSSnevQw0Rizxfn5GzALOKdGfzKUUvVadMQJtiRdhP3/dxGQBdGt6kHTiyCqNlGLSHMRaVnyHrgYWOnvwJRSdcxHH5EycycRnYAwQOzPiG4RpAypB00vgsiTh4lRwCyxA6Y0Bt41xnzu16iUUnXLG2/A7beTWFwMFwxm3D2r2LZyOx16deDZsc+S2LsePuELoGoTtTFmI9AnALEopeoaY+CZZ+CRR+zn8eNJnDCBRB0Jz6e0eZ5SqmaKi2HcOHj+eduZZdIkCEATvIZIE7VSyntHjsCtt8L06RAeDtOmwXXXBTuqeksTtVLKOwcOwDXX2IGVmjeHWbPgoouCHVW9polaKeW53bvhd7+D77+Hdu3gs8/g7LODHVW9p4laKeWZ3FwYNgxWr7YDacydCz16BDuqBkGHOVVKVe/nn+20WatXQ8+e8N13mqQDSBO1UqpqP/wAAwfC5s02WZfM0KICRhO1UqpSx//4IwweDLt22brpL7/06fRZyjOaqJVS7s2YQe/HHrMTD44ebVt3VByZXwWEJmql1LEmTYJRo2hUVGQ7tbz5pm0vrYJCE7VSqowx8PjjdiJaYMMdd8A//gGNNFUEkzbPU0pZRUVw113w6qsQFgavvcbm2FhODnZcSkvUSing0CHb2/DVV6FpU1sfPWZMsKNSDk3USjV0eXlwySU2ObdubVt2XH55jQ+Xmgqxsba2JDbWfla1o1UfSjVk27bBpZfCsmXQsSN88QX07l3jw6WmQlKSbSgCkJ1tP0P9nHQ2ULRErVRDtWEDnH++TdKnnGJ7G9YiSQMkJ5cl6RIFBXa5qjlN1Eo1RMuW2SS9cSPExcG339p6ilrKyfFuufKMJmqlGpqvv4YLL4Tt22HIEEhPhxNO8MmhoyuZw7ay5cozmqiVakhmzbIj4O3bZ1t5fPoptGzps8OnpBzbeTEy0i5XNaeJWql6KnVFKrEvxNJoYiNiX4gl9cXb4Oqr4fBh2156xgyIiPDpORMTYepUiImxs3PFxNjP+iCxdrTVh1L1UOqKVJJmJ1FwxD7Zy87LJqnwdegJiSMnwPjxNpP6QWKiJmZf0xK1UvVQ8rxkm6Q3A/8DNkNBE0i+5nh44gm/JWnlHx6XqEUkDMgEfjXGXOa/kJRStZWTl2OT9NtAERAG3Aw5XfYGNS5VM96UqO8H1vgrEKWU70RHdoQsbJI2zs8siG6lzS/qIo8StYh0Bn4HvObfcJRStTZtGinv7SSiE7YkLfZnRLcIUoZo84u6SIwx1W8k8j7wV6AlMM5d1YeIJAFJAFFRUXFpaWk1Cig/P58WLVrUaF9/0ri8o3F5xxdxNTp0iFMmTaLjnDkAvHL96Tzeegt71u3l+FOO5+74uxkaNTTgcflDfYwrISFhsTHmLLcrjTFVvoDLgFec9/HAJ9XtExcXZ2oqPT29xvv6k8blHY3LO7WOa9UqY3r2NAaMadrUmFdfNaa4OPhx+Ul9jAvINJXkVE8eJp4PjBCR4UBT4DgRmW6MubFGfzaUUr719ttw9912UI1TT4WZM2s9ZocKLdXWURtjHjXGdDbGxALXA/M1SSsVAg4csGNGjxljk/RNN8GiRZqk6yHt8KJUXbRyJVx7LaxZA82awcsv24St7aPrJa8StTEmA8jwSyRKqeoZYyeavfdeOHgQTjsN/vMf6Nkz2JEpP9KeiUrVFfn5MHo0jB1rk/SYMbaqQ5N0vadVH0rVBStW2NHu1q61w9G98grcfHOwo1IBoiVqpUKZMXbC2XPOsUm6Z09bitYk3aBoolYqVO3fDzfeaCcdPHTIVnn8+COcfnqwI1MBpolaqVD0009w1lnw7rvQvDlMmwavvXbsqPxV0NnA6w+to1YqlBgD//43PPCAHeC/d2/bgeXUU706jM4GXr9oiVqpULFvH9xwg5195fBhuP12+OEHr5M06Gzg9Y2WqJUKAS3WrbOJef16aNHCzl91ww01Pp7OBl6/aKJWKpiMgcmT6ffAA3DkCPTpY6s6unev1WGjo211h7vlqu7Rqg+lgiUvz3YDv+ceGh05AnfeCd9/X+skDTobeH2jiVopPzpmJvAVTtOLxYuhXz94/31o2ZJV48fD5MnQtKlPzquzgdcvWvWhlJ+4nQl8dhJ88QWJye9BYSH07QszZ7IjN9fn59fZwOsPLVEr5SelM4G7KDhSQPKv02ySvuce+O476NYtSBGqukJL1Er5SU6e08RiM3ai2VigC+S0wo54d/XVwQpN1TFaolbKT6JbRdsk/TYw3/m5GaJbnKhJWnlFE7VSfpISNYrwjUARYOzP8JxwUoY9E+TIVF2jiVopX9u6FUaNInHUX3ksGyQMEJDGwmOjHyOxtz7hU97ROmqlfKWoyI4T/fjjtjt406ZMGPtnhg0cSMaCBcTHx9O/f/9gR6nqIE3USvnCjz/aDitLl9rPl10GkyZB1670B/pfcEFQw1N1m1Z9KFUbe/bYQZTOO88m6S5dYNYs+Phj6No12NGpekITtVI1YQy88w706AFTpkBYGPzpT3ZW8Cuv1NnAlU9p1YdS3lq1Cu6+G775xn6+4AJbN62TzCo/qbZELSJNReRHEflJRFaJyMRABKZUyDlwAB55BM480ybp9u3h7bchI6PSJK2zrChf8KREfRgYbIzJF5Fw4FsRmWOM+d7PsSkVOj76CP7wBzugswjccQc8/TS0aVPpLjrLivKVakvUxsp3PoY7L+PXqJQKFVlZMGKErXfOybGDKC1caOulq0jSoLOsKN8RY6rPuSISBiwGugEvG2P+z802SUASQFRUVFxaWlqNAsrPz6dFixY12tefNC7v1PW45MgRusycScy0aYQdPszRyEg23XorW668EhMW5tG5Bg++EGOOfagoYpg//+saxRVoGpd3ahNXQkLCYmPMWW5XGmM8fgGtgXSgV1XbxcXFmZpKT0+v8b7+pHF5p07HNX++Maeeaoxt22HM9dcbs2WL1+eKiSk7hOsrJqaGcQWBxuWd2sQFZJpKcqpXzfOMMXuBDOCSGv3JUCqUbd8ON90EgwfDzz/bmVa+/BJmzICOHb0+nM6yonzFk1Yf7UWktfO+GTAU+NnPcSkVOCVdv3v0gOnT7Swrf/kLLF8OQ4fW+LA6y4ryFU9afXQE3nbqqRsBM40xn/g3LKV8K3VFKsnzksnJyyF6WTQpQ1Ls4EiZmbZnYWam3fDSS+Gll+Ckk3xyXp1lRflCtYnaGLMc6BuAWJTyC7dTYn18O7z1NonPf2Wrjjt3hhdfhN//XnsVqpCjXchVved2SqyjB0ku/tL2RHnoIdv1+6qrNEmrkKRdyFW9V+WUWEuXQu/ewQpNKY9oiVrVe9GRHd1PidU6WpO0qhM0Uav6a/16GDOGlOlbCd9A+SmxNoeTMuTpIAeolGc0Uav6x0nQnHoqvP02iasa8Vi7M5HGlE2JdZNOiaXqDk3Uqv5Yvx5uuaU0QQMwdiz88gsT/rOUBV9/x21jb2PB1wuYcOOEoIaqlDf0YaKq+zZsgKeegmnTbOeVsDC49VY7+pFLe+j+/ftz+PBhnbdQ1TlaolZ114YNNiH36AFvvWWX3Xor/PILvP56uSRdMi704MEX6rjQqs7RErWqezZssANmvPNOWQn6lltsCfrkk4/ZvPy40KLjQqs6R0vUqu7YuNHWOffoAW++aZfdcgusXQtvvOE2SYOOC63qPi1Rq9C3caMtQb/9tkcl6IpycrxbrlSo0RK1Cl0lJeju3W2JGWyzu59/rrIEXVF0tHfLlQo1mqhV6Nm0CW67zVZxVEzQb74J3bp5dTgdF1rVdZqoVUClrkgl9oVYGk1sROwLsaSucGl+UZKgu3e3rTaKi+Hmm2ucoEuUHxfa6LjQqs7ROmoVMG6HG52dBL/tIDFtlW1id/SoHdHu5pttHfQpp/jk3CXjQmdkfE18fLxPjqlUoGiiVgHjdrjRIwUkf/Igia9hE/To0fD44z5L0ErVB5qoVcBUOdyoJmilKqWJWgWGMURHnED2+u12mNEiIAy4GaJ7nFg2NodS6hj6MFH51969MGkS9OxJyrvb3Q83eskzwY1RqRCnJWrlFy3WrrUzes+YUdotMLFjR9addQpPLvgf5qjR4UaV8pAmauU7BQXw3nsweTJnLVpUtnzIEDvT94gRnDIznKg5C9m2LYOotvGcIjqSnVLV0UStam/tWpgyxTav27sXgCMtWxJ+221wxx224wqugyP1B/qzbZsOjqSUJ6pN1CLSBXgH6AAUA1ONMS/6OzAV4o4cgY8+gsmTYf78suXnnAN33cXCjh25YNiwcrtUNTiSJmqlKudJifoo8JAxZomItAQWi8iXxpjVfo5NhaLNm+HVV+G112DrVrssMhJGjbLVG/36AVCckXHMrjo4klI1U22iNsZsBbY67/eLyBqgE6CJuqEoLoYvv7Sl59mz7WeA006zyfmmm6B162oPEx0N2dnulyulKifGGM83FokFvgF6GWP2VViXBCQBREVFxaWlpdUooPz8fFq0aFGjff2pIcYVnpdHhzlzOHH2bJpt2QJAcePG7Bw0iF9HjCCvTx8Q8Tiur746gWef7cHhw2GlyyIiihg3bi1Dh/7ml+/gSVyhQOPyTn2MKyEhYbEx5iy3K40xHr2AFsBi4Krqto2LizM1lZ6eXuN9/ak+xTV9+XQT83yMkQliYp6PMdOXTy9bWVxszIIFxtx4ozEREcaAfUVHG5OSYsy2bbWKa/p0Y2JijBGxP6dPd7uZ39Sn32MgaFzeqU1cQKapJKd61OpDRMKBD4BUY8x/a/TnQoWESgdGOniIxCVHbPXG8uV2YxEYPtxWb1x6qR2wv5ZKBkdSSnnOk1YfArwOrDHGPOf/kJQ/lQ6M5DLeRkGXApLTbifxeacarH17O2B/UhJ07Rq8YJVSgGetPs4HbgJWiMgyZ9ljxpjP/BaV8pucvBybpCuMt5HT2cCgQbb0fNVVEBER3ECVUqU8afXxLeD+iZGqO/Ly4KOPiD7clOysg+XG2yALort3hG++CW6MSim3dFCm+mzfPjvexogRcMIJcPPNDP/kTOjUxJakBfuzcwTDW/4juLEqpSqlXcjrm/37bVvnmTPh88/h8GG7vFEjSEjgs58+Aj6Gq8bBzm3QrgMseZbPMhLhrqBGrpSqhCbq+mD/fvjkE5uc58wpS84icOGFcO21tt65QwdyGgG7E2FF+aYXOVq5pVTI0kRdV+Xnw6ef0vPll2HRIjh0yC4XsQ8Fr70WRo6Ejh3L7aa9A5WqezRR1yUHDsBnn9mS86efwsGDtC9Zd/75Zcm5U6dKD5GSUjKCXdmyyEi7XCkVmjRRh7qCAludMXOmrd5wzbD9+7O+Xz+6PfIIdO7s0eFKOpskJ9vBkKKjbZLWTihKhS5N1EGUuiKV5HnJ5OTlEN0qmpQhKXa2k4MH7YPAmTPtg8EDB8p2Ou88W3K++mro0oXcjAy6eZikS2jvQKXqFk3UQeK2K/essfDyyySmrrB10CXOOacsOcfEBClipVSwaKIOEvdduQ+T3GwhifnAWWeVJWftxq1Ug6aJOpCKi2HZMvjiC7IPZ0Mux3Tlzu4ssGE9nHRSUENVSoUOTdT+tmWLHXR/7lz7c8cOAMIe6ERR1q/HdOUOOz5ak7RSqhxN1L528CB8+y188YVNzitWlF/fpQsMG0bRvHjoOxbCDpeVqDtFUPRFCmhvbqWUC03UtWUMrFplk/IXX9iBjUo6n4BtpJyQABdfbF89eoAIMbGQvZTyXbmXPkvMPm2OoZQqTxN1TezcWVadMXeurd5w1bcvDBtmE/OAAW6HDLUdTxIpcOnKHRkJKVP9HbxSqq5p8Im60rbMrgoLabVsWVliXrLElqRLdOhgk/KwYTB0qB2prhra8UQp5akGnagrnZbKQGLTs8sSc3o6fV3bNUdEwAUXlFVn9O5d6SSvVdGOJ0opTzToRF3ptFTTbibx2aJy2x6IjaX5VVfZxDxokK2nUEqpAGiYidoYyMkhe29OJW2Zi6BtW7joIludcdFFLFq3jvj4+KCGrZRqmBpGoj58GJYuhe++s6+FC2HLlsrbMreOht822cH2S6xbF5zYlVINXv1M1Fu32mRckpQzM6GwsPw2bdpQNO9v0Pe2Y9syz30antVZypRSoaHuJ+qjR2H58rKk/N13kJV17HY9e0L//ra53IAB0L07MV2F7KWibZmVUiGt7iXqXbtsQi5Jyj/+WH6MZoCWLeHcc21C7t/fvj/++GMOpW2ZlVJ1QbWJWkTeAC4DfjPG9PJXIHdPTmXqxmSKmucQ9mk0SSel8ModN8Dq1WVJ+bvv4Jdfjt25W7eypDxggC09h4VVe05ty6yUqgs8KVG/BbwEvOOvIO6enMrkX5NgTwEshaLYbCbnJsHZr/LKkq/Lb9y0KZx9dlkVxnnnedTBpDLallkpFeqqTdTGmG9EJNafQUzdmGyTdLlmcgVMvWA9r+zoUpaU+/eHPn2gSRN/hqOUUiFFjGtX6Mo2son6k6qqPkQkCUgCiIqKiktLS/M4iISMwfCtgfnYZnICDAYGCunx8z0+jj/l5+fTokWLYIdxDI3LOxqXdzQu79QmroSEhMXGmLPcrjTGVPsCYoGVnmxrjCEuLs54I2xcjGEshsYYxPk5FhM2Lsar4/hTenp6sENwS+PyjsblHY3LO7WJC8g0leTUkGgsnHRSCnSIhJuxJembgQ6RdrlSSjVwIZGoX7krkbs6TSXs+BgYKIQdH8Ndnabyyl36lE8ppTxpnjcDiAfaiUgu8IQx5nVfB/LKXYm8QiIZGRk6poZSSrnwpNXHDYEIRCmllHshUfWhlFKqcpqolVIqxGmiVkqpEKeJWimlQpxHPRO9PqjIDiC7hru3A3b6MBxf0bi8o3F5R+PyTn2MK8YY097dCr8k6toQkUxTWTfKINK4vKNxeUfj8k5Di0urPpRSKsRpolZKqRAXiok6VOdX0bi8o3F5R+PyToOKK+TqqJVSSpUXiiVqpZRSLjRRK6VUiAtYohaRS0RkrYisF5FH3KwXEZnkrF8uIv083dfPcSU68SwXke9EpI/LuiwRWSEiy0QkM8BxxYtInnPuZSIy3tN9/RzXwy4xrRSRIhFp46zz5/V6Q0R+E5GVlawP1v1VXVzBur+qiytY91d1cQXr/uoiIukiskZEVonI/W628d89VtmMAr58YWdB3ACcBDQBfgJOr7DNcGAOdiKu84AfPN3Xz3ENAI533l9aEpfzOQtoF6TrFY+dHs3rff0ZV4XtLwfm+/t6Oce+AOhHJTMRBeP+8jCugN9fHsYV8PvLk7iCeH91BPo571sCvwQyhwWqRH0OsN4Ys9EYUwikAVdU2OYK4B1jfQ+0FpGOHu7rt7iMMd8ZY/Y4H78HOvvo3LWKy0/7+vrYNwAzfHTuKhljvgF2V7FJMO6vauMK0v3lyfWqTFCvVwWBvL+2GmOWOO/3A2uAThU289s9FqhE3QnY7PI5l2O/ZGXbeLKvP+NyNRb7F7OEAeaKyGKxk/v6iqdx9ReRn0Rkjoj09HJff8aFiEQClwAfuCz21/XyRDDuL28F6v7yVKDvL48F8/4SO9l3X+CHCqv8do9VO3GAj4ibZRXbBVa2jSf71pTHxxaRBOw/pIEui883xmwRkROAL0XkZ6dEEIi4lmDHBsgXkeHAh8ApHu7rz7hKXA4sMMa4lo78db08EYz7y2MBvr88EYz7yxtBub9EpAX2j8MDxph9FVe72cUn91igStS5QBeXz52BLR5u48m+/owLETkDeA24whizq2S5MWaL8/M3YBb2vzgBicsYs88Yk++8/wwIF5F2nuzrz7hcXE+F/5b68Xp5Ihj3l0eCcH9VK0j3lzcCfn+JSDg2SacaY/7rZhP/3WP+qHh3UxHfGNgIdKWsMr1nhW1+R/mK+B893dfPcUUD64EBFZY3B1q6vP8OuCSAcXWgrMPSOUCOc+2Cer2c7Vph6xmbB+J6uZwjlsofjgX8/vIwroDfXx7GFfD7y5O4gnV/Od/9HeCFKrbx2z3ms4vrwRcdjn1SugFIdpbdCdzpciFedtavAM6qat8AxvUasAdY5rwyneUnORf8J2BVEOK61znvT9iHUAOq2jdQcTmfxwBpFfbz9/WaAWwFjmBLMGND5P6qLq5g3V/VxRWs+6vKuIJ4fw3EVlcsd/ldDQ/UPaZdyJVSKsRpz0SllApxmqiVUirEaaJWSqkQp4laKaVCnCZqpZQKcZqolVIqxGmiVkqpEPf/6s+mPuvhfMAAAAAASUVORK5CYII=\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": 3, "id": "df8e3774", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEICAYAAACqMQjAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt1UlEQVR4nO3deZgU1b3/8feXGVBQRDYFQRjgZ35sMgPMDxiXZAjRKDe4xoByXXlEfYJxuUZZopJ4CVxi9D5GlpBcLyqIcuPVoEFBtigBFRBCWCSyDDCKbLKoOAwM5/dH1bQ9TfdMz/RSTM/n9Tz9dFfVOXVO1dT0t+ucqlPmnENERCSaekFXQERETl0KEiIiEpOChIiIxKQgISIiMSlIiIhITAoSIiISk4KEZBQzG2pm84OuRzkza2hmb5jZITP7n6Drky5mNt3M/j1J6xprZjOSsS6pPgUJicrMbjKzlWb2lZntMrO3zOySoOtVFefcTOfc5UHXI8yPgXOB5s65G6IlMLOuZjbHDyRfmtliM7sobHmOmTn/b/GVme02szfN7LKI9RSZ2Tdh6b4ys2dTu3lgZreZ2dJUlyPBUJCQk5jZg8B/Ar/G+4JrB0wGrg6wWlUys+yg6xBFe+Cfzrnj0RaaWSfgb8A/gA7AecBrwHwzK4hIfrZz7kwgF3gHeM3MbotIM8g5d2bYa0QSt0XqIuecXnqFXkAT4CvghkrSnIYXRD7zX/8JnOYvKwSKgYeBPcAu4BpgIPBP4AtgdNi6xgJ/Al4BvgQ+AnLDlo8EtvjLNgDXhi27De8L9ml/vf/uz1vqLzd/2R7gELAW6B62nS8Ae4HtwC+AemHrXQo8CRwAtgFXVrI/ugBLgIPAeuAqf/4vgVLgmL9Ph0XJ+yIwN8r8KcC7/uccwAHZEWkeAnaH1bsI+EGcf+exwP8AM/x9+w/gO8Aof3/tBC6POC7+y/97furv6yx/20uAMn8bD/rppwOTgL/46/8A6BS2vouAFf7fZQVwUdiyDsBf/XzvAM8CM/xlp/t13u/v7xXAuUH/32TyK/AK6HVqvYArgOORX0gRaX4FvA+cA7QElgFP+MsK/fyPAfWBO/0v4peAxkA3/0ulo59+rP8l+mM//UP+l3J9f/kNeL+u6wGDga+B1v6y2/yy7gWygYZUDBI/BFYBZ+MFjC5heV8A/uzXKQcvgA0LW+8xv+5ZwD14wdCi7Iv6wGZgNNAA+L7/5fZ/w7ZvRiX78nPg9ijz+/tfvI2IHSQ6+vO7+NNFVC9IlPj7KNvfH9uAMWF/t21h6V8Hfg+c4f/dPwTuCttfSyPWPx0vcPfx1z8TeNlf1gwv+N7sL7vRn27uL18OPIX3Y+S7/v4sDxJ3AW/4+yUL6A2cFfT/TSa/Aq+AXqfWCxgKfF5Fmi3AwLDpHwJF/udC4Bsgy59u7H+R9Q1Lvwq4xv88Fng/bFk9vF+rl8Yoew1wtf/5NmBHxPLQF5b/hf1PoB/+r21/fhZwFOgaNu8uYEnYOjaHLWvkb0OrKPW5FO+LPnz9s4CxYdtXWZA4DlwRZX5nv8w2xA4Sp/vzL/ani/B/zYe97oxR7ljgnbDpQX7eyL/b2XhNjkeBhmHpbwQWR+7zsOXTgT+GTQ8EPvY/3wx8GJF+ub+edv4+OSNs2Ut8GyTuwPtR0iPo/5W68joV23AlWPuBFmaW7WK0o+P9st8eNr3dnxdah3OuzP/8jf++O2z5N8CZYdM7yz84506YWXH5+szsFuBBvC9K/HwtouWN5Jxb5HfcTgLamdlreGcqDfF+9UduQ5uw6c/D1nPEzMrLjnQesNM5d6KSdVVmH9A6yvzWwAm8X9jnxMhbXsYXYfOucc4tiLPsyL/Jvih/tzPxtrE+sMvfD+AF85j73vd52OcjfLv/Io8f+HafnQcccM59HbHsfP/zi/7nl83sbLympzHOuWNV1EVqSB3XEmk5XjPENZWk+QyvQ7ZcO39eTZV/AWBm9YC2wGdm1h74AzACrynibGAdXtNRuUqHMXbOPeOc643XzPUd4Od4X8zHomzDpzWo+2fA+X69a7KuBXhNapF+Aix3zh2pJO+1eP0Hm+Isq6Z24p1JtHDOne2/znLOdfOXV3co6cjjB77dZ7uApmZ2RsQyryDnjjnnfumc64rXr/Ej4JZqli/VoCAhFTjnDuH1J0wys2vMrJGZ1TezK81sop9sFvALM2tpZi389Ilcx97bzK7zr066H+8L6X289m+H16eBmd0OdI93pWb2/8ysr5nVx+vLKAHK/F/Ls4FxZtbYD0YP1nAbPvDX/bC/nwrxmm5ejjP/L4GLzGycmTXz63Mv3hffIzG261wzGwE8DoyKOItJOufcLmA+8FszO8vM6plZJzP7np9kN9DWzBrEucq5wHf8y6yzzWww0BV40zm3HVgJ/NLMGviXXQ8qz2hm/c3sQjPLAg7jBfuyKGVIkihIyEmcc0/hfWn+Au8Leifer/nX/ST/jvePvBbvqpiP/Hk19We8Tunyzszr/F+MG4Df4p3d7AYuxLuaKV5n4Z2JHMBrstiPd8USeJ3dXwNb8a5kegl4rroVd86VAlcBV+KdoUwGbnHOfRxn/k+AS/Auay3C+yV9PfBD51zkth40s6/x9vlAvCvQIuv8RsR9Eq9Vd5tiuAWviW4D3v78E982ky3Cu6rrczPbV9WKnHP78c4A/g3vb/Iw8CPnXHnem4C+eM1oj+N1qpdr5Zd9GNiIdxWUbrRLIfM7g0QCYWZjgf/jnPvXoOsiIifTmYSIiMSkICEiIjElJUiY2RVmtsnMNpvZyCjLzcye8ZevNbNeVeX1B/X61MzW+K+ByairnFqcc2PV1CRy6ko4SPhXGUzC67jrCtxoZl0jkl0JXOC/huMNORBP3qedc3n+a26idRURkepJxs10ffDuTt0KYGYv4w0EtyEszdXAC87rJX/fzM42s9Z4N0hVlTduLVq0cDk5OTXdDhGROmnVqlX7nHMtoy1LRpBoQ8U7L4vxLl+rKk2bOPKO8O+4XQn8m3PuQGThZjYc7+yEdu3asXLlyhpuhohI3WRmkXfAhySjT8KizIu8rjZWmsryTgE6AXl4147/Nlrhzrlpzrl851x+y5ZRA6GIiNRQMs4kigkbVgF/SIU40zSIldc5FxpXxsz+ALyZhLqKiEg1JONMYgVwgZl18G/LHwLMiUgzB7jFv8qpH3DIv9U/Zl6/z6LctXhj9oiISBolfCbhnDvujyMzD28I5uecc+vN7G5/+VS8sVoG4o27fwS4vbK8/qonmlkeXvNTEd5QziIikkYZNSxHfn6+U8e1iEj1mNkq51x+tGW641pEpDabORNycqBePe995sykrl4PHRIRqa1mzoThw+GI/9iR7du9aYChQ5NShM4kRERqqzFjvg0Q5Y4c8eYniYKEiEhttWNH9ebXgIKEiEht1a5d9ebXgIKEiEhtNW4cNGpUcV6jRt78JFGQEBGprYYOhWnToH17MPPep01LWqc16OomEZHabejQpAaFSDqTEBGRmBQkREQSleIb2oKk5iYRkUSk4Ya2IOlMQkQkEWm4oS1IChIiIolIww1tQVKQEBFJRBpuaAuSgoSISCLScENbkBQkREQSkYYb2oKkq5tERBKV4hvagqQzCRHJHBl8v0JQdCYhIpkhw+9XCIrOJEQkM2T4/QpBUZAQkcyQ4fcrBEVBQkQyQ4bfrxAUBQkRyQwZfr9CUBQkRCT5grjKKMPvVwiKrm4SkeQK8iqjDL5fISg6kxCR5NJVRhlFQUJEkktXGWUUBQkRSS5dZZRRFCREMlVQQ1ToKqOMoiAhkonKO4+3bwfnvu081lVGUk3mnAu6DkmTn5/vVq5cGXQ1RIKXk+MFhkjt20NRUbprI6c4M1vlnMuPtkxnEiKpFkSzjzqPJUkUJERSKahmH3UeS5IoSEjdEFQnblD3DKjzWJJEQUIyX5CduEE1+6jzWJJEQULSK4hf9EHeARxks8/QoV4n9YkT3rsChNSAgoSkT1C/6IPsxFWzj9RyChJBCqqdvK61zwf9a17NPlKLJSVImNkVZrbJzDab2cgoy83MnvGXrzWzXlXlNbNmZvaOmX3ivzdNRl1PMnEi/7h6DMXZOZywehRn5/CPq8fAxIkpKS7kxhspu21YhV/VZbcNgxtvzMxyIfp1+5XNT5aCAsqyT6swqyz7NCgoSG25ABMnsmDDeeRQRD1OkEMRCzacl/rja+JEFoxZXOG3wIIxi1NfbpBl17Vy01W2cy6hF5AFbAE6Ag2AvwNdI9IMBN4CDOgHfFBVXmAiMNL/PBL4j6rq0rt3b1dda68a7U6AW0Y/92tGumX0cyfArb1qdLXXVR1HzjrXOe9rusLryFnnZmS5QZb9zuhF7jCN3S7OdWWY28W57jCN3TujF6W03PKy99DCFbLIgXOFeNOpLjuocoMsu66Vm8yygZUuxvdqMp4n0QfY7JzbCmBmLwNXAxvC0lwNvOBX5n0zO9vMWgM5leS9Gij08z8PLAEeSUJ9K2j6l5m8Tz8GsJBSGtCAUkYzgBNvTGHxM+eelN7FeYd6Venc4d0YEJnKHd6NPf10nLWvvvJyiSjbHd6NPfVUysoFeO34AH7EK2RTFpp3jCzePD6Aa3/725SVO24KNGMoN/MCyxhAAct4kZv5YvIqxjRdlbJyAcZNhmZcz7/yI+pTQAHL6ce/8sXklYw+O77RAeI95sKNnwzNuI6h/Ihs+tGP9+nHUL6Y/CEjz/ow7vXUpOz/mAxNuYab+BFZ9KUfH9CXmzgw+X0ePvP9lJX7m8nQlKu5kR9Rjz705UP6MoQDk5fxUKNlca+numU/OQmaMoghfrl9+JA+DOagX66ZYeb918XzXp00Y6cYzRnCMAbRlX9lLK/yE2azbWZ/ipLU7ZWMINEG2Bk2XQz0jSNNmyrynuuc2wXgnNtlZudEK9zMhgPDAdrVoI35vLIdvMiNlNKAMrIpxfFPCmni3mff5s1R85T/gapSWbrDNOYsvvTSRc5PYYdqeLnhZR+mMWcVF6esXIClR85lL5cxgEWcTilHacACBrDpSCv6fvZZyso9cAAO0Ih59KCABbxDP7ZyJhz8nM8/T1mxABw8CAc5kwX04CIWsoACtnIWHNzDnj3xryfeY65iuU1YRC4XsYjFFLCVpnDwC774olqrqmHZzfkreVzCYv7KRWyjBRw8zOHDqS73HN6jJ5ewhKVczDZawcEjJ3WFJbPsQ4fgEOexjF5cyhLe4xKKaA0Hj/D116HWlLjeq5vGO7Ydr9GFJ/k9U3iUJfTHkvkVEusUI94XcAPwx7Dpm4HfRaT5C3BJ2PRCoHdleYGDEes4UFVdatLctDOrvVtGP9eQr10Wpa4hX7tl9HM7s9pXe13VcW/zGe4Ip1VsduE0d2/zGRlZrnPOtW//7enwL3k0dJrcvn1mlhtk2drmzC83mWVTSXNTMoJEATAvbHoUMCoize+BG8OmNwGtK8tbnsb/3BrYVFVdalOfRFDt5Gqfr51txrWl3CDLrmvlJrPsyoJEMq5uWgFcYGYdzKwBMASYE5FmDnCLf5VTP+CQ85qSKss7B7jV/3wr8Ock1PUkF17chHVXjeb8rF08wn9wftYu1l01mgsvbpKK4kJ+0GQFH4z+M/3af062naBf+8/5YPSf+UGTFRlZbnnZfx89m23t+2MG29r35++jZ6dlm4MoN8iytc2ZX266yk7KUOFmNhD4T7yrlZ5zzo0zs7sBnHNTzWvgexa4AjgC3O6cWxkrrz+/OTAbaAfsAG5wzlXamqqhwkVEqq+yocL1PAkRkTpOz5MQEZEaUZAQEZGYFCRERCQmBQkREYlJQUJERGJSkAjY8uUwfrz3LiJyqknG2E1SQ8uXw4ABUFoKDRrAwoXpGb1aRCReOpMI0JIlUHrUUVbmvS9Zkr6ydQYjIvHQmUSACpv/gwYnOlFKfRqcOEZh8y3AhSkvV2cwIhIvnUkEqGD/myysdzlP8BgL611Owf4301LukiVegCgr897TeQYjIrWLziSCVFhIwWlPUFD6vveTvvA36SqWBg2+PZMoLExLsSJSCylIBKmgwGvrWbLE+6ZOU5tPQMWKSC2kAf5EROo4DfAnJ9PlTSISBzU31UW6vElE4qQzibpIlzeJSJwUJOqi8subsrJ0eZOIVErNTXWRLm8SkTgpSNRVBQUKDiJSJTU3iYhITAoSIiISk4KEpJ/u0RCpNRQkJL2WL2d54SjGj/mK5YWjFChETnHquJa0Wv7CJwwonUspDWhQWsrCF/5EgTrQRU5ZOpOQtFrC9yilAWVkU0p9lvC9oKskIpVQkJC0KrylPQ1OM7KsjAan1aPwlvZBV0lEKqHmJkmrggJYuDhL9/GJ1BIKEpJ2Qd3Ht3y5bjIXqS4FCakTNPCtSM2oT0LqBA18K1IzChJSJ2jgW5GaUXOT1Aka+FakZhQkpM4oYDkFLAEKAUUJkXgoSEjdoJ5rkRpRn4TUDeq5FqkRBQmpG9RzLVIjam6SuiHgnmvdyCe1lYKE1B0B3eqt7hCpzRJqbjKzZmb2jpl94r83jZHuCjPbZGabzWxkVfnNLMfMvjGzNf5raiL1FAmSukOkNku0T2IksNA5dwGw0J+uwMyygEnAlUBX4EYz6xpH/i3OuTz/dXeC9RQJjLpDpDZLNEhcDTzvf34euCZKmj7AZufcVudcKfCyny/e/CK1Wnl3yBNPqKlJap9E+yTOdc7tAnDO7TKzc6KkaQPsDJsuBvrGkb+Dma0GDgO/cM69F60CZjYcGA7Qrl27hDZGJFWCGvlWJFFVBgkzWwC0irJoTJxlWJR5roo8u4B2zrn9ZtYbeN3MujnnDp+0IuemAdMA8vPzq1qviIhUQ5VBwjn3g1jLzGy3mbX2zwJaA3uiJCsGzg+bbgt85n+Omt85dxQ46n9eZWZbgO8AK+PZKBERSY5E+yTmALf6n28F/hwlzQrgAjPrYGYNgCF+vpj5zayl3+GNmXUELgC2JlhXkeAsXw7jx3vvIrVIon0SE4DZZjYM2AHcAGBm5wF/dM4NdM4dN7MRwDwgC3jOObe+svzAd4FfmdlxoAy42zn3RYJ1FQmGbpSQWiyhIOGc2w8MiDL/M2Bg2PRcYG418r8KvJpI3UROGdFulFCQkFpCYzeJpFqAN0qolUsSpWE5RFItoHGj1MolyaAgIZIOAdwooVYuSQY1N4lkKA0HIsmgMwmRDKXneksyKEiIZDANByKJUnOTiIjEpCAhIiIxKUiIiEhMChIiIhKTgoSIiMSkICEiIjEpSIiISEwKEiKZTCP8SYJ0M51Ipgp4hL/ly3W3dyZQkBDJVAGO8KcRaDOHmptEMlWAI/xFi09SO+lMQiRTBTjCX3l8Kj+T0Ai0tZeChEgmC2iEP41AmzkyPkgcO3aM4uJiSkpKgq5KrXD66afTtm1b6tevH3RVpJbTCLSZIeODRHFxMY0bNyYnJwczC7o6pzTnHPv376e4uJgOHToEXR0ROQVkfMd1SUkJzZs3V4CIg5nRvHlznXWJSEjGBwlAAaIatK9EJFydCBIiIlIzChIRZs6EnByoV897nzkz8XVmZWWRl5dH9+7dGTRoEAcPHgSgqKiI7t27h9L94Q9/oFevXhw4cCA078knn8TM2LdvXyhPw4YNycvLIy8vj7vvvjvxCoqIxJDxHdfVMXMmDB8OR45409u3e9MAQ4fWfL0NGzZkzZo1ANx6661MmjSJMWPGVEjz4osv8rvf/Y5FixbRtGlTAHbu3Mk777xDu3btKqTt1KlTaH0iIqmkM4kwY8Z8GyDKHTnizU+WgoICPv300wrzZs+ezYQJE5g/fz4tWrQIzX/ggQeYOHGi+glEJDAKEmF27Kje/OoqKytj4cKFXHXVVaF527dvZ8SIEcyfP59WrVqF5s+ZM4c2bdqQm5t70nq2bdtGz549+d73vsd7772XnMqJiEShIBEmolWnyvnx+uabb8jLy6N58+Z88cUXXHbZZaFlLVu2pF27dsyePTs078iRI4wbN45f/epXJ62rdevW7Nixg9WrV/PUU09x0003cfjw4cQqKJJBNDp6cilIhBk3Dho1qjivUSNvfiLK+yS2b99OaWkpkyZNClt/I9566y2mTp3KTL+XfMuWLWzbto3c3FxycnIoLi6mV69efP7555x22mk0b94cgN69e9OpUyf++c9/JlZBkQxRPvrso4967woUiVOQCDN0KEybBu3bg5n3Pm1aYp3W4Zo0acIzzzzDk08+ybFjx0LzW7Zsydtvv83o0aOZN28eF154IXv27KGoqIiioiLatm3LRx99RKtWrdi7dy9lZWUAbN26lU8++YSOHTsmp4IiyRTAT3qNPpt8uropwtChyQsK0fTs2ZPc3FxefvllLr300tD8Dh06MGfOHAYOHMj//u//0rdv36j53333XR577DGys7PJyspi6tSpNGvWLHUVFqmJgB4oodFnk09BIg2++uqrCtNvvPFG6PO6detCn3Nzc0+68gm8eyPKXX/99Vx//fXJr6RIMgX0wCONPpt8ChIiknwB/qTX6LPJpSAhIsmnn/QZQ0FCRFJDP+kzgq5uEhGRmBIKEmbWzMzeMbNP/PemMdJdYWabzGyzmY0Mm3+Dma03sxNmlh+RZ5SffpOZ/TCReoqISM0keiYxEljonLsAWOhPV2BmWcAk4EqgK3CjmXX1F68DrgPejcjTFRgCdAOuACb76xERkTRKNEhcDTzvf34euCZKmj7AZufcVudcKfCynw/n3Ebn3KYY633ZOXfUObcN2OyvJ/VSMFZ4TYYK//nPf07nzp3p0aMH1157bSjPhx9+GBomPDc3l9deey3h+omIxJJokDjXObcLwH8/J0qaNsDOsOlif15l4s5jZsPNbKWZrdy7d2/cFY+qfKzw7dvBuW/HCk8wUJQPy7Fu3TqaNWtWYViOcuVDhc+fP5+mTZty2WWXsW7dOtauXct3vvMdxo8fD0D37t1ZuXIla9as4e233+auu+7i+PHjCdVPRCSWKoOEmS0ws3VRXlfHWUa0ca5dsvI456Y55/Kdc/ktW7aMs0oxpGGs8HiHCr/88svJzvYuPuvXrx/FxcWAN9ZT+fySkhINIy4iKVXlJbDOuR/EWmZmu82stXNul5m1BvZESVYMnB823Rb4rIpia5IncSkeK7x8qPBhw4aF5pUPFb569eoKQ4WHe+655xg8eHBo+oMPPuCOO+5g+/btvPjii6GgISKSbIk2N80BbvU/3wr8OUqaFcAFZtbBzBrgdUjPiWO9Q8zsNDPrAFwAfJhgXauWorHCqztUeLhx48aRnZ3N0LABpfr27cv69etZsWIF48ePp6SkJKH6iYjEkmiQmABcZmafAJf505jZeWY2F8A5dxwYAcwDNgKznXPr/XTXmlkxUAD8xczm+XnWA7OBDcDbwE+dc2UJ1rVqKRorvLpDhZd7/vnnefPNN5k5c2bUZqUuXbpwxhlnVBj/SUQkqZxzGfPq3bu3i7Rhw4aT5lVqxgzn2rd3zsx7nzGjevmjOOOMM0KfP/roI3f++ee70tJSt23bNtetWzfnnHNbt2517dq1c2+//bZzzrm33nrLdenSxe3Zs6fCurZu3eqOHTvmnHOuqKjItW7d2u3duzfhOoar9j4TkVoNWOlifK+qMTtSiscKj3eo8BEjRnD06NFQ01S/fv2YOnUqS5cuZcKECdSvX5969eoxefLkCs/FFhFJJvOCSGbIz893K1eurDBv48aNdOnSJaAa1U7aZyLVt3x57R3P0MxWOefyoy3TmYSISIICesZSWmiAPxGRBGXyY1MVJEREElT+jKWsrMx7bKqam0QkswTQOZDJz1hSkBCRzBFg50CmPmNJzU0ikjkyuXMgIAoS4SZOhMWLK85bvNibn4DyocK7detGbm4uTz31FCdOnAgtX7p0KX369KFz58507tyZadOmhZaNHTuWNm3ahIYanzPHG9Fkx44d9O/fn549e9KjRw/mzp17Unl5eXlcddVVofnDhg0jNzeXHj168OMf/5ivvvoqoe0SOeVkcudAUGLdZVcbXwnfcb1okXMtWnjv0aZrKPyO6927d7sBAwa4xx57zDnn3K5du9z555/vVq1a5Zxzbu/eva5Xr17uzTffdM459/jjj7vf/OY3oW1p3ry5Kysrc3feeaebPHmyc8659evXu/bt20ctL9yhQ4dCnx944AE3fvz4qOl0x7XUasuWOffrX3vvEhcqueNaZxLh+veH2bPhJz+Bxx7z3mfP9uYnyTnnnMO0adN49tlncc4xadIkbrvtNnr16gVAixYtmDhxIhMmTDgpb5cuXcjOzmbfvn2YGYcPHwbg0KFDnHfeeVWWfdZZZwHeD4NvvvlGw4xLZioogFGjMrODIAAKEpH694d77oEnnvDekxggynXs2JETJ06wZ88e1q9fT+/evSssz8/PZ/369Sfl++CDD6hXrx4tW7Zk7NixzJgxg7Zt2zJw4EB+97vfhdKVlJSQn59Pv379eP311yus4/bbb6dVq1Z8/PHH3HvvvUnfNhHJLAoSkRYvhilT4NFHvffIPookcf5wKM65qL/ow+c9/fTT5OXl8dBDD/HKK69gZsyaNYvbbruN4uJi5s6dy8033xzq59ixYwcrV67kpZde4v7772fLli2hdf33f/83n332GV26dOGVV15JybaJSOZQkAi3ePG3TUy/+tW3TU9JDhRbt24lKyuLc845h27duhE53tSqVavo2rVraPqBBx5gzZo1vPfee6FBAf/rv/6Ln/zkJ4D3tLuSkhL27dsHEGp66tixI4WFhaxevbrC+rOyshg8eDCvvvpqUrdLRDKPgkS4FSsq9kGU91GsWJG0Ivbu3cvdd9/NiBEjMDN++tOfMn36dNasWQPA/v37eeSRR3j44YcrXU+7du1YuHAh4A3IV1JSQsuWLTlw4ABHjx4FYN++ffztb3+ja9euOOfYvHkz4J29vPHGG3Tu3Dlp2yUimUk304WL9sXcv3/C/RLlT6Y7duwY2dnZ3HzzzTz44IMAtG7dmhkzZnDnnXfy5Zdf4pzj/vvvZ9CgQZWu87e//S133nknTz/9NGbG9OnTMTM2btzIXXfdRb169Thx4gQjR46ka9eunDhxgltvvZXDhw/jnCM3N5cpU6YktF0ikvk0VLicRPtMpG6pbKhwNTeJiEhMChIiIhKTgoSIiMSkICEiIjEpSIiISEwKEiIiEpOCRBqUD93dvXt3Bg0axMGDB6uVv7Cw8KS7skVE0kFBIg0aNmzImjVrWLduHc2aNWPSpElBV0lEMsjy5TB+vPeebAoSaVZQUMCnn34KwIcffshFF11Ez549ueiii9i0aRPg3aE9ZMgQevToweDBg/nmm29C+e+55x7y8/Pp1q0bjz/+eGh+Tk5OaOymlStXUqiHrYjUCeVPbH30Ue892YGizg3LkYpnKMR713pZWRkLFy5k2LBhAHTu3Jl3332X7OxsFixYwOjRo3n11VeZMmUKjRo1Yu3ataxduzb0rAmAcePG0axZM8rKyhgwYABr166lR48eSd8mEakdoj2xNZmP0qhzQSKIYUjKx24qKiqid+/eXHbZZYD3sKBbb72VTz75BDPj2LFjALz77rv87Gc/A6BHjx4VgsDs2bOZNm0ax48fZ9euXWzYsEFBQqQOK39ia2lpap7YquamNCjvk9i+fTulpaWhPolHH32U/v37s27dOt544w1KSkpCeaKd8Wzbto0nn3yShQsXsnbtWv7lX/4llCc7Ozv0PInw9YhIZisogIULveekLVyY/AfyKUikUZMmTXjmmWd48sknOXbsGIcOHaJNmzYATJ8+PZTuu9/9LjNnzgRg3bp1rF27FoDDhw9zxhln0KRJE3bv3s1bb70VypOTk8OqVasA9JwIkSCksve4Cql8YquCRJr17NmT3NxcXn75ZR5++GFGjRrFxRdfTFlZWSjNPffcw1dffUWPHj2YOHEiffr0ASA3N5eePXvSrVs37rjjDi6++OJQnscff5z77ruPSy+9lKysrLRvl0idlure4wBpqHA5ifaZSDWNH+8FiLIyyMry2n5GjQq6VnHTUOEiIqlU3nuclZWa3uMA1bmrm0REkq6893jJEi9ApKJzICAKEiIiyVBQkFHBoZyam0REJCYFCRERiUlBQkREYkooSJhZMzN7x8w+8d+bxkh3hZltMrPNZjYybP4NZrbezE6YWX7Y/Bwz+8bM1vivqYnUM2jlQ4Xn5ubSq1cvli1bBkBRURENGzakZ8+edOnShT59+vD888+H8k2fPp2WLVuSl5dHXl4et9xyS1CbICJ1VKId1yOBhc65Cf6X/0jgkfAEZpYFTAIuA4qBFWY2xzm3AVgHXAf8Psq6tzjn8hKsX80sX57UqxTKh+UAmDdvHqNGjeKvf/0rAJ06dWL16tUAbN26leuuu44TJ05w++23AzB48GCeffbZhOsgIlITiTY3XQ2U//R9HrgmSpo+wGbn3FbnXCnwsp8P59xG59ymBOuQXCm+c/Lw4cM0bRr1hIuOHTvy1FNP8cwzzyS1TBGRmkr0TOJc59wuAOfcLjM7J0qaNsDOsOlioG8c6+5gZquBw8AvnHPvRUtkZsOB4QDt2rWrTt2jS8G4u+WjwJaUlLBr1y4WLVoUM22vXr34+OOPQ9OvvPIKS5cuBeC+++4LnWGIiKRDlUHCzBYAraIsGhNnGdEe4FDVWCC7gHbOuf1m1ht43cy6OecOn7Qi56YB08AbliPOOsWWgnF3w5ubli9fzi233MK6deuipo0cJkXNTSISpCqDhHPuB7GWmdluM2vtn0W0BvZESVYMnB823Rb4rIoyjwJH/c+rzGwL8B0g9Q96TvGdkwUFBezbt4+9e/dGXb569WqNmyQip4xEm5vmALcCE/z3P0dJswK4wMw6AJ8CQ4CbKlupmbUEvnDOlZlZR+ACYGuCdY1fCu+c/PjjjykrK6N58+YcOXKkwrKioiIeeugh7r333pSULSJSXYkGiQnAbDMbBuwAbgAws/OAPzrnBjrnjpvZCGAekAU855xb76e7Fvgd0BL4i5mtcc79EPgu8CszOw6UAXc7575IsK6BKe+TAK856fnnnw8N571lyxZ69uxJSUkJjRs35t5771W/g4icMjRUuJxE+0ykbtFQ4SIiUiMKEiIiEpOChIiIxKQgISIiMSlIiIhITAoSIiISk4JEGpx55pknzdu0aROFhYXk5eXRpUsXhg8fHlq2dOlS+vTpQ+fOnencuTPTpk0LLRs7dixt2rQhLy+PCy64gOuuu44NGzZUWPfq1asxM+bNm5e6jRKROkFBIorly2H8+KQPAFvBz372Mx544AHWrFnDxo0bQ3dZf/7559x0001MnTqVjz/+mKVLl/L73/+ev/zlL6G85fk++eQTBg8ezPe///0Kw3zMmjWLSy65hFmzZqVuA0SkTlCQiJDikcJDdu3aRdu2bUPTF154IQCTJk3itttuo1evXgC0aNGCiRMnMmHChKjrGTx4MJdffjkvvfQS4N3R/ac//Ynp06czf/58SkpKUrMBIlInKEhEiDZSeCo88MADfP/73+fKK6/k6aef5uDBgwCsX7+e3r17V0ibn5/P+vXrY64rfHjxv/3tb3To0IFOnTpRWFjI3LlzU7MBIlInKEhEKB8pPCsraSOFR3X77bezceNGbrjhBpYsWUK/fv04evQozjnMTh5dPdq8cuFDq8yaNYshQ4YAMGTIEDU5iUhCEh3gL+OkeKTwCs477zzuuOMO7rjjDrp37866devo1q0bK1eu5KqrrgqlW7VqFV27do25ntWrV5Ofn09ZWRmvvvoqc+bMYdy4cTjn2L9/P19++SWNGzdO3YaISMbSmUQUBQUwalRqA8Tbb7/NsWPHAK+zev/+/bRp04af/vSnTJ8+PfSQov379/PII4/w8MMPR13Pq6++yvz587nxxhtZsGABubm57Ny5k6KiIrZv387111/P66+/nroNEZHgpfBqG51JpMGRI0cqdFI/+OCDFBcXc99993H66acD8Jvf/IZWrbwHAM6YMYM777yTL7/8Eucc999/P4MGDQrlf/rpp5kxYwZff/013bt3Z9GiRbRs2ZJZs2Zx7bXXVij7+uuvZ8qUKdx8881p2FIRSbvyq23Kn6a5cGFSf+FqqHA5ifaZSC0yfrx3OWZZmdeZ+sQTXlNINWiocBGRTJXiq23U3CQiUpul+GqbOhEkYl1WKifLpOZHkTqjoCBlV9pkfHPT6aefzv79+/XlF4fyS2bLO9NFRDL+TKJt27YUFxdXGNtIYjv99NMrXIklInVbxgeJ+vXr06FDh6CrISJSK2V8c5OIiNScgoSIiMSkICEiIjFl1B3XZrYX2J7AKloA+5JUnWRSvapH9aoe1at6MrFe7Z1zLaMtyKggkSgzWxnr1vQgqV7Vo3pVj+pVPXWtXmpuEhGRmBQkREQkJgWJiqYFXYEYVK/qUb2qR/WqnjpVL/VJiIhITDqTEBGRmBQkREQkpjoRJMzsCjPbZGabzWxklOVmZs/4y9eaWa9486a4XkP9+qw1s2Vmlhu2rMjM/mFma8xsZWTeFNer0MwO+WWvMbPH4s2b4nr9PKxO68yszMya+ctSub+eM7M9ZrYuxvKgjq+q6hXU8VVVvYI6vqqqV9qPLzM738wWm9lGM1tvZvdFSZPa48s5l9EvIAvYAnQEGgB/B7pGpBkIvAUY0A/4IN68Ka7XRUBT//OV5fXyp4uAFgHtr0LgzZrkTWW9ItIPAhalen/56/4u0AtYF2N52o+vOOuV9uMrznql/fiKp15BHF9Aa6CX/7kx8M90f3/VhTOJPsBm59xW51wp8DJwdUSaq4EXnOd94Gwzax1n3pTVyzm3zDl3wJ98H0jHGN6JbHOg+yvCjcCsJJVdKefcu8AXlSQJ4viqsl4BHV/x7K9YAt1fEdJyfDnndjnnPvI/fwlsBNpEJEvp8VUXgkQbYGfYdDEn7+RYaeLJm8p6hRuG92uhnAPmm9kqMxuepDpVp14FZvZ3M3vLzLpVM28q64WZNQKuAF4Nm52q/RWPII6v6krX8RWvdB9fcQvq+DKzHKAn8EHEopQeXxn/PAm8U7BIkdf9xkoTT96ainvdZtYf75/4krDZFzvnPjOzc4B3zOxj/5dQOur1Ed5YL1+Z2UDgdeCCOPOmsl7lBgF/c86F/ypM1f6KRxDHV9zSfHzFI4jjqzrSfnyZ2Zl4Qel+59zhyMVRsiTt+KoLZxLFwPlh022Bz+JME0/eVNYLM+sB/BG42jm3v3y+c+4z/30P8BreqWVa6uWcO+yc+8r/PBeob2Yt4smbynqFGUJEU0AK91c8gji+4hLA8VWlgI6v6kjr8WVm9fECxEzn3P9GSZLa4yvZHS2n2gvvbGkr0IFvO2+6RaT5Fyp2/HwYb94U16sdsBm4KGL+GUDjsM/LgCvSWK9WfHsjZh9gh7/vAt1ffromeO3KZ6Rjf4WVkUPsjti0H19x1ivtx1ec9Ur78RVPvYI4vvztfgH4z0rSpPT4yvjmJufccTMbAczD6+1/zjm33szu9pdPBebiXSGwGTgC3F5Z3jTW6zGgOTDZzACOO2+Ux3OB1/x52cBLzrm301ivHwP3mNlx4BtgiPOOyqD3F8C1wHzn3Ndh2VO2vwDMbBbeFTktzKwYeByoH1avtB9fcdYr7cdXnPVK+/EVZ70g/cfXxcDNwD/MbI0/bzRegE/L8aVhOUREJKa60CchIiI1pCAhIiIxKUiIiEhMChIiIhKTgoSIiMSkICEiIjEpSIiISEz/Hy72ijt5QnkRAAAAAElFTkSuQmCC\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", "plt.show()" ] }, { "cell_type": "code", "execution_count": 4, "id": "17efc931", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkGElEQVR4nO3de3xU1b338c+PgMV4LRCrgCTg87QIlnAJ1ku1pB5PLS2tpw9iezy2VB6pedV66bFU4ahUD8+xlNbW1uJJb5wWvKRYPWq1YjU9qLVIsJGCeE+ACEKgIqKgCL/nj70TJiGXmcye2dnJ9/16zSsz+7LW2nv2/LJm7TVrmbsjIiLJ1SfuAoiISHYUyEVEEk6BXEQk4RTIRUQSToFcRCThFMhFRBJOgVzyzswuMLNlcZejiZkdamb3m9mbZvbbuMuTL2a2yMz+PaK05prZ4ijSkswpkCeYmf2zmdWY2S4z22xmD5nZx+MuV2fcfYm7/2Pc5UgxFfgQMNDdz2trAzMbZWb3hcH+LTOrNrPTUtaXmJmH78UuM9tiZg+Y2dmt0qk3s90p2+0ys5/k9vDAzKab2RO5zkfioUCeUGb2TeCHwP8jCELDgJ8Cn4+xWJ0ys75xl6ENxcCL7v5+WyvN7ATgSeBvwHBgMHAPsMzMTm21+dHufjhQCjwC3GNm01ttM8XdD095XBrhsUhv5O56JOwBHAXsAs7rYJsPEAT6TeHjh8AHwnWTgAZgFrAV2AycC0wGXgT+DsxOSWsusBS4C3gLeAYoTVl/NfBKuO454J9S1k0nCII3h+n+e7jsiXC9heu2Am8Cq4GTUo7z10AjsB74N6BPSrpPAAuAN4A64NMdnI8TgT8BO4C1wOfC5d8B3gP2hud0Rhv7/gZ4sI3lC4Hl4fMSwIG+rba5CtiSUu564B/SfJ/nAr8FFofn9m/Ah4FrwvO1EfjHVtfFL8L387XwXBeEx74H2Bce445w+0XArcDvw/RXACekpHcasDJ8X1YCp6WsGw78T7jfI8BPgMXhuv5hmbeH53sl8KG4Pzc9+RF7AfTowpsG5wDvtw4arba5AfgLcAxQBPwZuDFcNync/zqgH3BxGCxvB44ARocf/BHh9nPDQDc13P6qMHD2C9efR1BL7QOcD7wNHBeumx7m9Q2gL3AoLQP5p4BVwNEEQf3ElH1/Dfx3WKYSgn8yM1LS3RuWvQCoIPiHZW2ci37Ay8Bs4BDgk2EA+kjK8S3u4Fy+Dny1jeXlYXAspP1APiJcfmL4up7MAvme8Bz1Dc9HHTAn5X2rS9n+XuA/gcPC9/1p4Gsp5+uJVukvIvjnenKY/hLgznDdAIJ/kBeG674Uvh4Yrn8K+AFBheHM8Hw2BfKvAfeH56UAmAAcGffnpic/4ssYfklQq1gTUXrDgGXAOoJaYUncJzeH5+4C4PVOtnkFmJzy+lNAffh8ErAbKAhfHxEGm4+lbL8KODd8Phf4S8q6PgS1vjPaybsW+Hz4fDqwodX65qBCEFRfBE4hrLWGywuAd4FRKcu+BvwpJY2XU9YVhsdwbBvlOYMgGKemfwcwN+X4Ogrk7wPntLF8ZJjnENoP5P3D5aeHr+sJa8Upj4vbyXcu8EjK6ynhvq3ft6MJmtfeBQ5N2f5LQHXrc56yfhHw85TXk4Hnw+cXAk+32v6pMJ1h4Tk5LGXd7RwI5BcRVBzGxP1Z6S2PONvIFxHULKPya+B77n4iQQ1ja4RpdzfbgUGdtDcPJmiOaLI+XNachrvvC5/vDv9uSVm/Gzg85fXGpifuvp+gaWYwgJl92cxqzWyHme0ATgIGtbVva+7+GMHX8luBLWZWaWZHhvsf0sYxDEl5/XpKOu+ET1PL3GQwsDEsd3tpdWQbcFwby48D9hPUVNvTlMffU5ad6+5Hpzx+1sH+rd+TbW28b4cTtPP3AzanvA//SVAz78jrKc/f4cD5a339wIFzNhh4w93fbrWuyW+Ah4E7zWyTmc03s36dlEOyEFsgd/fltLy4MbMTzOwPZrbKzB43s5HppGVmowhqQo+Eae9K+WD3RE8RfOU+t4NtNhF8uJsMC5d11fFNT8ysDzAU2GRmxcDPgEsJvnYfDawhaCZp0uEQm+5+i7tPIGjS+TDwLYLgubeNY3itC2XfBBwflrsraf2RoPmotWnAU51ca/9EUKl4Ic28umojQY18UMo/iCPdfXS4PtNhTltfP3DgnG0GPmhmh7VaF2Tkvtfdv+Puowja2T8LfDnD/CUD3a3XSiXwjfBDfRVBL4x0fBjYYWa/M7O/mtn3zKwgZ6WMmbu/SdC+fauZnWtmhWbWz8w+bWbzw83uAP7NzIrMbFC4fTb9fCeY2RfCbwFXEASNvxC0xzpBGztm9lWCGnlazGyimX0srLG9TXhTLqx1VgHzzOyI8B/GN7t4DCvCtGeF52kSQTPFnWnu/x3gNDObZ2YDwvJ8gyA4fbud4/qQmV0KXA9c0+rbQOTcfTNB0+L3zexIM+sTVow+EW6yBRhqZoekmeSDwIfDLq59zex8YBTwgLuvB2qA75jZIWGX1ylNO5pZuZl9NPwM7iT4h7yvjTwkIt0mkJvZ4QT/vX9rZrUEXwuPC9d9wczWtPF4ONy9L0E76FXARIIbTNPzfQz55O4/IAhs/0YQRDcS1IrvDTf5d4IP22qC3g7PhMu66r8JbmQ23QD7Qljzeg74PsG3hC3ARwl6qaTrSIIa/RsEX8+3E/REgeAG6dvAqwQ9VG4nuLeSEXd/D/gc8GmCmv5PgS+7+/Np7v8S8HGCLoX1BDXS/wN8yt1bH+sOM3ub4JxPJuhZ1LrM97fqR35PpsfUji8TNEc9R3A+l3KgSegxgt46r5vZts4ScvftBDXpfyV4T2YBn3X3pn3/GfgYwbfq6wmaNpscG+a9k+Ce1f+QXSVCOmHu8U0sYWYlBP/hTwrbRV9w97baIjtL5xTgJnefFL6+EDjF3b8eZXl7KzObC/wvd/+XuMsiIgfrNjVyd98J1JnZeQAWKE1z95UEbXZF4etPEtRKRER6vNgCuZndQfB1/CNm1mBmMwi61c0ws2cJvgam9SvFsD31KuBRM/sbwY22jnoCiIj0GLE2rYiISPa6TdOKiIh0TSwDGA0aNMhLSkriyFpEJLFWrVq1zd2LWi+PJZCXlJRQU1MTR9YiIollZq1/bQuoaUVEJPEUyEVEEk6BXEQk4brNbC179+6loaGBPXv2xF2UROjfvz9Dhw6lXz8NKifS23WbQN7Q0MARRxxBSUkJZtb5Dr2Yu7N9+3YaGhoYPnx43MURkZh1m6aVPXv2MHDgQAXxNJgZAwcO1LcXkQzMf3I+1XXVLZZV11Uz/8n57eyRHN0mkAMK4hnQuRLJzMTBE5m2dFpzMK+uq2ba0mlMHDwx5pJlr9s0rYiI5FL58HKqplYxbek0KsoqWFizkKqpVZQPL4+7aFnrVjXyuBUUFDB27FhOOukkpkyZwo4dOwCor6/npJMOzJXws5/9jPHjx/PGGwdm+FqwYAFmxrZt25r3OfTQQxk7dixjx47lkksuyeuxiMjByoeXU1FWwY3Lb6SirKJHBHFIcCBfsgRKSqBPn+DvkiXZp3nooYdSW1vLmjVrGDBgALfeeutB2/zmN7/hxz/+McuWLeODH/wgABs3buSRRx5h2LBhLbY94YQTqK2tpba2lttuuy37AopIVqrrqllYs5Brz7yWhTULD2ozT6pEBvIlS2DmTFi/HtyDvzNnRhPMm5x66qm89lrLKR2rqqq46aabWLZsGYMGHZhb+Morr2T+/PlqtxbpxpraxKumVnFD+Q3NzSw9IZgnMpDPmQPvtJru9p13guVR2LdvH48++iif+9znmpetX7+eSy+9lGXLlnHsscc2L7/vvvsYMmQIpaUHz4FRV1fHuHHj+MQnPsHjjz8eTeFEpEtWblrZok28qc185aaVMZcse4m82blhQ2bL07V7927Gjh1LfX09EyZM4Oyzz25eV1RUxIABA6iqquLKK68E4J133mHevHksW7bsoLSOO+44NmzYwMCBA1m1ahXnnnsua9eu5cgjj8yukCLSJbNOn3XQsvLh5T2inTyRNfJWTdGdLk9XUxv5+vXree+991q0kRcWFvLQQw9x2223sSRsw3nllVeoq6ujtLSUkpISGhoaGD9+PK+//jof+MAHGDhwIAATJkzghBNO4MUXX8yugCIibUhkIJ83DwoLWy4rLAyWR+Goo47illtuYcGCBezdu7d5eVFREX/4wx+YPXs2Dz/8MB/96EfZunUr9fX11NfXM3ToUJ555hmOPfZYGhsb2bdvHwCvvvoqL730EiNGjIimgCIiKRIZyC+4ACorobgYzIK/lZXB8qiMGzeO0tJS7rzzzhbLhw8fzn333cdFF13EihUr2t1/+fLljBkzhtLSUqZOncptt93GgAEDoiugiEgoljk7y8rKvPXEEuvWrePEE0/Me1mSTOdMpHcxs1XuXtZ6eSQ1cjM72syWmtnzZrbOzE6NIl0REelcVL1WfgT8wd2nmtkhQGFnO4iISDSyDuRmdiRwJjAdwN3fA97LNl0REUlPFE0rI4BG4Fdm9lcz+7mZHRZBuiIikoYoAnlfYDyw0N3HAW8DV7feyMxmmlmNmdU0NjZGkK2IiEA0gbwBaHD3pr54SwkCewvuXunuZe5eVlRUFEG2IiICEQRyd38d2GhmHwkXnQU8l226cWgaxnb06NGUlpbygx/8gP379zevf+KJJzj55JMZOXIkI0eOpLKysnnd3LlzGTJkSPMwuPfddx8AGzZsoLy8nHHjxjFmzBgefPDBg/IbO3Zsi3FdZsyYQWlpKWPGjGHq1Kns2rUrD0cvIonl7lk/gLFADbAauBf4YEfbT5gwwVt77rnnDlrWocWL3YuL3c2Cv4sXZ7Z/Gw477LDm51u2bPGzzjrLr7vuOnd337x5sx9//PG+atUqd3dvbGz08ePH+wMPPODu7tdff71/73vfaz6WgQMH+r59+/ziiy/2n/70p+7uvnbtWi8uLm4zv1Rvvvlm8/Mrr7zS/+M//qPN7TI+ZyKSaECNtxFTI+lH7u61HjSbjHH3c939jc73ykIexrE95phjqKys5Cc/+Qnuzq233sr06dMZPz5oNRo0aBDz58/npptuOmjfE088kb59+7Jt2zbMjJ07dwLw5ptvMnjw4E7zbhpYy93ZvXu3hscVkQ4l8if6OR/HNjRixAj279/P1q1bWbt2LRMmTGixvqysjLVr1x6034oVK+jTpw9FRUXMnTuXxYsXM3ToUCZPnsyPf/zj5u327NlDWVkZp5xyCvfee2+LNL761a9y7LHH8vzzz/ONb3wj0uMSkZ4lmYE8V+PYtsHDIQzcvc2aceqym2++mbFjx3LVVVdx1113YWbccccdTJ8+nYaGBh588EEuvPDC5nb3DRs2UFNTw+23384VV1zBK6+80pzWr371KzZt2sSJJ57IXXfdFflxiUjPkcxAnqtxbFt59dVXKSgo4JhjjmH06NG0Hh9m1apVjBo1qvn1lVdeSW1tLY8//jhnnHEGAL/4xS+YNm0aEMw6tGfPnuZ5PZuaWUaMGMGkSZP461//2iL9goICzj//fO6+++5Ij0tEepZkBvJcj2MLNDY2cskll3DppZdiZnz9619n0aJF1NbWArB9+3a+/e1vM2vWwYPVpxo2bBiPPvooEAxytWfPHoqKinjjjTd49913Adi2bRtPPvkko0aNwt15+eWXgeBbwP3338/IkSMjOy4R6XkSOUNQ83i1c+YEzSnDhgVBPMtxbJtmCNq7dy99+/blwgsv5Jvf/CYQzPizePFiLr74Yt566y3cnSuuuIIpU6Z0mOb3v/99Lr74Ym6++WbMjEWLFmFmrFu3jq997Wv06dOH/fv3c/XVVzNq1Cj279/PV77yFXbu3Im7U1paysKFC7M6LhHp2TSMbYLpnIn0LjkdxlZEROKjQC4iknAK5CIiCadALiJ5Nf/J+VTXVbdYVl1Xzfwn58dUouRTIBeRvJo4eCLTlk5rDubVddVMWzqNiYMnxlyy5Epm90MRSazy4eVUTa1i2tJpVJRVsLBmIVVTqygfXh530RJLNfIUTcPKnnTSSUyZMoUdO3ZktP+kSZMO+vWniBysfHg5FWUV3Lj8RirKKhTEs5TIQJ6rNrZDDz2U2tpa1qxZw4ABA7j11luzSk9E2lZdV83CmoVce+a1LKxZeNDnWTKTyECejza2U089lddeew2Ap59+mtNOO41x48Zx2mmn8cILLwDBL0G/+MUvMmbMGM4//3x2797dvH9FRQVlZWWMHj2a66+/vnl5SUlJ81grNTU1TJo0KbIyiyRB0+e1amoVN5Tf0NzMomDedYlsI891G9u+fft49NFHmTFjBgAjR45k+fLl9O3blz/+8Y/Mnj2bu+++m4ULF1JYWMjq1atZvXp181jlAPPmzWPAgAHs27ePs846i9WrVzNmzJhIyieSZCs3rWzxeW36PK/ctFJNLF2UyEAOLdvYrj3z2kgugKaxVurr65kwYQJnn302EEwI8ZWvfIWXXnoJM2Pv3r0ALF++nMsuuwyAMWPGtAjUVVVVVFZW8v7777N582aee+45BXIRYNbpBw80Vz68XEE8C4lsWoHctLE1tZGvX7+e9957r7mN/Nprr6W8vJw1a9Zw//33s2fPnuZ92hqjvK6ujgULFvDoo4+yevVqPvOZzzTv07dv3+bxyFPTERHpqkgCuZnVm9nfzKzWzHLebSPXbWxHHXUUt9xyCwsWLGDv3r28+eabDBkyBIBFixY1b3fmmWeyJJxebs2aNaxevRqAnTt3cthhh3HUUUexZcsWHnrooeZ9SkpKWLVqFYDGGReRSERZIy9397FtjcwVtY7a2KIybtw4SktLufPOO5k1axbXXHMNp59+Ovv27WvepqKigl27djFmzBjmz5/PySefDEBpaSnjxo1j9OjRXHTRRZx++unN+1x//fVcfvnlnHHGGRQUFERWXhHpvSIZxtbM6oEyd9+WzvYaxjYaOmcivUuuh7F1YJmZrTKzme0UYKaZ1ZhZTWNjY0TZiohIVIH8dHcfD3wa+LqZndl6A3evdPcydy8rKiqKKFsREYkkkLv7pvDvVuAe4OQuphNFcXoFnSsRaZJ1IDezw8zsiKbnwD8CazJNp3///mzfvl0BKg3uzvbt2+nfv3/cRRGRbiCKHwR9CLgn7E/dF7jd3f+QaSJDhw6loaEBtZ+np3///gwdOjTuYohIN5B1IHf3V4HSbNPp168fw4cPzzYZEZFeJ7G/7BQRkYACuYhIwimQi4gknAK5iEjCKZCLiCScArmISMIpkIuIJJwCuUgvlKsJzCUeCuQivVA+JjCX/EnsnJ0i0nW5nsBc8ks1cpFeKnUC84qyCgXxBFMgF+mlcjGBucRDgVykF8r1BOaSXwrkIr1QPiYwl/yJZPLlTLU1+bKIiHQs15Mvi4hITBTIRUQSToFcRCThIgvkZlZgZn81sweiSlNERDoXZY38cmBdhOmJiEgaIgnkZjYU+Azw8yjSExGR9EVVI/8hMAvY394GZjbTzGrMrKaxsTGibEVEJOtAbmafBba6+6qOtnP3Sncvc/eyoqKibLMVEZFQFDXy04HPmVk9cCfwSTNbHEG6IiKShqwDubtf4+5D3b0E+CLwmLv/S9YlExGRtKgfuYhIwkU6sYS7/wn4U5RpiohIx1QjF4mR5s6UKCiQi8RIc2dKFDRnp0iMNHemREE1cpGYae5MyZYCuUjMNHemZEuBXCRGmjtToqBALhIjzZ0pUdCcnSIiCaE5O0VEeigFchGRhFMgFxFJOAVyEZGEUyAXEUk4BXIRkYRTIBcRSTgFchGRhFMgFxFJuKwDuZn1N7OnzexZM1trZt+JomAi+aLJHSTpoqiRvwt80t1LgbHAOWZ2SgTpiuSFJneQpMt6YgkPBmvZFb7sFz7yP4CLSBdpcgdJukjayM2swMxqga3AI+6+oo1tZppZjZnVNDY2RpGtSGQ0uYMkWSSB3N33uftYYChwspmd1MY2le5e5u5lRUVFUWQrEhlN7iBJFmmvFXffAfwJOCfKdEVySZM7SNJF0WulyMyODp8fCvwD8Hy26YrkiyZ3kKTLemIJMxsD/BdQQPCPocrdb+hoH00sISKSufYmloii18pqYFy26YiISNfol50iIgmnQC4iknAK5CIiCadALiKScArkIiIJp0Au3YZGIRTpGgVy6TY0CqFI12Tdj1wkKhqFUKRrVCOXbkWjEIpkToFcuhWNQiiSOQVy6TY0CqFI1yiQS7ehUQhFuibr0Q+7QqMfiohkrr3RD1UjFxFJOAVyEZGEUyAXEUk4BXIRkYRTIJcWNN6JSPJEMfny8WZWbWbrzGytmV0eRcEkHhrvRCR5ohhr5X3gX939GTM7AlhlZo+4+3MRpC15pvFORJIn6xq5u29292fC528B64Ah2aYr8dF4JyLJEmkbuZmVAOOAFW2sm2lmNWZW09jYGGW2EjGNdyKSLJEFcjM7HLgbuMLdd7Ze7+6V7l7m7mVFRUVRZSsR03gnIskTSSA3s34EQXyJu/8uijQlHhrvRCR5sh5rxcwM+C/g7+5+RTr7aKwVEZHM5XKsldOBC4FPmllt+JgcQboiIpKGKHqtPOHu5u5j3H1s+HgwisL1ZvphjoikS7/s7Kb0wxwRSZcmX+6m9MMcEUmXauTdmH6YIyLpUCDvxvTDHBFJhwJ5N6Uf5ohIuhTIuyn9MEdE0qXJl0VEEkKTL3eR+nOLSHenQN4J9ecWke5O/cg7of7cItLdqUaeBvXnFpHuTIE8DerPLSLdmQJ5J9SfW0S6OwXyTqg/t4h0d+pHLiKSEOpHLiLSQymQi4gknAK5iEjCRRLIzeyXZrbVzNZEkZ6IiKQvqhr5IuCciNISEZEMRBLI3X058Pco0hIRkczkrY3czGaaWY2Z1TQ2NuYrWxGRHi9vgdzdK929zN3LioqK8pWtiEiPp14rIiIJp0AuIpJwUXU/vAN4CviImTWY2Ywo0hURkc5FMrGEu38pinRERCRzaloREUk4BXIRkYRTIBcRSTgFchGRhFMgFxFJOAVyEZGEUyAXEUk4BXIRkYRTIBcRSTgFchGRhFMgFxFJOAVyEZGEUyAXEUk4BXIRkYRTIBcRSTgFchGRhFMgFxFJuKimejvHzF4ws5fN7Ooo0kw1/8n53HL5HBr6lrDf+tDQt4RbLp/D/CfnR51V25YsgZIS6NMn+LtkifLtaXnrmHXMSc7b3bN6AAXAK8AI4BDgWWBUR/tMmDDBM/Gjy2b7wG/hj5XgTvB34LfwH102O6N0umTxYvfCQnc48CgsDJYr356Rt45Zx5xrEeUN1HhbcbithZk8gFOBh1NeXwNc09E+mQbyjQXF/lgJPuhb+LXlwd/HSvCNBcUZpdMlxcUtT37TozjHefe2fOPMW8esY861iPJuL5BbsK7rzGwqcI67/9/w9YXAx9z90lbbzQRmAgwbNmzC+vXr085jv/WhD8515XDjJ+Da/4EbqmE/Rh/fn1X5O+PWB+Pgc+QYlsu8+/QJ3urWzGB/D8w3zrx1zPnLN868e8Axm9kqdy87KPmsChem3cayg0rs7pXuXubuZUVFRRllsKlgGNUlsLAsCOILy6C6JFiea6+1k0d7y6Oya0Db6be3PDLD2km/veU9IW8dc/7yjTPvHnzMUQTyBuD4lNdDgU0RpNvsd1+/gPPOg6rfBjXxqt/CeecFy3Pt2/vm8TaFLZa9TSHf3jcvp/nOpu18Z5PbfJ+Y3Ha+T0zObb4AzJsHhS3zprAwWN4T840zbx1z/vLNR95ttbdk8gD6Aq8Cwzlws3N0R/tk2kb+3Se+6z+6bLZvLCj2fZhvLCj2H10227/7xHczSqcriovdv8RiryPIu45i/xKLc96sZtZ2vma5zTeu423yeMXiFu/z4xV5uBHlHtx0Ki4OTnxxcX5ugMWdt445ccdMrtrIAcxsMvBDgh4sv3T3Dv/NlJWVeU1NTdb55sOSJTBzJrzzzoFlhYVQWQkX5PALQUkJtHUbobgY6utzl2+czYhxnWuRpMhlGznu/qC7f9jdT+gsiCfNBRcEgaS4OAhmxcX5CSxxfQuMsxlxzpyWQRyC13Pm5D7vOLsXi2StrWp6rh+ZNq30VnF8C4yzq61Z2z20ct2cFOcxi2SCXDatZCpJTSu90ZIlQS14w4agJj5vXn6aNuJqToorX5FM5bRpRXqWCy4IAtj+/cHffLVPx9WctGFDZsujpmYdyZYCuXQbcd2PiPO+QNMN3vXrg0ad9euD1wrmkgk1rUivF2dvGTXrSCbUtCLSjri+CUD8zTrSMyiQixDffYG4m3XUNt8zKJCLxCiuG7xqm+9ZFMhFYhRXs06cP77qjeY/OZ/quuoWy6rrqiObHEeBXCRmcTTrqG0+vyYOnsi0pdOag3l1XTXTlk5j4uCJkaSvQC7SC8XZNt8blQ8vp2pqFdOWTuO66uuYtnQaVVOrKB9eHkn6CuQivVCcI7pC77zRWj68nIqyCm5cfiMVZRWRBXFQIBfpleLsctlbb7RW11WzsGYh1555LQtrFh7UZp4N/SBIRPKqN/4Iqrqums8vnkbh76vY+nQ5x5xczTufmcZ//0tmzSv6QZCIdAu98UZr5e9X8t7tVWxZUY47bFlRznu3V1H5+5WRpK9ALiJ51RtvtD61YBbvPt+y5v3u8+U8tWBWJOkrkItIXsV9ozUOuf4WokAuInkV943WOHrL5PpbSFaB3MzOM7O1ZrbfzA5qgBcRaUscP4KKs7dMrr+FZFsjXwN8AVgeQVlERHImzmEJcv0tpG82O7v7OgAzi6Y0IiI5EndvmQsuyN03j7y1kZvZTDOrMbOaxsbGfGUrIgL07N4ynQZyM/ujma1p4/H5TDJy90p3L3P3sqKioq6XWESkC3pyb5lOm1bc/R/yURARkVxqataYMydoThk2LAji+ZpEJJeyaiMXEUmSXLZTxynb7of/ZGYNwKnA783s4WiKJSIi6cq218o9wD0RlUVERLpAv+wUEUk4BXIRkYRTIBcRSbhYJpYws0agjaHl0zII2BZhcaKicmVG5cqMypWZ7louyK5sxe5+0A9xYgnk2TCzmrZmyIibypUZlSszKldmumu5IDdlU9OKiEjCKZCLiCRcEgN5ZdwFaIfKlRmVKzMqV2a6a7kgB2VLXBu5iIi0lMQauYiIpFAgFxFJuG4VyM3sHDN7wcxeNrOr21hvZnZLuH61mY1Pd98cl+uCsDyrzezPZlaasq7ezP5mZrVmVpPnck0yszfDvGvN7Lp0981xub6VUqY1ZrbPzAaE63Jyvszsl2a21czWtLM+rmurs3LFdW11Vq64rq3OypX3aytM+3gzqzazdRbMY3x5G9vk7hpz927xAAqAV4ARwCHAs8CoVttMBh4CDDgFWJHuvjku12nAB8Pnn24qV/i6HhgU0/maBDzQlX1zWa5W208BHsvD+ToTGA+saWd93q+tNMuV92srzXLl/dpKp1xxXFth2scB48PnRwAv5jN+daca+cnAy+7+qru/B9wJtJ6F6PPArz3wF+BoMzsuzX1zVi53/7O7vxG+/AswNKK8sypXjvaNOu0vAXdElHe73H058PcONonj2uq0XDFdW+mcr/bEer5aycu1BeDum939mfD5W8A6YEirzXJ2jXWnQD4E2JjyuoGDT0R726Szby7LlWoGwX/dJg4sM7NVZjYzojJlUq5TzexZM3vIzEZnuG8uy4WZFQLnAHenLM7V+epMHNdWpvJ1baUr39dW2uK8tsysBBgHrGi1KmfXWHeaIcjaWNa6b2R726Szb1elnbaZlRN82D6esvh0d99kZscAj5jZ82GtIh/leoZgbIZdZjYZuBf432num8tyNZkCPOnuqTWsXJ2vzsRxbaUtz9dWOuK4tjIRy7VlZocT/PO4wt13tl7dxi6RXGPdqUbeAByf8noosCnNbdLZN5flwszGAD8HPu/u25uWu/um8O9Wgkk4Ts5Xudx9p7vvCp8/CPQzs0Hp7JvLcqX4Iq2++ubwfHUmjmsrLTFcW52K6drKRN6vLTPrRxDEl7j779rYJHfXWC4a/rt4s6Av8CownAMN/qNbbfMZWt4seDrdfXNcrmHAy8BprZYfBhyR8vzPwDl5LNexHPjR18nAhvDcxXq+wu2OImjrPCwf5ytMs4T2b97l/dpKs1x5v7bSLFfer610yhXjtWXAr4EfdrBNzq6xyE5uRCdjMsHd3leAOeGyS4BLUk7WreH6vwFlHe2bx3L9HHgDqA0fNeHyEeGb8iywNoZyXRrm+yzBjbLTOto3X+UKX08H7my1X87OF0HtbDOwl6AGNKObXFudlSuua6uzcsV1bXVYrjiurTD9jxM0h6xOea8m5+sa00/0RUQSrju1kYuISBcokIuIJJwCuYhIwimQi4gknAK5iEjCKZCLiCScArmISML9fz8ITEuLs2DgAAAAAElFTkSuQmCC\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", "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": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABkC0lEQVR4nO2deXwV1dn4v8+9IRC2sC8hJGERgSQk7AiVpaJVUKi7SFFb69bat9WfTS20aq1QTVv1rbUufd2luNQFFNygIBUUQYQsrGEJhLAbAiGEJPee3x8zk0wuN8lN7jL35s7385lP7pw558wzZybzzDnnOc8jSilsbGxsbKIXh9UC2NjY2NhYi60IbGxsbKIcWxHY2NjYRDm2IrCxsbGJcmxFYGNjYxPl2IrAxsbGJsqxFYGN34jIbBH51Go5DEQkTkQ+EJFSEXnbanlChYi8LCKPBKiuh0Tk9UDUFWhE5CMRubkJ+e8QkSd9yNdaRLaJSA+/BIxAbEUQRojIjSKyQUTKROSg/sB/z2q5GkMptVApdYnVcpi4BugJdFVKXestg4gMFZElurI4JSIrRWS86XiKiCj9XpSJyGER+VBELvaoZ6+InDHlKxORvwf38kBEbhGRL4J9nnBEKXWZUuoVX/KKSCzwO+DPPtR7FngR+I1/EkYetiIIE0TkXuBJYAHaSywJ+Acw00KxGkVEYqyWwQvJwA6lVLW3gyIyAFgD5AL9gATgPeBTEbnAI3snpVR7IAP4DHhPRG7xyHOFUqq9abs7gNdi4x8zgW1KqQM+5v8XcLOItA6iTOGHUsreLN6AeKAMuLaBPK3RFEWxvj0JtNaPTQaKgCzgCHAQ+CEwDdgBfAfMNdX1EPBv4E3gFLARyDAdvx/YpR/bAlxpOnYL2kv0Cb3eR/S0L/Tjoh87ApQCOUCa6TpfBY4ChWhfag5TvV8AfwFKgD3AZQ20xxBgFXACyAdm6Ol/ACqBKr1Nb/VS9jVgmZf0Z4DV+u8UQAExHnnuAw6b5N4LTPXxPj8EvA28rrdtLjAI+K3eXvuBSzyeixf0+3lAb2unfu0VgEu/xhN6/peBp4Glev3rgAGm+sYD6/X7sh4YbzrWD/hcL/cZ8Hfgdf1YG13m43p7rwd6NuW++CKfRz31nlOv/6e+PDdoX/i/M+1fD+wGOur7lwGHgO6mPDuBSVa/F0K5WS6AvSmAS4Fqz5eOR56Hga+AHkB3YC3wR/3YZL38A0Ar4Da0l+2/gA5Aqv7i6K/nfwjtRXmNnv8+/R+olX78WrSvZIf+j3Ma6K0fu0U/1y+AGCCOuorgB8A3QCc0pTDEVPZVYLEuUwqakrrVVG+VLrsTuAtN4YmXtmgFFABzgVjg+/qL5XzT9b3eQFseAn7sJX0K2su1LfUrgv56+hB9fy9NUwQVehvF6O2xB5hnum97TPnfB54D2un3/WvgDlN7feFR/8toynmMXv9C4A39WBe0F+Uc/dgsfb+rfvxL4HG0D46JensaiuAO4AO9XZzASPQXaRPvS73yeamr3nNyriKo97lBUyDXetS9UJelq573co/jS4D/sfq9EMrNcgHsTQHMBg41kmcXMM20/wNgr/57MnAGcOr7HfSX1VhT/m+AH+q/HwK+Mh1zoH11XljPuTcBM/XftwD7PI7XvJT0f/4dwDj0r2Y93QmcBYaa0u4AVpnqKDAda6tfQy8v8lyI9jI3178IeMh0fQ0pgmrgUi/pg/Vz9qF+RdBGT5+g7+9F/yo3bbfVc96HgM9M+1foZT3vWye04cGzQJwp/yxgpWebm46/DPyfaX8a2rAIaArga4/8X+r1JOlt0s507F/UKoKfoH14DGvkGW3svtQrn5e66j0n5yqCep8btK/7Sz3KdwL2ofXInvNS/0LggUD8b0fKZs8RhAfHgW6NjLcnoA2nGBTqaTV1KKVc+u8z+t/DpuNngPam/f3GD6WUG21oKQFARG4SkU0ickJETgBpQDdvZT1RSv0HbVjhaeCwiDwvIh318rFerqGPaf+QqZ5y/adZZoMEYL8ud311NcQxoLeX9N6AG+1LuT6Mc3xnSvuhUqqTaftnA+U978kxL/etPdo8RyvgoOk+PIfWM2iIQ6bf5dS2n+fzA7VtlgCUKKVOexwzeA34BHhDRIpFJFtEWnk5ty/3pT75PPH1nHXq9PLclKApWEx5TqAN0aUBf/VSXwc0hR412IogPPgSbcjghw3kKUZ7ORgk6WnNpa/xQ0QcQCJQLCLJwD+Bu9GGDToBeWjDPAaqoYqVUn9TSo1EG5IaBPwa7eVb5eUafJ3EM1MM9NXlbk5dy9GGvzy5DvjS9DLxxpVo4/nbfTxXc9mP1iPoZlIwHZVSqfrxBu+BFzyfH6hts4NAZxFp53FMO5FSVUqpPyilhqLNM1wO3FTPOfy5LzU04ZyNkYP2DNYgIploPY5FwN+8lBkCbG7GuSIWWxGEAUqpUrTx/adF5Ici0lZEWonIZSKSrWdbBPxORLqLSDc9vz923iNF5Cq9F/IrtJfOV2jj0QptjgER+THal5NPiMhoERmrf72dRp/U1L963wLmi0gHXeHc28xrWKfXnaW302S0YZY3fCz/B2C8iMwXkS66PL9Ae9F4NR0UkZ4icjfwIPBbj6/egKOUOgh8CvxVRDqKiENEBojIJD3LYSBRN4/0hWXAIN1EOUZErgeGAh8qpQqBDcAfRCRWN1m+wigoIlNEJF1EnMBJNIXu8nIOf+9LDU04Z2MsA4w2Q0SMSei5wI+BPiLyM9PxPmjzKV8141wRi60IwgSl1ONoL8bfob2E96N9lb+vZ3kE7Z81B21sc6Oe1lwWo00EGxOIV+lfYVvQustfor1s0tGshHylI1qPogRtWOA4mkUHaBPMp9GsNr5AG4d+samCK6UqgRloFh/H0Mxsb1JKbfOx/E7ge2gmoXvRvoivBn6glPK81hMichqtzaehTTx6yvyBxzqC95p6TfVwE9pw2ha09vw3tUNa/0GzyjkkIscaq0gpdRztq/r/od2TLLRJUqPsjcBYtCGvB9Emsg166ec+CWxFsy46R4H7e1888OmcPvABMFhEjGHUPwFFSqlnlLZu4EfAIyJynn78RuAV/VjUYMys20QRIvIQMFAp9SOrZbGxCTYicjuakcKvGsnXGm1IaKJS6kgoZAsXwnExkI2NjU3AUEo972O+s2iWY1GHPTRkY2NjE+XYQ0M2NjY2UY7dI7CxsbGJciJyjqBbt24qJSWlWWVPnz5Nu3btGs9oIeEuoy2f/4S7jLZ8/hOOMn7zzTfHlFLdzzlg9dLm5mwjR45UzWXlypXNLhsqwl1GWz7/CXcZbfn8JxxlBDYo28WEjY2NjY0ntiKwsbGxiXJsRWBjY2MT5UTkZLGNjU10UVVVRVFRERUVFQDEx8ezdetWi6VqGCtlbNOmDYmJibRqVZ/D1rrYisDGxibsKSoqokOHDqSkpCAinDp1ig4dOjRe0EKsklEpxfHjxykqKqJfv34+lQnI0JCIvCgiR0Qkr57jIiJ/E5ECEckRkRGmY5eKyHb92P2BkMcysrPJnTmPopgU3OKgKCaF3JnzIDu78bI2gSc7m+XzVpKSAg4HpKTA8nkr7fsRgVRUVNC1a1dEpPHMUY6I0LVr15reky8Eao7gZbRwi/VxGXCevt2OFhsW3cXs0/rxocAsERkaIJlCTu6aUtKWLCDRVYgDRaKrkLQlC8hdU2q1aFHJ8tLRjF0wk28Ku1CthG8KuzB2wUyWl462WrTopc7HklAq8ZRKPG6RRj+cbCXgO01tq4AoAqXUaupGbPJkJvCqbsr6FdBJRHqjxS4tUErtVpoL2zf0vJGD6cFOW7IAQXPmbzjuEKDz0oXWyRfFLHmumFgq6EIJAnShhDjKyFhwrd1js4i6H0sQz0niOYkD7A8nCwnVHEEf6oY3LNLTvKWP9VaB7kr2doCePXuyatWqZglSVlbW7LLeUIt3MXnt8whavLwdaOGUYtCUgQB9XIXsdyZRMO4yktJg/6xZIZUx0ISzfH0XLaJ6czX719/E/7r3U40We9KJFpR3LYphHK958fRZsoCVx25HxqwKqZzh3IYQPPkGfPgaAJVoX6FGpJlW+r4AnT58jVWrLq5TLj4+nlOnTtXsu1yuOvuhoFOnTqSmplJdXU1ycjLPP/88nTp1orCwkOuuu45169YB8PLLL/PCCy/w/vvv15T929/+xu9+9zv27NlD165dKSwsZPTo0Zx3nhYGYfTo0Tz55JMBlbeiosL3e+htlVlzNrRg33n1HFsKfM+0vwIYiRYu0BzMeg7wVGPnCqeVxfudyUqBOguqSt/coFygToNaA0rpmxtUzoy5IZcx0ISzfDkz5iq33t4nQM0HVaG3fTWoF/T74jbdl/3O5JDLGc5tqFTw5HMh6j/6/4nxP2HcC1fNXzmn3JYtW+rsnzx5ssHzvP66UsnJSolof19/3X/Z27VrV/P7pptuUo888ohSSqk9e/ao1NRUpZRSr776qkpPT1dHjx6tkXHfvn3qkksuUUlJSero0aPnlAkWnm2mlPUri4swxchFj4/bQHrEkODax3a0r5kYfRN9vwQYTu1Xjz1MFHw6L12IoPUC2qEFS45Fa3snWqDaNvq+MXyX4NpngaRRhGn4dAmKJ6kdkxZq70UFWigyQfk1bLdwIdx+OxQWapqlsFDbXxjAf70LLriAAwfqhmJ+6623ePTRR/n000/p1q1bTfo999xDdnZ2WM9xhEoRLAFu0q2HxgGlSovJuh44T0T66bFXb9DzRgy7HInMo25kd9Ae7D5oLx0n9ksnVBjtq9CUcivOvTeCNkxUre8XO5OwCR7GvMBZVyF5aLEiPV88gvYy+C1Qhn/zBfPmQXl53bTyci09ELhcLlasWMGMGTNq0goLC7n77rv59NNP6dWrV036kiVL6NOnDxkZGefUs2fPHoYPH86kSZP473//GxjhmkmgzEcXocW4PV9EikTkVhG5U0Tu1LMsQ4tTW4AWz/ZnAEqparS4vJ+gfQy8pZTKD4RMQcX0hfOKez8OtJe9gQJO0hHQHnAX2pgo2C+dYFPsTGIXWnubI20ooLJNx5p5m5NoUdX3AAmuQnviOIh0XroQF5CAFjHeMAtUQCkdKUW7L5OAt4G2+vHm9qD31fOtVV+6r5w5c4bMzEy6du3Kd999x8UX185jdO/enaSkJN56662atPLycubPn8/DDz98Tl29e/dm3759fPvttzz++OPceOONnDx50j8B/SBQVkOzlFK9lVKtlFKJSqkXlFLPKqWe1Y8rpdTPlVIDlFLpSqkNprLLlFKD9GPzAyFPsDG+cKpdhbRBs39VaMMRRc5k8mbMpXzIyJoXUTVwFVqk8ZLps60ROkoomT6bucC71PYEFJA3Yy6xf5hH3oy5FDmT6QRciHbvBNtiJZgkuPaxjdoJYQOFEK9KiVelKP1uxeB/Dzqpnm+t+tJ9JS4ujk2bNlFYWEhlZSVPP/10zbG2bdvy0Ucf8eyzz7JQH4Pas2cPe/bsISMjg5SUFIqKihgxYgSHDh2idevWdO3aFYCRI0cyYMAAduzY4Z+AfmD7GmoGxjh0L7QvnO5oL5NiZzKJ1XtJXzyf3rdcWvPSaYWQIZ24r89Q0ifEWyp7S0cGnOGzVnGMdvTFjdQo5vQJ8ZCVRfri+SRW76XYmUwW8EdsU99gU+xM4rl60j1/Gz3oKi95fGX+fGjbtm5a27ZaeiCIj4/nb3/7G3/5y1+oqqqqSe/evTsff/wxc+fO5ZNPPiE1NZUjR46wd+9e9u7dS2JiIhs3bqRXr14cPXoUl0ubPdy9ezc7d+6kf//+gRGwGdiKoBkkuPZxBO3rxeGRXoPppeNQbub+4V6+PrqX//z2SduGPdCYhuoefuIJ7qqO5ezlc3A89miNYiYrq06RBNc+WgGtaeAe2gSEDRMv4d+cO1Rn7h2XTJ9dc7wcuBxtYVJzetCzZ8Pzz0NyMohof59/XksPFMOHDycjI4M33nijTnq/fv1YsmQJP/nJT1i/fn295VevXs2wYcPIyMjgmmuu4dlnn6VLly6BE7CJ2L6GmkGxM4kPXYX8xEt6Yj1l9myo4NeV5ayinO9Ta8Oey1zSgyxvS8cYqtsDZAD3qlLaNNK2xc4kEl2FCNrXp2Ht1dA9tGkeyyvymDloAod3FZHg2kexM4mS6bPr9I7TJ8STy1w6L11Igmsf7SSOR4eOJ7uZPejZswP74gdtbYWZDz74oOZ3Xl6td52MjAwOHDhwzjqHvXv31vy++uqrufrqqwMroB/YPYJmcHzajfwJKDSleX7heNJ56UL+B7gfeygi0BhDdYloVidxNN62xheoAO+hTVI2dg9tfMTk40nkJM+s28bkH/wPiQt+hkO5vffSPHrQv77zZhYWbyX56ftwOKCoCE4eOAmHDll2WS0ZWxE0gxMJ39G6Q3finMnnjkPXQ4JrH23RhiLEI93GPxJc+3CjPcy+DvOkT4ivmcPpCMwnttF7aOMbZh9Pa4nnHncZ05/6aZN8PJ3udA3dSg7yz33xVCuht2s/7Q4WcNLVtvHCNk3GHhpqBq+7XPz0d/eRqH/RJOpbQxhDEQrNiijGnB48UaOCYmcSR12FDKHuA91g22Zl6cNG8+lVXc2Rvn1p9dgcGDw4yNK2fJY8V8yFVNCeKkYDY6iiWk+f6uOE7QfPH+QBhEmU6ebZbgTFiaNVdOwTROGjFLtH0EQqKip45513uPHGG5tUzhiKcAHz0VYd20MRgaFk+mxeBv5lSmtK28bExHDjjTfy+uuvB0G66OPe4/Nordv9GAYVrani3uO+r+i69/g8rsBVZ32OA0XP6gP1lrFpPrYi8AWTVcriuDhSS85Q8vNnmmTxYwxFHHUmkwe8IF3soYgAkTa+I2/FdSTZ0Rs3wn5H36a1bXY2474p5ZX5j1Ltgztkm4ZJQhuSq0b78BGPdF/riNXLu03psTVLM20CiT005AOGVYoAZ4H7qdAWHzXF4sc0FHH1G2/w6quvct/iiFg/F/Z8c9FFxL/0Et/fuhURYdeqVUyePNnn8rlrSrn68xdYjOYKt59t0eUXZzv2IO7kYe4Afgf0N6c3sY5P0VxOZOrpytnqHJchNv5j9wh8wLBKUcCNwDT8s/iZPn06X3zxBaWl9irWQPDOO+9w1VVXNdupV+elC3EALwHG8iXboqv5xE2bQrEzlnep9SjpimlN3LQpTarDFdOakWj+aBSgRHDEWxee0ul0kpmZSVpaGldccQUnTpwANLPQtLS0mnz//Oc/GTFiBCUlJfz6179m8ODBDBs2jCuvvLKmzNdff01mZiaZmZlkZGTw3nvvWXBFtdiKwAc8HZmJR3pT6dChA5MmTWLp0qUBkS+aUUrVKILmYtzHQLg3sAEWLeKzW29hatu2tNJXdDlffgEWLWpSHc6XXyAhOZl+wFmHA0lJAV9X3y5cSJ0YpQFwPWq4mMjLy6NLly51XEwYvPbaazz11FN8+umndO7cmYsvvpi8vDxycnIYNGgQf/rTnwBIS0tjw4YNbNq0iY8//pg77riD6urqc+oLFbYi8AFjmbsR5MQzvclkZzO2CF6ffZu9yri56PM2/4npQ8XOnfQce3Wz29B8H6uoVQS2g8Dms+zECaY99RS43bB3b/NWd82eDXv3Mm3ePM60bw+6b55GCYEfal/dUF9yySXExGgj8OPGjaOoqAjQfBMZ6RUVFZa7qLYVgQ8YFj83Akf0NH8sfnLXlHLnpg+ppJxqO7ZxszDmbQrdB/kH0Ne9r9lt6Lm4bDG2RZc/VFdX8+mnn3LppQ2FMfedadOmcebMGd8LBNkPdVPcUJt58cUXueyyy2r2161bR2pqKunp6Tz77LM1isEKbEXgA+kT4vn4op+xGgfdwacFZA3ReelCugEfUjtbb49JNw1j3uZGYLqe1tw2NC8uKwdeoK1t0eUHX375Jf369SMhISEg9Y0dOxaXy0VlpY8WQ0HyQ91UN9Rm5s+fT0xMDLNNPaOxY8eSn5/P+vXr+dOf/kRFRYVf8vmDrQh8ISuLAzcM5+Ibb8CpVL2OzHzFGHuOxXZ41lwSXPuooO6cjZHeZEzuDS4uKuLLLm1Ife+Pzb6/0c6yZcuYNm1awOpzOp3ExcX5blwRJD/UTXVDbfDKK6/w4YcfsnDhQq9DQEOGDKFdu3Z1/BWFGlsR+Mjy5cuZOnVqQOoyxp6NVcae6TaNU+xM4lvqtp+R7g99+vShR48ebNq0ya96og7TWptljz7K6D+9ELh5r+xsWlW5KSncj9qwgcoNOagD39XvdyjIfqh9dUMN8PHHH/PYY4+xZMkS2ppk2rNnT83kcGFhIdu3byclJSUg8jWHQEUou1REtotIgYjc7+X4r0Vkk77liYhLRLrox/aKSK5+bMO5tVuP2+1mxYoVXHTRRQGpz7zK+I/Aaewx6aZSMn02i4HXTGmBasOLLrqI5cuX+11PNGHM2cS6CrkGuNx9KGDzXrlrSmlfeQbDVCOWSjqUHaW8zOW9QAj8UPvqhvruu+/m1KlTXHzxxWRmZnLnnVrQxi+++IKMjAwyMzO58sor+cc//lEnznGo8Xt2QkScaIGeLkYLRr9eRJYopbYYeZRSfwb+rOe/ArhHKfWdqZopSqlj/soSLHJzc+ncuTNJ/oY40jG73P3cVchARw8yL/+pPSbdBNInxPPxf/twf2kVbvfROq6NV/lZ99SpU3nmmWfIsoeGfMaYs+mCFqzJaUrXnKr4V/epuy5koClNgJgTx9Eig3shCH6om+OGuqCgwGtdc+bMYc6cOQGVzx8CMU09BihQSu0GEJE3gJnAlnryzwKaYFBsPYEcFgLqrDKe9OCDbK2qYs4Ce5VxUyi/+252/uEPzDh1FEfbtnUd/61a5VfdkydPZs6cOZw9e5bWrVv7KWl0YMzNCHXjdwdi3ivBtY/tet3mEfZWtruJgBEIRdAHbWW+QREw1ltGEWkLXIoWsN5AAZ+KiAKeU0o9X0/Z24HbAXr27MmqZv6zl5WV+VS276JF7MuDgV99xHL3fn4oXVm52UVSGuyfNatZ5/ZGx44deemll7jkkkuaLKNVhIN8GzdupF+/fnz99dfnHPNHvpr7XlbJmjZtOM/Rl4JxlwX8vodDGzZEU+Ub4Eikr3s/lUAbal/YBxyJ7PLzOgc4an3IGma+AFXEctYj+Es44XK5zglOE0oqKip8v4dKKb824Frg/0z7c4Cn6sl7PfCBR1qC/rcHsBmY2Ng5R44cqZrLypUrfcqXM2OucoOqBPUgqFJQblA5M+Y2+9zeOH36tGrXrp06ffp0k2W0inCQ74EHHlD333+/12P+yGfc93+A+kRbjhSU+x4ObdgQTZUvZ8ZcdQTURaBcAW63nBlz1ZaPPlJV69erA+vXK7e+nd5Z5HfdweTkyZOWnn/Lli3npAEblJd3aiAmi4uodSkCWg+9uJ68N+AxLKSUKtb/HkFbzzMmADL5jTHm6UBznNWR4Nj6t23blmHDhvHVV18FtN6WzurVq5k4cWLA6zXu++3A9/U0e41H46RPiOeNMddSpfcH/F1r41l3ZVw8bmI5DpyiFafad6dte2ejZW18IxCKYD1wnoj0E5FYtJf9Es9MIhIPTEJbuGmktRORDsZv4BLAOmNaE+Yxzxgv6YFk0qRJrF69OuD1tlTOnj3L+vXrmTBhQsDrNu6vA+2+236HfCQriwNT+vP9h+6vPxylH3XH9uhE7KhhtO/albPJCUifLlDPCl6bpuO3IlBKVaON+X8CbAXeUkrli8idInKnKeuVwKdKqdOmtJ7AFyKyGfgaWKqU+thfmQKBYY9eRYD8CzXAxIkT+fzzzwNeb0tl/fr1DB48mI4dOwa8buP+CnV9S9lrPBpn7dq1jB8/PqjnaN++/TnWOzb+E5B1BEqpZUqpQUqpAUqp+Xras0qpZ015XlZK3eBRbrdSKkPfUo2y4UDJ9Nm4gauAk3paUGz9s7Pp9NQKvl71OWf0oChq3iLbAZ0npgVLn194IcM37giKoz5jjQfAUuyg9r5SWVnJxo0bGTvWq51IwGjXrp1liqA+N9S+MnnyZDZsCMulUvbK4vpInxDPZ1N/zgacdMB//0L1kbumlHEf/ZXrUBwGEl2FTF77vO2AzgNjwVKiq5DewBx1KiiO+sx+h74D/iW23yFf+Pbbbxk4cGBQemlm4uLiqK6ubthlc3Y2rFxZN23lSr8/GnxxQx2p2IqgPrKyKL1tIhfMvDwg/oXqw5icfA47KEpDmIMD3QRcSJDayeR3aNyWLeQm9wjKfW9phGJYCEBEaNeuXcMO2kaPhuuuq1UGK1dq+6NHB0wOsxvqr7/+mvHjxzN8+HDGjx/P9u3bAc1J3Q033MCwYcO4/vrr63hQveuuuxg1ahSpqak8+OCDNekpKSkcO6atrd2wYUOTIu35gx2qsgHWrVsX9K6uOSiK7YCufsztEezJe4Pzzz+fkpISjhw5Qo8ePYJ2npbA2rVr+eEPfxiSc7Vv375ht9RTpsBbb2kv/7vugmee0fan+B4hrSEMN9S33norAIMHD2b16tXExMSwfPly5s6dyzvvvMMLL7xA27ZtycnJIScnhxEjRtTUMX/+fLp06YLL5eKiiy4iJyeHYcOGBUS+5mD3CBogFIrAPAlpO6CrH6M9XPrmmR4MHA4Ho0eP9rpozaYWpVTIegTggyIA7aV/113wxz9qfwOgBOpzQ11aWsq1115LWloa99xzD/n5+QCsWbOGH/3oRwAMGzaszov+rbfeYsSIEQwfPpz8/Hy2bKnPEUNosBVBPVRVVbFp0yZGjRoV1POYg6IsAI5hT056w2in/wW+1NNC0U5jx45l3bp1QT1HxKJP4H8V0xd3cTExAycHP9LeoUM4Dpdy9swZ3Lon0vKCA+d6Il25UusJ/P732l/POYNmUJ8b6t///vdMmTKFvLw8PvjggzrDVt7cTu/Zs4e//OUvrFixgpycHKZPn15TJiYmBrdbs1ULZXwCWxHUQ25uLsnJyUGf/DImJw84k/kv8KGjO6vG325PTnpgtNNCYnERvMl7T2xFUD/GBP4p9wH+jH9R4nylvMxFu9LDdEbrQcdSSdyJg3U9kRpzAm+9BQ8/XDtMFABlAOe6oS4tLaVPH8353csvv1yTb8KECTWxCYy4xQAnT56kXbt2xMfHc/jwYT766KOaMikpKXzzzTcAvPPOOwGR1xfsOYJ6CMWwEFDHAd3ouXPZGxtLyuTJEKJJooghK4vBVVVs7/QkIw8fp3379rVO5oLImDFjWL9+PW63G4fD/m4yY0zgT8ZbpL3gWILHnDiOAMnU+hw6xxPp+vV15wSMOYP16wM2T2B2Q52VlcXNN9/M448/zve///2aPLfeeiv/8z//w7Bhw8jMzGTMGM1pQkZGBsOHDyc1NZX+/fvXWRj54IMPcuutt7JgwYLQvH90bEVQD+vWreOCCy4I6TlHjRrFiy++GDJLgUgjPz+flJQU2rdvH7Jz9uzZk/j4eHbu3Mn5558fsvNGAuZV2KEydDA8jjboidSbhdeUKX4rgYbcUO/YsaPm9x//+EdAG0ryjFdgYO45mLnwwgvr1BUq7E+cevj6669DqpFBUwQbNmwwnPHZeLBhw4agz9l4wx4e8k6xM4lqoJJaVxxGerCoIrbmt6on3abp2IrAjD75tcWZxP6tW+k0cmbwJ79M9O3bF7fbXWNHbFMXSxRBdjZJ+cdZccvPcYuDopiUkD4T4UzJ9NlsA26h9us82BP41Z261iiAYjQLMqWn2zQfWxGYMCa/Trv38xcgJQSTX2ZEhFGjRtUsSLGpixWKIHdNKTPzV+BUZThQJLoKQ/pMhDPpE+JZMvwKKqQt7gB7HK2Ptu2dnOnUmypiOQmcpBVnOvW2PZH6iT1HYMKY/MoEjKUfwZ788sRWBN45e/YsW7ZsISMjI6Tn7bx0IYPQngfDzDfUz0TYkpXFsUOH+N4N38ORlVU3Slyw6NWLtsCpUx1p+913VMXF0dle7Oc3do/AhHnyK9Dh9nxl9OjRtiLwQk5ODueddx5t27YN6XkTXPtoDcRSd3LSXvmtsXHjRoYPH27Judu2bUt5ebkl525p2IrAhDHJFcrJL09GjhzJ9u3b7QljD6yaKDavaA62O/JIw+128+2339qKoAVgKwITJdNnUwJcTq0iCPUq34SEBGJjYyksLAzZOSMBqxSBsaL5fcAwFrRXfmvs2bOHTp060a1bN0vOHxcXR0VFRc1K3GBjuKHOyMhgxIgRrF27FoC9e/cSFxfH8OHDGTJkCGPGjOGVV16pKffyyy/TvXt3MjMzyczM5KabbgqJvE3BniMwkT4hnpeO/4jSNW8DlRQ5kyiZPjt0q3yzs8ldU8rwY2V83a8fMc7k2vNHo/dLvT06L13IBlchV768jNyPi0LaHukT4sllLo4PXyDbfZiR5nsS5Vg5LASaL6jWrVtz5swZ2rVrF/TzGS4mAD755BN++9vf1gSUGjBgAN9++y0Au3fv5qqrrqK8vJy77roLgOuvv56///3vQZexuQREEYjIpWhuYJxogewf9Tg+GS1E5R496V2l1MO+lA0pWVmcavM3Rqe3x/HMM6GZ/DJhWC3dgBYjOdFVSJ8lC8hlrr76OLow2uMsMBO4xH2IVqFuD33ld2LJffwkKYmE0t0k2iuMAU0RmD1qWoExPBQKRWDm5MmTdO7c2eux/v378/jjj3PPPffUKIJwx+8nWkScwNPAZcBQYJaIDPWS9b9KqUx9e7iJZUPGpk2byMzMtOTchtXSzcDFelo0xyYw2qMV8AC1E7ZWtEfnzp3p2rUru3fvDvm5w5VwUgShwPA+OnjwYH7605/y+9//vt68I0aMqLNC+M0336wZGnrppZdCIW6TCESPYAxQoJTaDSAib6B9wPniV9WfskFh06ZN3HHHHZac22y1ZARhEaLXQsW4bsE6Ky4zGRkZbNq0iYEDB1py/nBCKWXpRHEwnEE2ZqBhHhr68ssvuemmm8jLy/OprmgYGuoD7DftFwHefDNcoAepLwbuU0rlN6EsInI7cDto/l9WrVrVLGHLysrqLVtdXc2WLVsoKSlpdv3+MMCRSF+31hxVaC8/AQ44EtllgTz10VAbBhKjPaqpG7insfYIlnydOnVi8eLFAZkcDVUbNpfG5Dt69CjV1dXs2LGDnTt3Bl2e+Ph4Tp06VbNfUlKC0+nE7XZTUFDAeeed59Xlc1Mw199YnrS0NI4ePcqePXsoLy/H7XbXKb9mzRoGDRrEqVOnqKiooLKy0qf6A0lFRYXvz5hSyq8NuBZtbN/YnwM85ZGnI9Be/z0N2OlrWW/byJEjVXNZuXJlvcdycnLU4MGDm123v+TMmKvcoBSof4L6EpQbVM6MuZbJ5I2G2jCQGO1xO6gCvV18aY9gyffuu++q6dOnB6SuULVhc2lMvsWLF6tLL700NMIopbZs2VJn/+TJkzW/c3JyVHl5edBlaNeuXc3vrVu3qq5du6rq6mq1Z88elZqaWnNsz549avjw4eof//iHUkqpl156Sf385z8PunyeeLaZUkoBG5SXd2ogegRFQF/TfiLaV79Z2Zw0/V4mIv8QkW6+lA0lmzZtCvnKVTOGhUqnD18jz72fPdKJdlf8LGotVNInxJOjfsu/PniUP6IosthiJzMzk82bN1ty7rBBt+T6/IOnyVSlFMWkWGfZdugQ5WUu2pyt4nR+Pk5iqe7UVXM30atXwE9nzBGA9gH9yiuv4HRqg5a7du1i+PDhVFRU0KFDB37xi19wzTXXBFyGYBEIRbAeOE9E+gEHgBuAG80ZRKQXcFgppURkDFov/zhworGyocTKiWKgxkJl1aqLydizhxUrVpD+ehS7McjKouN1e+n4zSv00AOFh9KKy5OUlBROnTrFsWPHLLOdtxrDksuJ9o9qpWVbeZmLuBMH6YI2jBpLJa1OHKSc3gRj/bnL5fKanpKS4jV0pjEUdMstt3DLLbcEQaLA4bfVkFKqGrgb+ATYCryllMoXkTtF5E492zVAnj5H8DfgBr2n4rWsvzI1F8sVgYlhw4bVRDSKZqwO6m1GRBg2bFhU9woMS64/AGl6mlWWXEaQmi6A0UesDVJj0xQCso5AKbUMWOaR9qzp998Br1Pm3spagVKKzZs3h40iGDp0KDt37qSyspLY2Oj1tZ6Tk2PpcJ0nmZmZbNq0iYsuushqUSwhwbWPM2gmvVb7XjIHozFb2dUJUmPjE/bKGJ0DBw4QExNDryCMLTaHuLg4+vXrx7Zt26wWxVLCqUcA9jxBsTOJHWjxgj3TQ405GI27nnQb37AVgY7VE8XesIeHwlMRGLbk0UjJ9NmsAR43pVnle8kIUiNoE46nsYPUNBfb15DJCiLDaisID6JdEZSXl1NYWBg+sYKzs3GvPs6O3DzKRfguCn1BpU+I55GcEQwt3INbnaA41P64TLRt76Sc3sScOM4ZKqnEiaNTDztITTOIekVgWEHEAtcRXv59hg0bxlNPPWWxFNaRn5/P+eefT6tWrawWBdCelZFLs7kGKCG8npWQkZVF8QcfcMcL/8bx/e+H3B9XHfQgNdCHuKNHKSsro22/PlZJE9FE/dCQYQXxO8AYgAgX/z7R3iMIt4li41l5Aeitp4XLsxIqlFLk5uaSnh5eqi8uLs6rCWcgMdxQp6amkpGRweOPP17HBfYXX3zBmDFjGDx4MIMHD67jU+ihhx6iT58+ZGZmkpaWxpIlSwDYt28fU6ZMYfjw4QwbNoxly5adc77MzExmzJhRk37rrbeSkZHBsGHDuOaaaygrK/P72qK+R5Dg2kcl4WEF4Unfvn2pqKjgyJEj9IjCcHzhNj9gPBOtqPsFFQ7PSqgoKiqiTZs2dO/e3WpR6mDEJlBK+e1qoqFzGPNDR44c4cYbb6S0tJQ//OEPHDp0iBtvvJH333+fESNGcOzYMS6++GIGDBjA9OnTAbjnnnu477772Lp1KxdeeCFHjhzhkUce4brrruOuu+5iy5YtTJs2jb17955zPjNPPPFEja+le++9l7///e/cf//9fl1b1PcIip1J7Ebz7eOZbjWG3Xpubq7VolhCuCkC45lQ1LWaCYdnJVSEY28AtK/nVq1acfbs2drEL7+EP/1J+xtgevTowfPPP8/f//53lFI8/fTT3HLLLTXeWLt168bDDz/Mo4+e61V/yJAhxMTEcOzYMUSEkyc1xwulpaUkJCQ0em5DCSilOHPmTEAUX9QrgpLps/kKMN+ucIpAFa3DQ0qpsFMERrSyI8Cf9bRwelZCQV5eXlgqAtC+oGtcUn/5JVx0Efz+99rfICiD/v3743a7OXLkCPn5+YwcObLO8eHDh5Off+762HXr1uFwOOjevTsPPfQQr7/+OomJiUybNq3OnGBFRQWjRo1i3LhxvP/++3Xq+PGPf0yvXr3Ytm0bv/jFL/y+lqhXBOkT4lk1YCxnpRNuhCJnMnkz5lrv3yc7m9yZ8+j79L/YfO+9FMWkkDtzHmRnWytXsNGve0NMX2K++46qPmPD5rrTJ8STN2Mu1Y4kFgB5jsTweFZCSLj2CMBjnmDVKqisBJdL+xskT69Kdzdd35CUOe2JJ54gMzOT++67jzfffBMRYdGiRdxyyy0UFRWxbNky5syZUzPvsG/fPjZs2MC//vUvfvWrX7Fr166aul566SWKi4sZMmQIb775pt/XEfWKgKwsDp/XmQmLX8Wh3CRW7yV98XzLzQENa6bp6jtGoVmopC1ZQO6aUkvlCjbGdVe5D7CAMLvurCzSF8+nr6uQ1LFjOfH5orB4VkJJxCiCyZMhNhacTu3v5MkBP9/u3btxOp306NGD1NRUNmzYUOf4pk2bGDq0Ns7WPffcw6ZNm/jvf//LhRdeCMALL7zAddddB8AFF1xARUUFx44dA6gZJurfvz+TJ0+uCYVp4HQ6uf7663nnnXf8vhZbEaB1d9PS0hrPGEIMC5WhwJ3ULp9v6RYqxnWPBm7R08LxutPS0uoNStJSqaqqYvv27XVebuFEHUVwwQWwYgX88Y/a3wsuCOi5jh49yp133sndd9+NiPDzn/+cl19+uWZy9/jx4zzwwANkNfKRkJSUxIoVKwDYunUrFRUVdO/enZKSkpr5jmPHjrFmzRqGDh2KUoqCggJA64V88MEHDB482O/riXqrodLSUk6cOEFycrLVotTBHJ0rBk0RmNNbKuYobeEQlaw+olER7Ny5k8TERNq2DYZvT/9p3bo1VVVVuFwuzT30BRcEVAEYbqirqqqIiYlhzpw53HvvvQD07t2b119/ndtuu41Tp06hlOLOO+/kiiuuaLDOv/71r9x222088cQTiAgvv/wyIsLWrVu54447cDgcuN1u7r//foYOHYrb7ebmm2/m5MmTKKXIyMjgmWee8fvaol4R5OfnM3ToUBxhFpC82JlEoqsQqLVQiTHSLZMq+BjXXUVtyM6adAvl8iQtLe2cCbyWTjgPCwE4HA5at25NRUVFUILZ1+eG2mDixImsX7++Zt8ckeyhhx7yWmbo0KGsWbPmnPTx48d7tRZ0OBxe8/tLeL39LCAch4Wg1kIFYDnwKtFhoVIyfTbVwGVQ40MyHK87NTWVvLy8RuPctiTCXRFAaBaWtUSivkcQrorAiFbWeelCjrgKeV/aMvqKX7V4C5X0CfEsLbuTgv/8k1a4KbLQl01D9OrVC6UUR44coWfPnlaLEzx0X1ydly4kx1XIdEc3cjdVh59/JT1aWesTpZz57jsq9xYHNVpZSyMgikBELgX+F21Y9/+UUo96HJ8N/EbfLQPuUkpt1o/tBU4BLqBaKTUqEDL5Sl5eXqPjeJagRyuD+aR/+y1/vukmzUKlpZOVRfX5ixnWZh+OpUut9WXTACJSM0/QkhWBYcUlwIXApe5jJFnkX6mhVcNGtDI3cIbgRysLd5raU/V7aEhEnMDTaL35ocAsEfE0K9gDTFJKDQP+CDzvcXyKUioz1EoANEWQmpoa6tM2icGDB1NQUEBVlef655ZJfn5+2N8TiI4JY8OKSwH3AElYY8XVpk0bjh8/Xu8LzohW1g4wgohGa7QypRTHjx+nTZs2PpcJRI9gDFCglNoNICJvADOBLSbB1pryf0WYfOQdOXKE6upqevfu3XhmC4mLiyMxMZGCggKGDBlitThBJy8vjx/84AdWi9EoaWlpbNy40WoxgorZWiumnvRQkJiYSFFREUePHgW0VbfmF506Voy3voICZKs13ms9ZQwlbdq0ITHR99dsIBRBH2C/ab8IGNtA/luBj0z7CvhURBTwnFLKs7cAgIjcDtwO0LNnT1Y1c6VgWVlZTdmNGzfSt29fPv/882bVFSzMMhr07NmTt956i0mTJlkjlAlv8gWSdevWMWnSpIDc42BSWVnJ2rVrm3WuUMnYXAz5BjgS6evejwvtC9sw6T3gSGSXhfKXlZXRvn37mv0B02+nr1t7DVVTK+t+R192rXg1LGQMNYWFhb5nVkr5tQHXos0LGPtzgKfqyTsFLUh9V1Nagv63B7AZmNjYOUeOHKmay8qVK2t+/+///q+66667ml1XsDDLaDB37lz10EMPhV4YL3iTL1BUVVWpNm3aqLKysmbXEUz5zBw7dkx16NBBud3uJpcNlYzNxZAvZ8Zc5Qb1FKiVoBQoN6icGXPDQj4DQ04F6h1Q74aBnOF4j4ENyss7NRDmo0VAX9N+IlDsmUlEhgH/B8xUStUM3CmlivW/R4D30IaagkbfRYvInTlP893zy1/S59k3w8aXTUOkpqZ6dWDV0ti1axe9e/cOih14QMnO5tvHcygvb4/DsZ+UFFg+b2XYP0dNxfCv9B5xlED4+OLywJCzyJnMXmCJdAhLOcOVQAwNrQfOE5F+wAHgBuBGcwYRSQLeBeYopXaY0tsBDqXUKf33JcDDAZCpXvblweS1mhXESGCs+k7zZRPmUaZSU1OZP7/lWw1FykTx8tLRjF0wk9s5wy9IpkdhZ2IXVLN87mKmWi1cINGt13b3+xfpn31G4sCB4THB54nJyi71k09Ymp0dHVZ2AcLvHoFSqhq4G/gEbdjnLaVUvojcKSJ36tkeALoC/xCRTSJieGfqCXwhIpuBr4GlSqmP/ZWpIQZ+9VGNFcRPgUzC05eNJ+effz67d++msrKy8cwRTH5+fliu6/BkyXPFxFLBk1RzPtCVEmKpYMlz53SGI56ysjIOHz5Mv379rBbFJ6Kl9xxIArKOQCm1DFjmkfas6fdP0d67nuV2AyGNRdjHXVTz20oriKbSpk0bkpKS2LlzZ0R8MTeX/Px8Lr/8cqvFaJR7j8+jNVW4qf2aak0V9x6fB4TXKmh/2bZtG+eff77mvycC6NOnD2fOnOH48eN07drVanEigqhzMXHAoXVs3URelKlo+NKJlKGhJLQPB89oZUZ6S8LwxxUpiAhDhw5ly5YtjWe2AaJQERSMuwwFvE2tDWs4+rLxRktXBFVVVRQUFATErW6wOdtRiyF9Bm2FpNsjvSWxZcuWiFDOZmxF0DSiThEkpUHejLkslvbsJ3ytILzR0hVBQUEBiYmJxMXFWS1Ko8RNm4IrpjXtgReAfYArpjVx06ZYLFngibQeAWj/K7Yi8J2oUwT7Z80iffF8Cselk/b552ETkaxRsrNp9dwKNr2zGLc4WlboSj085eq0iQwtKIiMa1u0COfLL0ByMkOBLT16aPuLFlktWcCJ1B5BS/5oCjRRpwhAW0QXaQ937ppSpq/6P/pQjRsVXiEc/cRwbBbjPsbPCLPwlA0xezbs3cvQX/6SLb/+tbbfwjh9+jSHDh2if//+VovSJOyhoaYRlYrg4MGDtG7dOqIsCjovXUgb4GNql/lHgtmrLxiOzX4EXKynRdK1teSXzrZt2xg0aFDEWAwZ9O3bl7KyMkpKSqwWJSKISkWwZcuWiBvzNMxbY6GOc61wN3v1BeManNR9ICPl2lqyIojE/xXQLIeGDBnSYu9LoLEVQYRgmLe69M0zPZIpdiZRjRaRTHmkRwLGC0e1wGhlkWLO642WblwRSKJSEUSiFYQRunIV8IaeFilmr41RMn02BcDN1PZ2IunaunbtStu2bTlw4IDVogScSPxoMmjJPbVAE5WKIBIfbsOpFo7eLCCyzF4bI31CPB+PvpoTxOFGIvLaWtRLJzsbNW8RRTEp5H/wAV2uuTv8rbg8yc6mckM1zz2Xj8NBi3UKGCiiLmaxUioiewSGU62BZ37Hrs6d6XlqJ4mtrAm4EXCysjhdNZ+R3x+I49FHwzY8ZUMYiuCSSy6xWhS/yV1TyuS1z3MWzXvkBHcxzghwzGhmeeloZr55BWspZwlCSUt1Chggoq5HUFJSgojQo0dkrgCNi4ujb9++FBQUWC1KQInEXpqZlrSAybDiigUeRPtajCQrLtCcAvbnLG+hEFq2U8BAEHWKoLCwkNTU1HqDYEcCLWoYQifSFUFLWsAULuEp/eHe4/OIo7qOlV2tU0AbT6JOEezduzeiXzigvXS2bt1qtRgBw+VysX379ojwMVQfhnJuCZZDhrVWNZFroWY4/3NR6wfKnG5TF1sRRCAtzT5679699OjRw9L4rv7SvXt3YmJiOHTokNWi+I1hofZ7tAAjEFlWXFDr/M/sXNKcblOXqFMEhYWFEa8IWtrQUKQPCxm0lPuSPiGeVeNv521iECLTQs1wCtgJeEpPa6lOAQNBQBSBiFwqIttFpEBE7vdyXETkb/rxHBEZ4WvZgKA7NSuKSaF482Y6/eCWyDOHMzF48GB27NiBy+VqPHMEEPGKQH++kv67kfypUyPDaV5DZGVR9eDVHGwTw6DKyshxzGhGdwo4NCGBLQDJyZHtFND0DguG00m/FYGIOIGngcuAocAsEfH8r74MOE/fbgeeaUJZvzGcmrVzFfJzYKT7QGQ4NauH9u3b06NHD/bs2WO1KAFh69atEa0IjOdrhiojiQhymtcARUVF9O/fn1aRbKI8ezbJ+/ZxvG1bTuXmRrRTQOMZS3QV4giC08lA9AjGAAVKqd1KqUq0ha8zPfLMBF5VGl8BnUSkt49l/cYwh+sIzEOzIog0czhPWsowBGg9giFDhlgtRrMxnq+rACPIZqQ/Xy1hLg3A6XQyaNAgtm3bZrUofmE8Y1XUumEJ5DMWiAVlfYD9pv0iYKwPefr4WBYAEbkdrTdBz549WbVqlc8CTtTN3hyc67CtKfWEirKyskbl6tChAx9++CEdO3YMjVAmfJHPV5RS5OXlcfz48YDVGUj5fKE5z1eoZWwqO3fuJC4uLmxlbEr7de3alXfeeYfTp08HVygPAnmPJ7r2UUHtmg6DgL3DlFJ+bcC1wP+Z9ucAT3nkWQp8z7S/AhjpS1lv28iRI1VT2O9MVgrO2fY7k5tUT6hYuXJlo3leeOEFNWfOnOAL4wVf5POVwsJClZCQELD6lAqsfL5gPF9uUJX638aer1DL2FQmTZqkFi1aZLUY9dKU9nvkkUdUVlZW8ISph0De4/3OZJUH6oyf7zBgg/LyTg3E0FAR0Ne0nwh4Lt+rL48vZf3GMIczE2nmcJ60lLUEET9RTO3zJWhhKzcR+c9XSxkagpYxjFoyfTbrAPPUcCCfsUAogvXAeSLST0RigRuAJR55lgA36dZD44BSpdRBH8v6jeGwrciZHLFOzTwZMmQIW7duxe12N545jGkJisD8fH0DfCRdIvr5qqqq4uDBgwwaNMhqUQJCS1AE6RPiWT1gLGekU1DeYX7PESilqkXkbuATtNgiLyql8kXkTv34s8AyYBpQAJQDP26orL8ynYPusA3ms2rVKiZPnhxxTs3qkJ3NvjWltD9dQaHTSStnMiXTZ2sPRaSY+GVnk7umlPUf/I2JqoyipxdH3jUYmJ6v1CefZNeuXaQ/Nd9ioZpPQUEBPXr0oE2bNlaLEhAGDBhAcXEx5eXltG3b1mpxmkdWFkc//5yrn3gVxxVXBNwxY0C8jyqllqG97M1pz5p+K+Dnvpa1aRjDlOw6tGAu/VyF9Ikw75DGNXQALkUzuYy0a/DG0KFDWbIk4J3akLJlyxaSk5OtFiNgxMTEMHDgQLZv387w4cOtFqfZBNNrctStLG4JGKZkf0FbmAGRZ65oyPokkKKnRdo1eKMlDEO0NEUAkX9fysrKOHLkCCkpKUGp31YEEYjhBTKGyIzxC5qsp2h5MZj79OlDeXk5x48ft1qUZrNly5agvXCsItIVwbZt2zj//PNxOp1Bqd9WBBGI4QVSoXmI9EyPBIqdSWyjrvxGeiQjIhFv0WX3CMKPYBtV2IogAjHMFU8Dj6AphEgzVyyZPpu1aP5FDCLtGuojkl86LpeLnTt3kpQU2QrZk0i+J2ArAhsvGOaKJ53JPAesd/SJOHPF9AnxfJGcSYV0aTEmvQYR+dLRnZqtjU2i55kzpE6/PbId55nJzubs/a9RuGMHZ0Qi0ilgsBVB1MUsbhGYzBWHXnQRJVlZjPnBDywWqolkZXFoyRK+98q7OCZNisg4xfUxdOhQPv30U6vFaBKGFdcXaOEp+7r3k9gCrLhAu7bhHz7GD4FSItNCze4R2DRIampqRIZIVEoF1RzOSiKxR2BYoo1D8/MCLcOKC2qv7RXACEsTSdd25swZDhw4wIABA4J2DlsRRDiR+NIBOHToEK1ataJ79+5WixJwkpKSOHHiBKWlkeOGOsHkOM/pJT2SiXQru23btjFw4EBiYoI3gGMrgggnUoOmtwTXEvXhcDgYPHhwRFkOGdZaZjfH5vRIJtKt7LZs2UJqampQz2ErgggnNTU1IoOmt9RhIQPjvkQKJdNnU40WIapST2spVlyGld1e4G96WiRdWyg+mmxFEOF07dqVNm3acODAAatFaRKh+MqxkkjrqaVPiOeT799JAU5aIex39G0xVlyGlZ04kpgL7HIkRdS1hUIR2FZDLQDj6zMxMXLsbvLz87n++uutFiM4ZGfT9t1v2Pj1B7gff4JiZ1L4O9TLyqL6/MUMa7MPx9Kl7NKdM7YITFZ2KYMHU/7226Snh7m9kO6UsfPShWxxFdJlyTpyL/82aM+Q3SNoAUTa16dhMdRSewS5a0q57Ou36URFUOLLBouWfE8MIsXKzjDn7eEq5DrgQvfBoD5DtiJoAUTaePThw4dxOBwt0mIINLPEfsBCghNfNljYiiB8MExeY9DWdbQiuM+QrQhaAJHWIzDmB0Sk8cwRSIJrH04iz6GerQjCB/OzElNPeiCxFUELINIsh1q6xZBhlugC3F7SwxGXy8WOHTsYMmSI1aIElUhRBMazUo32HHmmBxq/FIGIdBGRz0Rkp/63s5c8fUVkpYhsFZF8Efml6dhDInJARDbp2zR/5IlWunXrRmxsLAcPHrRaFJ9o6RZDhrni28BHelq4myvu2rWLXr160a5dO6tFCSqDBg2isLCQiooKq0VpEOMZ+h1grEYJ5jPkb4/gfmCFUuo8YIW+70k18P+UUkPQVrD/XETMn4NPKKUy9c2OVNZUdGdhA46eJLdPn/B1qKXLWRSTQv6zz9L9l38MTzkDgGGu6HJ050mICId6eXl5LVo5G8TGxtK/f3+2bdtmtSgNYjxDbxKDk+A/Q/6aj84EJuu/XwFWAb8xZ9CD1B/Uf58Ska1AHyByZjfDGMO64GqgPeHrUMuQE+BiYKr7CF3CUM6AoJsrdiy8nd9ecAGJxXvD3qFeNMwPGKSlpZGfn09mZqbVotRPVhb9T5/maPcnOO/kGWJiYoL6DIk/48oickIp1cm0X6KUOmd4yHQ8BVgNpCmlTorIQ8AtwElgA1rPoaSesrcDtwP07Nlz5BtvvNEsmcvKymjfvn2zyoaKpsg44KKb6OveTxV1/cTsd/Rl14pXLZfPwJBToY15OtEmUoMhZ7jcY6UUl19+OYsWLaJjx451joWLjAYPP/ww48aN45JLLgHCTz5P/JHvlVdeobKykttuuy3AUtXF3zbcvn07f/nLX/jnP/8ZMJmmTJnyjVJq1DkHlFINbsByIM/LNhM44ZG3pIF62gPfAFeZ0nqivRMcwHzgxcbkUUoxcuRI1VxWrlzZ7LKhoikyuhClQLlBKdPmQsJCPgNDTlcI5Aynezxu3Di1evXqc9LDSUallEpLS1PffPNNzX64yeeJP/K9/fbbasaMGYETph78bcOXXnpJzZ49OzDC6AAblJd3aqNzBEqpqUqpNC/bYuCwiPQG0P8e8VaHiLQC3gEWKqXeNdV9WCnlUkq5gX8CYxqTx6YuZisCs8OwcLNQMVvShMIKIlxIS0sjLy/PajEapKqqioKCAgYPHmy1KCEhUiyH8vLySEtLC8m5/J0sXgLcrP++GVjsmUE0Y/EXgK1Kqcc9jvU27V6J1tOwaQKGdYEAfwYOEZ4WKoacfwbW62nhKGegiQRFUFBQQGJiIm3btrValJAwcOBADhw4QHl5udWiNEgkKYJHgYtFZCfaHOCjACKSICKGBdAEtFgX3/diJpotIrkikgNMAe7xU56ow7AuKHImswr4xNE9LC1UDDnfJpZKIsOSJhCkpqaGryLQLblWp01kaEFB+FqcBZLsbLZd8xDJFS62tGsX1tccSkXgl9WQUuo4cJGX9GJgmv77C+ousDTnm+Mt3aYJmBxqpf+//8fBbt1I/+1vLRbKC1lZpLrd7Oz4vwwrOkynTp3C3pImEBgWKkqpsFtJbVhyfQ38jLoWZy11kNa45mvRJifD1crOCGyUlBSaoVN7ZXELIj09PXy/PoE9e/bQpUsXOnXqZLUoIaNnz56A5l8p3DD82fwIrTsPkeETyR+Ma34AGKanheM1G+a8DkdoXtG2ImhBpKenk5uba7UY9ZKXlxf+7n8DjIiE7TyB4bfGMNvzTG+JmK85nENyhnJYCGxF0KIYOnQoO3fupKqqympRvJKbmxt1igDCd8K42JnEWbSIZC0tPGV9mK8tnMNWhnqlt60IWhBxcXEkJSWxfft2q0Xxiq0IwouS6bPJR1ulacxetHRLLsN6DeAR4BThec12j8DGL8J5eCg3NzekD3e4EK526+kT4lmWeTmnpS1uJCosuQzrtWJnMp8Anzp6htc1L1yISk4md9Uq0n78Y1gYmrkLO1RlCyNcvz7Pnj3Lnj17ombRUg3Z2cjKQ+R9tQ6XCAedyZRMn03fXsfA6lCQWVmUHD7MuBu+h+M3vyERai25Vq2yTq5gYrKyG3bHHRxMT+fqu++2WCidWbNw/fs9jlWfRQG9DhzAdcutOD/8EBYtCuqp7R5BCyNcewRbt26lf//+tG7d2mpRQkrumlImfPwEU1CchJqwlfvCRFdv3ryZjIwMq8WwhGHDhpGTk2O1GDWcWbYSZ/VZ9gF/Qhuuc1af5cyylUE/t60IWhjhqgiidX7AMFd8CzAGHwQY+NVH9RcKEUopNm/ezLBhwxrP3AIJN0XQ+qTmoScTuNVLejCxFUELY8CAARw5coRTp05ZLUodolURGGaJraj7z9bHXWSJPGYOHz6MUorevXs3nrkFYqy7cbvdjWcOAfvQLJfMXoTN6cHEVgQtDKfTyZAhQ8JuniAa1xBArVmi4X7b4IDD+nXVOTk5DBs2LOxWPIeKTp060a1bN3bv3m21KAA83nU+Z2lVx3nkWVrxeNf5QT+3rQhaIOE4PBStFkOGueI+4Ek9TQEF4y6zTCYDQxFEM+E0PDTjjgRO05pL0dZ2HKczlbRhxh0JQT+3rQhaEroTscRXFpN7xx3WO9TS5clz9uVEUREx500JWwdfwcIwV3Q4kngA2O7oS96MuSSFgU60FYGmCDZv3my1GABMjV/POz99kjUxg4kTxcjk71g3dzFT49c3XthPbEXQgjAcal2uShhGrYVK7ppSS+WpdhfxGJDk3mepPJaQlUX64vmkuAoZOmoU333xJumL57N/1iyrJbMVAeHVIyArizYTY7n66gzcbti7F6bOnwJZWUE/ta0IWhCGhcpo4MdQE6fAKodahjypwB16Wjg6+AoVGRkZbNq0yWoxAC0YzY4dO6ImTnF9ZGRkhI8iADZt2mRJLGVbEbQgDAsVoe5KQascahnn9bSCCDcHX6EiMzMzbBTBtm3bSE5OJi4uzmpRLGXgwIEcOnQobKzsIlIRiEgXEflMRHbqf70GrheRvXoAmk0isqGp5W18w7BQETSHWi6PdKvkMVtBWCmP1YSTIrCHhTScTidDhw4NCys7Y11HxCkC4H5ghVLqPGCFvl8fU5RSmUqpUc0sb9MIZodaK4CXsdahVsn02ZwFLkVTBlgsj9UMGzaM/Px8XC5X45mDgT55XxSTwuYf/Yjktz+Nusn7OujtMeCbrWwaP95y44oDBw7gdDrp1atXyM/tryKYCbyi/34F+GGIy9uYMIetPAm8SZylDrXSJ8SzZNKtFNOKmChxatYQHTt2pFevXuzcudOS8xuT94muQoYAM9SJ6Ju8N2G0x5XqNAlYb1xh1bAQgCilGs9VX2GRE0qpTqb9EqXUOcM7IrIHKEH7IHxOKfV8U8rrx25H85hLz549R77xxhvNkrmsrIz27ds3q2yoCISMxcXF/PKXv+Ttt98OkFS1NEW+jz76iI0bNzJv3ryAy1Ef4XyPH3jgASZPnsyYMWNCLuOAi26ir3s/oA0dOvRtv6Mvu1a8WidvOLchBEY+oz1caC8mY17NW3s0h6bK+Prrr3P69GnuuOOOxjM3kylTpnzjMSqjoZRqcAOWA3letpnACY+8JfXUkaD/7QFsBibq+z6V99xGjhypmsvKlSubXTZUBEJGt9ut4uPj1eHDh/0XyIOmyHf33Xerv/71rwGXoSHC+R4//PDD6je/+Y0lMroQpUC5oc5fF3JO3nBuQ6UCI19T2qM5NFXGa665Ri1cuDAg564PYIPy8k5tdGhIKTVVKZXmZVsMHBaR3gD6X6/ekZQWzB6l1BHgPWpDY/tU3qbpiAgjRozg22+/tVSOjRs3Mnz4cEtlCCcyMzMtW8BkTNK70XoE4pEebXgaVyiP9FBj5dCQv3MES4Cb9d83A4s9M4hIOxHpYPwGLkHrUfhU3qb5DB8+nI0bN1p2fpfLRU5Ojq0ITFhpOWQYE7wKfKanRfPkvdm44k20oQ+r2uPUqVMUFxczaNCgkJ8b/FcEjwIXi8hO4GJ9HxFJEJFlep6ewBcishn4GliqlPq4ofI2gWHEiBGWKoKdO3fSo0cPOnXqZJkMYUV2NiU/f4aKQ0cYNGVKyK1UDGOCN4jjJET95L3ZuGI/8LZ0CH176JZLKzqdR2p5OYfaDLTEcsmvCGVKqePARV7Si4Fp+u/dgNfIF/WVtwkMI0aM4MEHH7Ts/PawUF1y15SS/sGfmA2UAYNchfRZsoBc5hISv6x6dK78xFcY89//ktivH9b7QLUQU7Sy0StWsOwPfyB9cfA9fZoxLJd2oNnOJ4b6mdCxVxa3YAYNGsShQ4c4ceKEJef/9ttvGTFihCXnDkcMlxuPAwP1tFC73Dh48CBnzpwhJSUlZOeMBIz5tFCv8TCeiSuotZ23wg2LrQhaME6nk2HDhlk2Jm33COpiuNaIoe4/XihdbnzzzTeMGjUqamMQ1Efnzp3p1asX27dvD+l5zW5YrHomjPPbtGCsshxSStk9Ag/M1ijV9aQHmw0bNjBq1Llm5DYwatQoNmzY0HjGAFLsTKICLf6AlW5YbEXQksnOpvsXe/ji//0etzhCMzmpT359FdOXuJISqvqMjW43BiYMKxUBFgDHCb2Viq0I6scKRVAyfTbrgbupNee1wnLJVgQtmNw1pVyxeRmt1WkcqJAsoTcmv866D7AA65fthxOGlcoBZzJrgPcd3UNqpaKUshVBA1ihCNInxLN46EW4pQNuC92w+GU1ZBPedF66kJ7A/1Ebm6B2Iio41hHG5Nd44EI9LdjnjBhMViq9br6ZXX36cOuC0LXJgQMHUEqRmBjVtkL1Mnz4cDZv3kx1dTUxMSF6NWZlsf+bb7j8/qdxzJlDIlhiyWX3CFowCa59tAJiqe12GunBPCdo8QfsGAT1M2TIENatWxfScxq9AXui2DsdO3YkKSmJ/Pz8kJ533bp1jB07NqTn9MRWBC0Ys0sBl5f0YJ3TDZzFjkHQEEOGDGHDhg0hNVe0h4UaJ9TDQ4cPH6a0tJSBAwc2njmI2IqgBWNMTq4CFulpwZ6IKpk+m3w0fyFWTn6FO/Hx8XTv3j345oqmGATr588n5ZFn7Mn7+sjOJmHTIVbd9quQGVesW7eOMWPG4HBY+yq2FUELxpic7Ojow33AfkdS0Cei0ifE8+9hl4G0s3TyKxIYO3Zs0IeHjMn7Pq5Cvgdc7j5qT97XQ+6aUq7IW04rVRYy44pwGBYCe7K4ZaNPTir1CM7ERCpX/4f0AQOCfs7CLVuY+rO/4rjjDssmvyKBMWPGsG7dOn784x8H7RzG5L0b+C21//D25P25dF66kIHACEJnXLFu3TruueeeoNTdFOweQRQgIkyYMIG1a9eG5Hxr165l/PjxITlXJBOKHoF5kj6mnnQbjQTXPuIInXGF2+1mw4YNjBkzpvHMQcZWBFHC+PHjQ6IIjhw5wpEjRxg6dGjQzxXpZGZmsmPHDsrLy4N2DmOSvhqtV+CZblNLqI0rtm/fTteuXenevXtQ6m8KtiKIEkKlCL788kvGjRuH0+lsPHMU03fRInZe/0fOK69mfbt2QZuYNAwGfgIU62n25L13jLb6L2C4fAtKW+kT+B+lXsjo3btD7o7cG7YiiBIyMzPZtWsXJ0+eDOp57GEh39iXB2lLFjCHSoTgrcBOnxDPf3/wK5Yg9MCOQdAQhnFFZ0civwL2Bcm4wpjAr1bH+QXhsfreniyOEmJjYxk5ciTr1q3j4osvDtp51q5da2kMhEhh4FcfIcAvqR2PDsrEZFYW3w16nwmyjdiPPrIn7xvCtPK7+/nnU/LWW2RkeA2l4hfGBP7/AK31NKtX3/vVIxCRLiLymYjs1P929pLnfBHZZNpOisiv9GMPicgB07Fp/shj0zDjx49nzZo1Qav/7NmzfPvtt2FhDhfu9HEXAbUrsI3Fd8GYmPz888+ZOHFiwOttyUycOJHVq1cHpe4E1z4q0L7CQ7XivzH8HRq6H1ihlDoPWKHv10EptV0plamUygRGAuVoAewNnjCOK6WWeZa3CRDZ2ST8Zxf/+cNjgV8so495ftI2hfNOn6a0c7rlY57hzgFH7Xd5sAOnr1692lYETWTSpElBUwTFziTWU9cVuZFuFf4qgpnAK/rvV6gNslMfFwG7lFKFfp7Xponkrinlhq/fJpYKCPBiGWPM87T7EA8RHmOe4U7BuMtqbNUXA/8mOBOTpaWl7Nixg9GjRwe03paO0SNQSjWeuYmUTJ/NW8DbpjSrJ/D9nSPoqZQ6CKCUOigiPRrJfwO13g4M7haRm4ANwP9TSpV4KygitwO3A/Ts2ZNVq1Y1S+CysrJmlw0VwZBxwIev0R34gLpj0p0+fI1Vq5o2Z+Ap34APX0OAq4BWftYdCCLhHvcaVMkqbmfgVx9x1r2fl2lD9/E3kdTrWEBl/+qrrzjvvPOabDEW7m0YCvkcDgevvfYaSUnN+1KvT8a+vY7xQVxnflARi1sd4YAjkYJxlwX83jcJpVSDG7AcyPOyzQROeOQtaaCeWOAYmvIw0nqiDZE60GZJXmxMHqUUI0eOVM1l5cqVzS4bKoIhowtRClQ1KGXaXIjf8rkQ5QZVFYC6A0Gk3ePDhw+r+Ph4VVVVFfDz/OY3v1EPPvhgk8uFexuGQr6bbrpJPfvss80uX5+MBw4cUF26dFHV1dXNrru5ABuUl3dqo0NDSqmpSqk0L9ti4LCI9AbQ/x5poKrLgI1KqcOmug8rpVxKKTfwT8D6JXYtFGP8URH4MInFziR2AFX1nNOmYXr06EFycnJgvV4uXAgpKXz+2GNMevZZbd+mSUxq3ZrV994LDgekpASsDVesWMHkyZPDaq2Nv3MES9AcTaL/XdxA3ll4DAsZSkTnSrSehk0QMBbLlAOPUOsmOhDjkiXTZ/Nv4DlTmtVjnpHGRRddxPLlywNT2axZuG65ldOFheQAYw8fxnXLrTBrVmDqjwZmzWLCi6/weXm5NnpRWBiwNlyxYgVTp04NgJCBw19F8ChwsYjsBC7W9xGRBBGpsQASkbb68Xc9ymeLSK6I5ABTAOu9L7VQjMUyJ53JLAf+7egRsMUy6RPieadLIt0dPWyPo81k6tSprFixIiB1nVm2Emf1WbYAfwXaAs7qs5xZtjIg9UcDZ5atZJCrkhGAMWkZiDZUSrF8+fKwUwR+TRYrpY6jWQJ5phcD00z75UBXL/nm+HN+myZgWiwzbcEC1h85wuwnA7N45buf/pSCRx7h6tNHcLRpYy9aagYTJ07kuuuuo7y8nLZt2/pVV+uT2ghtBponTc90m8ZpffIIgmbZ08oj3R927NiBw+GwPBCNJ7aLiShk2rRpLFsWuCUbn376KZMmTaJNmzYBqzPaaP+PfzCodWfe69DP73Ue+9DmZhzUDRdqpNs0jtFWraj7kmxWG5qCA302eDAXFH1H3g9/F1brbGxFEIVkZGRQVlbGzp07A1LfsmXLmDbNXhTuD7lrSpnxXREl7iN+B0V5vOt8CnBSRe1CtbO04vGudvwBX3m863zO0gqBmnZsbhsa62wSXYV0An6iTofdOhtbEUQhIsK0adP46KOP/K7L7Xbz8ccfc9lllwVAsuil89KFXAn81JRW63+macy4I4F/4eQvxKKA43SmkjbMuCMhQNK2fGbckUAlbfiOzrwOLKZds9vQ8C2k0BZSXULz722wsBVBlBKo4aFvvvmG7t27k5KS4r9QUUyCax9pnDtp1xz/M1Pj1/Ov3v34e4+lxIhiZPJ3rJu7mKnx6wMiazQwNX496+YuZmTyd/yEfzK77eXNbkPP4EDiJd1qbO+j0Uh2Nr1XHeGLTz7llAilzmRKps/WrHyyshot3nfRInKf+IzOSxey1FXIhdKR3JnzfC5vcy7FziQSXYWcRXtROM3pTaxr3w03cPSxxzh0aBKtamY6p+ibjU9kZTEV2DsfjhyZwaBB9/G9378MbZrehsa9raLWyWBNeuAk9gu7RxCF5K4pZdxHf+UWFEU03TeQ4Us/0VVID+A2dTLsxjwjDWOdRx7wtJ7W3LUY77//PldccQWtWrVqPLNNo/To0YPMzEw+++yzZpUvmT6bM8AVwCk9LdzW2diKIAoxxiyfAAbpaU0ZszR86Ss0508jmlje5lyMdR7dHUn8BfjM0bvpazH01cTv/vKXXL1smb2aOIBclZTEu7NmNWuVcfqEeJ4ddRWnaE3HMF1nYw8NRSHG2KQxXml4wfR1zNLwpa+wA6IHDNM6j9m//S2fuVxkZzfBQmXWLFz/fo/j1WfZBFx89CiuW27F+eGHsMjTz6NNk5g1ixlvv8vDrkqqgFb6KmOf2zYri9VfXslP7/oHjp/8JCzX2dg9gijE8AEkaH6Hqj3SG+OAI5Fqat1UeNZr4x9z5sxh4cKFuFyuxjPrGKuJvwKeAdpgryYOFGeWrSTFVclMwPCf35S2PX78OCtXruSaa64Jmoz+YiuCKMQYjwbYC/wYcOP7mGXBuMv4AC3UnmEBEW5jnpHM0A8/JP60m7diE3xeXGaseL0EuN5Luk3zMdrwWaCfl3SvZGej5i2iKCaFN7p1Y+LJKgrnPBZWi8jM2IogCjHGo4ucyfQDviKGVyb8yOcxy6Q0yO6WwjDpavsWCgK5a0q5rfQQO5uwuGwfSRxAG6rzeyWsTR2MNoyhbljRhto2d00pk9c+T6KrECfwG1Ue1gYVtiKIRrKySF88n8TqvTgfe4xr06bw7tp3cf/m/vq/Pk3L5M8+/zy7ju1j/LRbcTz2KInVe0lfPN82HQ0QnZcu5GbgPmpfOo1Nxj/edT5PIHWG6uzVxIHBvMrYhbbSuLG2NQwy3GiLBMcT3gYVtiKIcnLXlPKbvM+oVOWUNfD1aV4mvwt4GDcjlmaH7RdOJJPg2kcXoDW+BzcfMrOMRcAxOuPGXk0cSIxVxsfpzGngcqCQ2Abb1rhXhkFFOC4iM2Mrgiin89KFdEILYdleT/P25WJeJn8ZmtloOH/hRDLmyfwqtK9KczpQp4fmFgcfvHg3N3Qbwcud7rdXEwcY8yrj+TzGWecY/oyb7y/4fr096GJnEjlAJZFhUGErgijH+EIxvCwqj3TPfG7qjkOH6xdOJGNM5juAz4G5aPclwVVY8+I5+PLHNT20QyhiqOaRY99w+cRS3G7Yuxemzp9iD9cFgqwsps6fwt698KMZpbzh+pp9nKEMj8WYJuXc3VXIj4AviAyDCnsdQZRjLH83j39qLg4cIMIpOgLQEcXnaLFE21D7cIfTMvmWQvqEeHKZS+elCxnvKuRnwHTgQrQXT58lC4Dae9ATeAdNmWs9NHteIFh0XrqQ3mg9aMNVhABDljyGWuIiTd9fjebQ4yK0j6disxuXMMSvHoGIXCsi+SLiFpFRDeS7VES2i0iBiNxvSu8iIp+JyE79b2d/5LFpOmZTUgfaK+SngBMXDiCek8RzkkNo4ef2EBlfOBGNaTL/O2cyC4FR1J04NhS3C+2FZHzR2T204GLuQTvR2l+bB3DV3Bc32uTwE2j/U8XO5LA3qPB3aCgPuApNAXpFRJxo7lMuA4YCs0RkqH74fmCFUuo8YIW+bxNCzKakCuHXOHBSN8C9ArqjfQUN1fdtk9HQkODax2ggltqXTDVwGniSWl/55h6aTfAwz98otPafhxYLvALIJTKHT/0NVbkVNP/2DTAGKFBK7dbzvgHMBLbofyfr+V4BVgG/8UcmmyZicm0A0F4cPEnti8VY22p+UBRCYvVee0goBBhDdw60F4zx4m8N/ALtq9Szh2bfl+BRMn02fZYsqPn6F7SY0IL2P3I+tXNtkTR8Goo5gj7AftN+ETBW/91TKXUQQCl1UER61FeJiNyOZqxCz549WbVqVbOEKSsra3bZUGGljAMcifR170dBzYSlJwcciewK4zZsSfdYjf0BfdY+X/PSaU3dlwz6fpGjLwXjLiOp17GAXHu4t6FV8vXtdYxV429n4Fcf0cddhBMH83DV+7+igJ1jf0BBGLclAEqpBjdgOdoQkOc205RnFTCqnvLXAv9n2p8DPKX/PuGRt6QxeZRSjBw5UjWXlStXNrtsqLBSxpwZc5UblKpnc4P6z/jbLZPPF1rUPX7sMZUzY67a70xWLlAuj/tRQSv1i66vWyefRYSLfL/o+rqqoFWde2Lcp32OvipnxlylHnvMajFrADYoL+/URnsESqmpfuqaIqCvaT8RKNZ/HxaR3krrDfQGbMcoFmO2WElwFdZYDXXgZI3lQ1KvYxZLGUWYhu6Wz1vJ2AUzqSSGzpRQQmdiqbYXjVnIjDsSqFzQhjLa17kn6+YuJuZiYfLkyVaL6BOhWEewHjhPRPqJSCxa2M4l+rElwM3675uBxSGQx6YhTBYrDqWIV6XEq1IcStVYPuyfNctqKaMS88Ime9FYeNBS7olfcwQiciXwFJpRyVIR2aSU+oGIJKANB01TSlWLyN3AJ2hzWy8qpfL1Kh4F3hKRW4F9aMNINjY23jCFT6zFDkFpKQ3dk3CfFzDhr9XQe8B7XtKLgWmm/WXAOZHSlVLH0dZc2NjY2NhYhO1iwsbGxibKsRWBjY2NTZRjKwIbGxubKMdWBDY2NjZRjmhrDCILETlKbRzpptINCHdD+HCX0ZbPf8JdRls+/wlHGZOVUt09EyNSEfiDiGxQStXrKTUcCHcZbfn8J9xltOXzn0iQ0cAeGrKxsbGJcmxFYGNjYxPlRKMieN5qAXwg3GW05fOfcJfRls9/IkFGIArnCGxsbGxs6hKNPQIbGxsbGxO2IrCxsbGJclqsIhCRS0Vku4gUiMg5sZBF42/68RwRGRFC2fqKyEoR2Soi+SLySy95JotIqYhs0rcHQiWfSYa9IpKrn3+Dl+NWtuH5prbZJCInReRXHnlC3oYi8qKIHBGRPFNaFxH5TER26n8711O2wWc2iPL9WUS26ffwPRHpVE/ZBp+HIMr3kIgcMN3HafWUDXr7NSDjmyb59orIpnrKBr0Nm4W3aDWRvqG5u94F9EeL+70ZGOqRZxrwEVrUv3HAuhDK1xsYof/uAOzwIt9k4EOL23Ev0K2B45a1oZf7fQhtsYylbQhMBEYAeaa0bOB+/ff9wGP1XEODz2wQ5bsEiNF/P+ZNPl+ehyDK9xBwnw/PQNDbrz4ZPY7/FXjAqjZsztZSewRjgAKl1G6lVCXwBjDTI89M4FWl8RXQSY+SFnSUUgeVUhv136eArWixnSMNy9rQg4uAXUqp5q42DxhKqdXAdx7JM4FX9N+vAD/0UtSXZzYo8imlPlVKVeu7X4F1sdbraT9fCEn7QcMyiogA1wGLgnHuYNFSFUEfYL9pv4hzX7S+5Ak6IpICDAfWeTl8gYhsFpGPRCQ1tJIBWuztT0XkGxG53cvxsGhDtKh39f3jWd2GAD2VUgdB+wgAenjJEy5t+RO0Xp43Gnsegsnd+tDVi/UMrYVL+10IHFZK7aznuJVtWC8tVRGIlzRPO1lf8gQVEWkPvAP8Sil10uPwRrShjgy0KHDvh1I2nQlKqRHAZcDPRWSix/FwaMNYYAbwtpfD4dCGvhIObTkPqAYW1pOlsechWDwDDAAygYNoQy+eWN5+OrNouDdgVRs2SEtVBEVAX9N+IlDcjDxBQ0RaoSmBhUqpdz2PK6VOKqXK9N/LgFYi0i1U8unnLdb/HkGLRDfGI4ulbahzGbBRKXXY80A4tKHOYWPITP97xEseq5/Hm4HLgdlKH8z2xIfnISgopQ4rpVxKKTfwz3rOa/mzKCIxwFXAm/XlsaoNG6OlKoL1wHki0k//YrwBWOKRZwlwk275Mg4oNbrvwUYfR3wB2KqUeryePL30fIjIGLR7dTwU8unnbCciHYzfaBOKeR7ZLGtDE/V+gVndhiaWADfrv28GFnvJ48szGxRE5FLgN8AMpVR5PXl8eR6CJZ953unKes5rWfuZmApsU0oVeTtoZRs2itWz1cHa0CxadqBZEszT0+4E7tR/C/C0fjwXGBVC2b6H1m3NATbp2zQP+e4G8tGsH74Cxoe4/frr596syxFWbaifvy3aiz3elGZpG6IppYNAFdpX6q1AV2AFsFP/20XPmwAsa+iZDZF8BWjj68az+KynfPU9DyGS7zX9+cpBe7n3tqr96pNRT3/ZePZMeUPehs3ZbBcTNjY2NlFOSx0asrGxsbHxEVsR2NjY2EQ5tiKwsbGxiXJsRWBjY2MT5diKwMbGxibKsRWBjY2NTZRjKwIbGxubKOf/A6A70xHHSjgiAAAAAElFTkSuQmCC\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": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABEZ0lEQVR4nO2deZwU1bX4v4cZEJAtgKK4BDWLAmFxUDMmmpkQE8AnolGiibhGIk9enhpCJPkhLk9QFOIal+QZn0tE8hSCBo2KM08jxAAyIAMaN1SCcSEijHHY5vz+uFUzNU33TC/V3dXd5/v51KdruVX39K3uU7fOPfccUVUMwzCM4qdDvgUwDMMwcoMpfMMwjBLBFL5hGEaJYArfMAyjRDCFbxiGUSKYwjcMwygRTOEboSEiPxCRp/Ith4+IdBGRx0TkExH5fb7lMVoQkQYROTSF8g+JyLgkyg0RkaUZCVfEmMKPICLyfRFZ4f0p3hORJ0Tk6/mWqz1U9UFV/Xa+5QhwGtAP6KOqp8crICIDRWSR91DYJiI1InJs4PgAEVHvXjSIyPsi8riInBBznQ0i8lmgXIOI3Jbdr1e4qGo3VX0zmbIiMgQYCvwhieuuAbaIyEkZiliUmMKPGCJyGXATMBOnrA4GfgWcnEex2kVEyvMtQxw+D/xNVXfFOygihwEvAC8DhwD9gQXAUyJSGVO8l6p2wymep4EFInJuTJmTPEXmL5Mz/QLx2jXVto7ovUmFHwEPavKzRB/0zjFiUVVbIrIAPYEG4PQ2yuyFeyBs8pabgL28Y1XARmAq8AHwHjAOGAP8Dfgn8PPAta4E/hd4GNgGvAQMDRy/HHjDO7YOOCVw7Fycsvyld93/8vb92Tsu3rEPgE+ANcDgwPe8D/gQeBv4f0CHwHX/DNwIfAy8BYxuoz2OAGqBLUA9MNbbfxWwA9jptekFcc69H1gcZ/8dwHPe+gBAgfKYMlOA9wNybwC+leR97hBo283AfKB3TH0XAO8AzyVo6/baMLb8F4D/8+7FR8DDbcg31mvLLV7bHhE4tsH77mu8az0MdE5wnYR1et/xC976vcDtwB9xv7UXgcMCZd8Evh5zf/43sH09sAQQb/sA4DO8/4UtgXuSbwFsCdwMGAXsilUuMWWuBv4C7AvsAywFrvGOVXnnXwF0BC70FMLvgO7AIKARONQrfyVOIZ7mlZ+CU7AdveOn43q9HYDvAZ8C+3vHzvXq+g+gHOhCa4X/HWAl0Aun/I8InHsf7vW8O07B/Q1PIXvX2OnJXgZMwj3YJE5bdAReB34OdAK+6SmMLwe+3wNttOU/gPPi7K8GdgNdSazwD/X2H+FtbyB5hX+Jdw8PxD3A7wIe8o759d0H7B1o19i2bq8NY8s/BPzCu5edCSjQGNm+5N3nE7z2neq1cafA9/yr97voDawHLkpwrYR1sqfC/ydwtCfvg8A879jeXtl9Aud29b7vucBxuIfJgTF1bwWG5Ps/HbUl7wK0KyDcg+slrg3pegcDT3k/1HXAgHx/x4BsPwD+0U6ZN4Axge3vABu89Spcz6bM2+7u/VmOCZRfCYzz1q8E/hI41gH3VnBcgrrrgJO99XOBd2KOn0uLwv+m96f8Kl7P09tfBmwHBgb2/QioDVzj9cCxrt532C+OPMfhlHbw+g8BVwa+X1sKfxcwKs7+w706DyCxwu/s7f+at70B9yaxJbBcmKDe9cDIwPb+uIdceaC+Q2Pa9Z3AdjJtGHtv7gPuJkYxxpFtOjA/5jfxd6Aq8D3PChyfDdyZ4FoJ62RPhf+bwLExwCve+gFe2c4x5x+Ne0i8DZwZ5/p/B47P9n+20JZCsOHfi+v5hsV9wA2qegTuR/NBiNfOlM1A33Zsrv1xP3Kft719zddQ1d3e+mfe5/uB458B3QLb7/orqtqEMwn1BxCRs0WkTkS2iMgWYDDQN965sajqs8BtuFf190XkbhHp4Z3fKc53OCCw/Y/Adf7lrQZl9ukPvOvJnehabfERTtnGsj/QhDMpJcKv45+BfeNUtVdg+XWCcz+PGwPw23U97o2iX6BMbNsGt5Npw9jzp+LetP4qIvUicn4C2Vr9vry2fZcE9wf4F/HvTSp1tnXNLd5n92BhVf0rztQjOJNYLN0D5xoekVf4qvocrf9UiMhhIvKkiKwUkedF5PBkriUiA3E9tae9azcEFEoUWIYzuYxro8wmnMLwOdjbly4H+Ssi0gFnZtgkIp8Hfg1Mxnm59ALW4v5gPtrWhVX1FlWtwJmSvgT8FKdkd8b5Dn9PQ/ZNwEGe3Olc6xmc2SqW8cCydn4bp+A6C68mWVeQd3HjEsGHQ2dVDcod27bB7WTasNX5qvoPVb1QVfvj3gZ+JSJfiCNbq9+XiAjuN5Ly/Umhzrau8SnurfZLwf0icjHOHLYJ92AJHuuPeyCmc2+Kmsgr/ATcDfyHp0ym4LxYkuFLOJetR0VklYjcICJlWZMyRVT1E5z9/XYRGSciXUWko4iMFpHZXrGHgP8nIvuISF+v/AMZVFshIqd6bxWX4EwFf6HFdvohgIich+vhJ4WIHCUix4hIR5xNuBHY7b19zAeuFZHu3oPlsjS/w4vetad67VQFnATMS/L8q4BjReRaEentyfMfwNnAzxJ8r34iMhmYAUyLebtIljtx3//z3jX3EZGTkz05nTYUkdNF5EBv82Pcvd0dp+h84EQRGendu5/gfhMp+7anUGd7LAa+Ebjul3AD0WcBE3D3f1igfBXwrKpuT6OuoqbgFL6IdAOOBX4vInW4Aa/9vWOnisjaOMufvNPLcXbfKcBRuIG3c3P9HdpCVefi/rz/D6ds38X1shd6Rf4LWIHzkngZ51nzXxlU+QfcgOzHuD/Pqaq6U1XXAXNwbx3vA1/BeX4kSw/cG8LHOBPBZpznDbjBxE9xr+R/xg0q35Oq4Kq6A+dRMhrX6/0VcLaqvpLk+a8BX8e5Wm7AjV98F/iOqsZ+1y0i8imuzcfgPKliZX4sxg9/QYKqbwYW4dw/t+EesMckI3OAVNvwKOBFEWnw6v5PVX0rtpCqvopTpLfi2vQknLvpjhTlS7rOJLgb+IE4ynEPtutVdbV3D38O3C8ie3nlf4B7qBox+G5MkUZEBgCPq+pgzw78qqrGs722d52vAtepapW3PQH4qqpeHKa8hYKIXIkbODsr37IYRluIyO9wg8kL2yn3FeBuVY2dR2FQgD18Vd0KvCUip4OzMYrI0CRPXw58TkT28ba/ifPUMQwjwqjq99tT9l65l03ZJybyCl9EHsKZFb4sIhtF5ALcK9sFIrIaN0EkKfunZ/ucAiwRkZdxA5CJPCkMwzCKioIw6RiGYRiZE/kevmEYhhEOkQ6q1LdvXx0wYEBa53766afsvffe4QoUIlGXD6Ivo8mXOVGX0eRLnZUrV36kqvvEPZjvqb5tLRUVFZouNTU1aZ+bC6Iun2r0ZTT5MifqMpp8qQOs0AIOrWAYhmGEgCl8wzCMEiFjhS8iB4nLErTeC5D0n3HKiIjcIiKvi8gaETky03oNwzCM1Ahj0HYX8BNVfUlEugMrReRpdVPzfUYDX/SWY3AJDFKdSg7Azp072bhxI42NjW2W69mzJ+vXr0+nipyQC/k6d+7MgQceSMeOHbNaj2EYhUHGCl9V38PFIEFVt4nIelwo1aDCPxm4zxtQ+IuI9BKR/b1zU2Ljxo10796dAQMG4AL5xWfbtm1079494fF8k235VJXNmzezceNGDjnkkKzVYxhG4RCqDd+LeTMcF8UwyAG0js+9keRjlreisbGRPn36tKnsDRAR+vTp0+6bkGEYpUNofvheFMtHgEvUxbtpdTjOKXGn+IrIRGAiQL9+/aitrW11vGfPnjQ0NLQrz+7du9m2bVv7gueJXMnX2Ni4RxsmS0NDQ9rn5gKTL3NyKWOP+np61dWxZdgwtg4alNQ5UW/DqMu3B4n8NVNZcLkv/wRcluD4XQTSkOESE+zf3nXj+eGvW7cuKV/UrVu3JlUuX+RKvmTbKx5R9DEOYvJlTlZlXLpUdeZM97l0qWqXLqplZe5z6dL8yxcCUZSPbPrhexlx/htYry6WezwWAWd73jpfBT7RNOz3UaGsrIxhw4YxePBgTjrpJLZs2QLAhg0bGDy4JUfIr3/9a4488kg+/rglU96NN96IiLB58+bmc7p06cKwYcMYNmwYF110UU6/i2FkhWXLYORImD7dfd53H+zYAbt3u89C6hUXEWHY8L+GS5zxTXH5T+tEZIyIXCQivvZajEvU8DouOuW/h1Bvu8yeDTU1rffV1Lj9mdClSxfq6upYu3YtvXv35vbbb9+jzP3338+tt97KU089xec+9zkA3n33XZ5++mkOPvjgVmUPO+ww6urqqKur4847LW+DUQTU1rZW8ACdOkFZmfvs0wdmzXIPBiNnhOGl82fi2+iDZRTIeZKRo46C8eNh/nyornbK3t8Oi8rKStasWdNq3/z587nuuutYsmQJffu25Py+9NJLmT17NiefnHQ2O8MoTKqqnGLfscN9nn22W2prnbK/5JKWY0uWQKWFsM8FRT3TtrraKffx4+GKK1or/zDYvXs3S5YsYezYsc373n77bSZPnsxTTz3Ffvvt17x/0aJFHHDAAQwdumeulrfeeovhw4fzjW98g+effz4c4Qwjn1RWOkV+zTUtCr2yEqZNg82bzbyTJyIdLTMMqqth0iT3u5s+PRxl/9lnnzFs2DA2bNhARUUFJ5xwQvOxffbZh969ezN//nwuvfRSAP71r39x7bXX8tRTT+1xrf3335933nmHPn36sHLlSsaNG0d9fT09evTIXFDDyDXLljkFXlXVouRjie39V1XlVsYSpqh7+ODMOHfc4ZT9HXfsadNPB9+G//bbb7Njx45WNvyuXbvyxBNPcOedd/Lggw8C8MYbb/DWW28xdOhQBgwYwMaNGznuuOP4xz/+wV577UWfPn0AqKio4LDDDuNvf/tb5kIaRq6JHahNZJ+P1/s3ckJR9/CDNvvqareEadbp2bMnt9xyCyeffDKTJk1q3r/PPvvw5JNPUlVVRd++ffnOd77DBx980Hx8wIAB1NbWst9++/Hhhx/Su3dvysrKePPNN3nttdc49NBDMxfOMHJN7EBtbW1iZZ6o929klaLu4S9f3lq5+zb95cvDq2P48OEMHTqUefPmtdp/yCGHsGjRIs4//3xefDF24nELzz33HEOGDGHo0KGcdtpp3HnnnfTu3Ts8AQ0jV/imGt8Tx0w1kaOoe/hTp+65z+/pZ0LsTN/HHnuseX3t2rXN60OHDuXvf//7Hudv2LCheZbtd7/7Xb773e9mJpBhRAHfVBO04SdDrN3fyBpFrfANw8gxqZpqfLu/uWjmhKI26RiGEXHi2f2NrGEK3zCM/GF2/5xiJh3DMDIjExt8unZ/Iy1M4RuGkT5h2ODNRTNnmEnHMIz0MRt8QWEKPw0ShUdOlqqqKl566aXsCGcYucRs8AWFKfw0SCY8smGUBBYmoaAoDYW/bFnWYm9XVlY2T67661//yrHHHsvw4cM59thjefXVVwEXbO2MM85gyJAhfO973+Ozzz5rPn/SpEmMGDGCQYMGMWPGjOb9AwYM4KOPPgJgxYoVVFnPyYgqfhRMU/aRp/gHbbM4scMPj3zBBRcAcPjhh/Pcc89RXl7OM888w89//nMeeeQR7rjjDrp27cqaNWtYs2YNRx55ZPM1rr32Wnr37s3u3bsZOXIka9asYciQIaHIZxgFh826zSqhKHwRuQf4N+ADVR0c53gV8AfgLW/Xo6p6dRh1t0sqAZ2SJFF45E8++YRzzjmH1157DRFh586dgIuX8+Mf/xiAIUOGtFLo8+fP5+6772bXrl289957rFu3zhS+UZrE65wZoRKWSedeYFQ7ZZ5X1WHekhtlD1kZVEoUHnn69OlUV1ezdu1aHnvsMRobG5vPcal/W/PWW29x4403smTJEtasWcOJJ57YfE55eTlNTU0Ara5jGEWLefxknVAUvqo+B/wzjGuFThYHlfzwyDfeeCM7d+7kk08+4YADDgDg3nvvbS53/PHHN8fGX7t2bXNKxK1bt7L33nvTs2dP3n//fZ544onmcwYMGMDKlSsBeOSRR0KT2TAii3n8ZJ1c2vArRWQ1sAmYoqr18QqJyERgIkC/fv2ojXnK9+zZsznSZFvs3r27pdzgwW4BSOLcZPCv/YUvfIFBgwbx29/+losvvpiLLrqIG264geOPPx5VZdu2bZx11llMmjSJwYMH85WvfIWKigqampo49NBDGTx4MEcccQQDBgzgmGOOobGxkW3btvHTn/6Uiy++mH333ZcRI0a0/j4p0NjYuEcbJktDQ0Pa5+aCtuTrUV9Pr7o6tgwbxtZBg3IrmEfU2w+iJ2OPG25ouW/bt0dOvliiLt8eqGooCzAAWJvgWA+gm7c+BngtmWtWVFRoLOvWrdtjXzy2bt2aVLl8kSv5km2veNTU1IQnSBbYQ76lS1VnzlS96y7VLl1Uy8rc59Kl0ZAvgkRdRpMvdYAVmkCn5qSHr6pbA+uLReRXItJXVT/KRf1GCRAc8BOBpia3hDRQb8Rg3jQFSU4UvojsB7yvqioiR+PGDjbnom6jRAgO+HXo4OzAIi22YFNQ4WEx7AuWsNwyHwKqgL4ishGYAXQEUNU7gdOASSKyC/gMOMN79TCMzPAVeZ8+Tvn4Suimm2Dz5paBP1NQ4ZEFV2cjN4Si8FX1zHaO3wbcFkZdhuHTo74efvrT+Eo+qIBmzTIFFSa+N43f7uZNUzAU/0xbo2jpVVfXWpFv3uym+MdiCipcLIZ9wWIK3yhYtgwblpwiNwUVPhbDviApjeBpIeOHRx40aBBDhw5l7ty5zbNiAf785z9z9NFHc/jhh3P44Ydz9913Nx+78sorOeCAA/ja177G4MGDWbRoEQDvvPMO1dXVDB8+nCFDhrB48eI96hs2bBhjx45t3n/BBRcwdOhQhgwZwmmnnUZDQ0MOvn102DpoUPKT6oIBvrIYTM9om9mzoaam9XpNjVuH1utGFkjkrxmFJap++HvvvXfz+vvvv68jR47UK664QlVV33vvPT3ooIN05cqVqqr64Ycf6pFHHqmPP/64qqrOmDFDb7jhBt26dauuW7dO+/Tpo7t379YLL7xQf/WrX6mqan19vX7+85+PW1+QTz75pHn90ksv1VmzZu1RpqT88JNh6dKc+ehHvf1Ucy/js8+q9u3rPp99VrVnT9UePVq2/WP5ki9VoigfbfjhWw8/Q/bdd1/uvvtubrvtNlSV22+/nXPPPbc5Imbfvn2ZPXs211133R7nHnHEEZSXl/PRRx8hImzd6qYrfPLJJ/Tv37/dunv06AG4h/Znn30WN15PUeL10HvUx52s3TYWryWvVFfD/Pkwfrzrzas679lrr4Vx49yx6mpXtqYGXrzpQ3sbCxGz4YfAoYceSlNTEx988AH19fWcc845rY6PGDGC+jjK6cUXX6RDhw7ss88+XHnllXz729/m1ltv5dNPP+WZZ55pLtfY2MiIESMoLy/n8ssvZ9y4cc3HzjvvPBYvXszAgQOZM2dO1r5jZAj4gA8tL4cjj0zNlmwDuHmnuhomTXKWuOnT3b5rroGuXVvK1NTA7FOWsehfE+DxXeZOGxIF38MXkbhLjx49Eh5rb0kH9aYVqGrcawT3/fKXv+RrX/saU6ZM4eGHH0ZEeOihhzj33HPZuHEjixcvZsKECc3jAu+88w4rVqzgd7/7HZdccglvvPFG87V++9vfsmnTJo444ggefvjhtGQvKAI9dNm5M/UeumVoyjlBuz249Vtucc/tm29269OnQ3k5nHIKXHGFewO47fRaynfvtLexECl4hZ/IVrV169ZM4gKlxJtvvklZWRn77rsvgwYNYsWKFa2Or1y5koEDBzZvX3rppbzwwgs8//zzHHfccQD893//N+PHjwdcFq3GxsbmjFe+eefQQw+lqqqKVatWtbp+WVkZ3/ve90ojqmYgoqJ27JheD90yNOWUo45qMeHU1Dilrgpjxjhzjqrr9S9c6PT6Nde4N4DDzq+iqWNHi54ZIgWv8PPNhx9+yEUXXcTkyZMRES6++GLuvfde6urqANi8eTM/+9nPmDp1apvXOfjgg1niJXxYv349jY2N7LPPPnz88cds374dgI8++ogXXniBgQMHoqq8/vrrgHvoPfbYYxx++OHZ+6JRIdBDXz1njintAiBot7/2WqfgFy6EXbtgwQK3vny5K9upk+v5z5kDc5dVunvsvY3VNFaaB0+GmA0/DfyMVzt37qS8vJwJEyZw2WWXAbD//vvzwAMPcOGFF7Jt2zZUlUsuuYSTTjqpzWvOmTOHCy+8kF/+8peICPfeey8iwvr16/nRj35Ehw4daGpq4vLLL2fgwIE0NTVxzjnnNL/JDB06lDvuuCMXXz//eD7gW+0Vv2CItdtXV7cMzvqMH+8eANXVMHcuTJkCF130HX71qy9QU+OOz5+fH/mLBVP4abB79+42jx9//PEs97ssMVx55ZUAe8S2HzhwIC+88MIe5Y899lhefvnlPfZ36NAhbnkjRSyoWnJk2E41NXDHHU7Z33HHngp/+fLWHjpe/4lf/OIQ+vZ15wSPG+lhCt8oXSzqY3Kk2E6zZzu7fdC98pRT4Hvfg6uvdvv93rpfJp7F87LLoK7uXa65ZkDzW4GRGWbDN0oX88lPjhTbKThICzBvnrPbn3GG2/Zt+glegpupqYFFi/o3vxUEPX2M9ChIhZ+OJ00pUjTtlK1QCJZDNTlSbKfgIO0VV8Cjj7qB2WAPvbo6fq/ex7fZz5ixjquvbj1Zy0ifgjPpdO7cmc2bN9OnT5/SmVmaBqrK5s2b6dy5c75FyYxsml0sqFpypNFO8QZpU8G36Ytsab6e/1Zgpp30KTiFf+CBB7Jx40Y+/PDDNss1NjZGWtnlQr7OnTtz4IEHZrWOrJPtZBsW9TE5Umyn9gZp28Pv/dfWth4TCI4LLF/e9luCsSdhZby6B/g34ANVHRznuAA34xKY/ws4V1VfSqeujh07csghh7Rbrra2luHDh6dTRU6IunyRwUIhFBxBF0pfSccO0qaCPybgn28umukTlg3/XmBUG8dHA1/0lolAiTiMGxljoRAKjlgXy2QHaRMROyaQycOj1AkrxeFzIjKgjSInA/d5oTv/IiK9RGR/VX0vjPqNIsfMLgVFPDNLqiadeOdnMiZgOCQsTw5P4T+ewKTzOHCdqv7Z214C/ExVV8QpOxH3FkC/fv0q5s2bl5Y8DQ0NdOvWLa1zc0HU5YPoy2jyZU5YMj700EEcfvg2hg/f0rxv1apevPJKd848892M5Vu1qhdXXTWQsWM3sWhRf2bMWNeqrnwRxXtcXV29UlVHxD2YboCxOAHHBgBrExz7I/D1wPYSoKK9a8ZLgJIsUUxMECTq8qlGX0aTL3PCkjE2eUm8ZCbpUFNT03ytFbcuVZ05U1fcujSUa4dBFO8xbSRAyZWXzkbgoMD2gcCmHNVtGMlhYRbSJmhnnzQp3FAIy5fDkzOWUTHVuedWdOrEk7OXsGR5pZl2UiRXCn8RMFlE5gHHAJ+o2e+NKGFhFjImW3b2qVOBWbWt3HMrttVSMc3uT6qE4qUjIg8By4Avi8hGEblARC4SkYu8IouBN4HXgV8D/x5GvUYRk+tE4xZmIWNife9DnRVrs6JDISwvnTPbOa7AxWHUZZQA+ehtm79/RoTte78Hnnvu/11dS7cTq6gI/B5sElbyFGQsHaPIyUdv2/z9MyJs3/u4VFbSNHUao66qbH578B80Rx0VYj1FTMGFVjBKgHz1ts3fvzUpDGJnw/c+HtkcHC4FTOEb0cOCmuWfCA9i2ySs9DGTjhFNLNF4fknCrDZ79p4DszU1ZD3vbFYHh4scU/hGyZAvBVWQJOEVE5voJBf29ODgsMXJTx1T+EZRE1TyvoKaOxfGjHGfQQVlyj9AEoPY+QhqlpPB4SLGbPhGURMbWnfaNJgyBb71Lfd5440WcjchSQxi59qenqvB4WLFevhhkeuJQkZcYs02vpI/6STXC501C846C55+2n3OmuWU/7hxrXuO1ttPDrOnFxbWww+DWI+Gm26CzZvNwyQPxEuWMWsWnHqq64VOmABPPNGioEaPhvvvh65dW65hvf3kyPpkKyN0TOFngu+n/M47LR4N27fD5MnQ1GTKPw/E89OeNs0p/QkT4IEHnBnnssugVy9n1pkwAf7wBzjlFPjxj905T85YRsVfaqFzld23BLRlTzeFH01M4adLsFdfVgblXlOKOMXf1LSn8o+QL3MkSTNaZTDnKbjP0aNbevSzZrUoohtvdNvgPm+8EXbtgvPOgxNPdOfcfV5LZEa7b4nJtz099r6DhVloD1P4qRBUSEE/ZYALL4SDD4Y+feCSS9yxoPLPRgLuYiKDiT6xZpy5c11PfsIEeOQRp8SDimj4cLjhhj1t9p06wQEHwKbf1dK0cwcdmpwP+hv31PLI85WmRCKG5bpNHVP4yRLPTh+c/n/22S0K6itfcco9qPwtIFfbxJvok6TCD5pxRo9ubbY57zy3f/jw1m8Asb3C8eNhwQK3PfOkKhq3d6Jzhx1oeScm/76KqQtC/bZGCFiYhdQxhZ8ssQpp8+bE0/+D7my+8vfLeG8JPXr0sAdAkAzj5wTdAydMcMre39+eXTnWFs1jlZx80hLO2K+W//2oiqkLLNGGT2TMKN7/qLqqikmTKi3MQrIkSoUVhSVSKQ6XLlXt0kW1rMx9Ll2a0TV27bVXetfIITlP37bUpbBLtl2C8vlp8KZPDye13vTpquA+0yWK6e9iSVXGbKUyTERc+Vr9j7roqJ5LQ7vvociXZ2gjxWFYCVBGicirIvK6iFwe53iViHwiInXeckUY9eYE378eMg+fG3hLkJ07LclGLCnEz5k92yXJhhaTzLRp0K1b5tPtzbc8MfmYXbsHgf9R0/Yd3HZ6rYVZSJKMTToiUgbcDpyAy127XEQWqeq6mKLPq+q/ZVpfTok3kDhtWvrXC5gttLzcTDoZcNRRcMopAxk2zJkTfNdLX/mk6x4Y61v+/vtuUtbCha0HeEvZEyTv0Sq9/9Huxh106NSJw86vapbL3ELbJowe/tHA66r6pqruAOYBJ4dw3fwTdiKOQHyS1XPmmMdOBlRXw4wZ6xg/HhoaWit7/3g6CjnWnn/GGc7Zat48t20JNyLwBuT9j8quvYaymtZv2+ne91JBnMkngwuInAaMUtUfetsTgGNUdXKgTBXwCO4NYBMwRVXrE1xvIjARoF+/fhXz/H9aijQ0NNCtW7e0zvXpUV/P0J/8BNm5E+3YkdVz5rB10KCMrhmUr//bb9Orro4tw4aFdt0wCaMNs0lDQwPz5w/m/vsHMGHCBs4/f0NW6lm1qhdXXTWQsWM3sWhRf2bMWMfw4VuSki/K7Qepy+i3hd8Gsdv5li/XRFG+6urqlao6Iu7BRMb9ZBfgdOA3ge0JwK0xZXoA3bz1McBryVw7b4O2wcHDFAcSk2XlbbdlPgicZaI4IBVk7txVoQ7UtkU6g7hRbz/VODK283u//vo92/nZZ93+nMgXMaIoH20M2obhlrkROCiwfSCuFx98qGwNrC8WkV+JSF9V/SiE+sMlbLt9AnrV1aXtd16qBF0Ca2rgqqsGcsUVbqasP2CXjQHEWBNG0UZnTGLyW75n1xqZEYYNfznwRRE5REQ6AWcAi4IFRGQ/ERFv/Wiv3s0h1B0+OUqgvWXYsHYTTBitCSbcWL4cvv/9d5g1q+UhkI246CWVcCMfyeONnJJxD19Vd4nIZOBPQBlwj6rWi8hF3vE7gdOASSKyC/gMOMN79YgeOUqgvXXQIMvbCinFz4mdWfm73x3MggWJZ9CGQUkFCMtX8ngjZ4Qy01ZVFwOLY/bdGVi/DbgtjLqyTi4TaCeRYKKoSSN+TusZtZuorh6QVRGbTRjLlvF/V9fS7cQqqidXFqeLZoEnj4/MLOAIY6EV4lHqijhXpBE/J2hPv+WW/tTU5KCn7T2Yjtu+g8YnO7GSJVRMrizOYF0F/Nu3YGrtYxmvfPKdsSrf9eeDJBJlB4m1p/t++Fm3p3sPpg5Nu+ncYQdP/Kw2f7NMc0whJX6PxCzgiGMKH1pMC9Onu89cK918158v2kmUHats/Bm1/sDs8OFbcpPAOvBg6rBXJ/p9r4prrnGmpWJXJsGBcoj+xLOgya8U7k+qmMKH/Hsn5Lv+fNJG/JxYZXPUUTR75fjkZGZl4MG0cvYSfv5YZcnE2Sm0XnPeZwFHHLPhQ/69E/Jdf0SJVLzzykpqGitLModr3mPnJInl2G0f6+FDu6aFMAiaJ/xIjzU1MGYM1DRWsnL2Ev7vBFd/TWNlJG2k+SBKr+htuWgWM4XSay7V+5MSiabgRmHJemiFLIVNUN1zCvqzz6r27Kk6caJb33vvndqjh+qcOW5/jx7u2Jw5e8Ybz9a09faIwrTxtuLc51O+ZEIMRKH92qM9GXMd/z6WqLdhFOUj2/HwC5IsD5TG2p8BVOHhh/19ighs2eL2i8Bnn8GUKc6kHXQri+oAWbaJ8izXQhvMTBfrNRcXpWvDzyCHaiKCEz/8P8Ypp8CIEbB6tYupXlPjTxr6OwMGDGi2i0JLer6nr17GVx6vZebSKqbNbJ1er9gnkgTb0Fc2/v6pU6MzyzVS4wtZxGLnFBel28NP0Qc8GeL16nfscMMCkya5bd8W+uijB3DLLW795ptpXv9w0TIWNoykumY6T+4aydNXLyv6XmSQYBv6yib4naMU7zxK4wslTSnOYUmT0u3hhzSNvK1e/fLl7lkyZUqLUl+wwJ03d66gCr16OXOOqjv/hJW1lC3eQTm70d07+PnXajl1fCVDh7rrFXvmpULqOZdMFM0ok0Z4jlKmtHr4sT2BFHKoJqKtXv2uXU7BX321y5zkh4tbvhyuuWYtCxfCM8+4MgsXuqxKM1+oQjt2Yre4ST4zX6hi9OiW6/kUbG8/id5YIfScozy+UFKU8hyWdEg0mhuFJVQvnUCm+0wTjiTywBk50nnbBNdjy/leHPFG95uvG/AemjNHde+9nZeKf+1cJPxIJGNGJHkP2vLMyap8KVDMXjq5TnLSFu22YYj/63SI4j2mDS+dvCv1tpZQFf7Mme5HAe5z5sy0rx3PVa1LF3fprl2Tc2FL5ofSVj0TJuxZNuw/ZOg/5iTuQSpugFH8swWJunyqnowx7sn5dsXcQ772iONenauHVhTvcVsKv3RMOhkM0sbGdKmudpagk05y083HjXOXHDkSystbl8vEhS3WJQ5cPcOHwwMPwNy5bl/BmHeSuAfmBphbetTX7+GeXGjhFOKZZkvFbTZlEj0JorCEPvEqzYlWiXo8Eyak1qtvV74UZJgzR1XEyZCt3ldWei8h9sby1ruK+Q6J5J848fU8CJcab/zwhwnfutLJ4xs2mdzjZE2DmVBoPfxQFDMwCngVeB24PM5xAW7xjq8Bjkzmuqkq/OAfz78RK25dqrWjUlfybdnp+/ZtmRHbnq0+Ean+UOIpFf+Bk60/ZK5+zOmaEPLyZ4tjM04k/9y5q3IvX4qsvO22uDbwXCjLZMj0Hmf7oZWrTlEqZFXh49IavgEcCnQCVgMDY8qMAZ7wFP9XgReTuXaqCt//ka64dam+8cMf6vrL7tJP6aK7O6Q+oNOe/TxTG2dKP5Q4PwC/znQfOKHLmCHpKJi8KPwE4xDx5I9i7y+WorDhJ6Age/ghDEJnW+FXAn8KbE8DpsWUuQs4M7D9KrB/e9dOx6Sz4tal+ilddBcddDvl2iQd4r6uxiMZ75vp050JZ86c1uemqmST/qG006P0ZfSVfph/zlwrrFR7Y1Hp4fvEyl8wCj+GgvLSSUCuHlr5cGxoj7YUvrjj6SMipwGjVPWH3vYE4BhVnRwo8zhwnar+2dteAvxMVVfEud5EYCJAv379KubNm5eSPAc/+CCf/809lNOU7lcyDMPIO7v22ovVc+awddCglM6rrq5eqaoj4h5M9CRIdgFOB34T2J4A3BpT5o/A1wPbS4CK9q6dUQ9fyvRTuuj6y+5qc5DQf/LPmeP2BwdDk/WpT4dMevjx8HuXI0eG1zszG35yFLINPyptmIiSlC+LNvww3DI3AgcFtg8ENqVRJmNqamDUVZWsv3UJi0Zcwh/+YwnH3TeRH22YRk1jZXMuzqOOcuEPfvSjFhfLKVNg7Vo3CfSss+D++1tmyj7zjJsJG3TzyllMlyRi9Qen+C9f7r5b1N3Rgq6usUHSCs0VM5Er6SuvdM+vYEZhEkIEgIQkehIku+Di8bwJHELLoO2gmDIn0nrQ9q/JXDsTL525c1c199wnTtzTzh07a9X3ePEHZLM5GKoaXs8gXu8yrBm52exdhWFjLcneX8hEXUaTL3XIZg9fVXcBk4E/AeuB+apaLyIXichFXrHF3kPhdeDXwL9nWm88pk5t6WX5Ca5nzYJ+/Vpizvs93oUL4cc/dh3n0aPhiSdcaOIHHnAP17z26lMgXu9y4UIXvC3KsWgKbnJPCsRO1AOa3y6jSiHKbKROKDNtVXWxqn5JVQ9T1Wu9fXeq6p3euqrqxd7xr2icwdpsEAzC9Z//2aLgg6GKg0p+8GC48Ub3kKipKQzTQvAhF2T1ajdx8uabW/+R8/InThAwrRCCpKWDP8tz1apeQHTNakFsZmqJkKjrH4Ul05m2QT/coJnDX/dNNPlIK5itV8EwXTZDkbGNQedM/aSj+Drt49p+e94nLrVHsA2jMtkqSJTvsWo05aMUY+msWtWrlZkgGHM+GKp46lS47LLWPfkomm6aaSe8cNDEU13tBp1F4Npr82Q2CYSvbdreEr62psYNLp96anGGF66uhrFjNxXU20uxvnEZLRStwn/lle7Nym358paY88uXw113taz7RFrJ+ySRhzfWxFNd7UxZftatnP+JvYBpTR3KaGzqxMruVYCL/a/qHr6+nFE3n6VCTQ0sWtS/OTlKITzIYhO6FILMyWDjEwESdf2jsIQePC1CpCVfGrPwMgnBEFoben7FK25dGqrJIKr3ONYPP5+hCdrDb8MohVMIEsY9zuZ3i+JvkDZMOqWb4rAQ8cML++nc2gnxHMzKBM6EMm6ce7uB1seySmUlVFZSAUz6gObE7cVqMvDNaiJbgNZvL1H9zm2FpY6qzAlZtqxV6tJCSpuZbUzhFxIp5uGN/RMvWOCU/rXXOi+eXP/oSyUHrG8aDGbbi/p3jWfOjLrMcUmQ4zY4PlHMnY32KFobftGSwiy8fNvzg7ZT/21j2jTo1q2IBmmTyNEbSTy5e9TX51uScEmQ47ZYxydSxRR+CeH/6HPlnx/07V6+3Cn7WbPc/qIYpE1iED2SA4YBuYf+5CeF97BqizhZ1SzhfAum8EuE4I/+F79wrprjxrn92ZpkE7SdNjQ4ZR9rJ468Z1RbJOhNBonkhKaA3LJzZ1y5C5Y4sacsbWYLZsMvEXJlz/eD0wX/XKNHF6ntNIlB9EgOGHpy727cQVNZOR0CcvtvYwX9IPacBHyKZnwiBKyHX8ikYD9uy54/dOie5Wtq4KGHDtrzQDvE9mjnznWhKyZMKELbaRKRTCGCE5o8uTdccA1jOj1FTaOTOxJvH0Z2SeSvGYXF/PDbIMNUaPHCToQVz/3ZZ1VH9VyqDw+bqZUsbc4OFrZvd6Hc4yiGLPDxo8pGUTbVwrnHUYJSDK1Q9CRhP05E7CDWwoVu1uspp8C3vuVs+/Pnu4ijfvm2BhljByarOy9jYcNITq2bTk3ZSC6rdG8gpWg7jfqA4fDhW6L19mFkFVP4hUocb4RkaSuk8pIlLvGLTzKv+bFmnOf/q5ay3TsoZzcddu/gjXtqm8sW/EBtikR9wHDVql7mrlhC2KBtoZLiJKwgiRTu6tXuj3/zza63f9hhQ3j99T3jDt1wA/z0p27dH+CbNg1OOskFQ3t9cRX/17ETNO2gQ3knJv++iqnfL83eY5QHDGtq4KqrBrJgQYtMxZSXoKCImR2cLUzhFzIx3gjpEjQ7+H/8E0+El17qTdeuroyfFlIVZsxoWV+40J0/a5ZT9vffDxMmVNJxknsYlVVVMbWxsjCn6Bc5y5fDjBnrqK4eBhR4OIVCJsHs4GyQkUlHRHqLyNMi8pr3+bkE5TaIyMsiUiciOUl+YiRPrNkB3O/uyCP/SXl5S45c9bKGbdmyZwaxadNc1rDp091nTWPLjOBSM+MkImqTsKZObRmn8Sn2exW1ewBkNB6XKpna8C8HlqjqF4El3nYiqlV1mKqOyLBOI2SCLpu+Al+wAObMWcPChe43GJs1LLg+enTLpKooDkxGhUhOwioxInkPMhiPS5VMTTonA1Xe+v8AtcDPMrymkUeCvX2/o9GpExx7rLPti7TY+f31OXOc4i+KSItZJJKTsEqMSN6DDMbjUkXUT/2UzskiW1S1V2D7Y1Xdw6wjIm8BHwMK3KWqd7dxzYnARIB+/fpVzJs3Ly3ZGhoa6NatW1rn5oKw5etRX0+vujq2DBvG1kGDQrnm0qV7MXt2BTNmrANg+vTBgHL22W9z330DAOWaa1zwrauuGsiMGev2MBFkk0K9x/fcM4D77x/AhAkbOP/8DTmV6aGHDuLww7c136eGhgZee+1AXnmlO2ee+W5OZUmGbN3jsO5BmPLF3htwXlSp3pvq6uqVCS0piRz0/QV4BlgbZzkZ2BJT9uME1+jvfe4LrAaOb69etYlXyZPhJKxETJz4evNEnOuvb8mLO3p0y7qfQCUXeYBjKcR7nO9JWLGT3/yJV1GbcOWTjXsc5j0IU76wErXQxsSrjGbCAq8C+3vr+wOvJnHOlcCUZK5vCj9J0siElQwl1YZZIFa+qGSVCiq8nj23R1bZq4Z/j8O+B9mSL5OHUVsKP9NB20XAOd76OcAfYguIyN4i0t1fB77tvSEYYZHDQZ82KdTY8GGRQoJ5yN8krGBsn7FjN5XUGEJU7kEish13KdNB2+uA+SJyAfAOcDqAiPQHfqOqY4B+wAIR8ev7nao+mWG9RpAcDvokJIe+xJEk3vePIejuGIwqGvSQykWkymAykFtu6U9NTREPHMdMaIryRDjIfla4jBS+qm4GRsbZvwkY462/CcSJx2iESkiTsNImni9xKSn8FL+/7x7o9zZj8w+Hjqf4VnavYvxVlc31fu5z6xg/flj+PVWyQYF1QuJNgAx75rPF0jHCISpmpXyR4vcPugdecUWWQxoEMlwNuWwkT85Y1lzP8OFbImXSCJUcTmgKg1yYmyy0ghEOUTAr5ZN4378dBZOzxNoBxdeRHVRsqwVa7k+UTBqhkkSCmiiRC3OTKXwjPPJtVso3KX7/bNtrwY0VjOxeRUVA8a3sXsWS2cUdQgGwTkgcTOEXKzmKvmekRy7steDGCkaNr+TJ2Uuo2OZs+KM8G35JUOqdkBjMhl+MBGy2jBxZum6SESZX7oH+dUddVckVn01rVvZFacJJg0gGU8sipvCLkQIbrCpFYnMM+womaGYJS/FELqduhIhkMLUsYgq/GMmVx0ypT7QKkWwqntixAoti2kJOvaUigNnwi5FcDFYVmI9z1MlWFMdcjRUUMjnzlooA1sMvVipbEpBkBTMbhU5YppegXdofK/D3Ry2UQBQopTcgU/hGepT6RKssEJbiCZqH/DGBoHmo2LNapULwDSjnyXvyYBI1hV8KZOOH5ZuNrrnGzDkhEFQ83bq5l7NYm36yA7ilZpfOhLwFU/NMok2/mM7u6taedNn0EjIbfrGTTVu7+Ti3SY/6etf+SYyjxCoeP0+wr3hSjbNTSnbpTMhbMDXPJNpBd7Nz+w423FPLYZWVWY+pZAq/2Cn1oGb5Ytkyhv7kJ7BrV1IP2qDiiR3ATaeHnotZvEYGBMI+dCjvxOTfV3HU/tlPuWgmnWInbFu7uWImR20tHXbuTHtQO9UB3OBArd9LnDbNmYcsqXwECZhEy2qWcNSPK3MyT8IUfrETpq3dZvAmT1UVTR07pv2g9XvoI0e6hPFBZR3PxhscqF2+3Cn7p69exndfm0V152XmmRNFPE+6msbKnHkJZWTSEZHTcSkLjwCOVtUVCcqNAm4GynCJUa7LpF4jRcKytZt5KHkqK1k9Zw5Hbt2a8lyIWDvuKafAuHGwcKHb9o/FJlGZP9+VHTECOq5YxuONIyn77Q54sBPVS5ZQPdXuVdTI9TyJTHv4a4FTgecSFRCRMuB2YDQwEDhTRAZmWK+RLpmYZMwVMyW2DhqU1lyI4ABudTUsWAAiMHGiU/z+saOOcgr+Rz9qOXfHDvcid9mRtZTtsnkSUSfXXkKZZrxaD+ClL0zE0cDrXuYrRGQecDKwLpO6jTTI1GPHws3mhFjPkepq+PGPnVWua9fWx1Th4YehXz9n+unUCaZMgbm3VPHN8k6UURix4HNCOxFkg29MPtlOO5lrLyFxSc4zvIhILTAlnklHRE4DRqnqD73tCcAxqjo5wbUmAhMB+vXrVzFv3ry0ZGpoaKBbt25pnZsL8iHfwQ8+yCH33IM0NdHUoQMbzj+fd37wg4TlfRl71NfTq66OLcOGuV5rRCiVe7xqVS+uumogY8du4tFHDwCEU0/dyKJF/ZkxYx2rVvXi/vsHsNdeu5k162WGD9/CqlW9WDx9KzO+8Xu6jjk84X0rlTbsUV/P0J/8hA47d9LUsSOr58zZo038dp4xY11zGwa3sylfmFRXV69U1RFxD6pqmwvwDM50E7ucHChTC4xIcP7pOLu9vz0BuLW9elWViooKTZeampq0z80FeZFv6VLVLl1Uy8rc5113qc6c6fbHoaamZs9zEpTNB6Vwj599VrVvX/fpb3fpogqq06e3HB85UrVHj5Zyftnrr8++jNkkNPlmznS/YXCfM2fGLea35/Tprds96/KFCLBCE+jUdk06qvqtzJ43bAQOCmwfCGzK8JpGOgRNMn36wCWXtG/esYHavBJr4wV3u4491plwbrnF2fhjE6H7ZgHzvfdIMt1hsU9Yy4Vb5nLgiyJyiIh0As4AFuWgXiMeflC1zZv3VOTBAd1lyzj4wQfdg8EGavNGMG6+r9AXLIBnnoEzznA2fB8LjNYGSbonF3sgtUzdMk8BbgX2Af4oInWq+h0R6Y8z44xR1V0iMhn4E84t8x5Vrc9YciMzYns8ffq0DOiWlYEIh+zcCQ8+CDfd5B4QNlCbV2J7+3fd5ZT+8uWtvTyKrVcaGu24J+fKRTIfg8M+mXrpLAAWxNm/CRgT2F4MLM6kLiNkYj1ugqabpiYARNXt27zZvRUYeSVvcV9KhLZcJMPOMxx8kGQ7fk4Qi6VTysT2ePwev9fDb9q5kw5mxjFKhFw9ULOV7CYZTOEbjtgeP7Dhnns49PzzzYwTBu34gBulRb4Gh03hGy3E9Pjf2b6dQ005ZY6lgzRiyFc0UwueZhjZxtJBFhzB6KM+GScm8bzgVt62LG9ZtkzhG0a2yWUMIgtfHQrB6KPQMrDqp4lMmUCk2SGXjeTJGctyn2ULM+kYRvbJVQwiMx2FRugDq4G3vI7soGJbLdByb3Jl0jGFbxi5IBfpIG1WdKiEOrCa5EzfbGMmHcMoFix8daiEOus2zEREGWA9fMMoFix8dWiEOeu2ZWZty1termbWxmI9fMMoJvxYSabs26adwe0wE5OEPgCcAdbDNwyjtEhicDvMWbf5nFkbi/XwDcMoLVKcFxGGT35wAHjSpPzFPzKFbxi5xnzl80uKg9thmGSiEnbZTDqGkUvMVz7/pDi43ZZJJplJ07kKu5wM1sM3jFxiYRaiQYqD25mYZMIcAM4UU/iGkUuy4StvJqKsk6pJJmj397OWBe3+1dW5d8mEzDNenQ5cCRwBHK2qKxKU2wBsA3YDuzRRRnXDKHbC9pU3E1HWacskIxL/nHwmOWmLTHv4a4FTgeeSKFutqsNM2RslT5i+8mYiyjpBk4zfQw+aZOJ57ATt/ldckT+bfSwZKXxVXa+qr4YljGEYKWLhFLJOMJG833P3969a1avZYyfWfbO6Gv5jxDJ2XjOLmScty7uyBxANpr1P9yIitcCUNkw6bwEfAwrcpap3t3GticBEgH79+lXMmzcvLZkaGhro1q1bWufmgqjLB9GX0eRz9Kivp1ddHVuGDWProEEpnWttmDqrVvXiqqsGMnbsJubPP5Dzz9/A+PEbm/d///vvsHu38IUPX+LiBeewFzvYTicW/fgO9jvl81mXr7q6emVCS4qqtrkAz+BMN7HLyYEytcCINq7R3/vcF1gNHN9evapKRUWFpktNTU3a5+aCqMunGn0ZTb7MibqMUZVv+nRVUD3hhPe0b1/VZ591++fMURVRPeEE1WnM1N1Spgq6u0OZXtN1ZnO5bAKs0AQ6tV2Tjqp+S1UHx1n+kOwTR1U3eZ8fAAuAo5M91zCKGvOwKTiCHjt//Wtvpk1rsdXPmgVnnQVPPw2dR1XRobMzt3XYqxOjr6/KiytmkKxPvBKRvYEOqrrNW/82cHW26zWMyGMeNgVHrMfO5z63jpkzhzF6tPPRnzABnnjCPQxuvaOSE2cvcclOqqqoqKykIs/yZzRoKyKniMhGXOqWP4rIn7z9/UVksVesH/BnEVkN/BX4o6o+mUm9hlEUpOthY28F4ZNkm8ZOoho+fAvTpsGjjzpl/8ADzgHLz1U76qpKar4aneilGfXwVXUBzkQTu38TMMZbfxMYmkk9hlGUpJMFyd4KwieFNo2dLLVqVS9mzYLHHnMPgxtvdM+N4cNbz6iNgocOWCwdw8gf6UzCsjSG4ZNBm77ySvdWE7LAKXtfyecqV22ymMI3jHySaq7biORGLSoyaNMzz3yXqqrDWu2LmpIPYgrfMKLEsmVt9/gtjWH4lFCbmsI3jKjQli059kFQxEopL5RIm5rCN4yoEGtLvu8+t69PH7jkEhuoNTLGFL5hRIWgLbmsDH77W9i1y4VkbGpyiw3UGhlgCt8wokLQlvzOO/DrX7vefocO7gEgYgO1hUB74zB5xBS+YUQJ35a8bBn8z/+0mHFuugk2b46kEjECRHyehCl8w4giJeQ5UlREfJ6EKXzDiCol4jkSOTIxyUR8noQpfMMwDJ9MTTIRfzMzhW8YhuEThkkmwm9mmea0NQzDKB6KPGWk9fANwzB8Im6SyRRT+IZhGEHSMclE2Pc+SEYKX0RuAE4CdgBvAOep6pY45UYBNwNlwG9U9bpM6jUMw4gMEfe9D5KpDf9pYLCqDgH+BkyLLSAiZcDtwGhgIHCmiAzMsF7DMIxokG7msjyQkcJX1adUdZe3+RfgwDjFjgZeV9U3VXUHMA84OZN6DcMwckIyqQ8LaKA3TBv++cDDcfYfALwb2N4IHBNivYZhGOGTrKmmgAZ621X4IvIMsF+cQ79Q1T94ZX4B7AIejHeJOPu0jfomAhMB+vXrR22ar0cNDQ1pn5sLoi4fRF9Gky9zoi5jPuU7+MEHOWT7dqSpiabt29lwzz28s307AD3q6+lVV0f5l79Ms3SVlbB9e6RNOqhqRgtwDrAM6JrgeCXwp8D2NGBaMteuqKjQdKmpqUn73FwQdflUoy+jyZc5UZcxr/ItXarapYtqWZn7vOsu1Zkz3ae3f9dee7lyEQJYoQl0aqZeOqOAnwHfUNV/JSi2HPiiiBwC/B04A/h+JvUahmFknaCpJpiEJpCfQFQjFyCtLTL10rkN6A48LSJ1InIngIj0F5HFAOoGdScDfwLWA/NVtT7Deg3DMLJPZSVMm+ZCU/ueOE1NboC2rAzt2DHSg7SxZNTDV9UvJNi/CRgT2F4MLM6kLsMwjLwRGwXTy0+wukcPjiyQ3j3YTFvDMIz2SeCJszXKA7RxMIVvGIaRDBGOgpksFi3TMAyjRDCFbxiGUSKYwjcMwygRTOEbhmGUCKbwDcMwSgRT+IZhGCWCuNAL0UREPgTeTvP0vsBHIYoTNlGXD6Ivo8mXOVGX0eRLnc+r6j7xDkRa4WeCiKxQ1RH5liMRUZcPoi+jyZc5UZfR5AsXM+kYhmGUCKbwDcMwSoRiVvh351uAdoi6fBB9GU2+zIm6jCZfiBStDd8wDMNoTTH38A3DMIwApvANwzBKhIJW+CIySkReFZHXReTyOMdFRG7xjq8RkSNzLN9BIlIjIutFpF5E/jNOmSoR+cTLGFYnIlfkWMYNIvKyV/eKOMfz3YZfDrRNnYhsFZFLYsrktA1F5B4R+UBE1gb29RaRp0XkNe/zcwnObfM3m2UZbxCRV7z7uEBEeiU4t83fRBblu1JE/h64j2MSnJv1Nkwg38MB2TaISF2Cc7PefmmTKNlt1BegDHgDOBToBKwGBsaUGQM8AQjwVeDFHMu4P3Ckt94d+FscGauAx/PYjhuAvm0cz2sbxrnn/8BNLMlbGwLHA0cCawP7ZgOXe+uXA9cnkL/N32yWZfw2UO6tXx9PxmR+E1mU70pgShK/gay3YTz5Yo7PAa7IV/uluxRyD/9o4HVVfVNVdwDzgJNjypwM3KeOvwC9RGT/XAmoqu+p6kve+jZcTt8DclV/SOS1DWMYCbyhqunOvg4FVX0O+GfM7pOB//HW/wcYF+fUZH6zWZNRVZ9Sl2Ma4C/AgdmoOxkStGEy5KQN25JPRAQYDzwUdr3ZppAV/gHAu4HtjeypTJMpkxNEZAAwHHgxzuFKEVktIk+IyKDcSoYCT4nIShGZGOd4ZNoQOIPEf7J8tiFAP1V9D9yDHtg3TpkoteX5uDe3eLT3m8gmkz2T0z0JzGJRaMPjgPdV9bUEx/PZfm1SyApf4uyL9TFNpkzWEZFuwCPAJaq6NebwSzgTxVDgVmBhjsX7mqoeCYwGLhaR42OOR6UNOwFjgd/HOZzvNkyWqLTlL4BdwIMJirT3m8gWdwCHAcOA93Bmk1ii0IZn0nbvPl/t1y6FrPA3AgcFtg8ENqVRJquISEecsn9QVR+NPa6qW1W1wVtfDHQUkb65kk9VN3mfHwALcK/MQfLehh6jgZdU9f3YA/luQ4/3fVOX9/lBnDJ5b0sROQf4N+AH6hmcY0niN5EVVPV9Vd2tqk3ArxPUm9c2FJFy4FTg4URl8tV+yVDICn858EUROcTr/Z0BLIopswg42/M0+Srwif/anQs8W99/A+tVdW6CMvt55RCRo3H3ZHOO5NtbRLr767hBvbUxxfLahgES9qry2YYBFgHneOvnAH+IUyaZ32zWEJFRwM+Asar6rwRlkvlNZEu+4NjQKQnqzWsbAt8CXlHVjfEO5rP9kiLfo8aZLDgPkr/hRu1/4e27CLjIWxfgdu/4y8CIHMv3ddzr5hqgzlvGxMg4GajHeRv8BTg2h/Id6tW72pMhcm3oydAVp8B7BvblrQ1xD573gJ24HucFQB9gCfCa99nbK9sfWNzWbzaHMr6Os3/7v8U7Y2VM9JvIkXz3e7+xNTglvn++2jCefN7+e/3fXaBsztsv3cVCKxiGYZQIhWzSMQzDMFLAFL5hGEaJYArfMAyjRDCFbxiGUSKYwjcMwygRTOEbhmGUCKbwDcMwSoT/DxrmECXr9gUvAAAAAElFTkSuQmCC\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.9.5" } }, "nbformat": 4, "nbformat_minor": 5 }