{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

# Worked Solutions 5: Numerical methods for ODEs: (a) explicit methods

" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "plt.rcParams.update({'font.size': 14}) # Make the labels larger" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

### Question 1: Difference equations

" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0, 2)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEZCAYAAACNebLAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4lFX2wPHvmwAhhPQAoYYaCC0BAtgJrqjsyir+1o6r\nqKgUV111XXWVIHZREBVQVEARsKIUUUAYRIr0EhJKyqT3THqbzNzfHy8gYihJ3plJ4HyeJw+STN57\nYZc5ueWcoymlEEIIIU7n5uoJCCGEaJwkQAghhKiVBAghhBC1kgAhhBCiVhIghBBC1EoChBBCiFpJ\ngBBCCFErpwYITdPcNE2brmlaoqZpFcd/na5pmgQqIYRoZJo5ebz/AhOBfwIxwEBgEVAJvOzkuQgh\nhDgLZweIS4GVSqkfjv8+RdO0lcBwJ89DCCHEOTh7a+dXYKSmab0BNE3rC1wNrHbyPIQQQpyDU1cQ\nSqnXNU3zBmI1TbMB7sDLSqkPnDkPIYQQ5+bUAKFp2u3A3cDtQCwQAczWNC1JKbXAmXMRQghxdpoz\nq7lqmpYCvKGUeu+Uzz0H3KOUCq3l9VJqVggh6kEppTX0Gc4+g2gF2E/7nP1s81BKyYdBH1OnTnX5\nHC6UD/m7lL/PxvxhFGffYloJ/FfTNDNwCBgMPA4sdPI8hBBCnIOzA8QUYDrwPtAWyAQ+OP45IYQQ\njYizbzGVAf8+/iGcLCoqytVTuGDI36Wx5O+zcXLqIXVdaZqmGvP8hBCiMdI0DdUED6mFEEI0ERIg\nhBBC1EoChBBCiFpJgBBCCFErCRBCCCFqJQFCCCFErSRACCGEqJUECCGEELWSACGEEKJWEiCEEELU\nSgKEEEJcQKpzqw17lrOruQohhHAAi8lC7le5ZC/JNuyZsoIQQogLgQ1yv86l1+xehj1SqrkKIUQT\nl/lJJonPJNLvy374jfAzrJqrbDEJIUQTpeyKpOeTyPkih0G/DKJV71aGPl8ChBBCNEG2ShtHxh+h\nMqWSwdsG06JNC8PHcOoZhKZpSZqm2Wv5WOnMeQghRFNWnVfN/r/sRylF+M/hDgkO4PxD6kgg+JSP\nwYACvnDyPIQQokkqP1rO3kv34hflR98lfXFv6e6wsZzdkzr/1N9rmjYBKAK+cuY8hBCiKSr8pZBD\ntxyi+yvdaX9/e4eP5+oziPuAz5RSVS6ehxBCNGpZi7NI+HcCYUvCCLgmwCljuixAaJp2LdAVmO+q\nOQghRGOnlCL5xWQyF2QSsTECr35eThvblSuICcBOpVTM2V4UHR198r+joqKIiopy7KyEEKKRsFfZ\nOTLhCOWHyxm8fTAewR61vs5kMmEymQwf3yWJcpqmtQHSgIlKqU/O8jpJlBNCXJSsBVZibo6heUBz\nwhaH4d7q/A+jjUqUc1WpjfFAJbDMReMLIUSjZDFZqEioYM+le/CO9KbfV/3qFByM5KotpvuBpUqp\ncheNL4QQjVL2Z9kU/FBAyAshdJzY0aVzcXqA0DQtCugJ3OnssYUQojHLXpZNzrIc+n3dj8DRga6e\njhTrE0IIV7NssGCONlN6oBRbkY2QqSEA+EX54R/lX+fnSbE+IYS4ANSU1pD+XjooGH50OOlz0ukW\n3c3V0wKkH4QQQrhMZXIley/fSzO/ZoSvD6dFW8fUVKovCRBCCOECRVuK2HPJHoLvDab3x71x89Df\njv2i/Fw8s9/JGYQQQjhZ5sJMEv+TSJ9FfRxyGC1nEEII0cQomyLh6QTyvssjYlMEXmHOK5tRHxIg\nhBDCCWqKaoi9MxZ7pZ0hvw2heWBzh4yzwWIx7FkSIIQQwsEqEio4OOYgfiP96DmrJ27NjT/+3Wix\n8HpKCr8WFRn2TDmkFkIIB7JstLDn8j10nNKR0PdDHRIcKm02FmVlkVldTczQoYY9V1YQQgjhIBkf\nZJD0QhJ9l/bF/+q6J7yd1xhVVYyNiSGkZUu2Dh6Ml7txdZskQAghhMHsNXYSHk/Ast7CoC2DaNWz\nlUPG2V5UxD8OHWJSx44806ULmtbgi0t/IAFCCCEMYjFZaB3emthbY9GaaQzePphmvo55m12QmcnT\niYl80rs3NwQFOWQMCRBCCGGQ3K9zOfrgUQLHBNLjjR5o7sb+RA9gtdt5MiGBNQUFbIqIIMzLcVdl\nJUAIIYQB8lfnk7Ugi17v9qL9fe0dM4bVyq2HDtHCzY0dgwfj19wxV2VPkFtMQgjRAJYNFvaO3Muh\nOw5hL7dTmVJJUnQSFpNx+QgAB0pLGbp7N5He3qwaMMDhwQFkBSGEEPVWU1xD2uw0VJVi+OHhZHyY\n4ZBKrN/k5vLw0aPM7tmTO9q1M/z5ZyIBQggh6qHscBmHxh7CL8qPfl/2w62F8RsydqWINptZlJXF\njwMHMsTb2/AxzsbpW0yapgVrmrZQ07QcTdMqNE2L0TTtSmfPQwgh6ivv+zz2XbWPzk92JnRu6Mng\nYGQl1uKaGsbGxGAqLGTnkCFODw7g5BWEpmm+wBbgF2A0kAd0B3KcOQ8hhKgPZVeYo81kLcxiwKoB\n+Azz+cPX69P97XQmi4WOHh7cGBPDCD8/vurXjxZurjkudvYW09NAhlJq/CmfS3byHIQQos6shVbi\nxsVhK7YxZOcQWrRzTHOfjzIzWWex8GK3bjzUoYNDxjhfzg5LNwK/aZq2TNO0bE3T9mqaNtnJcxBC\niDopO1TGnqF78OzhSfjP4Q4JDkop3khJ4fu8PL7u18/lwQGc3DBI07QKQAEzgS+BCOA94Gml1Jxa\nXi8Ng4QQLpXzdQ7HJh6jx1s9CP5nsEPGWJWXx1MJCZTYbKRXVzM1JASAKD8/ovzrvm1lVMMgZweI\nKmCHUurKUz73MnCTUqpfLa+XACGEcAllUyT9L4nspdn0/6Y/3kMcc0gcU1rKzYcOMcrfn7d79uTV\n5GSiuzXsqmxT7SiXCcSd9rk44F9n+obo6OiT/x0VFUVUVJQj5iWEECdZC6zE3hGLqlH6eUMbx5w3\nLM3O5l/x8bzdowd3B9d/dWIymTCZTMZN7DhnryA+BzoppUac8rnpwFilVP9aXi8rCCGE01hMFpr7\nNyfm5hiCxgbR/bXuuDUz/qi2+ng9pR/y8/mmf3/CW7c++TWTxVKvbaVTNdUtpkj0a67TgC+AwcB8\n4L9KqXm1vF4ChBDCaWL+L4aiX4roObsn7e5wTMZyRlUVtxw6REDz5nzapw/+DiiZ0SQDBICmaaOB\nV4FQIAV4Vyn1/hleKwFCCOFwdqudxP8mkvlJJoNMg2gd3vrc31QPmwoLuSM2lsnH+ze4Gdy/4YQm\nGyDqQgKEEMLRcr7KIeGJBLQWGpUJlYRM1W8Q+UX5GZL4BvoV1rdSU5mRmsqnYWFcGxBgyHPPpKke\nUgshRKNh+dlC/KPxdJzckS7PdMH8otnwYnslNTWMP3yY5KoqfhsyhJCWLQ19viNJgBBCXHSUXZH8\ncjIZczMIWxzmsH7RsWVl3Hy8ZMbisDBaGtgv2hkkQAghLirVedXEjYvDXm5nyK4heHTwOPk1I4vt\nfZGTw5Rjx3ije3fGt3dMAyFHkzMIIcRFo2h7EbG3xdL29rZ0e7mb4VdYTRYLl/v68p/ExJMlMwa7\noAqrHFILIcR5UkqRPjud5FeS6T2/N0F/D3LIOE/Ex7OjpARvd3cWh4UR4ISub7WRQ2ohhDgPNcU1\nHLn/CBWJFQzePhjPbp4OGcdksfBhRgZPdunC8yEhDrvCek4rVhj2KOlJLYS4YJUeKGV35G6aBzVn\n0JZBDgkOP1ssjNy7lxsOHqTUbkcpxYtmMyaLsT2pz2ndOrjkEvjnPw17pKwghBAXpMwFmST+J5Ge\ns3rS7i7HZEVnVVXxSnIyCjg6fDgfZmQ0uNBevSQlwXPPQfv28MMPEBhoyGNlBSGEuKDYKmwcvv8w\nqW+kErEpwmHB4WeLhcG7d3O5ry/rw8Pp4OFx7m9yhG+/heHD4c474bvvwMAkPFlBCCEuCBaTBY+O\nHhz6xyG8+nsxeOdgmrU2/i3OdnwLaX5mJp/26cM1p7whR/kZd032nKqq4MknYfVq/WPoUMOHkAAh\nhLggZMzJoNBUSNcXu9LhoQ5oDjgkzqiq4q64ODRg95AhtD9t1dDQKqznLT4ebrsNunaFPXvAQYFJ\ntpiEEE2avcrOsceOUbC2gAE/DKDjwx0dEhzWFhQwZPduovz8WBce/qfg4DRffgmXXQbjx8PXXzss\nOICsIIQQTVjWZ1kk/CeBZn7NsBXZyF+VT/6qfEML7dXY7USbzSzIymJJWBgjnbVKOF1FBfz73/pt\npTVrYMgQhw8pAUII0SRlfZpFwhMJdJ2ubymZpxlfaC+9qoo7YmNp6ebGnshI2rVwTGe5czpyBG69\nFfr0gd27wdfXKcPKFpMQokmpKakh7u44Ul5LIXxDuMO2lH7Mz2fIrl1cFxDAjwMHui44fP45XHEF\nTJwIy5Y5LTiArCCEEE1Iye4SYm+PxW+kH0N2DcG91e/VUY0otHeiltLzSUkszs7mi379GOHMm0mn\nKi+HRx+FTZv0baWICKdPQQKEEKLRU3ZF2qw0Ul5Lodd7vWh7a9s/vcaIM4fv8/J4LikJb3d39kZG\n0sYVqwaTCdq107eUBg7Ut5RcUPAPnLzFpGnaVE3T7Kd9ZDhzDkKIpqU6t5qDYw6S82UOg38bXGtw\nMMLq/HzmZ2YyJjCQHwYOdE1wAHjnHbjqKn31sHixy4IDuGYFcRgYAZzYNLS5YA5CiCbAssFC3D/j\naDeuHd2md8OtufE/067Nz+d5s5nD5eWU2e1U2u28aDYT5efnvLwGgOJimDIFNm+GjRthwADnjX0G\nrggQNUqpXBeMK4RoIuw1dsxTzWQtyKLPoj4EjHJMD+eY0lKeSEwkrFUrfhw4kHfS0lxTS+n99+F/\n/4Pu3SE/H775Rv+IitI/XMQVAaK7pmnpQBXwG/CsUirJBfMQQjRClcmVxN4Zi3trdyL3RtKinfFb\nPUop3ktP58XkZN7o3p17g4MdchPqnGpq4JVX9ADxyScwdixER+sfjYCzA8R24F70baa2wPPAVk3T\n+iqlnFwbVwjRWFhMFvyj/Mn9NpejE4/S+cnOdH6iM5qb8W/a2dXV3Hf4MLlWK1sHDaJXq1Ynv+bU\nWkpmM4wbBy1b6uUyOnZ03tjnyakBQin106m/1zRtO5AE3APMqu17ok+JpFFRUUS5cLklhHAMyzoL\nuV/k6uUyVg7AZ5iPQ8b5IT+f+48cYXxwMNO6dqW52x/PNJx25rB0qX4I/dRT8MQTcOo86vEeZzKZ\nMJlMhk3vBJe3HNU0bQMQp5SaXMvXpOWoEBe40phS9l+9H/+/+BM6L5Rmvsb/3Fpps/F0YiLL8/L4\ntE8f5x4+n6q4GCZPhp07YckSGDzY8CFyy3Jp27qtIS1HXZpJrWlaS6APkOnKeQghnM+ywcL+6/az\nZ9gerLlWPEM9SZ2ZisVk7G5zTGkpw/bsIaO6mn2Rka4LDtu26cluXl56boPBwcFkNnHnN3fS/Z3u\nhj3TqVtMmqa9CawEUoB26GcQrYBFzpyHEMK1KtMqSX4lGXu5naEHh5L1WZbhdZSUUryfns60xnIQ\nPWcOzJsHN91k+BBl1WV8EfMFW1O3surOVUQ9F2XIc519SN0JWAIEAbnoh9aXKKVSnTwPIYSLZC/L\nJv5f8XR6tBOdn+6MWzPjNzJyqqsZf/gwObUcRDvV6QfRHToYPsSO9B2M+3Ycl3S6hP0P78e3pXG1\nmpy6xaSUukMp1Ukp1VIp1VkpdYtS6rAz5yCEcA1roZXYu2IxR5sZ8MMAQp4LORkcjKqjBLAmP5+I\nXbsIb93a+cHh1IPiJUtg2DB9xbB2reHBocZewzTTNMYsHcPLV7/Mp2M/NTQ4gNRiEkI4gWWjhcP3\nHiZwTCCReyL/UGQPjKmjtN5iYXleHsvz8lgSFuaaswaTST9bOHEQ/eOPDjmIPpZ/jHHLx+HX0o+9\nD+2lg7fxKxOQct9CCAeyV9lJeCqBuHFxhH4QSuh7oX8KDkaIKS1lfmYm6VVVrj2ITk39/SB6zx7D\ng4NSig92fcClH1/KuAHjWHPXGocFB5AVhBDCQUoPlhI3Lg7Pnp5E7o+kRZDxGdE/Wyy8lpzMlqIi\nKpSin5cXs9PSnFtHyWSCn3/Wy3Jv3qz3ig4Ohh07DC2TkV2azQMrHyC9OJ3N4zcT1ibMsGefiQQI\nIYShTpbmfjWF7m92J/gex9weSqioINpsxg04NGwYi7KyXFNHKSAAVq6Ezp31pLcZMwwf4vvD3/Pw\n6oe5L+I+vrn1G1q4O6fSrAQIIYRhKlMrOXzPYezVdgb/NhjP7p6Gj6GUYl5GBs8nJfFcSAiPduqE\nm6uur86YAW+9BW+8AffeC9OmGTpESVUJj//0OBuSNvD1LV9zeZfLDX3+ucgZhBCiQU4ktmUvyWb3\nkN34X+PPoE2DHBIc0ioruf7AAT7JymLzoEE83rnzyeDg1DpKR4/ClVfqt5N27YLx40HTDNlSMplN\nAGxN3cqgDwahlGLfw/ucHhygEZTaOBsptSFE45fwdAJVKVWU7islbHEY3kOMb3CjlOLz7Gz+nZDA\nIx078t8uXf5UR8kp7Ha98uq0aTB1qn5byeB5PL/heRSKj/Z8xNy/zWVs2Ng6P0PTNENKbcgWkxCi\n3grWFpAxL4PgfwYzZPcQh9xQyq2u5uGjRzlSXs6PAwcy2FUd1pKT4b779F7RW7dCaKjhQ8TlxvHx\n3o+JCI5g38P7CG4dbPgYdSFbTEKIOstbnceuwbs4dMshbMU2mgU2I+WNFMPrKH2Xm8vAXbvo6enJ\n7shI1wQHpWDBAoiMhFGj9JtKBgeHnxN/5ppPr2HIh0PILM1kaIehzNs17+R2k6vICkIIUScFaws4\nNukYAdcFEGGKIPXtVMPrKBVarTwaH8+WoiK+7tePy32NzRA+b1lZ8OCDkJKiX2UdONDwIeJy43h2\nw7N4Nffi0KRDLNq/iOioaMPHqQ9ZQQghzktNcQ1HJhzhyIQj9J7fm94f9qaZj/E/Y64rKGDgrl20\ndndnX2Sk64LDl19CeLgeFHbsMDw41NhreP3X17lywZXcG34v6/+5nm7+LrimexayghBCnFPB2gKO\nTDhCwHUBDD049A+Bwag6SkN9fPhPQgIr8/P5qHdvrg1wTB/qs0/EpAeCyZNh715YsQKGDzd8mLjc\nOO79/l68mnuxc8LOPwSGqK5Rho9XX3KLSQhxRjXFNSQ8mUDBTwX0nt+bgGsd86Z9/+HD/FJUxKU+\nPszu2RO/5s0dMs453XWXHiRuvRVefhkMLvRXY6/hra1v8ebWN5k+cjoPRT6Em2b8Ro7cYhJCONTJ\nVcO1AQw9MNQhnd7KbTamms18kZPDZ2FhjG3TxvAxzktRkZ4FvXo1fPedoSUyTojNjWX89+Np3aI1\nux7cRVe/roaPYbTz/l/8ePe3S4GugCd6P4c9SqkEx0xNCOEKzlo1zEpN5UWzmY4eHpTZ7ewvLWV/\naalz6yiB3szntdf0m0lFRfoKwmTSg4QBgeLEqmHGthlMHzmdB4c86JBVgyOcM0BomnY58CgwBmgO\nFAEVQADgoWlaIvAhME8pVeLAuQohHMwZq4aimhqeTkhgdUEBi8LCGBMURHRSkvPrKOXkwKOP6pnQ\nK1bowSA6Wv8wyKmrhp0TdjaJVcOpzhrGNE1bAXwJJAPXAt5KqcDjTX9aAb2Al4C/AEc1TRtVl8E1\nTXtG0zS7pmmz6zd9IURDWUwW/YbSg0c48sARen/Ym97zezskOKzOz6f/zp3YgZihQxkTFGT4GOek\nFHz+OQwYoBfY27/fsC2lE3kLJ24ojVg4gvsi7mP93eubXHCAc68gfgT+oZSqru2LSqlEIBFYpGla\nP+C8C5NrmnYJMAHYf77fI4QwXubHmRy+57C+ajjomFVDbnU1j8XHs724mEV9+nD1aVtITqujlJoK\nEyfqeQ2rV+vJb3+YSFSDHm8ym2jr1Zbx34/Hu4V3k1w1nOqsKwil1JwzBYdaXntIKbXufF6raZov\nsBgYDxSez/cIIYxlLbRy5MEj5K/Id9iqQSnFsuxsBuzcSXCLFhwcOvRPwQFw/JmD3Q7z5ukNfIYP\n17eVTg8O0KAAUWOv4deUX0+uGtbdva5JBweo2yF1IjBUKZV/2uf90A+ru9dh3A+BL5VSmxxRJ14I\ncWZKKZJfTCZlRgqterfCVmyjaFsRRduK8IvyM6T9J0B6VRWTjh4loaKC7wcMYLiPjyHPrbNjx2DC\nBKis1A+f+/Uz9PEms4klB5ew4sgKssuyeXT4o2SWZrIpeVOjymmoj7r8uNAVqK0SlwfQ8Xwfomna\nBKA7cEcdxhZCGKAytZJjU45RcbSCgWsG4neFH0nRSYaWylBK8VFmJs8mJTG5Qwe+6tePFq6ovFpT\nAzNnwuuvw//+B488Au7GFhMsrS5l5ZGVrDiygjdHvUl8QTzTRhrbE8KVzucW082n/PZvmqYVnfJ7\nd/QDavP5DKZpWijwMnC5Uspeh3kKIRpA2RTpc9MxR5vp9Egn+n3ZDzcP49+0EyoqmHDkCCU2GxvC\nwxnQurXhY5yXAwfg/vvB11cvk9G9Lhsc52fNsTVMXD2Rq0KuImZSDEGtgog2RRs+jiudzwri6+O/\nKuDj075mRQ8OT5zneJcCgUDsKVtL7sBVmqY9DHgppaynfkP0KVfOoqKiiHJAAosQF7LSg6UcffAo\nuMOgzYPwCvP6w9cbWirDZLFwpZ8fs9PSeDk5mWdCQni0Y0eaOXvVYDLBpZfqGdDz5um5DSca+Rgo\nuzSbx356jB3pO5g/Zj6jevx+edNVW0omkwmTyWT4c8+71IamaUnoZxB59R5M03yATqd9eiFwFHhZ\nKRV32uul1IYQ9WSrtJH8UjKZH2TS7aVutJ/QHs3N+DO/SUePsqekhJZubnzUuzc9DS5Pcd4eeAC2\nbYNevWDOHOhw3pcqz4tSigX7FvDMz88wPmI8L4x4gVbNXfRnPQenl9pQSjV4k1IpVQzEnvo5TdPK\ngILTg4MQov4sJgtHHzyK10AvIvdH4tHBw/AxKmw2XklJYWFWFjN79mRC+/au6Q1dXAzPPw/Llul9\nG/7xD8NXDUfzj/LQqocorS7lp3E/EREcYejzG6tzJcqN087zmpGmaSGapl1ZjznIEkEIg1gLrBx+\n4DCH7z5Mjzd70P/r/g4JDm+kpNBh61a+zc2lwm4ns6qKF81mTBZjGwadlVJ6689OnfQmPmVlcOiQ\n/jmDtluqbdW8/MvLXPbxZdzY+0a237/9ogkOcO4VxH3AC5qmLQRWAjGn7vlomhYAXAGMA6KOv75O\nlFJX1/V7hBB/pJQi98tc4h+PJ+jmIIYeGuqQXg2ZVVU8Hh/PjpISlvTty+jAQNeUyUhMhClT9IS3\nH36AK64wvEzG9rTtTFg5gS6+Xdj94G5C/EIMe3ZTca5EuauBx4ERwD6gVNO0JE3T4jRNy0Iv2PcB\nEA/0U0qtcvSEhRC/s5gsVKZUcnDMQczTzfT7ph+h74UaHhxsSvFeWhoDd+2ih6cnMUOHMjow0NAx\nzkt1tV5cb9gwGDEC9uzRg4NBTGYTxVXFTPlhCjd/cTPPX/U8q+5YdVEGBziPMwil1GpgtaZpQcCV\nQBf0aq55wF5gr1xZFcL5lE2R8loKJbtK6PRYJ/p/2x+3FsbfHNpTUsJDR4/Sys2NTRER9PX64y0o\np5XJ2LRJL5PRvbueCd216x+/bsANx3m75rEldQvX97ieQ5MO4e/pxKqyjVBdbjH9n1LqmzN87Wml\n1OuGzgy5xSTEmRTvKObopKNY86yE/xROq97G36Yprqnh+aQkluXk8Hr37twTHIxLKh/k5sJTT+k9\noWfPhptuMvwQOrUolcd+egyT2cQ3t37T5DOgjbrFVJcfNxZrmvaRpmkn/5+oaVonTdM2om9DCSEc\nzJpv5cANB9h39T48OntQlVxF9tJskqKTsJiMOSBWSvFVTg59d+yg1GYjdtgw7m3f3vnBwW6Hjz7S\nS2MEBkJsLIwda2hwqLZV8+DKB+n9Xm/yy/MpqCjAZDYRbYo+WZn1YlaXjcrhwBJgn6ZpdwI9gLnA\nb0C4A+YmhDhO2RWZH2eS9L8k2t7WlrDFYTT3a254mYzEigqmHDtGcmUlS/v25UpnbR+d7uBBePhh\nsNlg7VqIMP7m0PrE9Uz5YQo9A3oSMymG7v7diTZFEx0VbfhYTVVd8iAOaJoWCcwBtqFfT31SKSW9\nHIRwoJLdJRyddBTNXWPgTwPxjvA29Pkmi4XLfH2ZkZrKW6mpPNW5M//u39+59ZNOdHArK9OvqS5c\nCNOn60X2DJ5HWnEaT6x9gh3pO5h9/WzG9B5j6PMvJHW96hCOfqMpHugMDNM0zVs6yQlhPKvFStJz\nSeR+m0v3V7sTfE/wnzKhG1omA2BRVhaTjh2jW8uW7BoyhG6eng1+Zp2ZTHrC2yOPwJVX6iuIdu0M\nHaLaVs0729/h9S2vM2noJBbcuOBPmdBN/ezBaHU5pH4BeA54H/gv0A34HAgC7lZKbTZ8cnJILS5C\nyq7IWpRF4jOJtLm5Dd1e6kbzgOaGj5NdXc3TCQl8m5vLgrAwbg4Kcs0hdHIyXHednvg2Zw785S+G\nD7ExaSOTf5hMV7+uzB49m54BPQ0fozFxeqkN4GFgjFJq7fHfHzneFe4lYD162W8hRAOU7Cvh2ORj\nqBrFgFUD8Ik0vodCtd3OY/HxLMzMJKJ1a0rsdg6WlnKwtJQoPz/HN+854aef4NVX4bff9F4Nzz2n\nZ0S7uxvWAjSjJIMn1z7J1tStzLp+Fjf2vtE1QbCJqkuAGHh6oT6lVA3wX03TfjB2WkJcXKyFVswv\nmMn5IkcvrHe/YwrrrS0o4NH4eEI8PNgTGUkfLy/nZ0IrBd98A08+qSe8HT6s11AyMAvaarPy7o53\neWXzKzwc+TAf/f2jRltYrzGryyH1Gau4KqV+MWY6Qlw8LCYLfiP8yF6cTeLTiQTeEMjQQ0NpEdTC\n8LESKyr4d3w8B8vKmNWzJzcEBrrmJ+mYGPjXv/TchgULYORIwx5tMpuI6hrFJvMmJv8wmY4+Hdl6\n/1ZCA0N6NOEPAAAgAElEQVQNG+NiY3yxFiHEecn5IgfzC2ZsZTb6L++Pz3Djt5PKbDZeTU5mbkYG\nT3buzLK+fWl5Wlc1p2RCWywwdapecfWFF/QrrM1OefsxYEtp5ZGVfLTnI35J/oWZ183k5rCbZTup\ngVzQB1CIi5u1wMqxR4+R/Wk2be9oy5AdQwwPDkoplmVn02fHDpIqK9kfGckzISF/Cg6AY88cbDb4\n8EPo0wesVj3ZbcqUPwYHaFCAqLZVM3PbTObumktnn87ETY7j//r+nwQHA8gKQggnsVvtxD8RT9bH\nWbTq2wp7uZ3q7GrM0834RfnhH2XMG/X+0lL+dewYxTYbS8LCXJfstmWLfm3Vywt+/BEGDTL08Uop\nXt38Km9te4sAzwAqairwaObBm1vfJKprlFxZNYAECCGcIH9NPgn/TsCjkweDfxtM6/6tDc+Czrda\neSEpia9zc5nWtSsTOnTA3RU/Raenw9NP68X13ngDbr/d8NpJB7MP8u+1/ya9OJ3P/+9zru95vWRB\nO4BsMQnhQGWxZRwYfYD4R+Pp/mZ3Bq4dSOv+rQ17vsliwaYUc9PTCduxAw2IGzaMhzt2dF5wONGc\np6pKv7YaHq5XWo2LgzvuMDQ45JblMnHVRK757Bpu6n0T+x/ez/U9rzfs+eKPZAUhhANY862Yp5nJ\nWZpDl2e70H/yn0txG5UF/Vh8PH7NmrE+PJyBrY0LPudt40YoKYHHH4f+/fW8hh49DB2iqqaKd3e8\ny2u/vsbdA+/m8OTDfyrFLVtKxjvvTGpDBtO0ScBDQNfjnzoEvKSUqjWPQjKpRVNjt9rJmJtB8kvJ\ntLmlDV2ndXXYtdVnEhP5saCA+b17c0ubNq45lI2Nhb//XT90fucdPSPaQEopvj/yPU+ufZKwNmHM\nGDWD3kG9DR3jQuSKTGojpAL/AY6hb2/dC3ynadpgpVSMk+cihGGUUhSsKSD+3/G07NKS8A3hhm4l\nnVBgtfLwkSOszM/nUh8fim02YsvKmFZW5tws6G+/hZdf1pPcysvhf/+DbdvAw8OwLOj9Wft5/KfH\nySnLYc7f5nBtj2sNea6oA6WUSz+AfGDCGb6mhGiMCjYWnPzv0kOlat91+9T20O0qb1Westvtho9X\nabOpt1NSVJtff1UPHT6sMisrlVJKTU1MNHyssyotVWraNKUCApR64gmlCgqUmjq1wY/dmLTx5H9n\nlWSpCSsmqLZvtlVzdsxRVpu1wc+/2Bx/72zw+7PLDqk1TXPTNO12wAvY6qp5CFEfhaZCrPlWjj1y\njH0j9hFwfQBDDw4l8G/GZiirU5r3/GyxYIqIYF7v3gR7OLn0mc0Gn3wCoaH64fOuXTBjBhi0YjGZ\nTVTVVPHGljfoN6cf3i28OTLlCBOHTqSZmxyVuorT/+Y1TeuP3k+iJVACjFVKHXL2PISoL3u1neLt\nxex4fwdtbmvD0DjHlMfYWlTEkwkJVNrtzO/dm6treTN2Shb0Tz/pLT/9/GD5cr1+0h8mEdWgxyul\niM2NJez9MAa2G8i2+7fRK7BXg54pjOGK0HwYva+EL/AP4FNN00YopWJre3H0KQW8oqKiiDJof1OI\nurJssJD+fjqWDRZshTY6TOxA86DmlMWU0SLKuAARX17OM0lJ/FZczMvdunFXu3a4nWFV4tAzhwMH\n9MCQlKTnM9x4Y+1XVuv5b9JkNvHp/k9Zn7ie1OJU7h54N939u5Neki4Boo5MJhOmE9eNDeTUW0y1\nTkDT1gFmpdSEWr6mXD0/IQAK1hWQ+HQiWjON7q93p3BToaFJbqAnuk03m1mcnc0TnTvzWKdOeNZS\nGsPh0tPh+edh9Wr914cegubG9qOIzY3lmZ+fYV/WPl6MepEESwIvjnzR0DEuZkbdYmoMiXJuSC8J\n0UiV7C5h/6j9HJt8jC7PdmHwb4PxH2nQvrvFAkClzcaMlBT67NiBVSlihw3jmZAQ5wWHEz95lpTo\nAWHgQL2b29Gjet0kA4NDalEq939/P1ELoxgRMoIjU45wT8Q9uGmN4a1InM6pW0yapr0KrEa/7uoN\n3IXewvSvzpyHEOdSkVBB0v+SKNxUSMgLIbS/vz1uzX9/EzMiyW1jYSFZ1dU8k5TEQC8vNkdE0MfL\nq8HPrbMNG/TrqtOmwahRsHcvdOli6BAFFQW89utrfLz3Yx4a8hBHHzmKX8vf/w4lya1xcvYZRDDw\n2fFfi4ADwPVKqfVOnocQtarOqSZ5ejLZS7Pp9FgnQueH0qz1n/+ZNLSwnsli4aPMTNq3aMGC3r2d\nl79wKqVg1SqYOxcGDNC3lAYPNnSICmsF7+54lze3vsnNfW7m4MSDdPDu8KfXSYBonJwaIJRS4505\nnhDnq6akhrS300ibnUa7u9sxLG4YLdoYfzNpbno6b6amYrFaKbTZeKB9e0yFhYCDD5xPN3MmvP22\nXj8pLw+uvBJWrIDiYkMS3WrsNSzat4joTdEM6ziMzeM30yeoT8PnLZxKLhiLi5LFZME/yh+71U7m\n/EySpyfj9xc/huwagmc3T8PHiykt5XmzmZ3FxTzftSv3BQfzcnKyc1t9AuzYofd+TkqC117TK61O\nn96gdp8nOrnB76Uxnv35Wdp4teGrW77ikk6XGDN34XQSIMRFqXBjIdYcK0nPJdGyR0sG/DAA70He\nho+TUFHB1KQk1lksPN2lC0vCwlxzM+ngQf0Aevdu/dfx4w07fD4RIDYnb+bp9U9TWl3KjGtnMLrn\naGna08RJgBAXFaUUhRsKyZyfiUdHD0I/CMX/auO3dtKrqphuNvN1bi6PdurE3NBQvE/rouaUJLf4\neL3V588/6z0ali2Dli3/+JoGbinllOUwZukYDmYfZPrI6dw54E7c3VwQBIXhJECIi0bqrFTSZqZh\nK7VRU1BD+wntKfylENwafuh8Qm51Na+lpLAwK4sH2rfnyPDhBJ7hJ3WHnjmkpupbR99+C489BvPm\ngfcZVkj1CBAms4mvY79mU/ImYnJiuLb7tYwbOI7Ovp0lOFxAJECIC17hr4WYp5qpTK6k20vdaHtH\nW5JfSjY00a2opoa3U1N5Lz2dO9q2JWboUNo7u14SQE6O3rRn0SI9we3oUQgIMHSI+IJ4Ptn7CWvi\n1/DY8Mco7VXKq9e8augYonGQ7BRxwSraVsT+a/dz+O7DtBun30wKvjsYt2YN/7/9iSS3cpuNN1JS\n6PXbb6RUVbF7yBDeCw11XnA4keRWWKiX3A4L0wvrxcbqgcLA4JBoSeS+7+/jko8uoWdAT+Ifiee5\nq57Do5nkuV6oZAUhLjjFO4sxTzVTdqiMkP+FEHxPsOHd3H62WIgtL+fl5GQu8/VlU0QEYa5Iclu7\nVu/D8PbbeuOePXsgJMTQIcyFZl7+5WW+PfwtU4ZOIf5f8ZLkdpGQACEuGCV7SjBPNVO6r1Rv87m8\nP24eta8W6nvmUG23sygri3fT07nM15eVAwYw+Ex7+45UXg4ffgizZ8MNN8Cvv0JvYzutpRSl8Mrm\nV/gq9ismRk7k2CPHCPD884pEAsSFSwKEaPJK95dijjZTvKOYLs90oe9XfXFvaexBaaXNxn8TE1mQ\nlUWb5s0pstkY5u3Nirw8imtqnJfkVloKTz4JixdD585QVgZ9+sDSpfphswFJbunF6byy+RWWHVrG\ng4Mf5MiUIwS1Cmrwc0XTIwFCNDknktxKY0pJnpZM0a9FdH66M2FLwnD3NDYwVNhsfJiZyZspKQzy\n9mZdeDjDfHyITkpybpJbURG8957e93nkSH1bacAAPcGtAUlu8HseQ2ZJJq/9+hqfHfiM+wfdT9zk\nONp6tTVk+qJpkgAhmpzcr3LJmJdBoamQzk92ps/CPrh7GRsYSmtqmJeRwVtpaVzq48MKV20lWSx6\nUHjvPRg9GjZt0g+iDbTq6CpWHFnBwn0LuTfiXmInxxLcOtjQMUTTJAFCNBmlB0pJeTWFvBV5dH2h\nK70/6l1rIb2GKK6p4f30dGalpRHl58dPAwcysHXrP73O4UlueXl6vaR58/RGPdu2Qa9amug0YEsp\nsySTt7a9xfs73ufBIQ9yaNIh2nu3r/+cxQVHrrmKRq9oaxF7Lt3Dnsv2UJ1Xjb3cjq3CRuqMVCwm\niyFjFFqtvGg20+O334gpK2NjRARf9OtXa3AABya5ZWfrXdxCQ/UgsWuX3gu6tuAA9QoQiZZE/r70\n73R/pzubkzdTaavE39OfD3Z/gMlsatD0xYVFVhCiUVJKYVlrIfnVZKqSq+j8n86Ejw/HvaU7SdFJ\nDU5yM1ksRPn7k2+1Mistjbnp6dwQGMiWQYMIbdXKoD/FuSZh+v0NPj0d3nwTPv0U7rwT9u/XD6EN\ndDD7IK9teY2f4n/i4ciH+fjvH9PGqw3Rpmiio6INHUtcGCRAiEZF2RV5y/NIfiUZe6WdLs90oe3t\nbQ1JbjvV6vx8fiwoYH5mJje3acOOIUPo7ml8FdezMpmge3d4/XX9FtK990JMDHT4c7+EhtiWuo1X\nfn2FXRm7ePySx5n7t7n4ePgYOoa4MEmAEI2C3Won+/NsUl9Pxd3Hna4vdCVwTCCa25+rgTYkyS21\nspKZx1cM49u3Z09kJCGnF69zhvh4WLlSz2OYMEHv6NbWuBtDSinWJqzl1V9fJbkomacue4ov//El\nns3/HAQlj0GciaaUcvUczkjTNNWY5ycazlZuI/PjTFJnpOLZ05OQZ0Pwu9rP8DLRB0tLeSIhgc2F\nhUS0bs32khKmHs84jvLzc14ew5w5+o0ksxkqKvTzhlatDMthsCs7y+OW88qvr1BZU8kzVzzDbf1u\no7m7cX2lReOnaRpKqYb/I1JKOe0DeAbYgd5uNAdYAfQ7y+uVuLAUbCxQSillLbQq8ytm9Wu7X9WB\nGw+oou1Fho9lt9vVhoICdf3+/ar9li3qVbNZWaqrlVJKTU1MNHy8M7LZlFq5Uqkrr1QqJESpWbOU\nKilRaurUBj96Y9JGpZRSVTVV6pM9n6je7/ZWw+YPU9/FfadsdluDny+apuPvnQ1+z3b2FtNVwHvA\nLkADpgPrNU0LU0oVOnkuwgXyV+djWWsh48MMAkcHEr4+nNb9a78pVF81djvf5OXxZkoKZXY7T3bu\nzHf9++Ph5uRLe1VV8PnnMGOG3oPhqafgllugmXH/7NYlrONA9gFmbJ1BaGAoc/42h5FdR0qjHmEI\nZ/ekHn3q7zVNuxt9NXE5sNqZcxHOVX6knLRZaWQuzKT9+PYM2Wl8a88ym40FmZm8nZZGRw8PXuja\nlRsCA3Gr5c3SoXkMhYXwwQd6gtuAAfo5w1/+AqfPowFbSrlluczbNY9Zv83iuh7X8c2t3zC049CG\nzVuI07j6kNoHPRfDmMvsolFRSlG4sZDE5xIpO1CGd6Q3qlLRvG1zshZl4RflZ0ijntzqat5LT2du\nRgZX+PryeVgYl/r6nvV7HHLmkJoKs2bBwoXw17/CmjUQHn6WSUTVeYjY3Fj+s+4/rE9cT982fSm3\nljOw3UBWH1tNmbVMDpyFoVwdIN4B9gDbXDwPYSB7lZ2cZTmkvp2Ksio6Pd6JduPa4e5pbA5DfHk5\nb6elsTQnh1vbtOFXV+UwHDig5zCsXq33et63z9AcBnX8RtLM7TPZl7WPSUMn8cmNn9DWq63kMAiH\nclmA0DTtbeAy4PLjhyq1ij6lEFlUVBRRBtz0EI5RnVdNxrwMMt7PwGuAF91f707AtQG1XlVtiM+z\ns3k/IwNTYSEPtW/P4WHDaNeihaFjnNPGjXpjnjff1APEv/4F774LBm5dVVgr+Pzg58zaPgs3zY3H\nL3mc727/jpbNXHAtVzRqJpMJ04nmUQZyyTVXTdNmArcCUUqpY2d53dlih2gkyuLKSJuVRu6XuQTd\nHESnxzrRekDtB88nKrHWldVuZ3leHu+mp3OwtJRp3bpxf3AwrQ088D0vlZXwxRfw3/+Cv79eevuu\nu8DADnJZpVnM2TmHD3Z/wNAOQ3n8kse5utvVtR48n6jEKsSpjLrm6vRaTJqmvQPcBow8W3AQjdOJ\n2kdKKQrWFXDgrwfYF7WPFu1bMOzwMPp83OeMwQHq3qgnp7qal5OT6bB1K08lJND+eC8Gi9XKjNTU\nk60/HS4tDcaNg8BAmD4dsrL0G0kpKXohvXo4ve7RgewDjP9+PGHvh5Fblsumezex6s5V/KX7X854\nK0mCg3Akp/74pWna+8A44EagSNO0dse/VKqUKnPmXET9WNZbqEyqJG1mGsqu6Px4Z/p928/wBj27\nS0qYnZbGivx8/i8oiPUREYQfL5zntF4MSumd2mbPhp9/1gPEnj165zaD+jBcFXIVa46tYeb2mcTl\nxektPR+JJ7BVoCF/BCEawtlnEBMBBfx82uenAS86eS6iDqoyq8j8MJO0WWn4XelHjxk98B/lb+h9\n+2q7nW9yc3k3PZ2MqiomdezI2z17EtjcyVnAFRWwZIl+plBZCVOmwMcfg49x9YvKqsvYlbGLvu/3\nxauFF49f8ji39ruVFu5OPksR4iycnQch5cWbEKUUhaZCzNPMFG8vxqufF/YyO97DvSnaWoTWQjPk\nmmpWVRUfZmYyLyODPq1a8VTnzowJDKTZGRLbHJbDkJysl8L45BMYPhzeeAOuuQZqm0c9L0ss3LeQ\nuTvnciDnAJU1ldwbfi9dfLvQyaeTBAfR6EgtJvEn1kIr2Z9mkzE3A9ygw8QOBN8dTDPfZg2+pnri\niirAjuJiZqelsbqggNvatGFKx470P0P/BcOduKaqlH4j6d134Zdf4J57YPJk6NHDsKGqbdUsj1vO\n3F1zOZJ/hAcGPcCDQx7k470fyxVV4RBGHVK7Og9CNCIle0rImJtB7te5+F/nT+gHofhe6WvoNtJ6\ni4W0qireTU8n12plcseOvNurF/7O3kZau1avoPree3qQeOQR+OwzMDBApRSl8OHuD/l478eEBYUx\neehkbupzkxTOE02GBIiLnK3CRu6XuaTPTac6o5oOD3VgaNxQPIJrv7ZZ31LbiRUVfJSZyey0NC7z\n9eV/ISH8NTAQd2fXDIqNhfnz9Vae11+vH0CPHPnnMhj1ZFd2for/ibm75rIldQvjBoxjwz83ENbm\nz32k5QaSaOwkQFxETs1BKI8vJ2NeBlkLs/AZ6kPIsyEE/i0Qzf3sb5R1OXOostv5Pi+P11NSiCsv\nJ9zLizK7nct8fNhdUoK3u7tzymyXl8OXX+pF85KTISJCP3wOD9e3ldzc6nymcHr+QW5ZLp/s/YQP\ndn9AgGcAEyMnsvT/luLVwuuMz5AAIRo7CRAXEcsGCzWFNWTMyaB0bynB44MZ8tsQPHsYWzTvSHk5\n8zMy+Cw7m/5eXjzVuTNj27TBw83NeVdUAfbu1VcLy5bBpZfCSy/B3/4GzZs3+JqqyWxiRMgItqZu\nZe6uuaw6uoqxYWP54h9fSNE8ccGQAHERqEyuJGthFumz0vHq70WHiR3ov6K/obkLFTYbX+fmMj8z\nk2MVFdwbHMyWQYPo6azaSCcUF+vtO+fPh5wcuP9+w/s7F1UWsTN9J+HzwqmyVfHwkIeZPXo2AZ4B\nho0hRGMgAeICZauwkbc8j5S3UiiPLcervxe2Ehv+1/pTkVBB8fbiel1RPfUWEsCB0lLmZ2ayJDub\nYT4+PNapE2MCA2nuyCuqpxbKA/2Q+bff9KDw7bf6mcL06XDtteB+hiBYxy0lu7Izc9tMFu5fyLH8\nY1TZqrh74N108+vGoPaDJDiIC5IEiAuIUoqSnSVkfpJJ7pe5eA/1pstTXQi6KQj3lgZVUi0sJNLb\nmy9yc5mfkUF6dTX3BQefd29nQ84cTgSIggJYvFgPDBUV8MADEBcHwcHnfsZ5BohESyIL9y1k0f5F\nBHoG8uDgB7lzwJ28u+NduaIqLngSIC4A1dnVZH2WRdaCLOyVdoLHBxO5L5KWXYyr+qmUYndJCSvz\n85mdns6Vvr4837Ur1wcEOPcmklJ6P+dx42DVKhg9Wm/MExVVe0JbPZRVl/F17Ncs2LeAQ7mHuGvA\nXay4fQXhwWfp7SDEBUgCRBNxehVUu9VO/up8shZkUbipkKCbggidczxv4QzltetzRTW1spLpZjPL\n8/KoUYpCm41/d+qEt7s7Xm5uzgsOn3+urxQOHACLRd8+evBBvTFPPbKaT7+FpJRiS+oWFuxdwLeH\nv+XyzpfzyLBHGNN7TK0ZznIDSVwMJEA0EYWmQvyj/CmNKSVrQRbZi7Px7OVJ+/vaE7Y4jGbe5/6f\n8nzPHIpravgmN5fPsrPZX1rKP9q0YXn//lzu68s0s9l5t5Dy8/XrqZ99BvHxcPvtev+FVatg2rQG\nPfpEgEgrTuPT/Z+ycN9C3N3cGR8xnthJsbT3bn/W75cAIS4GEiCaAKvFSsnOEnYP201VehXB9wQz\naPMgWoUad0Ooxm5nrcXCZ9nZrMnPJ8rPj8kdO/K3gABanumg1xGqqvTObJ99Bhs26Mlszz4L112n\nX08F/esNUFlTSUxODNcvvp4d6Tu4pe8tfDr2U4Z3HG5o1rgQTZ0EiEbKVmYj+fVkcr/IpdJciapW\ntL2zLQHXB+B3tV+9gsPpN5CUUuwtLeWz7GyWZmfTtWVL7g4O5t2ePQk6Q4c2hxTKUwq2bNEPnL/6\nCgYOhLvv1ns719Zbuh5bSja7jVnbZ7H4wGIO5x2m0lbJ2D5jmRg5kVE9RnFJp0sa/McQ4kIjxfoa\nEXu1nYK1BeQszSF/dT4+l/jQ7o52BN0UROrM1AbfQDqRpJZaWcnn2dl8lp1Nhd3OuHbtGNeunXP6\nOZ96RfXYMX2lsHgxtGypB4W77oIuXQwZSinFtrRtLD24lK9iv6KTTyfu6H8Ht/W/jY/2fCS3kMQF\nS4r1XSCUTVH4SyE5S3PI/TaXVn1a0e6OdvSc2ZMWbY0r/1xcU8Pe0lKu3rfv5LnCB6GhXO5rbDG+\nc/rhBzh0SA8MZrN+rvDVVzB4sCH1kJRSHMg+wNKYpSyLWYZnc0/u6H8Hv4z/hdDA0IbPX4iLiAQI\nJzj9BpJSipJdJeQszSHnixyat2lOuzvbEbk7kpYhtV9Nrc8NpOKaGl5PSeGb3FySKiupVopb2rRh\nYocOXOPvzxWO6qtwuoIC+P57PRBs2AA33wxTp8KoUVCPntK19WGOL4hnWcwylsYspbS6lNv73c53\nt39HeLvwWgOgHDILcW5ODxCapl0JPAkMAToA9yqlPnX2PJzpxA2ksrgyPSgszQGg7R1tCV8fjlfY\nmQu6nVCXG0gr8/P5KieHDYWFXOXryzMhIfw9MJB30tKcewPpu+/0oLB5s75t1LevfggdGqpnPnt6\nNuiKakZJBl/EfMHSmKUkFyVzS99b+PCGD7m086W4aWfPiZAAIcS5uWIF0Ro4CCwCLujAAFCRWEHR\nliJ2RuzEmmul7W1tCVsShnekt2FbO0U1NazMy+Or3Fw2FhYyws+PW9q0YWGfPvg5s89CXt7vQWH7\ndn2FMH48fP31730WGlgkL688j90Zuxm5aCT7s/ZzY58beenql7i629U0c5MFsRBGcvq/KKXUGmAN\ngKZpi5w9vqMppSg7UEbqzFQK1hRgK7VhL7fT7p52tOzSEr+r/fAZWvfexqffQCqqqWHF8aBgOiUo\nLDpLUHBIHaTcXFi+XA8KO3bo11Hvv1+vieR17pXR+UgpSuHNLW+y8uhKMkszqbZVc2vfW5k0dBLX\ndL9GVgNCOIj8yGUAZVcUbysmd3kuecvzwA5BY4Po900/fC/1xTzdbEgNpIjWrVmRn89XublsKiwk\n6nhQ+CwsDN/z2Ms3rA5S376/B4WdO/VchYce0lcP5woK57GlpJQiLi+O5XHLWX54OeZCMzeE3sA7\n17/DtT2u5fUtr8sNJCGcQALEeTj9kBn0K6mFGwv1oPBdHi3atCDoZj0otA5vbdj2UXZ1Navz81mS\nk8PbaWmM9PPjtjZtWHyeQcEwaWmwciUsWgSzZuk1kCZN0oNDXa7HniFA2JWdnek7WX5YDwrl1nJu\n6n0Tb4x6g6tCrpLtIyFcoNH/q4s+Zb86KiqKqHocajbUiUNmW5mNgh8LyP02l4I1BbTq04qgsUEM\n+nUQrXqe+U2yLjeQlFIcLCtjZX4+n2dlkVhZSXdPT45VVPDfzp3xcHOjk4eH44OD3Q579uhBYckS\nPUD06qVfTX32WT2rOSCgbsGBP95AstqsbErexPK45Xx/5Hu8PbwZ22csi8cuJrJD5BmDrGwpCfFH\nJpMJk8lk+HNdmiinaVoJMPlMt5gaQ6KctcDK4fsOg4LCjYX4XOJD0Ngggm4MwqND7X2b66rKbsdU\nWMjKvDxW5eejaRpjAgMZExjICD8/WjirE1t5OaxfrweF1avBxwfGjNE/LrtMv5LawEPm535+jsgO\nkSw/vJzVx1bTw78HY/uMZWzYWPoE9THsjyLExUwS5RxE2RWl+0pJfz+dgnUFWLOtqGpF4E2BdJjY\ngYDRAXVutHP6ATNA7vGto5X5+ay3WOjn5cWYwEBWDRhAPy8vxySvnX7ADJCerhe/W7lS788cGakH\nhP/8R18xGCDRksiaY2tYE7+GdYnruKLLFYztM5ZX/vIKnXw6GTKGEMJ4rsiD8AJ6AhrgBnTRNC0c\nKFBKpTp7PqAXw7OstZC/Jp+CHwto5tuMgNEB9JnfB98RvqS8ltKgQ+YTt4wOHd86Wpmfz6GyMq7x\n92dMYCBzQ0Npe4baRycYdgPpqqt+3zpauRKSk/VzhHHj9JIX5xrnPLb4Kmsq2WTexJp4PSjklOXQ\nxbcLvQJ6UW2r5souV5JXnkd8QbwECCEaMadvMWmaNgLYCJw+8CKl1H2nvbbBW0y1HTCfWCUUrCkg\nf00+ZQfK8L3Sl8C/BhIwOgDP7p5/eH19O7GV2WxsKixkmtlMrtVKjVInt46i/PycVyW1uBg2boQX\nX4TMzNq3js5TbVnMAAkFCScDwubkzQxsN5DRPUczutdoIoIjTiauRZui5QaSEA7WZLeYlFKb0FcO\nTuUONq8AAA8PSURBVHHigNlqsWJZZzkZFJr5NCPgrwF0fb4rvlf54u555jfr8z1ktinFnpIS1lks\nfJmTQ2xZGe09PEipquKh9u1p16IFI/38jLluejZWq56TsG6dnqR29Ch06gRJSTBlCgQG6iuBq66q\n86NPBIgKawWbkjed3DoqqS7h+p7Xc0/4PSweuxh/Twf/GYUQDnfBVnNVNn2VkPBUAvZqO2X79VVC\nwF8DCBwdiGcPz3M/5DwkVVSwzmJhncXCBouFdi1aMMrfn1H+/ozw88O7WTPHHzArpQeBdev0j02b\noFs3PZN51Ci44gq9rEUDDpiVUsQXxPOvNf9C0zQ2p2wmIjhCXyX0HE14cPg5y1vAmVcgQgjjNNkV\nhKMou6IspozCjYXkfJVDya4S3Fu7U5NfQ9u72uI3wg//Uf51PmCGPx4yF1qtbCwsZK3FwrqCAkps\nNq7x9+dvAQHM7NGDTi2N6wP9x0mY/pzBvH69HhDWr9c/N2oU3HGH3pqzbdsGD5lcmMxG80aWxixl\ne+p27MpOqbWUf4T9g0mRkxjda3Sd3+wlOAjRdDTZAKGUovxIOYUbCincWEihqRB3X3f8R/rTcXJH\n/KL88GjvUe/zgxOsdjuLsrLYUFjIOouFmLIyLvPxYZS/Pw/3788ALy/cznHjyJAD5nXr9K2jE6uE\npCQYMUIPCk8/rRfAO9fNp3McMKcXp7PRvJGNSRvZaN5ImbWMqK5RjO0zltnXzyY0MJRpm6bJGYIQ\nF4kmEyCUUlQmVmLZaDkZFLQWGn4j/QgcE0iPt3vQsnPDf3qvsNn4rbiYzUVFbC4q4rfiYlq5uXFP\ncDAvdevG5T4+dT5crteZQ2kpbNumV0Ld/P/t3X9sldd9x/H3Fwhg8PWP8cs2TmwcwgBDhtmgmNaA\nK4rULKnUtVC1ZCkTlVYSKSNdOylqV4jWqH8kStOKRvsjWrZKndRI3bRIa7r1B2QOmOIE15iYyBDX\nFBvb2IDxNeBf+OyPc22Mc23fe33t+wCfl3Tky3PP89xzD8fP189zznNOhV9x7Z13fEA4dAg2bry9\nBGeMjhTCthH/butuuyMgXLl5ha2FWykvLOcbpd9g9aLVWoJT5D4W+ADR+m+tw0HBDTiyyrPI+nQW\ny763jLnL5k54Apuog7mzv5+jXV1UdHbyf9eucaq7mzXz51Mwdy4LZs3ib/PyeOnCBebOmEFFZycz\nSdKcRqO1t8O7794OCGfO+H6ErCw/VXZ/P2zf7vsb+vvjDg4Avzj7CzpudAwHhJbuFrYUbKG8sJyn\nNzzN2iVrNU22iAwLfCf16Z2nySrPIvvT2aStSIv7L9rRD6ld7O31VwednVRcu0ZDTw8bQyG2ZGVR\nlpnJJzIymD/qCmHSncyj+w+c888fDAWDigo//LS0FMrKfNqwwS/DOVyIg3F1MDvnqL9cT2VTJccu\nHKOyqZL6y/VsL9pOeWE55YXlrMtZx8wZ0zTUVkSmzX3TST1v9Tz62vrobell3p/GN++Pc46ft7fz\nUU/PcEDoHBjgU5mZlGVm8tc5OaxPT+eBGVM86vbwYVi48M6A0N9/Oxjs2wePPgrj3Lpq7GykcJyP\nCPeGOdF8gsqmSiqbKjnedJzQ7BBF2UXMmTmH0vxSTl86zYa8DXT3dRPuCys4iMi4Ah8g4ulgbu/r\noyocpioc5kRXF1XhMD2Dg1weGKAsM5NvPfQQq+bNm7BTebS4O5lbW/002FVV/nmEd96Bn/7UB4Pt\n2+GFF2D58rjWYD5SCHsir4eGnFY2VVJ5wQeEs1fOUpJTQml+KXtL9vL6E6+TG8q94xh5oTx1MItI\nzAIfIMYSHhjgZHf3cCCoCoe52t/PX4RCLJ49m/SZM9m9eDGvNjezIi2Ntr4+2vv6KE5gEZttNTVj\njwDq7IT33rsdEKqq4Pp1ePhh/+xBXh709PipLACKiuKe46i7r5sjhdBS8X2ONR3jeNNx0malUfpg\nKaX5pexZt4eS3BJmzxx/ug4RkXjcFQGid3CQU93dd1wdNPb08Gh6OhtCIT63YAH/tGwZj6Slfezq\nIHPWrMk/pDbUh3DjBvz+9/6qYCgYtLRASYnvM9i1C156yQeBkeVYuTLm/oOOGx1Ut1RT3VrNL8/9\nkpq2Grp6uxgYHGDT0k3kZ+Tz2mOvsbN4Z9xfQx3MIhKPwAeIje+/zwfXr7M8LY0NoRCbMjJ4dulS\n1syfP7V9B93dcPo01NTAW2/51dLq6/1qahs3+ltFzz8Pq1aN23cA0fsPnHM0dTVR3VpNdUs1J1tP\nUt1SzbXea6zLWcf6nPX+yiCnhJULV/JixYuTvj2kACEi8Qh8gFgzbx7bs7PZkZ2d0PDSbefP++Gi\nYxkYgHPn4NQpqK316dQpuHjRL4izeLEPEl/7Gjz+uA8McS5adLjA8cnL9T4QtJz0QaG1mhk2g/W5\n6ynJKeHJtU/yyo5XWJa9LKYpK0REplrgh7lOunxDw0Od853HI4NAbS18+KHvJ1i71o8kWrvWp+XL\nh68MGvfvofDVf53wowbdII2djZxpP0Ndex11HXWcaT9DdUs1uaHc4WCwPnc9Jbkl5KbnxjxsV3MY\niUis7pthrnFzzvcLnD3rbwm9/bZfCKe21r83FAjKyvyaysXFkJ4+7iFH3yLqv9XPR1c/oq7dB4C6\njjrq2uuov1zPgrQFrF60mvTZ6fQO9LJ60Wp+1/w79qzbA0BZQVlCJ3oFBxGZbndvgLhy5XYQqK+/\n/frsWf+Xf0aGv0VUUwO7d/uO5Mcfh/LymD8i3Bum4WoDhwscvz38Xc50+CuDhqsN5Gfks2rhKlYv\nWs2Ooh3s/8R+Vi5cSWhO6GPHeSjzIQ0vFZG7TvADRE3Nx4NAfT309fkJ6las8MNGn3ji9usRzy2M\nd3vo1uAtmsPNNFxt+Fj6sONDrvddJzstm/YZ7WxphEXzF/HcpufYvXY3aQ8kZ7pwEZGgCn6A2L37\n9om/rAz27vWvlyyJ6UGzc1fOca215s4A0Ol/nu88z4J5C3g4+2GKsosoyi7isUceG369ZP4SzGzS\nq6Dp9pCI3I1SEiDM7Gngm0Au8AGw3zn3btTMX/yi/7lt2x2jhwbdIJe622juauZi+CLN4eY7Xg/9\nXDerk0v/+aQ/6WcVsXLhyuEgUJhVOC1XAgoQInI3mvYAYWZfAl4Fvg4cBZ4B3jazVc65ptH5f/L5\nIn+y7/o5zT/70fCJv627jay5WeSF8liasZSloaXkhfLYlL+Jyzcv09jZSGh2iJdvvsyBVV8A/Ila\nHcQiIjFyzk1rAo4D/zxqWz3wYpS8bs1ra9zm1ze7Z/77Gffm6Tfd0T8edX+4+gfX09/jYvHVNw7E\nlG88hw+ndv976RhBKEMyjhGEMgTlGEEoQ1COEYQyOOecP7VP/nw9rU9kmdkDwJ8Dvxr11v8Cm6Pt\nU7uvlqN7j3LosUPsLN7J5gc3U5hVyJxZc2L6zMbGyZTYO3IktfvfS8cIQhmScYwglCEoxwhCGYJy\njCCUIZmm+5HdhcBMoG3U9jYgZyo+sPCONdTub42NR1JdhHuG6jK5VJ/BNK1PUptZLtAMbHEjOqXN\n7B+BrzjnVo3K7+DAiC3bQCf8STgYSTJ5B1FdJtNBVJ+TcSSShryQlCepp/sKogO4BSwZtX0J0Bpt\nB+cOjkjbcI640oED8eWfimMEoQzOwdatqS9HUOpisscIQl0G5Rj3StsMSl0kdoxtjDxXJsu0Bgjn\nXD/wPvCZUW99Bj+iSUREAmLaJ+szs13AT/DDW48C+4C/AYqdcxdG5Z3ewomI3COScYtp2p+DcM69\naWZ/Anwb/6DcaeCzo4NDJO+kv6CIiCQm0NN9i4hI6mhlGhERiUoBQkREokppgDCzp82swcxumtl7\nZvapCfKvMbMjZnbDzC5Enp+QiHjq08wKzGxwVLplZjums8xBZGZlZvZfZtYUqZenYthHbXMM8dan\n2ubYzOx5MzthZtfM7JKZvWVmxTHsl1D7TFmAGDFp3/eAdcAx/KR9+WPkD+Gn6GjBT9fxd8C3zOy5\n6SlxsMVbnxEO2IF/ij0HP2jgt1Nc1LtBOlALPAvcmCiz2uaE4qrPCLXN6LYAh4BSoBwYAH5tZllj\n7TCp9pmMCZ0SScQxaV/kvX1AJzB7xLZvAxdS9R2ClBKozwJgEFif6rIHOQFh4KkJ8qhtJrc+1TZj\nr8/5kSDxl+PkSbh9puQKIpFJ+4BNQIVzrm/Etv8B8sysIPmlvHskWJ9D/sPM2szsXTP7wpQU8N6n\ntjk11DYnloG/E3R1nDwJt89U3WJKZNK+nDHy2zj73C8Sqc9u4O+BXcBngd8APzOzr0xVIe9hapvJ\npbYZux8CJ4HKcfIk3D6Dv+SoTAnn3GXgByM2nTSzBcA/AP+emlKJqG3Gysxewd8h+KSL3DdKtlRd\nQcQ9aV9ke7T8bpx97heJ1Gc0J4BHklWo+4ja5tRT2xzBzH4AfAkod86dnyB7wu0zJQHCJTZpXyVQ\nZmazR2zbAVyMoYLuaQnWZzQl+JEOEh+1zamnthlhZj/kdnA4G8MuibfPFPa+7wJ6gL3ASvy9tC4g\nP/L+94Ffj8ifAVzEX2IWA38FXAP2p3okQRBSAvX5FPDlSN4VwDcj+z+b6u+S6oQfGfJn+OHC14Hv\nRP794Bh1qbaZ3PpU2xy7Ln8caVvb8FcBQ2n+iDxJa5+p/rJfBxqAm0AV/l7a0HtvAB+Nyl+MXxXj\nBn7hoe+k+j8sSCme+oz8En6AH3bYib+E/3Kqv0MQErAVP8zy1qj0L9HqMrJNbTNJ9am2OW5dRqvH\nW8B3R+RJWvvUZH0iIhKV5mISEZGoFCBERCQqBQgREYlKAUJERKJSgBARkagUIEREJCoFCBERiUoB\nQkREolKAEImTmR02sx+luhwiU00BQkREotJUGyJxMLM3gK/ip0q2yM9lzrk/prRgIlNAAUIkDmaW\nAbwNnAGexweJdqdfJLkHaUU5kTg457rMrA+44ZxrT3V5RKaS+iBERCQqBQgREYlKAUIkfn3AzFQX\nQmSqKUCIxK8R2GhmBWa2wMws1QUSmQoKECLxexl/FVEHXAIeTG1xRKaGhrmKiEhUuoIQEZGoFCBE\nRCQqBQgREYlKAUJERKJSgBARkagUIEREJCoFCBERiUoBQkREolKAEBGRqP4fwHlH9mZ/sRsAAAAA\nSUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def question1a(h = 0.1, x0 = 0.2):\n", " n = int(2/h)\n", " ts = np.zeros(n+1)\n", " xs = np.zeros(n+1)\n", " xs[0] = x0\n", " for k in range(n):\n", " ts[k+1] = ts[k] + h\n", " xs[k+1] = xs[k] + 2*h*xs[k]**0.5\n", " return ts, xs\n", "\n", "# 1(a): various initial conditions\n", "x0s = [0.0, 0.1, 0.2, 0.4, 0.6]\n", "for x0 in x0s:\n", " ts, xs = question1a(x0=x0)\n", " plt.plot(ts, xs, '-+');\n", "plt.xlabel('t'); plt.ylabel('x(t)'); plt.ylim(-0.1,8); plt.xlim(0,2);\n", "# The x0=0 solution is an unstable equilibrium." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Q1(b). Euler's method for an ODE $\\frac{dx}{dt} = f(x,t)$ yields a difference equation $x_{k+1} = x_{k} + h \\, f(x_k, t_k)$, where $h$ is the grid spacing.\n", "\n", "Therefore, in this case $f(x,t) = 2 \\sqrt{x}$. Thus, by separation of variables, $\\int dt = \\frac{1}{2} \\int x^{-1/2} dx$ $\\Rightarrow$ $t + c = x^{1/2}$ $\\Rightarrow$ $x(t) = (t + x_0)^2$. " ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEZCAYAAAAzL+qdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYFNW19/HvEgUUgxCMgBJ8QSXxEklexRww4JiAiUeP\n+MoREQ0mGBMjRIziIZ4ogibBY0wU7wgBQVG8BBFFHEVogzAIAhpRIxq5eUREkIvgiMJ6/9g1oWl6\nZnpmuqere36f5+lnpqt2VVcV7Sz33mvvbe6OiIhInOyT7wsQERFJpeAkIiKxo+AkIiKxo+AkIiKx\no+AkIiKxo+AkIiKxo+AkIiKxU+/Bycy6m9mTZva+me0yswEZHNPXzJaa2TYzW2FmQ9OUGWRmb5rZ\ndjN7y8x+nKZMHzN7w8zKzWyZmZ2drfsSEZHsyUfN6UDgdeByYHt1hc3sdGAycC9wLHAZ8Gszuyyp\nzC+BUcAI4Jjo511mdkZSma7AFOABoDPwEPCYmXXJxk2JiEj2WD5niDCzrcAgd59URZnJQFN375O0\nbTBwtbsfHr2fByxw96uSytwCnOTuPaL3U4CW7v7DpDLPAx+5+wVZvjUREamDQuhzagKUp2wrB9qZ\nWftqypxkZo2i912B51LKlALdsnitIiKSBYUQnEqB3mbWy4JOwJXRvrZJZQaa2YkA0c+Lgf2Ag6My\nbYB1KedeF20XEZEYiX1wcvexwB3ANGAHMB94ONq9K/p5IzADmGdmXwBPAPenlBERkQKxb74vIBPu\nfo2Z/TehlrMe6Bntei/aXw78zMx+AbQG1gK/ALa6+/qo7IfRvmSto+17MTNN1y4iUgvubnU9R+xr\nThU8WOvuXwL9gTJ335BSZqe7f+Ahy6Mf8FTS7jKgV8ppexFqYpV9pl5Zel1//fV5v4Zieul56lnG\n9ZUt9V5zMrNmwJGAEYJjezPrDGx09zVmNgro4u49o/KtgHOBBCHxYSDQB+iRdM6jgO8CC4CvEvqk\njgWSx1CNBl40s2GEJsJzgBLg5Fzdq4iI1E4+ak4nAkuBxUBTYCSwJPoJoemuQ8oxA4CFwEvA0cAp\n7r44aX8jQkB6lZAc0Rjo5u6rKwq4exmhNnUR8BpwIdDX3V/J5s2JiEjd1XvNyd1fpIqg6O4/TXm/\ngWrSvd39H8D/zeCzpwJTM7tSyaaSkpJ8X0JR0fPMHj3LeMrrINw4MzPXsxERqRkzwxtSQoSIiDQc\nCk4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4i\nIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIhI7Ck4iIlJ7\niUROTlvvwcnMupvZk2b2vpntMrMBGRzT18yWmtk2M1thZkPTlOmfVGatmT1gZq2T9l8Ufd7O6GfF\n742zfY8iIkUhXeBJ3VbxfssWePnlrH10PmpOBwKvA5cD26srbGanA5OBe4FjgcuAX5vZZUllTgYm\nAROAY4DewNHAgymn2wa0SXq1dfcddbwfEZHiUFngSTZzJsybB2PHwhVXwAMPwNe+BgcfDGefnbVL\n2TdrZ8qQu88EZgKY2cQMDrkQmO7uY6L3K81sFDAMuDva9m/AGne/PXq/yszuBG7f81S4u6+v0w2I\niBSqRAJKSjJ7v2kTrF4N990Hb7yx+7VhA0yfDk2bhqD03nswZAgcdBCcemp4ZUG9B6daaAKUp2wr\nB9qZWXt3Xw3MA35vZme6+9NmdjDQD5iRctz+ZrYSaAS8Clzn7q/m9vJFRPIkk2D03e/CW2/B66/D\nc8/BM8/AP/4B5eXwxRewZAnsuy+0aAH9+sFtt8F554XjS0qgWzcYMSLrl14IwakUuNXMegGzgKOA\nK6N9bYHV7r7AzM4HJpvZ/oT7eg74SdJ53gYGAq8BXwGuAOaZ2fHu/s96uRMRkVyqKhjt2hVqPVOn\nhkD097+H/X/4A7RsCYccAsuWhQDUtSv07g1/+9vegeegg/bclqOEiNgHJ3cfa2YdgWlAY2AzMBoY\nAewCMLNjgDuAkYSg1Ba4BbgPuCg6zwJgQcV5zawMWAr8ihCo9jIi6R+gpKSEkuR/dBGRfEoNRMnb\nNm0KwWfhQjjjDHj1VVi/PtSEnngC9tkHOnSAjRvh2muhUaNwXCKxZ+D529+qv4wWLUjkoOaEu+ft\nBWwFBmRY1ghBZ1/gR8BOoFW0bxLw15TyJxOC16FVnHM8MKOSfS4ikjdz5lT9/vrr3XfudH/3Xfe/\n/tX9uuvcO3Vyb93afb/93A87zB3c//3f3X/6U/enngrHpJ6jqvepn1nZtiTR3846x4fY15wqRDe9\nFkLaOFDm7hui3QcQglWyXYBTdUZiZ0LtSUQkv6rrH5o1KzSpLV0aXk8+CX/+c0hMaNEC2rSB5cth\n8GD46ldDYkJqTeiVV6q+htSaWLrWonpqQar34GRmzYAjCTWhfYD2ZtYZ2Ojua6JMvC7u3jMq3wo4\nF0gQkiMGAn2AHkmnfQq4z8wuJfRRHQrcCix29/ej8wwnNOu9AzQHhgDHAT/P6Q2LiKRTVTD69NOQ\nKXfnnSEhYenS0B80eTJ85SshEK1ZA1dfDQccEI4rKQmBqKr+oOqCT4y6LvJRczoRmEOo1UDoJxoJ\nTCQEnjZAh5RjBgA3EwJaGXCKuy+u2OnuE83sQGAQoa9pEzAb+E3SOVoAY6LzbybUmLonn0dEJCeq\n6h/aujUEn7Iy6NUr/L55M3z5Jbz0Uih71FGh3+iii8L7TLPkCigYpbLQWiapzMz1bESkVlKDUUWN\nZtu2EHwWLYIJE0Liwtq10Lo1/O//wllnQdu28J//GQJTcvBJrRWlvk/3uXlgZri71fU8BdPnJCIS\nW+ma6Lp1250xN20a/PWvoU/o4IPh0ENDOvell4aBrN///t79QxW1psrksT+oPig4iYjUVGowmjMn\nBJyXXw7BaOpUuOmmkMBw6KHw2mtwySVhep8f/KDo+odyQcFJRCRZVf1DFZ55BrZvD8FowYIwHmjM\nGGjVCg47DD74AK65Bho3Tj9+KJ0GHoxSKTiJSMNWXQr3l1/Cww/Dm2+GQFRWFjLpXnghZMq1axem\n+hk2LJRPl6yQbhaFBh58qqPgJCINS3XB6LPPYMaMEITmzw/JC02bhtTtRo1CRt0998B//EcoX1IC\n3/hGzZroKtsm/6LgJCLFrapg5B7mm/vNb8Kkp6tXh/ePPRbGEx1yCFx2Gdx8M5x0UjimpCRsr0t/\nkVRLwUlEiktVwejLL0PK9qBBMHt2CEbbt8Oxx0L79vDLX4Ya0g037HnO/fdXMKpnWqZdRApbVQvk\nffYZrFwJAwdCx47QrBmMGxcy6o45BsaPh+uvD7MvPPNMyKjbJ4M/iwpGOaeak4gUttRpf/75T7jw\nwpBB9+GHYSburl3DukXnnx+mA0odzJosk/4hBaOcU3ASkcKR2mRXEYwuuCAEo3XrQjDq3h1OOQX6\n9w/p3skBaMmSPc8Zo8lOZTcFJxGJr9Rg9PzzsHMn3H9/2FcRjHr0CLMs9O8fsuySg9HLL+95TtWC\nCoKCk4jER2owmj0bmjQJY4peeAHmzQuzMZx6aghQ3brB//zPnsGorGzPcyoYFSQFJxHJn3TTAB1y\nSFi7aNaskN790ENh/rkOHUIt6bTTQtlGjUIWXSoFo6Kg4CQi9SddmvfRR4dA9PzzYXLUMWPCfHQd\nO8Lnn4fkBgjHdepUfQKDglFRUHASkdyobI66rl1D81xpKUycCLfcEqYAOuKIkOBw5ZVgFo497ria\nTwMkRUHBSUSyo6rBr//8Jzz7bGiiu/nmsIz4EUeEhIbrrgtji0pKoEuXmk8DJEVJwUlEsiM5GJWX\nw7vvwpAhoalu0yY48kh45509lxZPna1bMy9IRMFJRGontaa0eXNoknv6aVi1CnbsCOndZ5wB550X\nfq/rGkbSYCg4iUhm0qV5v/lmWE5i+XL46CP41rfC65Zb9p6JIR0FI6mEgpOIpJeuD6lLl5BV9+ST\n8OijIaPuzDPDWKNnn91zwtTUmRhAwUgypolfRSS9iia3detg7Njd442uvDLMWbd9O/TpEwbJ7tix\n94SpmgZI6qDea05m1h0YCpwAHAr8xN0nVXNMX+AaoBPwEXCXu9+SUqY/cHVUZgswCxjq7uuSyvQB\nbgCOAN4FrnX3aVm6NZHCllxTevfdsNDet74VkhgqkhmGDQsL75WUhIlUq2q2UyCSOshHs96BwOvA\nRKDKoARgZqcDk4HBQClwNDDOzLa7+91RmZOjc10JPAm0Bu4GHgR6RWW6AlOA64AngD7AY2bWzd0X\nZfMGRWIvtcnOHR55JEwJNHMmbNsWXhdcAL17Q8+eNc+sE6mDem/Wc/eZ7n6tu08FPINDLgSmu/sY\nd1/p7jOBUcCwpDL/Bqxx99vdfZW7LwTuBL6bVGYIMNvdb3L3t939D0ACuCIb9yUSa+nWPHKHxYvh\nmmvCzAsPPxzGH02dGjLvrr8eHnwQfvc7NdFJvSuEPqcmQHnKtnKgnZm1j97PA9qa2ZkAZnYw0A+Y\nkXRMV+C5lPOUAt2yfsUicVMRnNxh0aKQ1HDEEXDWWTB3bphIdfNmaN487Js7d+9zKJlB6lEhZOuV\nAreaWS9CP9JRhOY7gLbAandfYGbnA5PNbH/CfT0H/CTpPG2AdexpXbRdpLgkN9u5w9q1YTmJmTND\n4sLGjfCLX0Dr1iEwlZSE+ezUhyQxEfvg5O5jzawjMA1oDGwGRgMjgF0AZnYMcAcwkhCU2gK3APcB\nF9X2s0ck/YdaUlJCif7jlLhKl/a9cmVYhnzZMvjkE/je98Jg2PPOgxdfrPkYJJE0EokEiXRzHtaR\nuWfS7ZMbZrYVGFRdtl5U1gi1nPVAT0KT3SHuvsHMJgHN3L1PUvmTgblAO3f/wMxWAbe7+5+SygyN\nPr9Dms/zfD4bkRqpmHlh1arQd3TrrdC4cQhE/frBU0/ByJF7l0+WbqJWkRoyM9zd6nqeQuhzAsCD\nte7+JdAfKHP3DdHuA4CdKYfsIiRcVNxjGVHmXpJewPwcXbJI7iT/n+qGDfDKKyHt++ij4fHHw2wN\nAwfCgQeGmb4t5W+FEhwk5vIxzqkZcCRghMDR3sw6AxvdfY2ZjQK6uHvPqHwr4FxCZl0TYCAhDbxH\n0mmfAu4zs0sJfVSHArcCi939/ajMaOBFMxtGaCI8BygBTs7d3YpkQboazaxZoWnu4YdhxYowCLZf\nPzj7bPjBD/ZO+06lQCQxl4+a04nAUmAx0JTQT7Qk+gmh6S61mW0AsBB4iTDO6RR3X1yx090nEpIk\nBhHGUD0C/AM4O6lMGSGD7yLgNUKKel93fyW7tyeSZcmZdmVlcOml8Kc/he3DhsH69SHt++GH4cYb\nVSuSolDvNSd3f5EqgqK7/zTl/QYySPd297uAu6opMxWYmtmViuRJak1pyxYYNQruvjtMGfTtb4cl\nKU45JfQxZTKHnUiBiX22nkiDk0iEZrqxY2Hp0rBQ3wknwA9/GGZsOPXU9AkNyRScpMApOInkW3JN\n6bXXwlikO+8MNaQbboA33oDf/77qcygYSZEpmGw9kaKROiaktBSGDoXDDoPu3WHhwlBD+t73wsDY\n/fbb+xwKRlLkVHMSqW8VNaUlS2DMGJg4EX70I7jvvvDzxhs1QFYaPAUnkVxLbrbbvj30I33zm2FK\noRNOgM8/D014ixZBs2Z7H69AJA2QgpNIriUSsGYN3HUX/P3v8NlncP75YY2k738fevSovqYk0sAo\nOIlkU3ItaedOePppeOCBMEvDxReHNZMmTNC6SCLVUEKESDYlEmGS1V/+Elq1gkGD4L334JJLwlx3\nK1bsfYyCkcheVHMSqYvkmtKbb4YJVkePhjPPhOeeg5NO0pgkkVpQcBKpizlzQur3fffBhx+Gpc2v\nuipMuLp9e/pjFIxEqqXgJFITFTWl8nKYPDlMKdS2Lfz2tyHJ4aab9q4lKRiJ1JiCk0hNzJwZEhwe\neQTatIGPP4bLLgtz3C1YkP4YBSeRGlNwEqlKRU1p9eowE/h994WlKRYsgOOOq74/SURqRdl6IlV5\n9NEw4eo3vhEGyZaXw+GHhwX90i1NrVqSSFao5iSSrKKmtHBhmGz1hRfgv/8bpkyBli1VUxKpJ6o5\niSSbOBFOPBF69YIvvgjZdzt2hPRw1ZRE6o1qTtKwJRJh0b7Zs2HkSFi2DP74R/jxj8OgWdWURPJC\nNSdpuNzhL3+B44+H886Dgw8OszusWQN/+INqSiJ5pJqTNBwV/UnuoaY0fDi8/TbcfnsITo0aqaYk\nEhOqOUnDkUjA3Lnwne9A375hnNKGDbB8eVhDSTUlkdhQzUmKV/K8d6+8EgbPTpwI118PF14I++6r\nmpJITNV7zcnMupvZk2b2vpntMrMBGRzT18yWmtk2M1thZkNT9k+IzrUz5efWpDIXpeyr+L1xLu5T\nYiCRgEmT4Jhj4NRTw+zgF14IK1fCSy+lP0Y1JZFYyEfN6UDgdWAiMKm6wmZ2OjAZGAyUAkcD48xs\nu7vfHRW7HBiWcuh8IJGybRvQEbCKDe6+o+a3ILGUXFN6/32YPj0kN1x9NQweDDffrHnvRApEvQcn\nd58JzAQws4kZHHIhMN3dx0TvV5rZKEIwujs651YguZZ0MiEIXbD3x/v6ut2BxFYiERb1u+kmWLw4\nzObwX/8VZgdfuDD9MQpOIrFUCH1OTYDylG3lQDsza+/uq9MccwmwzN1fTtm+v5mtBBoBrwLXufur\n2b5gyYMdO8J8d/fcE9ZSmjIFxo1Tf5JIgSqEbL1SoLeZ9bKgE3BltK9tamEzaw6cC9yXsuttYCBw\nFtCPEODmmdkRObtyya1EIqyndN55IfOutBTOOQe+/nV49929y6uWJFIwYl9zcvexZtYRmAY0BjYD\no4ERwK40h/yY0Kf0YMp5FgD/WtPAzMqApcCvgCvSffaIpP/rLikpoUR/3OLloYfCOKVPPglLWMyb\np5qSSD1LJBIk0g3DqCNz96yfNOMPD9l0g9w9k8QIA9oA64GewAzgEHffkFJuKfC6u2eSBTgeaO3u\nZ6TZ5/l8NpJGRcLDhx/CNdfAY4/Bn/8MF1+sAbQiMWFmuLtVX7JqhdCsB4RMBndf6+5fAv2BsjSB\n6SSgMzA2w9N2BtZm90olZ154AS69FDp2DANnt22DDz7YPYBWNVuRolHvzXpm1gw4ktD0tg/Q3sw6\nAxvdfU2UidfF3XtG5VsR+pAShOSIgUAfoEea0/8cWO7uc9N87nBCs947QHNgCHBcdIzE3cyZYUn0\nrl3h1VehUyfVlESKWI1qTmbWxMwGmtktZvZHM/uJmTWp4WeeSOjrWQw0BUYCS6KfEJruOqQcMwBY\nCLxEGOd0irsvTrm2A4G+VF5ragGMAd4kJFm0BbqnnkdioqINe8oU+OY34YILYOPGsJzFQw+ln2pI\nRIpGxn1OZnYMYXzSQYRBtADfIiQo/Mjd38rJFeaJ+pzy7Npr4YADwtLoV1wRBtLedNOeNSU15YnE\nTrb6nGrSrDeaMDbox+6+JbqI5oSsuNuAH9b1YqQBSw40c+bAvfeGJrxXXoEOqRXpiAKTSNGqSXA6\nmdAXtKVig7tvMbPfkpSiLVIriQRs3hwmZV2xArZsgRNOCBO1lpTsfolIg1CT4FRO6LdJdRB7z+Ag\nkjn3kORwzz2hb2nkyNCcp3nwRBqsmgSnp4CxZnYJu2tKXQlJBtOzfWFS5CoSGp54IkzQunIlXHIJ\nNG8e5sUTkQatJsFpCGEm8bnAzmjbPoTAlHaGBZFKzZ4NLVvC5MkwbBhs3Qo33JDvqxKRmMg4OLn7\nJsIcd0cB34w2v+XuaSYxE6nCm2/C+PFw5JFQVgZHHaUmPBHZQ16nL4ozpZJnWSIBO3fCqFEwfz58\n9hkMHw5muwORApJIwauXVHIzux24xt23Rb9Xyt0vr+vFSBF7/PGQFt6sWag53X+/ZncQkUpV16z3\nLWC/pN9FambXLrjtttCMd8stYW68fQpmSkcRyRM161VCzXp1lEjAunXwm9+E5rw1a8IYJlAznkgR\ny1azXk2mLxoO3OLu21O27w9c7e5FlWql4FQH7tCnD7z0Evz612Gp9BtvVDOeSAOQjyUzrgcOTLP9\ngGifCGzaBP37w9y58OyzYd2lRo3yfVUiUmBqMs7JgHRVie8AG7NzOVKwEglYtixM2NqpE3z8cRhc\nO326ph4SkRqrNjhFq9V69HrPzJIDVCPCshf35ubypCDs3BmmHHrrLZg0Cc46S2stiUidZFJzGkyo\nNY0HfktYIqPCDmClu5fl4NqkEHzwQWjGW7UqTDt02GH5viIRKQLVBid3nwhgZiuA+e7+Rc6vSuKt\nYnmLP/4x1Ji6dAlz442N1nlUM56I1FGtUsnNrA3QOHmbu6/O1kXFgbL1qjB8eBi/dP/98OCDIRCp\nGU9EyMNig9HCgncQlkJvnKaIUrIagrVrQ7/SUUeFZrzWrfN9RSJShGqSSv4noDNwNmH9pv7A1cD7\nwHnZvzSJjUQivAYODJl4q1aFVWrvuWf30hdqxhORLKpJKvnpwPnuPtfMdgKL3f0RM1sL/AJ4PCdX\nKPk3Zw589aswYwY8+ii8/LJmEReRnKpJcGoBrIp+3wy0At4FyoBxWb4uiYtPP4W//hUaN4YFC6BD\nhxCcRERyqCbNev8EOka/vwX0MzMDzqEGg3DNrLuZPWlm75vZLjMbkMExfc1sqZltM7MVZjY0Zf+E\n6Fw7U35uTSnXx8zeMLNyM1tmZmdnet0NTiIBQ4bAEUfAG2/A6afDxIm7M/VERHKoJjWn+4HjgQRw\nE/A0YQzUPoRVcjN1IPA6YVXdSdUVNrPTgcnRZ5UCRwPjzGy7u98dFbscGJZy6PzoWivO0xWYAlwH\nPAH0AR4zs27uvqgG11/8EokwP94jj8BvfwsbNoSUcRGRepJRKrmZ7Qe8BAxw97ejbe2BE4F33P31\nWn14qNkMcvdKg5SZTQaaunufpG2DCZPNHl7JMScTlpPv6u4vR9umAC3d/YdJ5Z4HPnL3C9Kco2Gm\nkrvDGWfAkiVhCfUf/EBp4iKSsXpNJXf3L8ysA0lz60XjmupjbFMTQnZgsnKgnZm1r2R81SXAsorA\nFOkKpC6YWAoMytqVFrovvoDBg8OigAsWQMeoFVfNeCJSz2rSrDeR8Ef/6hxdS2VKgVvNrBcwCzgK\nuDLa15aUABmNxzqXvZv52gDrUrati7Y3bIkEbN0a+pj23RfWrw9jmUCzPYhIXtQkODUDLoiCxGJg\nW/LOXC3T7u5jzawjMI0w+HczMBoYAexKc8iPCXMBPljXzx6R1JRVUlJCSbH+kZ46FZ57Ds4+O0xJ\npLWXRCRDiUSCRMV4xyyqyWKDc6rY7e7+/Rp/eAZ9TklljVDLWQ/0BGYAh7j7hpRyS4HX3X1AyvZV\nwO3u/qekbUOjz++Q5vMaRp/T3LkhE69iCXVQH5OI1Fq9T1/k7qfW9cPqIooUawHMrD9QliYwnUSY\nxSJdLa4M6EWY6aJCL0JWX8OTSMCdd8LMmbB9O3z4YQhIasYTkRioSbNeVphZM+BIQtPbPkB7M+sM\nbHT3NWY2Cuji7j2j8q0IfUgJQnLEQEIaeI80p/85sNzd56bZNxp40cyGEZoIzwFKgJOzd3cFYs6c\nkI338stQVhaa9VRTEpEYqckg3Gw5EVhK6LdqCowElkQ/ITTdpTazDQAWEtLZjwZOcffFyQXM7EDC\npLRj031otOZUP+Ai4DXgQqCvu79S91sqIDt3wm9+AxMmwPz5cPzx+b4iEZG91GrJjIagKPucysvh\nggtCrWnpUmjRImzXrA8ikiX13uckBW7GDBg0CJo1CwsD3nZb2K4+JhGJIdWcKlE0NadEAo45JmTk\nffe7cMcdShUXkZzJVs0pH31OUp+mTYPvfQ/OPBPuugsaaU1IEYk/NesVs3/8IyQ+/O538Ktf7d6u\nZjwRiTnVnIpRIhEG1HbpAlu2hFnFR4zQqrUiUjBUcyo2iURIenjiCRg/PqzFpP4lESkwqjkVm0mT\nwpIX48bBuefm+2pERGpFNadiMmdOWCBw2jTo1StsUxOeiBQgBadikEiExIfHHoPPPoN588JLY5hE\npEApOBW6iiSHmTPhmWfCe/UxiUiBU59Tobv/fujbFx59VLUkESkaqjkVsr/9LTTlPf307sCkACUi\nRUDBqRAlEjB5Mjz0UFiL6cUXw0t9TCJSJBScCk0iAS1bwvTpoSlv0SL1MYlI0VGfU6F5/PEwieud\nd4bxTCIiRUg1p0KyYgU88ADcfvvuAbZqxhORIqQlMyoRqyUzEomQ9DB+PHzyCVx/fdiuPiYRiZls\nLZmh4FSJWAWnLVtCEOrdG9zVxyQisaX1nBqK556Ds8+Grl1h+PB8X42ISL1QcIqznTvhiiugVavQ\nz2SmZjwRaRDUrFeJvDfrucNll8Gzz4ZFA5s0yd+1iIhkqGCb9cysu5k9aWbvm9kuMxuQwTF9zWyp\nmW0zsxVmNjRNmf3M7AYze8/Mys1spZkNTtp/UfR5O6OfFb83zvY91lkiEWYVf+IJWLkSRo3ac7FA\nEZEil49U8gOB14GJwKTqCpvZ6cBkYDBQChwNjDOz7e5+d1LRR4BDgZ8B7wKtgf1TTrcN6Aj8K6q7\n+45a30mufPghLF8OixfD2LFKgBCRBqfeg5O7zwRmApjZxAwOuRCY7u5jovcrzWwUMAy4OzrPacCp\nwBHuvjEqtzr9x/v6ulx/TiUS0KgRXH45zJoFhx2W7ysSEcmLQkiIaAKUp2wrB9qZWfvofW9gEXCV\nma0xs+VmNtrMmqUct3/U3LfGzJ4ys2/n+Npr5vHHw+DaBx+E448P25QAISINUCEEp1Kgt5n1sqAT\ncGW0r230syPQHTgeOAcYBPwImJB0nreBgcBZQD9CgJtnZkfk/hYy8PHHYSLX3/8eTjtt93YFJxFp\ngGI/fZG7jzWzjsA0oDGwGRgNjAB2RcX2iX4/390/BYiSIZ41s6+5+3p3XwAsqDivmZUBS4FfAVek\n++wRSX20adyGAAAM/0lEQVQ9JSUllOQiUCQS8MILYVqiTz6BNWtCH5NmfxCRApBIJEjkIFkrr6nk\nZrYVGOTumSRGGNAGWA/0BGYAh7j7BjO7H+jm7p2Syrcj9Dt1cffFlZxzPNDa3feaQbXeUsnd4Re/\nCEkQ3/kOjByZ+88UEcmRgk0lry0P1rr7l0B/oMzdN0S75wGHmtkBSYd8A3BgVRWn7QyszckFZyKR\ngDvugLKysD6T1fnfU0SkKNR7s16UpHAkIZ17H6C9mXUGNrr7migTr4u794zKtwLOBRKE5IiBQB+g\nR9JpHwKuBSaY2UigJXAb8Ji7fxydZzihWe8doDkwBDgO+HlOb7gqf/lLyMqbPx++8hU144mIRPJR\nczqR0NezGGgKjASWRD8hNN11SDlmALAQeIkwzumU5KY6d99GaOo7KCo3BZgDXJx0jhbAGOBNQpJF\nW6B7ZU1+Obd8OUydGhYM7BDdroKTiAig6YsqlbM+p0QCSkth3LiQoaflL0SkiGjJjBzLWXByh759\noUWLMMhWsz+ISBFpcAkRRSGRgD/+EVatCokQIiKSVuzHORWV8ePh+edh4UJo2lTNeCIilVCzXiWy\n3qy3ahUce2xYbl1BSUSKVLaa9VRzyrWKGSDGj4dt28L7REIJECIiVVDNqRJZrTn9+tfw3nvw7W9r\nBggRKWpKiCgUTzwB06bB/fdrBggRkQypWS+XHn4YhgwJ/UwtW6oZT0QkQ2rWq0Sdm/U+/zzM/DBs\nWAhQIiINgAbh5lidg9MVV4SZIN58U815ItJgKFsvrhIJGDs2NOVt2bI7AULZeSIiGVPNqRK1rjmt\nWxey8h55BGbP1vREItKgKFsvjmbPhp/8BH72M+jRo9riIiKSnmpOlahVzelHP4JNm2DuXNhvv92D\nbUVEGgglRORYjYPTa69B166wbBl07Ji7CxMRiTElRMRFIhFWsx0zBj77DCZNCtuVACEiUmuqOVWi\nRjWnq6+GFSvCxK6ankhEGjAlRMTF3LkweTLcc4/GM4mIZImCU118+mnIzrv3Xvja19SMJyKSJWrW\nq0RGzXpnnQWtWsGECfVzUSIiMaeEiHwrLYUXX4TVq/N9JSIiRafem/XMrLuZPWlm75vZLjMbkMEx\nfc1sqZltM7MVZjY0TZn9zOwGM3vPzMrNbKWZDU4p08fM3oj2LzOzs2t1E5s2hYG2vXvDQQfV6hQi\nIlK5fNScDgReByYCk6orbGanA5OBwUApcDQwzsy2u/vdSUUfAQ4Ffga8C7QG9k86T1dgCnAd8ATQ\nB3jMzLq5+6KMrz6RgKuugtat4YEHdo9pUuq4iEjW5LXPycy2AoPcvdIgZWaTgabu3idp22Dganc/\nPHp/GiE4HeHuGys5zxSgpbv/MGnb88BH7n5BmvLp+5wqpihatgz+/GfNnScikqQhpZI3AcpTtpUD\n7cysffS+N7AIuMrM1pjZcjMbbWbNko7pCjyXcp5SoFvGV7JtG1xySUgbb968RjchIiKZK4SEiFLg\nVjPrBcwCjgKujPa1BVYDHYHuwOfAOUAL4E5CM9+5Udk2wLqUc6+LtlcvkYCnngpTFJ1xRtimZjwR\nkZyIfXBy97Fm1hGYBjQGNgOjgRHArqjYPtHv57v7p/Cvpr9nzexr7r6+Np89IqnJrmTRIkoWLw7N\nef/aWFKb04qIFI1EIkEikcj6eWPf55RU1gi1nPVAT2AGcIi7bzCz+4Fu7t4pqXw7Qq2qi7svNrNV\nwO3u/qekMkOjz++Q5vN29zl9/jm0awd33AH9+tX+hkVEilyDG+cURYq1AGbWHyhz9w3R7nnAf5rZ\nAe6+Pdr2DcCBVdH7MqAX8KfdZ6UXML/SD00kwmvbNvj4Y3jrrZAAocw8EZGcqvfgFCUpHAkYoTmu\nvZl1Bja6+xozG0Wo7fSMyrci9BslCMkRAwlp4Mmr+T0EXAtMMLORQEvgNuAxd/84KjMaeNHMhhGa\nCM8BSoCTK73Y5CDUrJky80RE6kk+svVOBJYCi4GmwEhgSfQTQtNdajPbAGAh8BJhnNMp7r64Yqe7\nbyM09R0UlZsCzAEuTipTBvQDLgJeAy4E+rr7K9m9PRERqSvNrVeJvcY5aVVbEZFqaSXcHKvVMu0i\nIg1cQxqEKyIiDYyCk4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iI\nxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6C\nk4iIxE69Bycz625mT5rZ+2a2y8wGZHBMXzNbambbzGyFmQ1N2X9KdK7k104z65RU5qKk7cllGufi\nPkVEpPb2zcNnHgi8DkwEJlVX2MxOByYDg4FS4GhgnJltd/e7k4o6cAzwSdK29Smn2wZ0BOxfB7nv\nqMU9iIhIDtV7zcndZ7r7te4+lRBQqnMhMN3dx7j7SnefCYwChqUpu97dP0p6pZ7f3X2PMnW8HclQ\nIpHI9yUUFT3P7NGzjKdC6HNqApSnbCsH2plZ+6RtBrxiZh+Y2SwzK0lzrv3NbKWZrTGzp8zs2zm6\nZkmhPwDZpeeZPXqW8VQIwakU6G1mvSzoBFwZ7Wsb/VwLXAr0Af4f8DbwgpmdnHSet4GBwFlAP0KA\nm2dmR9TDPYiISA3ko8+pRtx9rJl1BKYBjYHNwGhgBLArKrMcWJ502Mtm9n+Aq4F5UZkFwIKKAmZW\nBiwFfgVckePbEBGRGrC9u2Xq8cPNtgKD3D2TxAgD2hCSHHoCM4BD3H1DJeWHA+e5+7FVnHM80Nrd\nz0izL38PRkSkgLm7VV+qarGvOVWIkhvWAphZf6CsssAU+U5F+Sp0JtSe0n1enR+uiIjUTr0HJzNr\nBhxJSGDYB2hvZp2Bje6+xsxGAV3cvWdUvhVwLpAgJEcMJPQt9Ug65xBgJfAGoenvx4S+pXOSygwn\nNOu9AzQHhgDHAT/P3d2KiEht5KPmdCIwh91p5COj10RC4GkDdEg5ZgBwMyGglQGnuPvipP2No/3t\ngM8IQerf3b00qUwLYEx0/s2EGlP3lPOIiEgM5LXPSUREJJ1CSCXPCTO7zMzeM7PPzOwVM/teNeWP\nM7OEmW2PxkldV1/XWghq8jzN7PBKpps6rT6vOY5qOb2XvpuVqOnz1HezcmZ2jZktNLPNZvaRmU03\ns0oTzpKOq9X3s0EGJzM7D7gN+B3wbWA+MNPM2lVS/ivA84QEixMI/VVXm9mv6+eK462mzzPiwGmE\nZtY2hDFrs3N8qYWgYnqvy4Ht1RXWd7NaNXqeEX030+sB3Al0BU4FvgRmmVmLyg6o0/fT3Rvci5AY\ncW/KtuXA7ysp/0tgE9A4adtvgTX5vpc4vGrxPA8njFH7v/m+9ji/gK3AgGrK6LuZ3eep72bmz7NZ\nFKDOqKJMrb+fDa7mZGb7ESL48ym7ngO6VXLYvwFzfc9JYkuBQ83s8OxfZeGo5fOsMNXM1pnZS2bW\nJycXWPz03cwNfTer15zQ+vZJFWVq/f1scMEJOBhoBKxL2b6OUIVPp00l5SsGBjdktXmenwJXAX2B\n04EXgEei8WtSM/puZpe+m5kbDSwhZFBXptbfz4IZhCvFw8Pg6VuTNi2JxrP9F/BQfq5KRN/NTJnZ\nnwktIyd71FaXbQ2x5vQxsBNonbK9NfBhJcd8WEl5r+KYhqI2zzOdhcBR2bqoBkTfzdzTdzOJmd0K\nnAec6u6rqile6+9ngwtO7v4FsBjolbKrF9EksWmUAd1TVs09Dfggg3+colbL55lOJtNNyd703cw9\nfTcjZjaa3YHpnQwOqf33M98ZH3nKMulLWDLjYuCbhLbTLUC7aP8oYFZS+ebAB4Rq/bGEaZE2A1fk\n+17i8KrF8xwAnB+V7QQMjY6/PN/3ku8XIQOqMyElfxtwbfT+65U8S303s/s89d2s/FneFX23Sgi1\nn4pXs6QyWft+5v2G8/igLwXeI0x3tIjQdlqxbwLwz5TyxxLm99sO/C9wbb7vIU6vmjzP6A/AG4TU\n3k2EZpPz830PcXgBpxBSmXemvMane5bRNn03s/Q89d2s8lmme447geFJZbL2/dT0RSIiEjsNrs9J\nRETiT8FJRERiR8FJRERiR8FJRERiR8FJRERiR8FJRERiR8FJRERiR8FJRERiR8FJpICY2Rwzuz3f\n1yGSawpOIiISO5q+SKRAmNkE4CLCcgMW/ezg7qvzemEiOaDgJFIgzKw5MBN4C7iGEKDWu/4jliKk\nlXBFCoS7bzGzHcB2d1+f7+sRySX1OYmISOwoOImISOwoOIkUlh1Ao3xfhEiuKTiJFJaVwElmdriZ\ntTIzy/cFieSCgpNIYbmFUHt6E/gI+Hp+L0ckN5RKLiIisaOak4iIxI6Ck4iIxI6Ck4iIxI6Ck4iI\nxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxI6Ck4iIxM7/B4p3oPZj9kYOAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Question 1(c): The ratio test.\n", "h0 = 0.025\n", "# Try changing h0. You should see that that the ratio gets closer to 2.\n", "t0s, x0s = question1a(h=h0)\n", "t1s, x1s = question1a(h=h0/2)\n", "t2s, x2s = question1a(h=h0/4)\n", "ratio = (x0s[1::1] - x1s[2::2])/(x1s[2::2] - x2s[4::4])\n", "plt.plot(t0s[1:], ratio, 'r+-');\n", "plt.xlabel(\"t\"); plt.ylabel(\"ratio\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Q1(d). The ratio takes the value of '2' because the ratio tends to $r \\sim 2^s$, where $s$ is the order of accuracy of the method. \n", "\n", "Euler's method is first-order accurate, therefore $s = 1$ and so $r \\sim 2^1 = 2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "