{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# A script to illustrate each of Euler's method, the improved Euler method, and RK4\n", "# for a system (also works on a scalar equation).\n", "# Key points:\n", "# The spatial variables are y[0],...,y[d-1] for a system with d dependent variables.\n", "# t0 is the initial time, T is the final time\n", "# The vector \"f\" is an expression involving t,y[0],...,y[d-1] that defines y'=f(t,y)\n", "# Define f as below (scalar and vector examples shown).\n", "# Solution times are returned in row vector \"tv\"\n", "# Solution values are return in matrix \"yv\". Rows correspond to times, columns\n", "# to variables y[0],...,y[d-1]." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# SYSTEM EXAMPLE: Define t and y, define f in y' = f(t,y). Initial condition is y(t0) = ic.\n", "var('t y');\n", "def f(t,y):\n", " return vector([-sin(t)*y[1],-y[0]+cos(y[1])])\n", "ic = vector(SR,[1.0,2.0]);\n", "\n", "#SCALAR EXAMPLE: Uncomment to test.\n", "#def f(t,y):\n", "# return vector([-y[0]])\n", "#ic = vector(SR,[1.0]);" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Set solution interval t0 <= t <= T, and stepsize.\n", "t0 = 0.0;\n", "T = 5.0;\n", "stepsize = 0.1;" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "#Call routine for Euler's method\n", "load('euler_sys.sage');\n", "[tv,yv] = euler_sys(f,ic,t0,T,stepsize)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "#Call routine for improved Euler's method.\n", "load('improved_euler_sys.sage');\n", "[tv2,yv2] = improved_euler_sys(f,ic,t0,T,stepsize)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "#Call routine for RK4 method.\n", "load('rk4_sys.sage');\n", "[tv3,yv3] = rk4_sys(f,ic,t0,T,stepsize)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Euler Method\n", "time 5.00000000000000\n", "y[ 0 ] = 1.84279309134410\n", "y[ 1 ] = 2.20217591665613\n", "Improved Euler Method\n", "time 5.00000000000000\n", "y[ 0 ] = 2.11492279933869\n", "y[ 1 ] = 2.30603856534005\n", "RK4 Method\n", "time 5.00000000000000\n", "y[ 0 ] = 2.12145709039255\n", "y[ 1 ] = 2.32454788276000\n" ] } ], "source": [ "#Print solution values at final time.\n", "N = len(tv)-1;\n", "d = len(ic);\n", "print('Euler Method');\n", "print('time',tv[N]);\n", "for j in range(d):\n", " print('y[',j,'] = ',yv[N,j]);\n", "\n", "print('Improved Euler Method');\n", "print('time',tv2[N]);\n", "for j in range(d):\n", " print('y[',j,'] = ',yv2[N,j]);\n", "\n", "print('RK4 Method');\n", "print('time',tv3[N]);\n", "for j in range(d):\n", " print('y[',j,'] = ',yv3[N,j]);" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAGSCAYAAAAsKv9nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAAByBElEQVR4nO3dd3xP1x/H8dfNRCRB7L1HzNqr9lajOlTVrFGrlKoquhW/VotK1ahRraKKqll716b2jlGrVkRC5v39cUmFhIR8v9+M97OP7yPNveee+/lKyCfnnvM5hmmaiIiIiEjcOTk6ABEREZGkRgmUiIiISDwpgRIRERGJJyVQIiIiIvGkBEpEREQknpRAiYiIiMSTEigRERGReFICJSIiIhJPSqBERERE4kkJlIiIiEg8KYESERERiackm0AZhpHZMIxscWw70DCMF2wdk4iIiKQMSTKBMgwjL+AH3I7jJV8DPQzDaGSzoERERCTFMEzTdMyNDcMZGAy8CtwA0gLLgM9N07z7mOs877V7zTTN8zGcrwy4m6a5/qHjGYCNQB3TNC8n2BsRERGRFMeRCdQEoBlQyTTNfwzD8ALWAxeBpmYsgRmGMRnYbprm5FjOXwLWm6bZOoZzw4BCpmm2T6j3ISIiIimPQx7hGYZRDXgLGGWa5j8ApmneAj4CGgOPJD/3rnsOqA1Mi+V8MSALsCWWW48HmhmGUfCZ3oCIiIikaI6aA9Xp3scFDx1fBoQAXWK57j1gmmma4bGcr37v44aYTpqmeQNYA7wZ91BFREREonNUAlULuP3wHCbTNMOA40DVe3Okotx7xNcSWPqYfp8HAoB9j2mzAXgx/iGLiIiIWBI0gTIMI4NhGFMMw1hgGMavhmG4PnR+lGEYvwF5gVuxdBMIpMZ6FPegekAE8PdDfb5mGMZOwzB2Am8AYcD2e8cqxdD/BqCIYRgP9y8iIiISJwk9AvUpMPTe62UgqmyAYRgG0BVIDzhjJUoxuX88w0PHqwI7TNOMePCgaZqzTdMsjzWqZABfmaZZ/t5rWwz9+9/7WDqub0pERETkQQmWQN2bwH3RNM1LWMkOwJUHmpTESp7W3Ps8tlIFYfc+ej90vBjwz2NCqH3v4/rHtAG4iTWSlfcJ7URERERilJAjUJmAmff+vz1wAtj+wPka9z7GOMH7AWnvfXw4wcqBlfzEphYQBOx8XOf3yiME8GiCJiIiIhInLgnVkWmaGwAMw8iNtRpu6EO1nGpgzXvaDAQTe/Lmee/jjYeOp8VKfGJTC9j8mBV6DwoB3OPQTkREROQRtliF1+rex3kPHa8BbLg3h+k01uO8mGTESnDOPHQ8AnCL6QLDMHIB+Xjy47v7MhD3bWBEREREorFFAlUBay7U0fsHDMMogrWq7v78p2VATsMwoiVRhmGkA3IDWx+eLI71+O7hieX3PTL/yTCMd+5NXI/GMAwPrNGnx82nEhEREYmVLRIoHx4dPap37+Paex8X3vv40kPtXrn38ccY+vUn9gSqAtYI1Q4AwzDyA3lj2Q7mfvmCozGcExEREXkiWyRQO4Hc9wthGoZRCvgMuM69ApemaW4GFgEfG4aR5167zMD7wDpgegz9/o31mC4m14GbpmmG3iu4OQL4OJa2VbDmYh2I5/sSERERAWywmbBhGGmA74BsWCvxbgM9gWWmab76QDt34EOgOdbk8EzAKmDwvX3xHu73eazkysc0zZsPncsAzAGuYZVB+Mw0zWOxxDcR8DJNs80zvVERERFJsRI0gbo35yiVaZp3HjjWEmvPu2amaS5+hr6dgAtAl6ft514fJ4Depmk+bksYERERkVgl9CO85cDlexO17ydU7wErniV5AjBNMxKYArz2DN28jFUratmzxCIiIiIpW0InUBWwimfeuTcH6pt792idQP2PBRoZhhHbXKgneR/oG8vkchEREZE4SehHePWB+kAaIDNWMjXWNM2wx14Yv3u8BVQ0TbNzPK8bjDV/6t2EikVERERSpgSfRG4PhmHMANaYpjkjju0bAm2AN2OoLyUiIiISL7YoY2APbwKV721gHBcXgU5KnkRERCQhJMkRKBERERFHSqojUCIiIiIOowRKREREJJ6UQImIiIjEkxIoERERkXhSAiUiIiIST0qgREREROJJCZSIiIhIPCmBEhEREYknJVAiIiIi8aQESkRERCSelECJiIiIxFOiTqAMi5dhGIajYxERERG5z5GbCcd6Yz8/P/z8/IiIiODYsWMEBATg5eVlz9hEREQkaXDIIEuiTKDuu3XrFt7e3kqgREREJDYOSaAS9SM8ERERkcRICZSIiIhIPCmBEhEREYknJVAiIiKSaJgmfDshnIAAR0fyeEqgREREJNH4cdFN3u7pwvgtux0dymNpFZ6IiIgkCqGhkLnEZYLyHeTi8tJkNHzicplW4YmIiEjKNfS7CwSczMiA0Rfimjw5jEagRERExOGuXTfJWjAQj1eX8u/3L+GKa1wv1QiUiIiIpEwdPzlNeDh8+0mG+CRPDmPTBMowjMGGYewwDCPQMIwrhmEsNAyjiC3vKSIiIknL/qOhLP4uFwWHzOWNLPUdHU6c2HoEqibgB1QG6gMuwJ+GYXjY+L4iIiKSRLQZeA5ynmd23yoYjnkiF28utuzcNM1GD35uGEYn4ApQzpb3FRERkaRh/uoADv5RgHpzJlMuVVdHhxNnNk2gYuB97+N1O99XREREEpmICOjWPxDnqof5+ZWWjg4nXuyWQBmGYQBfA5tM0zwQU5uQkBBCQkKiPr9165adohMRERF7Gz7tAtf+zkmfv9aT2ajs6HDixZ6r8MYDpYA2sTUYMWIE3t7eUa9cuXLZLzoRERGxm8BAGD40FWlf/50vK73s6HDizS4JlGEY3wLNgdqmaZ6Prd3gwYMJCAiIep07d84e4YmIiIiddR15gtCA1Hw9IhXuuDs6nHiz6SO8e4/tvgVeBGqZpnn6ce3d3d1xd096f4giIiISdyfOhDF3dC7yvPsrXXK3dXQ4T8XWc6D8gNeBFkCgYRhZ7x0PcGAFdBEREXGg1wafxkzvyaz3SyWZsgUPs3UC1ePex3UPHe9k4/uKiIhIIvTnX7fY9Uthnp8yg6ppOzg6nKemvfBERETELkwTclQ7zeXgQM7tykJ25ywJ0a32whMREZHka8zcC1zcmo/OXx9MqOTJYTQCJSIiIjZ39y5kKHoZSv/N9d+fJxWpEqprh4xA2bsSuYiIiKRA/Scc5875fIxZEZ6QyZPDaARKREREbCoo2CR9/ut4v7CJK1OaJ/TKO82BEhERkeTn3YnHCbvqzagPvOOWPP32GzywtVtipARKREREbCYo2OSHUT5k6rCUTvlrPvmC1avh5Zdh6VLbB/cMlECJiIiIzQyYdNQafRoSh9GniAjo3x+qVIGWLe0S39PSJHIRERGxieA7JlNHZSRT++V0zN/0yRdMnw5//w1bt4KRuCuUawRKREREbGLApCOE/ZuOUUO8njz6FBgIQ4ZAmzZQubJ9AnwGiTKB8vPzw9fXlwoVKjg6FBEREXkKwXdMfhiZiUztVtCxwPNPvmDUKAgIgJEjbR9cAlAZAxEREUlwb409zMQBhZh6ZCudCj4hgTp7FooUseY/DR8e31s55FmfEigRERFJUMF3TNIVuEq6Bju4PL3xkx/ftW1rrb47fhw8PeN7O1UiFxERkaSv/+TDhF0pzJdD47Dybts2mDULJk9+muTJYTQCJSIiIgnmzl0T7/xXSd9gB5eeNPpkmlC9Oty+Dbt3g7Pz09xSI1AiIiKStL0z+RBhl4vyvyHpnjz69OuvsGULrFz5tMmTw2gESkRERBJE1OhT/Z1cmtHo8QnU3btQrBiUKAF//PEst9UIlIiIiCRd90efvhwah9GnsWPh/HlYvtw+wSUwjUCJiIjIM7tz18S7wL+kr7ubSz82fHwCdeUKFCwIHTvCuHHPemuHjEAlykKaIiIikrS8M+UgYZd84jb69OGH1pynjz6yT3A2oEd4IiIi8kzu3DWZOiIzmV9fTbvC9R/f+MABq2TBV1+Bj499ArQBjUCJiIjIM3nnhwP3Rp/iUPfp3Xchf37o1cs+wdmIRqBERETkqUWNPrVZQ/siTxh9WrYMVqyABQvAzc0+AdqIRqBERETkqb0zdT9hFzPy5VDvxzcMD4cBA6BmTWjRwj7B2ZBGoEREROSphEeYTP/Kh0yvrKd90TqPbzxtGhw5Aj/9BIZDFs4lqEQ5AuXn54evry8VKlRwdCgiIiISixELjhByOgdD3nV/fMOQEPjsM3j1VShb1j7B2ZjqQImIiEi8mSb4VDlKSKoAAteVx+lxYzJ+fvD223DwIBQtmtChqA6UiIiIJA2zt5zlxrYidHk34PHJ05078MUX8PrrtkieHEZzoERERCTePvjqGs5FQxjZpMbjG06cCJcvW8UzkxGNQImIiEi8bD1+Ff/fS9O8/3FSOz1m/lNQEIwYAe3bQ6FC9gvQDjQCJSIiIvHS55tTkCkSv3ZVHt/wu+/g+nUYNsw+gdmRRqBEREQkzs5cDWLXtJJU7r2TbKnSx94wMBBGjYLOnSFfPvsFaCdKoERERCTOek04AIaJX48Sj2/47bdWEjVkiH0CszMlUCIiIhInQXcjWPZtQQp23EzZjLljbxgQYG0W3LUr5H5MuyRMCZSIiIjEyYCZe4m8mp7/vZPt8Q3HjIHgYPjgA7vE5QgqpCkiIiJPFBFpkra4P2mL/sO/C6rH3vDGDcib15r79M039ghNhTRFREQkcfpq6SHuHsnHoHedH99w9GgIC4P337dPYA6iESgRERF5ooy19hMcEsHtLaVxim0z4KtXrRV3PXrA//5nr9A0AiUiIiKJz7yd/lxbX5JO716PPXkC+PJL6+N779knMAdSIU0RERF5rPdHX8E5P3zZslrsjS5fhvHjoV8/yJjRbrE5SqIcgfLz88PX15cKFSo4OhQREZEUbeeZfzn5a1mavHOMNM6P2bZl1ChwcYEBA+wXnANpDpSIiIjEqnL/zWybXoxz5wxyesRSefzCBShQAAYNgo8/tmt8aA6UiIiIJCb/3Axi2+TSVOixK/bkCawNg1OlgnfesV9wDqYESkRERGLUa9JeCHVlfO+isTc6exYmTYJ33wVvb7vF5mhKoEREROQRwaHhLB6Xn3xvbKFitlyxNxw5Ejw94e237RdcIqAESkRERB4xaM4uIv7Jxsj+mWNvdOkSTJ1qrbzz9LRbbImBEigRERGJxjRh6tfpyNB4G68WLx57wzFjwNUVevWyW2yJhRIoERERieaHLYcJ3luEPn0jY2908yZ8951VdTz9YyaYJ1NKoERERCSa4X43cSl4miH1K8be6LvvIDQ0Ra28e5ASKBEREYly4NJV/OeVo1GvU7g6xbJxcHCw9fiuY0fIls2e4SUaSqBEREQkyjuTD4FrGGM7Phd7o6lT4do1GDjQfoElMkqgREREBIA7YeGs/b4wRd7YQf50GWJuFBZmbRrcurVVfTyF0mbCIiIiAsCwhbuIuFCJT3rdiL3R7NlW8cw//rBfYImQ9sITERERANLX2ktkJARsKBNzg8hIKFkS8uWDxYvtGttjaC88ERERcYwF+09wc30ZOva+HXujP/6AQ4dg8GD7BZZIJcoRKD8/P/z8/IiIiODYsWMagRIREbGx4j02cOT3wtw6kwEPV7dHG5gmVKkCbm6wYYP9A4ydQ0agEmUCdZ8e4YmIiNje2YAA8uRwpubAHaz7qHbMjdauhTp1YOlSaNzYvgE+nh7hiYiIiP31m7EHQtz5plux2BuNGAGlS0OjRvYLLBHTKjwREZEULDwyksV+ucn10g6ey1Y15ka7dsHKlfDLL2A4ZMAn0dEIlIiISAr2v1W7CTuWn0G90sbeaMQIq+bTyy/bL7BETiNQIiIiKdg4v3BSlTpGj+olY25w9CjMnw/ffw8uShvu0wiUiIhICrXB/yyX/6jIq72v4BTbo7n//Q+yZoUOHewbXCKnBEpERCSFevf7ExhegXz1etmYG5w/DzNnQv/+4O5u3+ASOSVQIiIiKdD1u8HsnFKasp33kckjTcyNRo+GtGmhe3f7BpcEKIESERFJgd6dswPzmg+jeuSNucHVqzBpEvTuDZ6edo0tKVACJSIiksJEmiazv81Ipka7qFsod8yNvv3W+vj22/YLLAmxaQJlGEYNwzD+MAzjgmEYpmEYLW15PxEREXmyydv3cWdXcXr3imXieHAwjB8PXbpAxoz2DS6JsPUIlAewD+ht4/uIiIhIHI3wu4VrvnMMblwm5gYzZsDNm9Cvnx2jSlpsWtDBNM1lwDIAQ5VLRUREHG7flYucmVOJFsO34+qc69EGkZHwzTfw0kuQL5/9A0wiElVFrJCQEEJCQqI+HzN5LgBt+04mh6dBVrcQsjuHk8sljDxGBHkNSGPeS8xcXCB3bqtSaoECkD07OGmKl4iIyIP6/3AQnKoxpnPpmBssXgzHj8OPP9o3sCTGME3TPjcyDBN40TTNhfcOPXLjjz/+mE8++STq83xV3+T0lh8g1UW4mzXmjr1v4ux9HXfvf8mZfh9lnHZT/9ouXvznKD5Zc/2XUN1/FSxofVRyJSIiKUxweCie+a9QuP5ZDv8Qy753NWtCeDhs3mzf4J6eQx5xJaoE6uERqFu3bpErVy4CAgKIcHXD/0YAZ64Hcu5aMBeu3+Xy9TD+vR7B9evw7wUXzu/NyJ2D+SHCBVzCcC96mJy5DlDGdS91L2yk1dE9ZAkMAR8fqFULate2XsWKaXNEERFJ9t5fvIlRzaoze8cJWpcv+GiDnTuhQgWYN896hJc0KIF62K1bt/D29iYgIAAvL6843efWnTB+//s0K3ZdZfduOLMrI8EH8kG4KziH4+Z7lGIld/NG+Ep6Lp5DmuBQyJIlekJVqJASKhERSXYyNdvKnQvpub2raMwNXn8dtm2DY8fA2dm+wT09JVAPe5oEKia374azeP8ZVuz6l23bTY4tLUjE5UyQ8Sp56m6nZc49vLtnKTnXb4OICGv+VO3a0LQpNG8OHh5PfW8REZHEYOM5f2rkzUXHCduZ1q3Kow3OnoX8+a0J5H362D/Ap5f8EijDMNIC98cI9wD9gbXAddM0zzzp+oRKoB4WEWkye/spJi68zI7fs3H3SD5IdQeferupW/EEA0L2UXHZeti920qeWraEtm2hXj1wdU2wOEREROyl2scr2TK6MpcvuJDZM/WjDQYOhClT4Nw5a/uWpCNZJlC1sBKmh80wTfOJ2zrbKoF62MZjlxnz+2nW/u7JjS3FAEhdZR+1m53mS/bh++OvcPiwVUzs1VetZKpKFT3mExGRJCEoPASvfFfxbXKG/RNjmDx+6xbkygU9esDIkfYP8NkkvwTqCez2CC8+Tv97m9GLj/L7AoPzS0uBWyi5Wm/lnabX6bN9Gy6/zLF2p86b13pW/PrrULy4XWITERF5GgP/2MBXzWvw687TvFwuhtpOY8ZYI1D+/pAjh73De1ZKoB7miATqQUcv3mLg1MMsn5yTsDM5cC11mIbdTvNVkTCKzF0Cv/5qVWotV87aK6h1a3B3t3ucIiIij5Ox2RZCLmYgcGcMk8fDw60SP88/DzNn2j+4Z+eQBErFkB6jSDYvFg2pRPDJ7Hy99AjZ8t1l8dsNKdqiHnnNVxi39HciFv4GmTNDhw6QJw98+ilcueLo0EVERABYe+4E15ZWonW3gJgbzJ8PZ87AgAH2DSyJ0whUPB3+5xYDfzjCn5NzE3Y+K25l99Os1zm+r5iFjN/9YO0fFB5uzZPq2xdKx1LpVURExA4qf7SCbV9X498LrmT0fOgpiWlC5crWpPHVqx0T4LPTCFRSUCyHF4s/rEjQ6Sx8uegombNF8tubTcjcJDONfZtw/vjf8NlnsHIllCljlUP4/XerPIKIiIgd3Q6/y/YfSlLq9f2PJk8AW7bA9u3Qv7/9g0vilEA9JVcXg3ebFeHc4tKsOXCFIjWusLxvY3KXTUtNp2Kc+HsLzJkDISFWGYTChWHcOAgOdnToIiKSQny49C/Mf7LzcffsMTcYPRqKFoXGje0bWDKgBCoB1C6emcM/leOvozcp9cI5NgxuRKGCaah0LDV/L50Nf/0FlSpZGX6BAlYidfeuo8MWEZFkbvokNzzLHeXFsnkePXnyJCxcCO+8o/1hn4L+xBJQpYI+7J1Snn0ngqnU5jTbP69P6TzelP7jCn+N+xSOHoWGDa1v1oIFYcIEa4RKREQkga06e4wbyyrxWvdbMTcYM8aqb9iunV3jSi6UQNlAqTze/DW+PCdOR1Kz63H2f1OXKnmyUez7Q+z8+kOrKGetWtCrl/Vob/JkCAtzdNgiIpKMDPnhFEaaO4x4LYbFTDduwNSp0LMnpI6hKrk8UaJMoPz8/PD19aVChQqODuWZFMiWhnVflefsGScavXOEo9/XpkKB9NT44xAXfhgLBw5YFc27dYMiRWDaNGsFn4iIyDMIDL/DjimlKNX2AD6ebo82mDjRWtzUs6f9g0smEmUC1atXLw4dOsSOHTscHUqCyJkxFcs+L8fJE05UbHOKjYNeIGexQNruP8GdX6bD339D2bLQuTMUK2YVMtOqPREReUrDlm7FvJCdT7rFUFU8NBS+/dZ6dJc5s/2DSyYSZQKVXOXLkoZt35Vj4/4A8pQIZNZrzUlX5TAfBVwict5c2LPH2hamfXsoXx7WrXN0yCIikgTNmJgKr/JHaVE216Mn58yBCxes+bjy1JRAOUD1Yj6cXlSS6WvOkjrMi0+fr4/PS+v42SPUWhGxdau1JUzt2tCqlbVSQkREJA7+PHuEm8sq06bb7UdPmiZ8/bVVtsDX1/7BJSNKoByoQ+3cXN9RgCEzjxC8w5c3fJ+jQL9FbC2UwUqifv4Zdu60HusNHAgBsZThFxERuWfIlNMYaYMY0abUoyfXr4e9e1U4MwFoK5dEIuhOJG+O2c/cEfkxnSKo/fka5veoT7oQZ6vQ2ciR4OFh7bXXpQu4uDg6ZBERSWQCwoNInyeAMs3PsXtCpUcbtGwJJ07A/v1gOGQHFFvQVi4pmUdqJ2YPLs2ZE25UeNWftX1akbnyKb468hfmsKFw7Bg0aQI9esBzz1lbxYiIiDxg6BJr8vhn3WKY+3TqFCxaZO3TmnySJ4dRApXI5MrszvZJZZi3+SKpQrwYWKEm+fv+wd+eYTB9OuzYAenSQYMG0KwZHD/u6JBFRCSRmDkpNd4VjtL0uRi2bhk/HtKntza7l2emBCqReqlqNq7uykvH/x3mzJR6lC7qzuu/LiCkXCnYsAHmzrWGYEuWtDYvVkVzEZEUbcmZAwQsq0LbbjHsuRoYCD/8YNUdTJPG/sElQ0qgEjE3V4NpA0pw8LBJwUrX+OXVF/Fpso1fTm+HV16BQ4esZaiffgqlS1uTA0VEJEX6cMoZjLRBDH+t5KMnp0+HoCBrBwxJEEqgkoBiuT04vqAEY3/3J+xQIV4vXornhs/jnPNtGDHCqh+VMaO1PUynTnD1qqNDFhERO7oZfps9P5SlbNvDpEv70CKjyEircObLL0POnI4JMBlSApWEvN08L1cPZaZhn+Ps/aglecvcZPCmJZgliluP9SZNsupIFS1qbQvjuBWWIiJiR0OXbMG8mI1PY5o8vmyZNV+2b1/7B5aMqYxBErV+/3Ve6X6Nf/8qQN5+v7NyeEUKps4Bly/DgAFWDamaNeH7762ESkREkq30zTYSeTErATsLPXqyfn2rjuC2bcl19Z3KGEjc1SyZgYsbC9Hpf0fw/64xRcoE89FfyzCzZIaffrLKHPzzD5QqBcOGwZ07jg5ZRERsYOU/h7i5tCptusZQefzgQVi1Cvr1S67Jk8MogUrCnJ1h6ru+bNlzl/TeTnxarQGF35+P/91LUK+etUrv/fdh1ChrkvmWLY4OWUREEtiH009BqhCGtynx6Mlx4yBbNmv+kySoRJlA+fn54evrS4UKFRwdSpJQpVg6Lm0pwBufH+HE180oWO4mw3euwEzlbq3Q27cPMmSA6tWtLWE0GiUikiwERd5h+w8lKPHqIXy8XKOfvHYNfvwRevYENzfHBJiMJcoEqlevXhw6dIgdO3Y4OpQkw8UFZg4uzvpdQXi5uzO0cl2KDZvH+dAr1l56mzdb28F8+61VyfyvvxwdsoiIPKPP124m8nRehnTJ8ujJyZOtxUTdu9s/sBQgUSZQ8vRqlEzP5W35ePXDoxwd2ZK8Fa7w5d4/red9770Hu3eDlxdUqwaDBsHdu44OWUREntIPUwzSFD1D66oPrb4LCwM/P3j9dciUyTHBJXNKoJIhV1eY82FxVm0PxMP05L0KtSk1fDbXIm6Cr681F2r4cBgzBsqWhe3bHR2yiIjE07ZrJ/h3fnVadrn26Pzw+fPh/HmVLrAhJVDJWN3nMnBlZ25avneM/cNeJUfdIyw8v9N63vf++7Brl1XSv0oVGDxY28GIiCQhH8w8DKbByPa+j54cO9Yqrly6tN3jSimUQCVz7m4GC4YX55c1VzBP5OfF0vl4feFcIoiAEiVg61Zrovno0VCuHOzc6eiQRUTkCULMUNZPKUTBlgfJlSlV9JM7dlj/tmv0yaaUQKUQr9XKytl9GSjy/BV+efFVcvVaxPE7563nfUOGWImTmxtUrmw93ouIcHTIIiISi2+2bSbiYFHe7ZLu0ZNjx0K+fNCsmd3jSkmUQKUgWXxcOLygGP38jnPxh8YUrXiLMQdXWidLlbKq1L7/vlV4s04dOHvWsQGLiEiM/KaE4JbnIl3r5Yt+4sIFmDMH+vSxFg+JzSiBSmEMA77pWYiNO0LwiPTknfLVqfH9zwSZwdZo1Oefw7p1cPq09ex87lxHhywiIg84FHiO87Or07jzBZwe/ik+YQKkSgWdOzsktpRECVQKVb2kNxd35KRGx9Ns7NGWbC9tZuP1g9bJGjWs4psNGkDr1tCxIwQGOjReERGxvD9nDwSnYWSnh/Y5vXvX2v+0Y0fw9nZIbCmJEqgUzCONwfoJvoz57RzB6ypSo7Q3/TbMx8SE9Olh9myYPh1++80qvrltm6NDFhFJ0SKIYMWUnORqdJCiuTyin5w1C65etR7fic0pgRL6tsrFob3uZM93l7G1W1Dmf7O4ZQZaz/s6dIA9e8DHxyq+qQnmIiIOM3H/FkK3laVXl4dW3pmmNXm8aVMoXNgxwaUwSqAEgMK5U3FmTUFavnecvwe1JUerbWwLOGKdLFgQNm36b4J57dqaYC4i4gBf/xCAS+Zr9HuhYPQT69fD33+rdIEdKYGSKC4usGBEUfwW/kPw2opUKe/CqL+XWScfnGB+5oy1am/BAofGKyKSkvjfvcTJmVWo1dEfd7eHSo+PG2fte1qvnmOCS4GUQMkjerbIwZ6dLninceX9yjVpMHMmoYRaJ+9PMK9XD1q1gn79IDTUofGKiKQEgxdsh+s+fPHmQ6NP/v7w++/w9ts8uqeL2IoSKIlRqYJpOL81N5VeOcfK9u3I1WMxJ0POWyfTpYNff7V+4/nuO3j+eWtUSkREbMLE5PcpGclc4wgVCj+0ws7Pz9okvl07xwSXQiXKBMrPzw9fX18qVKjg6FBSNI80BlunF+H9iae4MrUpRWtcZtbZjdZJw7BWemzaBJcvW6v0/vjDsQGLiCRTs05u486aqnTp8tCJoCCYMgW6dAEPjxivFdswTNN01L2feONbt27h7e1NQEAAXl5e9ohJYrFq502avXyXu7dd6DDrT6Y2eA2n+/n3jRtW3ZFFi2DgQGulnqurQ+MVEUlOSgxZyGG/Oty64IlHmgce033/PfTqBSdPQt68DovPwRzy3DJRjkBJ4lOvfDrO7MpEgfI3mNHoNYp9NovrkTetk+nTw8KF8NVX8M031g7g5887MFoRkeTjSvh1Dk6rSKW2J6MnT6ZpTaVo0SIlJ08OowRK4iyzjzNHlxSi7UcnOfbR6+R5eTu7bx+zThoGDBgAGzZYJQ7KlIHlyx0ar4hIcjB02Wa4mJ1Pu+SOfmLVKjh8WKULHEQJlMSLszP89FEhJiy8RPDKalSoFsoM//X/NahSBfbuhUqVoHFjGDIEwsMdFq+ISFJmYjJ7SlrSlTtJved8op8cN84qKVOjhmOCS+GUQMlTeat5djb/FUmq2z50rOBLzw2zrS1gwKpa/scfMHIkjBoFjRpZ2wuIiEi8LLmwh8Alz9Ouy93oJ06cgCVLrNEnlS5wCCVQ8tQqF/fk9PbM5C5xiwl1X6La5Onc5d5fcicnGDTIGmL++28oXx5273ZswCIiScwnM05juIXxcZuHNg4ePx4yZIA2bRwTmCiBkmeT2ceZE38WoF6302zt1om8ff7gbNjF/xrUqgW7dkGmTNZeej/+6LBYRUSSkoDIQHZNeY7Srx4lg7fzfydu3YKpU6FbN0id2nEBpnBKoOSZubrCSr/CDJrgz+XvW1Ko0UlWXX9gtClXLti40fpNqUMHq1puWJjjAhYRSQI+XbcB81R+hnbJFv3EjBkQHAw9ezomMAFUB0oS2G/rrvHay85EpLvGyEWHec/3hf9OmqZVs+Ttt63J5r/+ClmyOC5YEZFEzOf1FYTs9iXwcK7/pjlFRkLRolbx4jlzHBpfIqI6UJL0vVTLh7+3p8Y7tRuDKtfglcUziCDCOmkY0KMHrF0Lx49DuXKwbZtjAxYRSYTWXzvA9d9q8XKXgOhzxJcvt/79VOkCh1MCJQmuWH53zmzJSfE6V5jXvB1lvv6RIDP4vwbVq1vzonLntpbfTpniuGBFRBKhYT8fAdNgePsi0U+MG2f98lmlimMCkyhKoMQmvDwN/p5fkJcHneLAgE7k77OYf8Iv/9cge3ZrJKpTJ+jaFbp3h5AQxwUsIpJIBJt32DzZlyItjpIj8wPbYh05AitWqHRBIqEESmzGyQl+HVGQIRPPcOX7VhR68QA7gw7/18Dd3ZoTNXkyTJ8OderAlSsOi1dEJDH4csc6Ig/4MqjLQ4Uzv/3Wmjf66quOCUyiUQIlNvd5tzz8uPgGIesqU7nmXX69uCl6gy5dYP16azPMChWsulEiIinUxCnhuOe+TPt62f87ePOmtfrurbesXz7F4RJlAuXn54evry8VKlRwdCiSQNo1ysSGTSYul3LyauVcfHHw9+gNKleGHTuswnBVq8KiRY4JVETEgfbcPs7FX2rxQucrOD9Q+ompUyE01EqgJEEZhlHZMIya8b5OZQzEnk6fD6N800tcP+NJu/kLmVGnA8aDK1CDgqBdO1i40NoKZuBAPesXkRSj0Q+/sqLrSxz1D6Vw7lTWwYgIKFjQWoAzc6ZjA0ycnumHhGEYl4D1pmm2js91iXIESpKvfDldObUxJ4UqXWdmw7ZUnzGZUEL/a+DhAfPmwQcfWFvBdOqkyeUikiKEEsrqKfnI0/Dof8kTwOLF4O9v1dCTBGUYRjEgC7AlvtcqgRK78/YyOLg4P7U7nmVLx24U/ngW180b/zVwcoLPP4effoLZs6FuXU0uF5Fkz+/AesL/Kk/fLh7RT4wbZ5Ut0LQWW6h+7+OG+F6oBEocwtUVVk8qQLcv/DnzSUfydVzP0dDT0Ru1bQvr1lm7jlesCPv3OyRWERF7GPvDbVwzX6dXs9z/Hdy/H9as0eiT7TwPBAD74nuhEihxGMOAiYPz8vWsi9ya3YRSTc6xJfChJKlyZdi+HdKlsyaX//GHQ2IVEbGlYyFnOPNjDep2OI+b2wMnxo2z6ua99JLDYktuDMN4zTCMnYZh7ATeAMKA7feOVYprP0qgxOHeaZONBX8GEb6jLM/XDuf3fx96FJ07N2zaBPXqQYsW8NVX1r56IiLJxAcLt8N1H4a/WeC/g//+a00a793bGraXBGGa5mzTNMsDL2JNQP/KNM3y915x3l9MCZQkCi1rpmfDegOX83l4sXpGJvqviN4gbVr47Td4/31rZV7PnhAe7phgRUQSUAQRLJ6SlazPH6dskQfmP02aZA3Vd+vmuOCSt9r3Pq5/mouVQEmiUa2MB3s2e5Am3Iu3qpXgkwO/RW/g5ARffGHtnTd5MrRsCbdvOyRWEZGEMuP0BkJWPU/3Lg8UfgoNBT8/aN8efHxiv1ieRS0gCNj5NBcrgZJExbeAO0c3ZyJj5kg+fr4O3TbPwHy4ZNibb8LSpbBhA9SsCRcvOiZYEZEE8L8fruLkFcjAl/P9d/DXX61/2/r2dVxgyV8tYLNpmk/1OEMJlCQ6ObI6c3xdTvKWDmBy/VdovmQikURGb9SgAWzcCJcvWxPNDx50TLAiIs/gfPgljk6rSrW2/nikuVcP0jRhzBjr3zlfX4fGl1wZhpELyMdTPr4DJVCSSKXzNji8PC+lG15mcYsuVP3x++gFNwFKl4a//rJW6FWrBmvXOiRWEZGnNXT5JriQg0+75Pnv4JYtsHMn9OvnsLhSgEfmPxmG8Y5hxH3rCyVQkmilSgU7f81HvU7n2NahJyVHT+c2D815ypnTGomqWBEaNrSKb4qIJAEmJvOmeJOh7GlqlX1gu7IxY6BIEevfNLGVCkAEsAPAMIz8QF4zHvvbKYGSRM3FBf6clI+2H/hz7N1uFBo0nyvmv9EbeXnBkiXwxhvWPnqff64yByKS6M2/uJWgxbXp0OWBKThnzsD8+dbcJyf9iLah68BN0zRDDcPwAkYAH8enAxdbRPUkhmEYAQEBjxwPCQkh5IF9zwIDAwFrU2FJ2b4blIEMXof59v2WFPznd7Z8W57czjmiN/rmG6vg3LBhcOyY9blqp4hIIvXxxCPgWoC+TTP/93Nu9Gjw9LRWGetnX5x4e3t7AYHxGT0CxgJVDcOYjVVIc5hpPrin2JMZ8btfwriX7T2aQYmIiIjEn7dpmnbNOB2VQBkBAQGRDx9/eATq4sWLVKxYkUOHDpEjR46Hmz+1ChUqsGPHjkTbny36TMj+bt26Ra5cuTh37hxeXl5PviCO4hrjtD+u0q+jF+71NrJ6Rm5Kpir0aKP16wls0QLPMmWs5cCZMtk1Rkf1Z4uvTUr7/rZFn47+O5Nc+kvoPh35dRm48XcmvdCCX5bdoEnV9NbBSZNg0CD4+2/IlSte/dkiRkf2F5+vjbe3tzfxH4F6Zg55hBffN+np6Zmg39zOzs6Juj9b9GmLGL28vBwSY9+2XmTyuc4bLzam9ht/sXHhGSqlKRm9UbNmdMiblwUXL0LjxrBiBeTLF3OHNojRUf3dl5Bfm5T6/Z2c/s4kl/5s1ae9vy4mJnNmpSdt4UBaN8yDYQCRkVYC9fLLULx4vPqzRYyO7u++uHxt7D3ydF+KnKHWq1evRN2fLfq0RYwJLT4xvt4oAwuWhRCxpSLVGwWy5tajhWTrDRgAmzdbE8qrVbN+q7NjjI7ozxZS6vd3SvzaJPb+bNVnQntSjMuubSdgXn3adA0iatH8smVw/HiMpQuSwp9jUvi6JDSHPMK754k3Pn/+fNQQXs6cOe0Rk8TBrVu38Pb2JiAgwCa/ccTHqr9u06iRiVn4GL8vD+GFDFUfbXTlijUKdfIkLFoENWrYP1A7SUxfG/mPvi6Jk6O+LuW/nsWuwS9x+R9XMme8N45Rv741afyvvyDupYiSrXh+bRzyB5aoR6Dc3d2jfZTEwd3dnY8++ihRfF3qVU7LhrUuOJ8uQPPaaZl7Zd2jjTJntopsli9vVfZduNDeYdpNYvrayH/0dUmcHPF1uWZeZ/ekcpRudfK/5OnAAVi1yhp9UvIEJI2/M4l6BEq/tUlc7TkUQpW6QYR4/8sPq/3pnCOGAnQhIdbGnPPmwcSJ0KWL/QMVkRTt7Q2/8m3NV5i/9gYv1ro3ebxrV2t/T39/lV55OhqBEnlaz/m6s2tDWlIH+/BmjYKM81/0aCN3d5g1C3r0sP7BGj5cBTdFxG5MTGZMTE3aQhdoWfNe8vTvvzBzJvTureQpiVECJclG8UJu7N+YDk8jLX2ff47hx359tJGzM3z7LXz2GQwdalX7jXykooaISIJbcu0vbs2rT9tuD0wenzTJemzXrZtDY5P4UwIlyUqBPC4c2pCJ9J4uDK3xPB8dmvNoI8OwkqeJE8HPD15/HUJDH20nIpKAPppxCjD4tGMB60BoqPVvUPv24OPj0Ngk/pRASbKTM7sTR9ZnxSdzJJ/WqcmHh2bH3LBbN2s+1MKF0KIFBAfbNU4RSTmumtfYM6l89Mnjv/4KFy9aI+GS5CiBknjZsGEDzZo1I3v27BiGwcJEuqItcyaDQ6uz4ZM5ks/q1Io9iXrxRWvy5saN0KhRktx7asSIEVSoUAFPT08yZ85My5YtOXr0qKPDEmDChAmUKlUqqhhglSpVWLZsmaPDkoeMGDECwzDoF0MNpoTy0cZVmEeL8HH3bNYB04QxY6yVwb6+NrtvUvLxxx9jGEa0V9asWR0dVqyUQEm8BAUFUbp0acaPH+/oUJ7ofhKV8X4SdfiXmBvWqWMtId6/H+rWhWvX7BvoM1q/fj29evXir7/+YuXKlYSHh9OgQQOCgoIcHVqKlzNnTkaOHMnOnTvZuXMnderUoUWLFhw8eNDRock9O3bsYNKkSZQqVcpm9zAxmTkxDZ6FLtKiZjrr4JYtsHNnjIUzU7LixYtz8eLFqNf+/fsdHVLsTNN01OuJAgICTMAMCAiIS3OxM8BcsGCBo8N4ostXIs2MJf8xyXLRHHZoVuwN9+wxzUyZTLN4cdO8cMFu8SW0K1eumIC5fv16R4ciMUifPr05ZcoUR4chpmkGBgaahQoVMleuXGnWrFnT7Nu3r03us/DqRhO3u+ZbX5747+BLL5lmkSKmGRFhk3smRR999JFZunTpp7nUIXmMRqAk2cucyeDg6mz4ZDL5rHZthh2eFXPDMmWsR3kBAfD881ZNliQoICAAgAwZMjg4EnlQREQEs2fPJigoiCpVqjg6HMHafqRp06bUq1fPpvf5ZIY/BgafdsxvHTh5EhYsgHfeASf9GH7Q8ePHyZ49O/ny5eO1117j1KlTjg4pVonyK+fn54evry8VKlRwdCiSTGTOZHBoTVZ8Mpl8XrtO7ElUkSJWEgVQvTocOWK/IBOAaZr079+f6tWrU6JECUeHI8D+/ftJmzYt7u7uvPXWWyxYsABfzXlxuNmzZ7N7925GjBhh0/v8a15lz6QKlG51ikwZ79UuGDMGMmSwVt9JlEqVKvHjjz+yYsUKJk+ezKVLl6hatSrXEuu0CkcNfcVlTE6P8BI3ksgjvAddvhJp+pS4YJLlojn00M+xN7xwwXqUlymTae7ebb8An1HPnj3NPHnymOfOnXN0KHJPSEiIefz4cXPHjh3m+++/b2bMmNE8ePCgo8NK0c6ePWtmzpzZ3Lt3b9QxWz3Ce2vdLyaY5sK1N60DV6+aZpo0pvnxxwl+r+Tm9u3bZpYsWczRo0c/qake4YnY2v2RqIz3RqKGHvkp5obZssH69ZAnD9SubU34TOT69OnDokWLWLt2rTbfTkTc3NwoWLAg5cuXZ8SIEZQuXZqxY8c6OqwUbdeuXVy5coVy5crh4uKCi4sL69evZ9y4cbi4uBAREZEg9zEx+XmSNXm8eU1v6+D331vFe3v1SpB7JGceHh6ULFmS48ePOzqUGCmBkhQncyaDA6utJGp4rXqxJ1E+PrB6NZQqZe2UvmqVfQONI9M06d27N/Pnz2fNmjXky5fP0SHJY5imSUhIiKPDSNHq1q3L/v372bt3b9SrfPnytG3blr179+Ls7Jwg91l4bSOB8xryRrc7VuXxu3etnRA6dYKMGRPkHslZSEgIhw8fJlu2bI4OJUYujg5Akpbbt29z4sSJqM9Pnz7N3r17yZAhA7lz53ZgZPGTJbOVRJWoe4nhterhvH4WnxR5/dGGXl6wfDm89BI0bWoVvmve3P4BP0avXr2YNWsWv//+O56enly6dAkAb29vUqdO7eDoUrYPPviAxo0bkytXLgIDA5k9ezbr1q1j+fLljg4tRfP09HxkjqCHhwc+Pj4JOnfQmjxemU863vul5qef4MoVa/K4POLdd9+lWbNm5M6dmytXrvD5559z69YtOnTo4OjQYuaoZ4dxef6pOVCJz9q1a03gkVeHDh0cHdpTuXQ50szge9Ekxzlz5Mm5sTcMCTHNVq1M08XFNOfNs1+AcRDT1wMwp02b5ujQUrzOnTubefLkMd3c3MxMmTKZdevWNf/8809HhyUxSOg5UJciL5sUOWI+99oR60BEhGkWLWqaL76YYPdIblq3bm1my5bNdHV1NbNnz262atUqrvMFHZLHGKbpsN3on3jjW7du4e3tTUBAAF5eXvaISVKgCxdNfGteISAsmG83/E3vXC1ibhgeDu3aWaNQP/8MrVvbN1ARSTK6r5/FpFqv8/vaWzSv5QWLF0OzZrB5M1St6ujwkhvjyU1scFMlUCJw5lwkxZ+/TpDbdaZsOM6bWZvG3DA83Jq/MGsWzJgBb7xh30BFJNGLJBLvtn9g7KhIwNFs1vynWrWszYOTwIKUJMghCZQmkYsAeXI5sWdNelIHZaRLvbz8cnVlzA1dXGD6dOjQwarhMm2aXeMUkcRv4bWN3J7XiHbd7lrJ044d1qred991dGiSgJRAidxTKL8z21Z74nYlJ20bZGThzXUxN3R2hilToGtX6NwZJk2ya5wikrjdrzz+cce81oHRo6FAAWgRy/QASZKUQIk8oGRRVzavSo2Lf0FeapyaFYGxDLc7OVn1XHr1gu7dwc/PvoGKSKJ00bzE35Mq8Vyr01blcX9/a95k//7WL1+SbCiBEnlI+VJurPnTFeNQcZo2D2dD8M6YGxqGVdOlXz/o3dvankFEUrQPN6yCo0X5uHt268CYMZA+PXTs6MiwxAaUQInEoHr5VCxb6oS5vQJ1W91kW8jemBsaBnz9Nbz3nlXb5csv7RqniCQeEUTwyyRPvApd4oWannDjhvW4v2dPSJPG0eFJAlMCJRKL+tXSsOCPCCLWPU+N1hfYG3Yw5oaGASNHwtChViI1fLh9AxWRRGHOv2sImteI9t1CrMnjEydaK3e1bUuypARK5DGa10nLL/NDCVtanyrtjnM44ljMDQ0DPvsMPvnESqQ++cS+gYqIw3085TyGk8nHnXNDSAiMG2et1s2SxdGhiQ0kyq1c/Pz88PPzS7ANHUWeResmngTPDqDzqy9QweM3Dk1JTW4jV8yNP/wQXF3hgw+spOrDD+0brIg4xNGIExz/vg7V25zDJ0MhmDYLLl6EAQMcHZrYiAppisTRuJ+u07ddBjIMmMrhL5uR2cgUe+MvvoAhQ6xRqaFD7RekiDhEi99/YFHLN9m08y7VyrpDiRJQsCD8/rujQ0sJHFJIM1GOQIkkRm+/kYF/r1/h876dec5nPIcHt8eLWBL7Dz6AyEgYNswqefDBB/YNVkTsJogglo7PT7ZKZ6lWLjcsWwaHDlmlTiTZUgIlEg+fvZ2Zf69fYuIHvSmfYTT7uvckNaljbjx0qJVEDRliJVHvv2/fYEXELr4+uoTwVa/y7o//Wge++goqVoTq1R0bmNiUEiiReJrwUVauXr/Abz3eoXq6r/ir9Tu44hpz4w8/tJKowYOtJOq99+wbrIjYlInJ2O/CcMt0k56vZILdu2HNGpg715oHKcmWEiiReDIMmDsmO/Vu/MPadv1olO4bVjZ8F6fYFrV+9JGVRA0aZCVR2g9LJNlYeXsr16a/QOte10iVKp21bUu+fPDii44OTWxMZQxEnoKTE6yYmoOyDf9lTatetN4yFjO2dRGGYZU1GDIEBg60Cm+KSLLwwc8H4XZaRr6VF06dgjlzrKK6LhqfSO6UQIk8JVdX2DQ3B4XKBzCvaUd6/v2YCaP360QNHmwta9a2LyJJ3gXzIrv8KlOimT95czvB//4HGTLAm286OjSxAyVQIs8gdWrYsSg72fOF8H3Dlgw7OT32xoZhVSm/v+3LuHF2i1NEEt6Qjcthf0k+7ZUF/vkHpk2zNg3Wti0pghIokWfk7Q27l2chvZcLn9evwdgLv8be+P62L+++C337wvjx9gtURBJMGGHM9vPBu/AlWtRNa819SpPG2vdOUgQlUCIJIEtmg91/ZsQjLB39Ghblp+tLY29sGNZQf//+0KePtV+WiCQpUy8u4+78xnTvGYnT9avW3+M+fUBFn1MMJVAiCSRvHoO//vTG7WIe2r+Qnj+DN8Xe2DCsWjG9e8Nbb8GMGfYLVESe2ReTruLkFs7gDtmtOY2GYY0qS4qhBEokAZUo5syaZalw+rsMTVvfZk/4/tgbGwaMHQtdu0LnzjB7tv0CFZGntivsb85ObESddhdIZwRYj+Lfegt8fBwdmtiREiiRBFatghtzf4skfHldqnc/gL95JvbGTk7Wdg9t28Ibb8CCBfYLVESeyvsLt8LF7IzqlRv8/ODuXW0anAIpgRKxgVYNPRg3NYjgqW0oN2wpV7kae2MnJ5g6FV56CVq3hiVL7BeoiMTLTW6yZrwveZ4/Q9n8ofDNN9YIcrZsjg5N7CxRJlB+fn74+vpSoUIFR4ci8tT6tEvHe//7l+vDe1DebxpBBMXe2MUFfvoJmjSxEqmVK+0XqIjE2Rf7/yByw/MM7uUNkyfDjRvaoimFMkwzlurJtvfEG9+6dQtvb28CAgLw0soGSYJME94YcIlZYzLz3NyRbHt5YOz75gGEhFhbQKxbZ+3oXrOm3WIVkceLJJL0PWYTurAxAcfS4FY0P9SvD9OnOzq0lM4hmw4myhEokeTCMGDmV1mp9dpl9rQdQLN1o2Pf8gXA3R1++w2qVoWmTWHLFvsFKyKPtTBgLbdmNqdNt0DcfpkBFy9auwtIiqQESsTGnJxgxfRslKhxnRUtevDm32Mff0Hq1PD771C2LDRuDDt32idQEXmsD388AXdT8WnnbDBqFLz8MhQp4uiwxEGUQInYgZsbbP4tGzkK3GVao1f50H/q4y/w8LAmkxcrBg0awL599glURGJ0yjzNQb+alHvxLDk3zrE2Dv7gA0eHJQ6kBErETry8YNeyLKRL7c5nDasx8epvj7/A0xOWL4d8+aBePTh82D6Bisgj3l+9Eo4W5fMe2WDECOsRe5kyjg5LHEgJlIgdZckC21dkINWNbPR4ISeLglY//oJ06eDPPyFrViuJOnXKLnGKyH/ucIeFfjnwKX6JhjeWwaFDMGSIo8MSB1MCJWJnhQoarF2aBqcDpWj1ajg7wvc8/gIfH6usQZo0ULcunD9vn0BFBIBvzswnbFEj3u7phPHFcKhdG6pUcXRY4mBKoEQcoHJ5F+bPN4j8sy41exzitOn/+AuyZoXVqyEy0hqJunLFLnGKpHSRRPLVuBBcve7QP9s+2L1bo08CKIEScZjmDVLxzeQg7kxpS8Xhv3Od64+/IHduK4kKCLBqz1x/QnsReWZzA1ZwY/LLtO5+k7RffwqVKkGdOo4OSxIBJVAiDtS3ozdvf3qVq8P6UnnGBO5y9/EXFCxoPc775x+rxEFgoH0CFUmhhvxwCuNOGkZWOA+bNlmjT4ZD6jZKIqNK5CIOZprQsttlFk3PwPNLR7Ku/hCcnvS7za5d1m/Bzz0HS5da86NEJEFtC99F5QIZqVET1l/uBpcuwd69SqASH1UiF0mJDAPmfZeFsvWvs/GlvrTfN/rJF5UrZyVOO3ZYe+eFhNg+UJEUpv+8LXA2D1/Xv2Gthv3gAyVPEkUjUCKJxO3b4FvzCucuhTFs6zI+zd3lyRetWmXVo3nhBZgzx9qUWESe2RnzLHkrXqawV1aOGp2s0ad9+8DZ2dGhyaM0AiWSkqVNC9uXZMbbLTWfNanMtJsLnnxRvXrw66/W1i+dO1ur9ETkmb23aRHsrMAXDS9Zizc+/1zJk0STKEeg/Pz88PPzIyIigmPHjmkESlKUQ0ciea5qMGGld/LncmfquT//5Itmz4bXX4fu3eG77/SYQeQZBBBAxpYb8T5WgSverXCKCINt2/T3KvHSCNR9vXr14tChQ+zYscPRoYjYnW9RJ5YvcsfYWpUmnS5zIPLQky967TWYMgW+/x4GDbJmpovIUxlxfB7hi5rwfqOTOP21xRp9UvIkD0mUCZRISle7uitTZ4YRNrsV1T5YxwUuPPmizp1h7Fj48kv44gvbBymSDIURht8YF1JlvE3vVf2gRg2r7prIQzTjVCSR6vCKB6fO3+TT/j2pmPsLDvXsjRdPeJT99ttWoc2hQ63NiN9+2z7BiiQT064t4va0V+jx0h5S/bQDNm7U6JPESAmUSCL2yTvpOHXmKj/1GUTNnJ+yvflQXHF9/EVDh8KtW9C3L3h5QceOdolVJKkzMfn4+0s4mc58svVdq1ht9eqODksSKSVQIonc9NEZ8T93hU1t3uPl9SNZWH4oxuPmTBoG/O9/VhL15pvW8r6XX7ZfwCJJ1MqQDVwc34rGdXaTaelfMHeXo0OSRExzoEQSOWdnWDEzM/lLBrPoha4M9Pd78kWGYa3Ga93aWp23fLntAxVJ4t77ZQ9cysbovUOsArVlyzo6JEnEEmUZg/tUSFPkP1euQNEq17nhfomJWw7QLd2rT74oLMz6QbBqFaxYAc/HoSSCSAp02DyCb+lQSqWJYN/2cnDgAPj6OjosiRuVMRCR2GXODFuWpsftUm7eapWJlaHrn3yRqyvMnQtVqlgVy3futH2gIknQuyuXw/5S/O/05/DGG0qe5Ik0AiWSxKzeEEb9+iYur/3G3ull8DWKPfmi27etpdjHj8OGDfrhIPKAK1whW6O9ZDmVn39OFcM4dhTy53d0WBJ3GoESkSerW8OVydPDCfuxDVU/WcElLj35orRprc2Hc+Swtn85dcr2gYokER8fmEfkigYMC/gSo2sXJU8SJzZNoAzDGGIYxhbDMIINw7hpy3uJpCRvtknDe1/cJOCTflSZMYEggp58Ufr01o7ynp5Qty7884/tAxVJ5O5wh6lfp8Mj02W6BMyyyoCIxIGtR6DcgF+BCTa+j0iKM/L9dLzY5Rr+XYZSb/UXhBP+5IuyZIGVK61Nh+vVg3//tX2gIonY+EvzCPn5JXq6fItrr27WKK1IHNg0gTJN8yPTNL8B9tvyPiIpkWHAnO98KFsngL9eGki7gyMxnzy1EHLntlbl3bgBDRtalctFUqBIIhnlF4izczgf3JoO77/v6JAkCdEcKJEkzNUV1v6akZx5IpjdpB0fXfw+bhcWKmSNRPn7W6vzguLwCFAkmZkfvJxrE17lRZ9ppHunE2TK5OiQJAlJVAlUSEgIt27divYSkcfz8oKtS3xIG56Oz16oyMzb8+N2YcmSsGwZ7NsHL74IISG2DVQkETExGTDtANxIz6hbk2HAAEeHJElMvBMowzA+NgzDfMKr/NMEM2LECLy9vaNeuXLleppuRFKcnDlh/RJPXI750rGNOxsjtsTtwkqV4I8/rA1TX3sNwuMwj0okGfgj5E/OjmxDnSyzyT+4DaRL5+iQJImJdx0owzAyAhmf0MzfNM27D1zTERhjmma6B9o8cuOQkBBCHvgt+NatW+TKlUt1oETiaNGyMFo0M0jVYxp/j6tNIaNg3C5csgRatrSSqBkzwClRDU6LJCgTkwKTRnH6rfc46FMTX//l4OHh6LDk6TmkDlS8NxM2TfMqcNUGseDu7o67u7stuhZJEZo3dmW0XxAD3upK1QLDOdyvOxmf+PsO1jyon3+GNm2sMgd+ftYsdZFkaHnYak5/8Ro1sszFd8SbSp7kqdi6DlRuwzDKALkBZ8Mwytx7pbXlfUVSsv7dPeg68CZX+w/m+YWjucvdJ18E8OqrMHkyTJhgrUZy3C4FIjZjYvL2j9vhTF6+zbAA2rd3dEiSRNl0KxfDMKYDHWI4Vds0zbVPul5buYg8nchIqNf6GmuXpKbeuuGsqPgZTnH9fWnMGHjnHRg+HD74wKZxitjbqvB11C+SiyrBe9kyP7u1T6QkdclvKxfTNDuapmnE8Fpny/uKpHROTrDkRx+KlLnDqmZv08v/y7hf3K8ffPopDBkC335rsxhFHKHPT5vgVAH8ntul5EmeiTYTFknG/v0XilW+yTX3C3yzZTv90nWM24WmCQMHwujRMG0adIzjdSKJ2NqIDdTxzUKF28fYvu05a/mqJAfJbwRKRBwrUybYtNQb90t5eKdVHhaFLo/bhYYBX34J3brBm2/C3Lm2DVTEDvrMXA3HijC++SUlT/LMNAIlkgKs3RBBvfqROLWZy/ZpxXnOKBO3CyMioEMHmDMH5s+HZs1sGqeIrWyM3EyNEul57vYFdh+tBqlTOzokSTgagRIR26hdw5kp0yIIn9GWmp+v5Dzn43ahszNMnw7Nm8PLL1vbv4gkQX2mLYLDvnzbHSVPkiA0AiWSggz6PJD/DfMk18wPOPDG+3gRx79XoaHWdi9r18Ly5VCjhm0DFUlAW0M3U7W8JyWDAvn7RFXVOEt+NAIlIrY1cognLTvd4Fznj6m7/lPCCIvbhW5uMG+etWqpaVPYvt2mcYokpN7TZsP+UowbklHJkyQYjUCJpDBhYVC58XV27zJ4acvX/FrsU4y4/gJ3+zY0bAiHDsG6dVC6tE1jFXlWO26somLdDBQLduLQkTKODkdsQyNQImJ7rq6w5rcM5MwJvzXpzODLY+J+cdq0sHQpFCgA9evD4cM2i1MkIfT+6RfYU5ZxI3I7OhRJZjQCJZJCnT0LJSoHEpjzMJPXnqSLR5u4X3ztGtSqBdevw4YNVkIlksjsPjaXcm/kpXCwF0f2F9XTu+RLI1D3+fn54evrS4UKFRwdikiylTs3rFuSFpdDpej2elpWRqyJ+8U+PtaKvLRpoW5dKxsTSUxMkz6LfoMdFRn7v7xKniTBaQRKJIVbtDScls0NXHtOYefYapQ0SsT94nPnrBV5rq7WSFTWrLYLVCQe9q36hjIfViL/7eyc2KcEKpnTCJSI2F/zJi587RdK6LfdqTHmNy5wIe4X58oFq1dDcDDUq2ftHSPiaHfv8va21bC1Kt+MzKnkSWxCCZSI0K97anoOCuTmgGFU++1rAgmM+8X588OqVXD1qpVEXbtmu0BF4uDgj4PY8OdA8pS6QLPGLo4OR5IpJVAiAsC3X3jSqPUt/N/4jPpbPySc8LhfXLSoNRJ18aKVRF2/brtARR7nyBF6XfobNtRk9GeZNPokNqM5UCIS5e5dqFD/BgeOhPPqX98wu8DwuNeIAti/H+rUsWaor1oF6dPbLliRh0VEsPmtMlQ/8h25A4vhvyejEqiUQXOgRMSxUqWCdQvTkzWDG3Mbd2TYtTHx66BkSStx8veHBg3g5k0bRCkSM3PcWDoXKAibnmfCqPRKnsSmlECJSDQ+PrB5mTceN3MwvEUFpt39JX4dlC5tJVEnT0KjRnDrlm0CFXnQiRPM3j6MY1NHUa7BNZo0dHZ0RJLMKYESkUfkzw+r/kiD866KvNnBhZWRq+PXwXPPWXWijhyxkqjAeExKF4mvyEhC3upEn2Jd4ERBpn7p4+iIJAVQAiUiMapcyeCXWU6Yv75E04GH2Me++HVQrpyVRB08CE2aWPvoidjChAmMLvw318YNo2XHAEqVcnRAkhIogRKRWL3yogujvw0j7Os+PP/NfPzxj18HFSrAn3/Cvn3QtCkEBdkkTknB/P25PvI9PkszFJdgL8Z/poULYh9KoETksfr3cqfP+0EE9v+EqnPGco141nmqVAmWL4fdu6FZM6vopkhCME3o2pVBvfJwd3wfeg8IJUcORwclKYUSKBF5orFfeND8jVtcbD+SGus+5A534tdB1aqwbBls3w7Nm8OdeF4vEpMffuDUyVX8cGAIHt7hfPpeWkdHJCmIEigReSLDgF9/8KJijbscajmcpgcGxa/QJkD16rB0KWzdCi+8oMd58mzOn4cBA3jr4xcwf27LF5+44enp6KAkJVEhTRGJs1u34LmaAZz6N5C2f41nZs4R8Su0CbBxozWpvEwZWLIE9Hdb4ss04YUX2OK0nWpBc8h2sSxn96fDRbu2pFQqpHmfn58fvr6+VKhQwdGhiMgDvLxg4xJvfFy8+blxW4be/Cr+nTz/vLU6b/9+qF8fbtxI+EAleZs5E3PpUt5s+xKsrcOELz2VPIndaQRKROLt8GEoV+0Od0ptY+KKM3Rz7xD/TnbvthKo3LmthCpjxoQPVJKfixfB15e5w56j9ZTxlM6SlT1rMqjqeMqmESgRSRqKFYM//0iF87aqvNUhFYsjl8a/k7JlYd06uHABatWCy5cTOkxJbkwTevYkNK0bvdNUhsO+/PCVkidxDCVQIvJUqlczmD3LGXPuK7w48ATb2R7/TkqWhPXrrcd4NWvCP/8kfKCSfMyZAwsX8s2cl/j347dp2vYm5co5OihJqRLlI7zQ0FDCw8O5ffs2NWrUYMOGDaRNq+WpKY2Liwtubm6ODkOe4OvxoQzo44bH10PZ/U57ClM4/p2cPAl16oCLC6xZA3nyJHygkrSdPg3lynGjRQ2y56pG2P/6cfKoq75VBBz0CC/RJVChoaEcPHiQyMhIe8cjiYyTkxPFixdXEpUE9H3/DuNGpcbn57fZ+/p75CRn/Ds5c8ZKosLDYfVqKFgw4QOVpOnuXahWDW7e5K0/2zCx5Af06hPJ+FH6xVoAJVCW4OBgDh8+TN68eUmdOrW9Y5JE4s6dO/j7+1OsWDHSpEnj6HDkCUwTWncO4tef3Mj1ex92N/mcjDzFpPB//oG6da16CWvWQNGiCR+sJD1du8JPP3F61zwKjr5Cqt9b88+JNKRL5+jAJJFwSAKVaBd+pk6dWj84hYCAAO7evYu7uzseHh6ODkdiYRgwa7IH167fZs3LX/P8yrfZXu0bPIlnZcMcOaw5UfXqWXOiVq5EO8OmcFOnwpQpMHUqb4WvI3LaKD4ZG6bkSRwu0Y5AaeQhZbv/fbB7926CgoLw8vKidevWSqISubt3oXqjQHbti6Di+oFsKDUed9zj39HVq9CwIZw6BYsWWbWjJOXZs8faBqhdO5ZOfJGmDZ3JcroyZw96oSf78gCVMRB5mIeHB25ubty6dYuQkBBHhyNPkCoVrFnkSeH8zmxv+CktTr0T/y1fwKoJtXatVeqgfn1YuDDBY5VE7sYNeOklKF6cwHHDaTdrGaxswORvPJU8SaKgBEoSNXd3d9zdn2IEQxzGyws2LvMku6cnK+oPoP2l9zCfXDc35o6WLoUWLawfpBMnJnywkjhFRkL79nDzJsybR//Ab7je90NeaB1EsxdU9EkSByVQIpLgMmeGLSvTkj4kK7807MDbNz98uiTK3R1++QV69YK33oKPP7ZmrEvyNmKEtU/izz+zNe9FpvQvRprItEwZq0f4kngogRIRm8iTBzb96UGa84UZ36wBnwWPfrqOnJxg7Fj44gv45BMrkYqISNhgJfFYuRKGDYNhwwhpXIfWy6fBT+349ms3smRxdHAi/1ECJSI24+sLa5amxnVPRT56tSgTwiY/XUeGAYMHWyuyfvgBXn4Z7txJ2GDF8c6dg9dfhwYN4MMP+eT2aM51H0LleoF06qAfV5K46DvyMdatW4dhGNy8edPRoTwTW74PwzBYqAm+8hiVKsHiBW44/dmInp1SMyfy16fvrFMna0L5ihXWKr0bNxIsTnGwkBArMU6TBn7+mYPORxg1NC0u/2bj54me2u9OEp1km0B17NgRwzAeeTVq1MjRodlUrVq16Nevn6PDEImmQX2DWT87wazXeb3fZZaZy5++sxdesCqVHzwINWrA+fMJF6g4Tv/+sHcvzJtHhE86Wm/7mshxvfn8Myfy53d0cCKPSrYJFECjRo24ePFitNcvv/xi1xhCQ0Ptej+RxKr1K074TTCJ/LY3zT7cwypWPX1nVarApk0QEGDVCTp8OOECFfv76Sf47jsYNw4qVGBs6AQOdulHkbLBDOibaOs9SwqXrBMod3d3smbNGu2VPn16APz9/TEMg71790a1v3nzJoZhsG7dulj73LJlCzVq1CB16tTkypWLt99+m6CgoKjzefPm5fPPP6djx454e3vTtWvXGPupVasWffr0oV+/fqRPn54sWbIwadIkgoKC6NSpE56enhQoUIBly5ZFu+7QoUM0adKEtGnTkiVLFtq1a8fVq1cBa9Rt/fr1jB07NmrEzd/fP+raXbt2Ub58edKkSUPVqlU5evRotL4nTJhAgQIFcHNzo0iRIsycOTPa+ePHj1OjRg1SpUqFr68vK1eujPXPSSQmPbs788X/won4fDCNP9vCOtY9fWfFisHWreDtbSVRf/6ZYHGKHW3eDN26QYcO0K0bZzjD+/+7jnHYl9lT0uKi/EkSqaSRQAUHw+7dj38FB9s8jP3799OwYUNatWrF33//zZw5c9i0aRO9e/eO1u7LL7+kRIkS7Nq1i2HDhsXa34wZM8iYMSPbt2+nT58+9OjRg1deeYWqVauye/duGjZsSLt27Qi+994uXrxIzZo1KVOmDDt37mT58uVcvnyZV199FYCxY8dSpUoVunbtGjXilitXrqj7DRkyhNGjR7Nz505cXFzo3Llz1LkFCxbQt29fBgwYwIEDB+jevTudOnVi7dq1AERGRtKqVSucnZ3566+/+P777xk0aFCC/dlKyjF4oAsffx5G+Icf0mDUajay8ek7y5HDGomqUgUaN4avv1aZg6Rk1y5o0sSaKDdhAqYBbxz5nLDPBtFvYDhlyjg6QJHHME3TUa8YBQUFmTt37jSDgoL+O7hrl2la/yzG/tq1K1o/HTp0MJ2dnU0PD49or08//dQ0TdM8ffq0CZh79uyJuubGjRsmYK5du9Y0TdNcu3atCZg3btwwTdM027VrZ3br1i3afTZu3Gg6OTmZd+7cMU3TNPPkyWO2bNkytrcXpWbNmmb16tWjPg8PDzc9PDzMdu3aRR27ePGiCZhbt241TdM0hw0bZjZo0CBaP+fOnTMB8+jRo1H99u3bN1qb++9j1apVUceWLFliAlFxV61a1ezatWu061555RWzSZMmpmma5ooVK0xnZ2fz3LlzUeeXLVtmAuaCBQue+H7j6/73wbx588wpU6aY33zzjXnt2rUEv484zuAPQ00wTbevB5mbzc3P1ll4uGm+9571b0GHDqZ57/taErEDB0zTx8c0K1UyzVu3TNM0zRkRM02qbzCzFwo0g4MdHJ8kJQ7JY5LG4GjRotZvKk9q85DatWszYcKEaMcyZMjw1GHs2rWLEydO8PPPP0cdM02TyMhITp8+TbFixQAoX758nPor9cAmqc7Ozvj4+FCyZMmoY1nuFT25cuVK1P3Xrl1L2rRpH+nr5MmTFC5cOM73y5YtW1TfuXPn5vDhw3Tr1i1a+2rVqjF27FgADh8+TO7cucmZM2fU+SpVqsTpfYrEZPjHroSGhjK6/0jquA1gfS9nKlHp6TpzdoZRo6BkSejSBY4cgQUL4N73uSQyJ05YG0bnzAnLloGnJ//yLz0n7YVNbzBrHaRO7eggRR4vaSRQadJYe2LFk4eHBwULFozxnJOT9fTSfGC4Pyws7LH9RUZG0r17d95+++1HzuXOnTvafePC1dU12ueGYUQ7ZtxbtxsZGRn1sVmzZowaNeqRvrLF4QfF4/p+8Nh9pmlGHXvwzym29iLxYRjw5RduhIaF8m3v0dR27cuGbs6UJ26/gMTojTegSBFo2RLKl7dKHlSokFAhS0I4exbq1rXmrv35J9ybl9rtn08Iem8Eb3S9Q82ayp4k8Usac6BsIFOmTIA1r+i+ByeUx6Rs2bIcPHiQggULPvJys8Pulvfvnzdv3kfufz9pc3NzI+IpqjQXK1aMTZs2RTu2ZcuWqFE1X19fzp49y4ULF6LOb9269RnejYiVRI390o3ufUK4030sNadNZw97nq3TChVgxw5rdOP55+GBEWNxsEuXrJEnJydYtcra8wdYbC5hYc/6eHs48+3/lDxJ0pCsE6iQkBAuXboU7XV/xVrq1KmpXLkyI0eO5NChQ2zYsIGhQ4c+tr9BgwaxdetWevXqxd69ezl+/DiLFi2iT58+9ng79OrVi+vXr9OmTRu2b9/OqVOn+PPPP+ncuXNU0pQ3b162bduGv78/V69ejTbC9DgDBw5k+vTpfP/99xw/fpyvv/6a+fPn8+677wJQr149ihQpQvv27dm3bx8bN25kyJAhNnuvknIYBkwY607nt0IIfnMcNX6axH72P1un2bPD+vXQurU1KjVokLZ/cbRr16B+fQgKsup43ZsOcIMbdJy3GBa14Ae/1KRL59gwReIqWSdQy5cvJ1u2bNFe1atXjzo/depUwsLCKF++PH379uXzzz9/bH+lSpVi/fr1HD9+nOeff57nnnuOYcOGxenxWULInj07mzdvJiIigoYNG1KiRAn69u2Lt7d31CPJd999F2dnZ3x9fcmUKRNnz56NU98tW7Zk7NixfPnllxQvXpyJEycybdo0atWqBViPPBcsWEBISAgVK1akS5cuDB8+3FZvVVIYw4DJfu607RTG7Q7jqTZ3LAc5+GydpkoF06fD6NHw1VfQvLlVN0rs79YtaNTIGoFatYr7lTEjiaT1lT5c7/MxDV8M5qVWmhYgSYcR09wWO4nxxsHBwRw+fJhixYqRJk0ae8ckicT97wN/f39u3rxJYGAg7du3f6ZFAJL4RUTA651CmDvLGe9fu/LXi4MoyqMLROJtxQprNCprVpgzB0qXfvY+JW6Cg63kaf9+WLuWB2sTfB4+kmENKpLuUDUO7XHXnH95Wg7JvJP1CJSIJC3OzjBrmjsvvhLOrdaTqLRoOHvZ++wdN2wI27dbo1IVK8I330AcH2/LMwgJgRdftGr1LVsWLXlazWqGfWhirK/FgtlKniTpUQIlIomKszPMnZmKZi0jCXxpGlV/+ZbNbH72jgsXhm3boFcva9+1xo3hgUUkksDCwuC116y5aIsWQeXKUafOc55Wf0yDEYMZ/gXcmykgkqQogUpmOnbsSMuWLW1+n7x58zJmzBib30dSJhcX+G2WO6+/EcmdtpOpM+kX/iQBtmpxd7eqla9YAX//bdWN+v33Z+9Xort6FRo0gCVL4LffoE6dqFOhhNL8VD8C24+nUYsQ3n9PP4YkaUq237n2SiSSmnXr1kXtk/fw69KlS44OTySKiwv8+IMbPXpHEtp9PI2/XMNv/JYwnTdoYM3JqVbNqhn11lt22Q4qRdi716rBdfCgNWG8adNop9+5O5g9Lw8hu08qfpnujsrJSVKVbBMoR4mIiIhz6QBHOnr0aNR+efdfme/VZLGHJxUtFQGrXJDfWBcGD40g8r2RvDzsMFPNaQnTecaMVqHN77+HH3+0ivXuecYaVCnd3LlWUurjAzt3Qo0a0U7PYQ7f9SmK6+GSLJ6XSiULJElLMQlUrVq16NOnD/369SN9+vRkyZKFSZMmERQURKdOnfD09KRAgQIsW7Ys6pr7ozVLliyhdOnSpEqVikqVKrF//381aqZPn066dOlYvHgxvr6+uLu7c+bMGW7cuEH79u1Jnz49adKkoXHjxhw/fhyAgIAAUqdOzfLly6PFOH/+fDw8PLh9+zYA//zzD61btyZ9+vT4+PjQokUL/P39o9pHRETQv39/0qVLh4+PD++9916MFcNjkjlzZrJmzRrtdb8UQq1atejXr1+09i1btqRjx46x9hcQEEC3bt3InDkzXl5e1KlTh3379kWd//jjjylTpgxTp04lf/78uLu7xzlWSdkMA774zJmR/4uEz4fyZt9Avokcm3Cdd+9uTXJOk8ba1ParrzTBPL4iI2HIEGulY4sWsHEjPLA7A8BhDtN++hqY0pUJfs7aKFiSvBSTQAHMmDGDjBkzsn37dvr06UOPHj145ZVXqFq1Krt376Zhw4a0a9eO4IeG8gcOHMhXX33Fjh07yJw5M82bN482ghIcHMyIESOYMmUKBw8eJHPmzHTs2JGdO3eyaNEitm7dimmaNGnShLCwMLy9vWnatGm0PfUAZs2aRYsWLUibNi3BwcHUrl2btGnTsmHDBjZt2kTatGlp1KgRoaGhAIwePZqpU6fyww8/sGnTJq5fv86CBQts/wf5ENM0adq0KZcuXWLp0qXs2rWLsmXLUrduXa5fvx7V7sSJE8ydO5fffvvtiVXfRR42aKATE743YXxv+r/pxYfhn2LGXA0l/ooWha1boV8/GDjQesR36lTC9J3cBQRYSdOIEdZ+hD//bCWjDwgkkCb7BhPWYywd3gzjzc56bifJgKN2MY5tS+WgoCBz586dZlBQ0H/HzCBz1xP+CzKDovXToUMHs0WLFlGf16xZ06xevXrU5+Hh4aaHh4fZrl27qGMXL140AXPr1q2maZrm2rVrTcCcPXt2VJtr166ZqVOnNufMmWOapmlOmzbNBMy9e/dGtTl27JgJmJs3/7fD/NWrV83UqVObc+fONU3TNOfPn2+mTZs26n0GBASYqVKlMpcsWWKapmn+8MMPZpEiRczIyMioPkJCQszUqVObK1asME3TNLNly2aOHDky6nxYWJiZM2fOaO/7Yfffk4eHR7RX4cKFo/1Z9e3bN9p1LVq0MDt06BD1eZ48ecxvvvnGNE3TXL16tenl5WXevXs32jUFChQwJ06caJqmaX700Uemq6ureeXKlVhje9D974N58+aZU6ZMMb/55hvz2rVrcbpWkrefZ0WaTs4RJi/9ava+O8CMMCMS9garVplmrlym6e5umkOHmubt2wnbf3Jy5IhpFilimt7eprl0aYxNIs1I88WbHU2j4HGzWJk7ZnCwfUOUFMEheUyS2Ez4CEcoR7nHttnFLsry+A2HS5UqFfX/zs7O+Pj4ULJkyahjWbJkAeDKlSvRrqtSpUrU/2fIkIEiRYpw+PDhqGNubm7R+j58+DAuLi5UqvTfzvI+Pj7RrmvatCkuLi4sWrSI1157jd9++w1PT08aNGhgvZ9duzhx4gSenp7RYrl79y4nT54kICCAixcvRovNxcWF8uXLx+nR2MaNG6P17eLy9N8Ku3bt4vbt2/j4+EQ7fufOHU6ePBn1eZ48eaL2IBR5Wq+3MfBMa9DqlZaMb5GWgPk9mZpmPC4J9c9Z3bpw+DCMHAlffmlVM//qK3j1VTTj+QFLl0KbNta2Odu3W2UiYjDO/JYFHZuT5t/cLF7hRmptdSfJRJJIoIpSlF3semKbJ3F1dY32uWEY0Y4Z9/5xjMskcOOBf0hTp04d7fPYEhjTNKPaubm58fLLLzNr1ixee+01Zs2aRevWraMSmcjISMqVK/fIYz4gQZKQfPnykS6WGZxOTk6PvIfHTfqOjIwkW7ZsrFu37pFzD97j/obHIs+qWTNYsdSFJs3rMrNhWm4s7sgc70mkIYF2L/DwgM8+g06dYMAAq57Rd9/BuHGqYm6a1qO6Dz6wVtj99BN4e8fYdCtb6f/VP7DwbX75PWoHF5FkIUnMgUpDGso+4b8E+4czBn/99VfU/9+4cYNjx45RtGjsCZuvry/h4eFs27Yt6ti1a9c4duwYxYoVizrWtm1bli9fzsGDB1m7di1t27aNOle2bFmOHz9O5syZKViwYLSXt7c33t7eZMuWLVps4eHh7Nr1+EQzLjJlysTFBwoMRkREcODAgVjbly1blkuXLuHi4vJIrBkzZnzmeERiUqcOrFvliseBiiyp+R6Vz7/MP/yTsDfJnx8WLLDqRl25Yq3U69XL2hg3JTp40KrqPniwlUD9/nusydMlLtF8w1dEDh7OwPcjaN7czrGK2FiSSKAc7dNPP2X16tUcOHCAjh07kjFjxsfWmCpUqBAtWrSga9eubNq0iX379vHGG2+QI0cOWrRoEdWuZs2aZMmShbZt25I3b14qP1Cpt23btmTMmJEWLVqwceNGTp8+zfr16+nbty/nz58HoG/fvowcOZIFCxZw5MgRevbsyc2bN+P0nq5cucKlS5eive6PMtWpU4clS5awZMmSOPVbr149qlSpQsuWLVmxYgX+/v5s2bKFoUOHsnPnzjjFI/I0KleGrRvcyHKjKAcrTKP0ji7sxAbfcw0aWIU3v/rKGnEpXBgmTLA270sJrl61EsfSpa3J9YsXw+efW3UmYvAv/1LjnzZcb/0dVWtE8MVnznYOWMT2lEDFwciRI+nbty/lypXj4sWLLFq0CDc3t8deM23aNMqVK8cLL7xAlSpVME2TpUuXPvLIsE2bNuzbty/a6BNAmjRp2LBhA7lz56ZVq1YUK1aMzp07c+fOHby8vAAYMGAA7du3p2PHjlSpUgVPT09efPHFOL2nIkWKkC1btmiv+6NXnTt3pkOHDrRv356aNWuSL18+ateuHWtfhmGwdOlSatSoQefOnSlcuDCvvfYa/v7+UfPKRGylZEnYu92N5/Km53qNBVSdO4a5zE34G7m6wjvvwLFjVvHNnj2hXDmr2nZyLXsQFgZjx0KhQlbiOGqUNQr1UHHMB13jGjUvteZUnclkdcvA/F/ceYYpliKJlhGXCcc2EuONg4ODOXz4MMWKFSNNGts9louLdevWUbt2bW7cuBHrfCGxjfvfB/7+/ty8eZPAwEDat29PhgwZHB2aJFJ370LHN8OZM8sFPvmQj4Y585HxIYatNmrfvt1KqLZssUak+vaFDh2s+VNJnWlak8QHDIDjx6FrV/j0U3hCsd2b3KTGvy9xqPZ4fG4UZMsGVwoUsFPMkpI5ZHWHRqBEJFlIlQp++cmFzz434aNP+eT1wrxypz13uGObG1asCJs2webN1qOtPn0gZ04YNAjOnrXNPe3h0CFro+UXXrBW2O3ZY1Vrf0LydItb1Lv+KofqjyHd1YJsWKPkSZI3myVQhmHkNQzjB8MwThuGcccwjJOGYXxiGMbjn32JiDwlw4ChQwx+/RXcfn+F+bX7UPVSKy5wwXY3rFrV2sLk1Cno0gUmTrQmn7duDQ8s8kj0rl61ksBSpeDECWvy/OrV1udPcJvbNAx4lb0NR5L2fDHWrXKlSBE7xCziQLYcgSp6r//uQHHgHeAt4Asb3jNB1apVC9M09fhOJIl5+WXYvMGFjOfKsL/iFJ7b14Hd7LbtTfPksepGnT8PY8ZY28NUqWLNdJ8925pPlNjcumXNbWrWzBptmjHDqih+8KA1zysOda+CCaZJYGu2N/6ENCdKsW6VCyVK2D50EUezWQJlmuZy0zQ7mab5p2map0zTXAR8BbSy1T0f1LFjRwzDwDAMXFxcyJ07Nz169ODGjRtRbfLmzcuYMWMejJkBAwbg6enJmjVrHumze/fuGIYR7RoRSZzKl4c9293wzZSZq9V+p8qiEcxjnu1vnDYt9O4NR4/CH39Yc6LatIEsWazMbuJEx24TExRkjZi99JL1WK5dO6ssw1dfwcmT1lY27u5x6uoud2ke/BqbX3if1AfLsfpPF+1xJymGvedAeQPXYzsZEhIStZFuQmjUqBEXL17E39+fKVOm8Mcff9CzZ88Y20ZERPDmm2/y448/smbNGurUqRPt/MKFC9m2bRvZs2dPsPhExLZy5ICtG1xp1tCd0JZzeGXUDnqZvW03L+pBTk7WPKLVq60SCH37wqVLVjmAAgWsV48eMH8+xLH8yFO7excWLrQSucyZrceLZ89apQjOnLEmwr/9NsSjSG8oobx4tw1rW/TFbVcV/lzmQoUKtnsLIomN3RaXGoZRAOgDDIitzYgRI5g9e3aM1befhru7O1mzZgUgZ86ctG7dmunTpz/SLiQkhDZt2rBjxw42bNgQrdglwD///EPv3r1ZsWIFTR+zfFdEEh8PD5j/qzNDh5mMeH8U329YxqqpDZmTZRxlKGOfIEqWtF4ffWRtvrtuHaxcCX/+aU3QdnKCChWgfn14/nnIls1KZjJmJF41AIKDrYTo/svf35rP9Oef1uO6UqVgyBBrW5qCBZ/67YQRxishb7DipW64ba7F8mXOVK361N2JJEnxTqAMw/gY+OgJzSqYphlVzc4wjOzAcuBX0zSnxHbR4MGD6dmzJ+fOnYtvWE906tQpli9f/sh2Lrdv36Zp06acO3eOzZs3kzt37mjnIyMjadeuHQMHDqR48eIJHpeI2J6TE3wx3KB6NWjXqT6nSlWkwrTOjGxSg3d4Byd7DsZ7e0OLFtYLrETnfjL13XfWqNCDMmSwkqlMmazRo/v/7+0NFy9aSdL9hOnq1f+uc3a2huDy5LHKLbRuDQ/9cvg0wgmnbVhH/mj9Bi6rG7B4sTM1az5ztyJJztOMQI0HZj+hjf/9/7mXPK0FtgLdHneRu7s7adOmfeR4cDAcOfL4GxYtCg+XjVq8eDFp06YlIiKCu3fvAvD1119Ha/PZZ5/h6enJoUOHyBzDMt1Ro0bh4uLC22+//fgARCTRa9IEDv3tQsfO6Vje9Hfe7TOOxaNe4KfUk8lBDscElSePtXqvSxersrm/P/z7r7V1zL///ve6/7m/v/X/AQGQNSvkzQvPPWdN+s6T579XjhzxG72Kg+tc543wjixv2w7npS+wcKET9eol6C1Ekox4/+0yTfMqcPWJDQHDMHJgJU+7gE6maT5Vud4jR6yCv4+za5e1TdWDateuzYQJEwgODmbKlCkcO3aMPn36RGvToEEDVq1axRdffPHI5PBdu3YxduxYdu/eHW2zYBFJurJkgaWLnfHzg/7v9mLTmgYUm9WaaaXe4SVecmxwzs7/zY9KZHawg1bXu3C53ZcYK+rx6zwnmjRxdFQijmPLOlDZgXXAOeBdIJNhGFkNw8ga376KFrUSpMe9Ytrb18PDg4IFC1KqVCnGjRtHSEgIn3zySbQ2devWZdGiRUyaNOmR5Grjxo1cuXKF3Llz4+LigouLC2fOnGHAgAHkzZs3vm9DRBIJw7AWyu3e6Uxhp4IEVVzDy2M30DmyC7dJuIUsyYGJyXjGU2V7X66UXYbHX3VZstiJx2wHKpIi2HISeQOg4L3X+YfOxWs4J02aR0eXnsZHH31E48aN6dGjR7TVdPXr12fx4sU0a9aMyMhIxo8fj2EYtGvXjnoPjU83bNiQdu3a0alTp2cPSEQcqkQJ2LXdhfcHm4ztN5Yfl61kzbT6zMk2hkpUcnR4DhdIIF3Mrsz1y4hT/w2UKevErxuceGiqqEiKZMs6UNNN0zRietnqnk9Sq1YtihcvzhdfPFrLs06dOixZsoQZM2bQq1cvTNPEx8eHEiVKRHu5urqSNWtWiqjMrkiykCoVjPnGYNkySL+3Nv+UWkrVP0bwNm9zjWuODs9h9rOfsoE1mf/aK9BnPH16urBRyZNIlBS3F17//v2ZPHlyjCv9atWqxdKlS5k5cyY9evTAgRsti4idNWoEh/a70KCyN5HNFzLh1VrkO1WX0YwmhBBHh2dXM5hBhf2dOVt+Pu7LWjJ3rlVc3U0bcYlEMRyYJMR44+DgYA4fPkyxYsVI8/CyOkkx7n8f+Pv7c/PmTQIDA2nfvj0ZMmRwdGiSzJkmzJwJgwZHcOVqJGbfseQZMpP/eQ/lZV7GcMzG73Zxhzv0oQ8/zAjDucckihVy5bdfnShc2NGRiTyWQ/5SprgRKBGRxzEMaN8eThxz5sMPXEnl158LBTfw6oQ1VAuvyTa2OTpEmzjKUSrdqcX0LtWg4wzav+bOtq1KnkRiowRKRCQGHh5W4fDjx5x4/QVvjF7fsbvUNCov+5g2tMH/v3J3SZo//nSlK8WPvcjRKtNx+bkDU6fC1KmP1tYTkf8ogRIReYwcOWDaNNi506BS5vzQZBnzG3Wj8IFWvM/7XI99e89EzR9/utGNgv71+KlbDczi+8kVVIRtfzmhRcYiT6YESkQkDsqWhXVrDRYsgFwnaxJeeidfvVWQbOcq8jqvs5rVRPJUtYLt6gxn6E53Cp2px8zu1aHwUdIufINRI5zZt9eJ0qUdHaFI0qAESkQkjgzD2jHl0EEnvh7tRNo5nQnLe4xFTbtRb8F4CoQV5TM+4/wjpe8c7yxneYu3KHi2Dj++VQWz0FE85rdjxHBn/E8bvPuu9dhSROJGCZSISDy5uUG/fnD2jBMTv3ei+NWa0GoBl3Lu4NNBHuQ+Vo8mNOE3fiOUUIfGepaz9KAHBc7VYkbPipgFj+HxWwe++NxKnAYOVOIk8jSUQImIPCUvL+jaFbZtM9i3D7q95k3aye9gFjnCppqf8/LMheQILsS7vMs2tnGHOzaP6RrXmM98+tCHEpQgz7lqTOtVHrPgMdLM7cTwz6zE6b33IIa920UkjlQHShIl1YGSpOruXViwAKZMgTVrwM3rDkbbXwh5ZSZOZffh652DspTlOZ6jLGUpQxm88Hrq+93kJhvYwNp7/+27fg7W18B7XQuc1tbjxv6cpM9gMvBdg969wdMzAd+sSOLgkDpQttwLT0QkxUmVCtq0sV4nT8LUqamZNq0TFyd0JhI4V+gyF8ru4+ey64go+yk8t4eCPhmikqpCFCKSSMIII5TQaB8f/P+b3GQTm9h98xTmhup4rmuO69p+GPtyY5oGPvmhVi2o/R60aGEocRJJYIl2BCpv3rykTp3a3jFJInHnzh38/f01AiXJQkQEHD4Mu3f/99qzx+T2besXZ88813Ere4DAcusILXQAIp0gzDXq5RSWCpewVDiFpcIpzB3nsFQ438yA06Ya3NiTFzPSidy5TWrXNqhd20qc8uRx6FsWsSeHjEAlugQqNDSUgwcPEhmZ+JcDi22ZpsmpU6cIDAxUAiXJTmQknDgRPanavdvkxo3oPwsMw8TNzcDVlWivNGmgYkWoXdt65ctnrRIUSYGUQN0XGhpKeHg4Fy9epH///rRq1Yq0mu2Y4kRERBAeHs6dO3eUQEmKYJpw8ya4uPyXKDk7OzoqkURPc6Duc3Nzw83NDRcXF86ePcuNGzcIDw93dFjiICEhIY4OQcQuDAPSp3d0FCISF4kygbrP3d2dkJAQQkNDCQwMdHQ44kBeXl64u7s7OgwREREgkT7Cu+/WrVv4+Phw4sQJPLWEJEVzd3fHQ9X+RETkUXqEF5Pw8HDSp0+Pl9fT10kRERERSUiqRC4iIiIST458hPdEhmF4AQGAt2matxwdj4iIiAgk/gTKADyBQDMxByoiIiIpSqJOoEREREQSI82BEhEREYknJVAiIiIi8aQESkRERCSelECJiIiIxJMSKBEREZF4UgIlIiIiEk9KoERERETi6f+vnmc6HDBtUQAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 3 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Plot first component of solution versus time, as simple example.\n", "yp = yv.column(0);\n", "p1 = line(list(zip(tv,yp)),rgbcolor=[1,0,0],legend_label='Euler method')\n", "yp = yv2.column(0);\n", "p2 = line(list(zip(tv2,yp)),rgbcolor=[0,1,0],legend_label='Improved Euler')\n", "yp = yv3.column(0);\n", "p3 = line(list(zip(tv3,yp)),rgbcolor=[0,0,1],legend_label='RK4')\n", "p = p1+p2+p3;\n", "p.set_legend_options(loc='lower left');\n", "p.axes_labels(['$t$','$y0(t)$'])\n", "show(p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.2", "language": "sage", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }