{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Foundations of Computational Economics #13\n",
"\n",
"by Fedor Iskhakov, ANU\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Two very important algorithms for solving equations\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"\n",
"[https://youtu.be/gUdYS5PmvWo](https://youtu.be/gUdYS5PmvWo)\n",
"\n",
"Description: Bisections and Newton-Raphson methods. Solving equations of one variable. Accuracy of solution. Rates of convergence."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Classic algorithm for equation solving\n",
"\n",
"1. Bisection method \n",
"1. Newton-Raphson method \n",
"\n",
"\n",
"Solve equations of the form $ f(x) = 0 $\n",
"\n",
"Focus on the scalar case today."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"#### Bisection method for solving equations\n",
"\n",
"Solve equation $ f(x)=0 $, conditional on $ x \\in [a,b] \\subset \\mathbb{R} $ such that $ f(a)f(b)<0 $\n",
"\n",
"Algorithm: similar to binary search, but in **continuous space**."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"def bisection(f,a=0,b=1,tol=1e-6,maxiter=100,callback=None):\n",
" '''Bisection method for solving equation f(x)=0\n",
" on the interval [a,b], with given tolerance and number of iterations.\n",
" Callback function is invoked at each iteration if given.\n",
" '''\n",
" pass"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solution is x=-1.000, f(x)=0.000000834465\n"
]
}
],
"source": [
"f = lambda x: -4*x**3+5*x+1\n",
"a,b = -3,-.5 # upper and lower limits\n",
"x = bisection(f,a,b)\n",
"print('Solution is x=%1.3f, f(x)=%1.12f' % (x,f(x)))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Converged in 22 steps\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAHSCAYAAADlm6P3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5RV1cGG8WfP0MWoFMUSC0rEFg1SNGpUUIK9t9gbYhdjwW4sgC1ixYJiwwJK7FjQGFtEsWBFxA42LCgaEXHO98cePlBBZpg7s++95/mtddZt586845kLL8d99g5ZliFJkiTlQUXqAJIkSVJDsfxKkiQpNyy/kiRJyg3LryRJknLD8itJkqTcsPxKkiQpNxo15Ddr06ZNtvzyyzfkt5QkSVLOPP/8859nWdZ2bq81aPldfvnlGTt2bEN+S0mSJOVMCOH9eb3msAdJkiTlhuVXkiRJuWH5lSRJUm5YfiVJkpQbll9JkiTlhuVXkiRJuWH5lSRJUm5YfiVJkpQbll9JkiTlhuVXkiRJuWH5lSRJUm5YfiVJkpQbll9JkiTlhuVXkiRJuWH5lSRJUm5YfiVJkpQb5V9+swzeeSd1CkmSJBWBRqkD1LcHTjsNxo2j15Ah0LZt6jiSVCsPPPAAAL169UqcROWokL9fNflas/YZM2YMAN26dfv/1+Z8X7H+3hdrrmJVrP+9yr78ftKuHUyaBKecAldckTqOJNXKJ598kjqCylghf79q8rV+uc+83lOsv/fFmqtYFet/r/If9tCiBSy1FFx9NYwblzqNJEmSEir/8guw/PKw2GJw1FFxDLAkSZJyKR/lt1EjOOsseOwxGDkydRpJkiQlko/yC3DggfDHP8Ixx8D06anTSJIkKYH8lN/KShg0CN57D/75z9RpJEmSlEB+yi/AxhvD9ttD//4weXLqNJIkSWpg+Sq/AOedBzNnwgknpE4iSZKkBpa/8tu+Pfz973DjjVA9ybYkSZLyIX/lF+JZ3yWXhCOPhKqq1GkkSZLUQPJZflu2hIED45nfYcNSp5EkSVIDyWf5BdhjD+jaFfr1g2+/TZ1GkiRJDSC/5beiAi66CD76KJ4FliRJUtnLb/kFWGedeAb4/PPh3XdTp5EkSVI9y3f5hXjWt7ISjj02dRJJkiTVM8vv0kvDiSfCHXfAY4+lTiNJkqR6ZPkFOPpoWH55OOKIuACGJEmSypLlF6B5c7jgAnjlFbjiitRpJEmSVE8sv7Nstx1suimccgp89lnqNJIkSaoHlt9ZQoCLL45z/p54Yuo0kiRJqgeW3zl17Ah9+8I118Czz6ZOI0mSpAKz/P7SKafAkkvCYYdBVVXqNJIkSSogy+8vLbwwnHcePPccDB2aOo0kSZIKyPI7N3/7G6y/PvTrB199lTqNJEmSCsTyOzchwKWXwpdfwqmnpk4jSZKkArH8zsuaa8Ihh8Dll8O4canTSJIkqQAsv7/ljDOgVat48VuWpU4jSZKkOrL8/pbFFoMBA+DJJ+Hmm1OnkSRJUh1Zfudnv/2gSxc49liYNi11GkmSJNWB5Xd+KirixW8ffwxnnpk6jSRJkurA8lsTXbvC/vvDhRfC+PGp00iSJGkBWX5rqn9/WGghOPxwL36TJEkqUZbfmlp8cTjrLBg9Gm6/PXUaSZIkLQDLb20cfDD86U9w1FFe/CZJklSCLL+1UVkJgwfHi99OPz11GkmSJNWS5be2unWD3r3hoovg5ZdTp5EkSVItWH4XRP/+cQGMgw+GqqrUaSRJklRDlt8F0aoVnHcePP00XHdd6jSSJEmqIcvvgtprL1h/fTjuOPjii9RpJEmSVAOW3wVVUQGXXw5Tp8IJJ6ROI0mSpBqw/NbFGmvEac+uvhqeeSZ1GkmSJM2H5beuTjsNll46Xvw2c2bqNJIkSfoNlt+6WnhhGDQIXnopDoOQJElS0bL8FsIOO8Bf/wonnxwXwJAkSVJRsvwWQghw6aUwYwb8/e+p00iSJGkeLL+FstJK0K8f3HILPPJI6jSSJEmaC8tvIR1/PKy4Ihx6KPzwQ+o0kiRJ+gXLbyE1bx6HP7z5JpxzTuo0kiRJ+gXLb6H16gW77AJnnw0TJqROI0mSpDlYfuvDoEHxLHCfPpBlqdNIkiSpmuW3PrRrF4c9/PvfcMMNqdNIkiSpmuW3vhx4IKy7bpz67PPPU6eRJEkSlt/6U1EBV10FX38Nxx6bOo0kSZKw/Nav1VePxfe66+Cxx1KnkSRJyr0ald8QQt8QwmshhFdDCLeEEJqFEFYIIYwJIbwVQrgthNCkvsOWpJNPhvbt4aCDYPr01GkkSZJybb7lN4SwNHAE0DnLstWBSmBX4BzgwizLOgBfAfvXZ9CS1aIFDB4cpz0bODB1GkmSpFyr6bCHRkDzEEIjoAXwMdAduL369euBbQsfr0z07Al/+xsMGADjx6dOI0mSlFvzLb9Zlk0Gzgc+IJber4HngalZls2s3m0SsHR9hSwL//xnPAvs3L+SJEnJ1GTYw2LANsAKwFLAQsBmc9l1ro0uhNA7hDA2hDB2ypQpdcla2pZYAs49F/7zn3gBnCRJkhpcTYY9bAK8m2XZlCzLfgRGAn8GFq0eBgGwDPDR3N6cZdlVWZZ1zrKsc9u2bQsSumTtvz+svz4ccwzk+R8CkiRJidSk/H4ArBNCaBFCCEAP4HXg38CO1fvsDdxVPxHLSEUFXHklTJsWC7AkSZIaVE3G/I4hXtj2AvBK9XuuAo4Hjg4hTARaA9fUY87yseqqcNxxcdnj0aNTp5EkScqVRvPfBbIsOw047RdPvwN0LXiiPDj5ZBgxAnr3hldegYUWSp1IkiQpF1zhLYVmzeDqq+Hdd+HUU1OnkSRJyg3Lbyp/+Utc9W3QIHjuudRpJEmScsHym9I550C7dnEWiBkzUqeRJEkqe5bflBZZJC59/MorcQ5gSZIk1SvLb2pbbw077wxnnglvvJE6jSRJUlmz/BaDiy+OMz4ceCBUVaVOI0mSVLYsv8VgiSXgwgvhqafgiitSp5EkSSpblt9isddesOmmcPzx8OGHqdNIkiSVJctvsQghLn1cVQUHHwxZljqRJElS2bH8FpMVVoCzzoL77oNbb02dRpIkqexYfovNEUdA167x9vPPU6eRJEkqK5bfYlNZCUOGwNSp0Ldv6jSSJEllxfJbjNZYA/r1g5tuglGjUqeRJEkqG5bfYnXyybDqqtC7N3z9deo0kiRJZcHyW6yaNoWhQ+Gjj+CYY1KnkSRJKguW32LWtWssvkOGwEMPpU4jSZJU8iy/xe4f/4COHePSx998kzqNJElSSbP8FrtmzeLwh0mT4LjjUqeRJEkqaZbfUrDOOnHasyuvhEceSZ1GkiSpZFl+S8WZZ8If/gAHHADffps6jSRJUkmy/JaK5s3h2mvh/ffjHMCSJEmqNctvKVlvPTjySLjsMnjssdRpJEmSSo7lt9ScfTasuCLsvz98913qNJIkSSXF8ltqWrSIsz+8+y6ceGLqNJIkSSXF8luKNtgADjsMLr4YnngidRpJkqSSYfktVQMGQPv2sN9+8L//pU4jSZJUEiy/pWqhheCaa2DiRIc/SJIk1ZDlt5RttBEcfjhcdJGzP0iSJNWA5bfUDRwIHTrAPvvAtGmp00iSJBU1y2+pa9ECrr8ePvwQjj46dRpJkqSiZvktB+uuC8cdB0OGwP33p04jSZJUtCy/5eL002GNNeCAA+DLL1OnkSRJKkqW33LRtCnccANMmQKHHpo6jSRJUlGy/JaTtdaC006DW2+F4cNTp5EkSSo6lt9y068fdOkChxwCn3ySOo0kSVJRsfyWm0aN4vCH776DAw+ELEudSJIkqWhYfstRx47Qvz/cey9cd13qNJIkSUXD8luujjwSNtww3r7/fuo0kiRJRcHyW64qKmDo0DjsYb/9oKoqdSJJkqTkLL/lbIUV4J//hEcfhcsuS51GkiQpOctvuTvgANhss7gC3BtvpE4jSZKUlOW33IUA114LLVvC7rvDjBmpE0mSJCVj+c2Ddu3g6qvhxRfjIhiSJEk5ZfnNi223hf33h3POgSeeSJ1GkiQpCctvngwaBO3bw557wtdfp04jSZLU4Cy/edKyJdx0E0yaBEcckTqNJElSg7P85s0668BJJ8UlkEeMSJ1GkiSpQVl+8+jkk6FrVzjoIJg8OXUaSZKkBmP5zaPGjePwhx9+gH32cfU3SZKUG5bfvOrQAS68EEaPhosvTp1GkiSpQVh+8+zAA2GrraBfP3j11dRpJEmS6p3lN89CgCFDYJFF4upvP/yQOpEkSVK9svzm3eKLwzXXwMsvxwvhJEmSypjlV7DllnDwwXD++fDww6nTSJIk1RvLr6Lzz4dVV4W99oIpU1KnkSRJqheWX0UtWsCtt8JXX8G++0KWpU4kSZJUcJZfzbbGGvEM8H33wSWXpE4jSZJUcJZf/dyhh8bpz449FsaNS51GkiSpoCy/+rkQ4NproXVr2HVX+N//UieSJEkqGMuvfq1NG7jxRnjzTejbN3UaSZKkgrH8au569IDjjoOrroKRI1OnkSRJKgjLr+btzDOhSxc44AD48MPUaSRJkurM8qt5a9wYbrkFfvwR9tgDfvopdSJJkqQ6sfzqt624Ilx+OTz+OPTvnzqNJElSnVh+NX977gm77w7/+Ac8/XTqNJIkSQvM8quaufxyWG452G23uAqcJElSCbL8qmZ+97u4/PHHH8N++7n8sSRJKkmWX9Vcly5wzjlw550ufyxJkkqS5Ve1c9RRsPXWcMwxMHZs6jSSJEm1YvlV7YQAQ4dCu3awyy7w9depE0mSJNWY5Ve116pVHP/7/vtxAQzH/0qSpBJh+dWC+fOf4eyz4fbb4YorUqeRJEmqEcuvFtyxx8Jmm0HfvvDSS6nTSJIkzZflVwuuogKuvx5at4add4Zp01InkiRJ+k01Kr8hhEVDCLeHEMaHEN4IIawbQmgVQng4hPBW9e1i9R1WRahtW7jlFnj7bejTx/G/kiSpqNX0zO9FwANZlnUE1gTeAPoBj2RZ1gF4pPqx8ugvf4lLH998M1xzTeo0kiRJ8zTf8htC+B3wF+AagCzLZmRZNhXYBri+erfrgW3rK6RKwAknQI8ecPjh8OqrqdNIkiTNVU3O/LYHpgBDQwgvhhCGhBAWApbIsuxjgOrbxesxp4pdZSXcdBMssgjstBN8+23qRJIkSb9Sk/LbCOgEDM6y7E/Ad9RiiEMIoXcIYWwIYeyUKVMWMKZKQrt2MGwYTJgAvXs7/leSJBWdmpTfScCkLMvGVD++nViGPw0hLAlQffvZ3N6cZdlVWZZ1zrKsc9u2bQuRWcWsRw8444x4EdzgwanTSJIk/cx8y2+WZZ8AH4YQVq5+qgfwOnA3sHf1c3sDd9VLQpWeE06AzTeHo46CZ59NnUaSJOn/NarhfocDw0IITYB3gH2JxXl4CGF/4ANgp/qJqJJTUQE33gidOsXxvy+8EOcCliRJSqxG5TfLspeAznN5qUdh46hstGoFI0bA+uvDnnvCvffGUixJkpSQbUT1p0sXGDQIRo2C/v1Tp5EkSbL8qp716QO77w6nngqjR6dOI0mScs7yq/oVAlx5JayyCuy2G0yalDqRJEnKMcuv6t9CC8Edd8D06bDzzvDjj6kTSZKknLL8qmF07AhDhsB//wvHHZc6jSRJyinLrxrOLrvA4YfHi+Buvz11GkmSlEOWXzWs88+HddaB/faD8eNTp5EkSTlj+VXDatIEhg+HZs1gu+1g2rTUiSRJUo5YftXwfv97uO02eOst2GcfyLLUiSRJUk5YfpXGxhvDuefCyJEwcGDqNJIkKScsv0qnb1/YdVc46SR46KHUaSRJUg5YfpVOCHH6szXWiCX43XdTJ5IkSWXO8qu0FlooDn3IsngB3P/+lzqRJEkqY5ZfpbfiinDzzfDyy9C7txfASZKkemP5VXHYbDM44wwYNgwuuSR1GkmSVKYsvyoeJ54I22wDRx8Njz+eOo0kSSpDll8Vj4oKuOGGOAxip51g8uTUiSRJUpmx/Kq4/O53cOed8cK3HXaAH35InUiSJJURy6+KzyqrwPXXw5gxcMghXgAnSZIKxvKr4rT99nDKKXDttV4AJ0mSCsbyq+J1+umw7bbxArhHHkmdRpIklQHLr4rXrAvgOnaMF8C9/XbqRJIkqcRZflXcFl4Y7r47LoW89dbwzTepE0mSpBJm+VXxa98eRoyAN9+EPfaAqqrUiSRJUomy/Ko0dO8OF14I99wDp56aOo0kSSpRjVIHkGrssMNg3Dg4+2xYYw3YZZfUiSRJUonxzK9KRwhw2WWw3nqw777wwgupE0mSpBJj+VVpadoU7rgD2rSJ06B9+mnqRJIkqYRYflV6llgiLoH8+edxCeQZM1InkiRJJcLyq9LUqRMMHQpPPQV9+rgEsiRJqhEveFPp2mUXeO01OPNMWGUVOPbY1IkkSVKRs/yqtJ1+epz/9/jjoUOHOA5YkiRpHhz2oNJWUQHXXQddusDuu8OLL6ZOJEmSipjlV6WveXO46y5o3Rq22go++ih1IkmSVKQsvyoP7drF1d+mToWtt4b//S91IkmSVIQsvyofa64Jt9wSF7/Yay+oqkqdSJIkFRnLr8rLVlvB+efHhTBOOSV1GkmSVGSc7UHlp29fGD8e+veHlVeOZ4ElSZLwzK/KUQhw2WXQvTsccAA88UTqRJIkqUhYflWeGjeG22+HFVaA7baDt99OnUiSJBUBy6/K12KLwb33xqWPt9gCvvwydSJJkpSY5VflrUMHuPNOePfduPrb9OmpE0mSpIQsvyp/G2wA118fx/7uu69ToEmSlGPO9qB82HVXeP996NcPll8eBgxInUiSJCVg+VV+HHdcHP4wcGC8EK5379SJJElSA7P8Kj9CgEsvhQ8/hEMOgWWWgc03T51KkiQ1IMf8Kl8aNYLbboM//hF23jkuhSxJknLD8qv8adkyToHWqlWcAu2DD1InkiRJDcTyq3xaaikYNQq+/z4OfZg6NXUiSZLUACy/yq/VVoORI2HCBNhhB5gxI3UiSZJUzyy/yrfu3WHIEHj0UTjggLganCRJKlvO9iDttVcc93vKKbDkknDOOakTSZKkemL5lQBOOgk+/hjOPTcW4KOOSp1IkiTVA8uvBHEO4Isvhk8+gb59oV27uCqcJEkqK475lWaprIRhw+Avf4lDIUaPTp1IkiQVmOVXmlOzZnDXXdCxI2y3nYtgSJJUZiy/0i8tumicA7hVK9hsM3j77dSJJElSgVh+pblZeml48EGYORN69YLPPkudSJIkFYDlV5qXjh3hvvtg8uS4Cty0aakTSZKkOrL8Sr9lnXVg+HB46SVXgZMkqQxYfqX52XJLuPpqePhh2HdfqKpKnUiSJC0g5/mVamLffeMcwCeeCIstBpdcEucGliRJJcXyK9VUv37wxRdwwQVxJogzzkidSJIk1ZLlV6qpEOC882DqVDjzzHgGuG/f1KkkSVItWH6l2ggBrrwSvv4ajj4aFlkE9tsvdSpJklRDll+ptior4aab4Jtv4MAD46IY22+fOpUkSaoBZ3uQFkTTpjByJHTrBrvtBqNHp04kSZJqwPIrLaiFFoqLYKy8Mmy7LTzzTOpEkiRpPiy/Ul0sthg89BC0axdXgXvlldSJJEnSb7D8SnXVrl0c9tC8OfTsCW+/nTqRJEmaB8uvVAjLLx9XgPvxR9hkE5g0KXUiSZI0F5ZfqVBWXRUeeAC+/BJ69IgrwkmSpKJi+ZUKqXNnuP9+mDw5ngH+/PPUiSRJ0hwsv1Khrbce3H13HPvbsyd89VXqRJIkqVqNy28IoTKE8GII4d7qxyuEEMaEEN4KIdwWQmhSfzGlEtO9O/zrX/Dqq7DZZjBtWupEkiSJ2p35PRJ4Y47H5wAXZlnWAfgK2L+QwaSS16sXDB8OY8fCFlvAd9+lTiRJUu7VqPyGEJYBtgCGVD8OQHfg9updrge2rY+AUknbdlsYNgyeeirenz49dSJJknKtpmd+BwHHAVXVj1sDU7Msm1n9eBKwdIGzSeVhl13g2mvjXMA77ggzZqROJElSbs23/IYQtgQ+y7Ls+Tmfnsuu2Tze3zuEMDaEMHbKlCkLGFMqcXvvDVdcEZdD3m03mDlz/u+RJEkFV5Mzv+sBW4cQ3gNuJQ53GAQsGkJoVL3PMsBHc3tzlmVXZVnWOcuyzm3bti1AZKlEHXQQDBoEI0fCXnvBTz+lTiRJUu7Mt/xmWXZClmXLZFm2PLAr8GiWZbsD/wZ2rN5tb+CueksplYsjj4QBA+CWW2DffS3AkiQ1sEbz32WejgduDSGcBbwIXFOYSFKZ69cvDns45ZT4eOhQqKxMm0mSpJyoVfnNsuwx4LHq++8AXQsfScqBk0+OtxZgSZIaVF3O/EqqCwuwJEkNzvIrpWQBliSpQVl+pdQswJIkNRjLr1QMTj4ZsgxOPTU+tgBLklQvLL9SsZh15tcCLElSvbH8SsXEAixJUr2y/ErFZs4C/NNPcP310MiPqiRJheDfqFIxOuUUaNwYTjgBpk+PK8I1aZI6lSRJJW++yxtLSqRfP7jwQhg5ErbfPpZgSZJUJ5ZfqZgddRQMHgz33QdbbQXffZc6kSRJJc3yKxW7Pn3ihW+PPgqbbQbTpqVOJElSybL8SqVgn31g2DB4+mno2ROmTk2dSJKkkmT5lUrFrrvCiBHw/PPQowd88UXqRJIklRzLr1RKttsO7rwTXnsNNtoIPv00dSJJkkqK5VcqNZtvHi+Ae+cd2HBDmDw5dSJJkkqG5VcqRT16wAMPwEcfwfrrw8SJqRNJklQSLL9SqdpggzgDxLRpsQC//HLqRJIkFT3Lr1TKOneGJ56Iq8FtuGGcDUKSJM2T5VcqdausAk8+CW3bwqabwoMPpk4kSVLRsvxK5WC55eIZ4A4d4kpwI0akTiRJUlGy/ErlYokl4LHHoGvXOCfwkCGpE0mSVHQsv1I5WXRReOihuArcgQfC+eenTiRJUlGx/ErlpkULuOsu2GUXOPZYOPFEyLLUqSRJKgqNUgeQVA+aNIFhw+KZ4AED4lLIl10GjfzIS5Lyzb8JpXJVWQmDB0Pr1tC/f1wK+ZZboHnz1MkkSUrGYQ9SOQsBzj4bLr0U7r4bNtkEvvwydSpJkpKx/Ep5cOihMHw4jB0bV4P74IPUiSRJSsLyK+XFjjvGmSA++gjWXRdeeSV1IkmSGpzlV8qTDTeMi2GEABtsAP/5T+pEkiQ1KMuvlDdrrAFPPw1LLRXnA7799tSJJElqMJZfKY+WXRaefBK6dIGdd44XxEmSlAOWXymvWrWChx+GrbeGww+Hfv2gqip1KkmS6pXlV8qz5s3jsIc+feCcc2C33WD69NSpJEmqNy5yIeVdo0Zw+eXQvj0cdxxMmgR33glt26ZOJklSwXnmV1Kc/eHYY2HECHjhhTgV2oQJqVNJklRwll9Js+24Izz6KHzzTSzATzyROpEkSQVl+ZX0c+uuC888E4c9bLIJ3Hxz6kSSJBWM5VfSr7VvH+cCXmcd2H13OOssyLLUqSRJqjPLr6S5a9UqLoe8xx5wyimw//4wY0bqVJIk1YmzPUiat6ZN4YYb4pngM86A99+PF8W1apU6mSRJC8Qzv5J+Wwjwj3/A9dfHVeG6dYPx41OnkiRpgVh+JdXMXnvFmSC+/jqOBX7oodSJJEmqNcuvpJpbbz147jlYbjnYbDO4+GIvhJMklRTLr6TaWW45eOop2GorOPJIOOggL4STJJUMy6+k2mvZEkaOhH794OqroWdP+OKL1KkkSZovy6+kBVNRAQMGwI03xkUxunaF119PnUqSpN9k+ZVUN3vsAY89Bt99Fy+Eu//+1IkkSZony6+kultnHXj2WVhxRdhySxg40AvhJElFyfIrqTCWXTbOA7zzznDCCfH2229Tp5Ik6Wcsv5IKZ6GF4JZb4Lzz4gVx664LEyemTiVJ0v+z/EoqrBDgmGPggQfgo4+gSxcYNSp1KkmSAMuvpPqy6aYwdmycF3iLLaB/f8cBS5KSs/xKqj8rrABPPw277gonnQQ77QTTpqVOJUnKMcuvpPrVogUMGwYXXAD/+lecGeKtt1KnkiTllOVXUv0LAY4+Gh56CD79NI4Dvvvu1KkkSTlk+ZXUcHr0iOOAV1oJttkmLo88c2bqVJKkHLH8SmpYyy8f5wPu0wfOOQc22QQ++SR1KklSTlh+JTW8Zs1g8GC44Ya4Mtyf/gT/+U/qVJKkHLD8Skpnzz1j+f3d76B793gm2OnQJEn1yPIrKa3VV4/jgHfcMY4B3nZbmDo1dSpJUpmy/EpKb+GF4dZb4aKL4P77oVMnePHF1KkkSWXI8iupOIQARxwBjz8OP/4I664bxwU7DEKSVECWX0nFZd1141nfjTeGQw6Jq8I5DEKSVCCWX0nFp00buO8+OO88uOsuWGsteOaZ1KkkSWXA8iupOFVUwDHHwFNPxfvrrx9ng6iqSp1MklTCLL+SilvXrnEYxPbbx9kgNtssLpEsSdICsPxKKn6LLAK33QZXXRUviFtzTRg9OnUqSVIJsvxKKg0hwIEHwnPPQevW0LMnnHhinBlCkqQasvxKKi2rrx4L8AEHwIABsMEGMHFi6lSSpBJh+ZVUelq0iEMghg+HCRPibBBDhjgnsCRpviy/kkrXTjvByy9Dt25xSMT228Pnn6dOJUkqYpZfSaVtmWXg4Yfh/PPj0shrrAEPPpg6lSSpSFl+JZW+igr4+9/h2WfjxXC9esGRR8L336dOJkkqMpZfSeVjzTXjxXBHHgkXXwxdusC4calTSZKKyHzLbwjh9yGEf4cQ3gghvBZCOLL6+VYhhIdDCG9V3y5W/3ElaT6aN4dBg+CBB+CLL+IiGeeeCz/9lDqZJKkI1OTM70zg71mWrQKsAxwaQlgV6Ac8kmVZB+CR6seSVBz++ld45RXYcks4/vg4JdqECalTSZISm2/5zbLs4yzLXqi+Pw14A1ga2Aa4vnq364Ft6yukJC2QNm3g9tth2DAYPz5OiXbRRVBVlTqZJCmRWo35DSEsD/wJGAMskWXZxxALMrB4ocNJUp2FAH/7G7z6KmSGviEAABSESURBVHTvDkcdFW/ffTd1MklSAjUuvyGElsAdwFFZln1Ti/f1DiGMDSGMnTJlyoJklKS6W2opuOceuPZaePHFOCXaFVe4MIYk5UyNym8IoTGx+A7Lsmxk9dOfhhCWrH59SeCzub03y7KrsizrnGVZ57Zt2xYisyQtmBBg333jWOA//xkOPjiODf7ww9TJJEkNpCazPQTgGuCNLMv+OcdLdwN7V9/fG7ir8PEkqR4su2xcCGPwYHj6aVh99XhG2LPAklT2anLmdz1gT6B7COGl6m1zYCCwaQjhLWDT6seSVBpCgD594vLIf/oT7L8/9OwJ77yTOpkkqR41mt8OWZY9CYR5vNyjsHEkqYG1bw+PPgpXXw3HHhvHAp91FhxxBFRWpk4nSSowV3iTpIoKOOggeP31OBPE0UfHMcGvvpo6mSSpwCy/kjTLMsvA3XfDLbfEqdA6dYLTToMffkidTJJUIJZfSZpTCLDrrvEs8K67whlnxDHB//1v6mSSpAKw/ErS3LRpAzfcAKNGwXffwXrrxXHA39R4mnNJUhGy/ErSb+nVK479PewwuPRS6NgRRoxwWjRJKlGWX0man4UXhosvhjFjoF072Hln2Hxzp0WTpBJk+ZWkmurSBZ59Fi66CJ56ClZbDc4+2wviJKmEWH4lqTYaNYpjf994A7baCk4+GdZaC/7zn9TJJEk1YPmVpAWx9NIwfDjcf38887vRRrDPPjBlSupkkqTfYPmVpLrYbLN4QdyJJ8LNN8PKK8Nll8HMmamTSZLmwvIrSXXVokUc+ztuXFwY47DDYO214fHHUyeTJP2C5VeSCmWVVeDhh+H222HqVNhwQ9htN5g0KXUySVI1y68kFVIIsMMO8YK4U0+Ff/0rDoUYMMBZISSpCFh+Jak+tGgB//hHLME9e8YxwautBvfemzqZJOWa5VeS6tMKK8Szvw8+GKdJ22or2GILGD8+dTJJyiXLryQ1hJ494eWX4fzz4YknYPXV4fDD4fPPUyeTpFyx/EpSQ2nSBP7+d5g4EXr3hssvh5VWggsucDywJDUQy68kNbTFF4/F9+WX4c9/hmOOgVVXjbNEZFnqdJJU1iy/kpTKaqvFFeIefDBeILfTTrDBBvDss6mTSVLZsvxKUmo9e8KLL8JVV8Fbb0G3brDHHvDee6mTSVLZsfxKUjFo1AgOPDCOBz7xRLjjjjg/8MSJ8OOPqdNJUtmw/EpSMVl44bhU8ltvwV57weTJMGZMnDN42rTU6SSp5Fl+JakYLbMMXH01dOkCiy0Gp58OK64Il1wCM2akTidJJcvyK0nFrEWLeGHcmDFxbuAjjoCOHWHYMKiqSp1OkkqO5VeSSkHXrvDII/DAA7DIIvGCuE6d4J57nB5NkmrB8itJpSIE+Otf4fnn4eab4dtvYeutYzEeNcoSLEk1YPmVpFJTUQG77QZvvAHXXBOXSN58c1hvPRg92hIsSb/B8itJpapxY9hvP3jzTbjySpg0CTbdFDbcEB57LHU6SSpKll9JKnVNmkDv3nF6tMsug7ffho03hu7d4cknU6eTpKJi+ZWkctG0KRxySCy/F10Er78el0vu0QP+/W+HQ0gSll9JKj/NmsUp0d55By64IJbg7t3jmOD777cES8o1y68klasWLeDoo+Hdd+Hyy+NqcVtsAWuvHZdPdp5gSTlk+ZWkctesGRx8MEycCEOHxinSdtwxLppx000wc2bqhJLUYCy/kpQXjRvDPvvEKdJuuQUqK2HPPWHlleNsEd9/nzqhJNU7y68k5U1lJey6K4wbB3feCa1bQ58+sNxycOaZ8MUXqRNKUr2x/EpSXlVUwDbbwJgxcTaILl3g1FPh97+Hww6LF8xJUpmx/EpS3oUAG20E990Hr74azwpfdRV06AA77wzPPZc6oSQVjOVXkjTbaqvBtdfCe+/BscfCQw9B165x1bh773WGCEklz/IrSfq1pZaCgQPhww/hn/+M06VttRWssgpceilMm5Y6oSQtEMuvJGneFl4Y+vaNq8bdfDO0agWHHw5LLw1HHRWnT5OkEmL5lSTNX+PGsNtu8N//xgvkttkmLpzxhz/AllvG4RGuHCepBFh+JUm107Ur3HgjfPABnHYajB0Lf/0rrLpqLMTffps6oSTNk+VXkrRg2rWL5ff99+NKcQsvDIceGodEHHYYvPJK6oSS9CuWX0lS3TRtCrvvDs8+C888Ey+MGzIE/vhHWG+9eJbY1eMkFQnLrySpcLp1i2eBJ0+GCy6Azz+HvfaKZ4P79oXx41MnlJRzll9JUuG1bg1HHx3L7qOPwqabwmWXxanSNtoIbr0VfvghdUpJOWT5lSTVnxBg443httvinMEDBsQL5XbbDZZZJk6XNm5c6pSScsTyK0lqGEssAf36xbmBH3ggngEePBjWWgs6dYJLLoEvvkidUlKZs/xKkhpWRUWcGm3ECPjoI7j44vj8EUfEleV22gnuvx9mzkybU1JZsvxKktJp3TquGPfCC/DSS3DwwfDYY7DFFrDccvFM8RtvpE4pqYxYfiVJxWHNNWHQoDhTxMiRsPbacP75cfGMTp3i7BGTJ6dOKanEWX4lScWlSRPYbju4+26YNCkW4kaN4Jhj4Pe/hx494NprYerU1EkllSDLrySpeLVrB0ceGRfQePNNOPXUOFvE/vvH13bYIZ4lnj49dVJJJcLyK0kqDX/4A5x+OkyYEMtwnz7w1FOxALdrB/vsA/fd5/zBkn6T5VeSVFpCgC5d4nCISZPgwQdh223hzjthyy3jlGp77QX33GMRlvQrll9JUulq1Ah69oTrroPPPotnfrfbLhbfrbeGxReHPfeEu+5yaIQkwPIrSSoXTZrA5pvD0KHw6acwahTsuGOcM3jbbWMR/tvfYPhw+Oab1GklJWL5lSSVnyZNoFcvuOYa+OSTODRil13goYfibZs28fXBg+PQCUm5YfmVJJW3xo3j0Iirr45nhB9/PK4m9/bbcMghcfq0zp3hzDPh5Zchy1InllSPLL+SpPyorIQNNoiLZ0yYAK+9Bv37x7HDp54aF9po3z6W41Gj4PvvUyeWVGCWX0lSPoUQV4874QR45hn4+GO46ipYbbV4lnjzzaFVq3h7ySUwcWLqxJIKwPIrSRLEuYIPPBDuvRe+/DKe+e3dO5beI46ADh3i5llhqaRZfiVJ+qXmzeMFcRddFIdHTJwYz/7+4Q8wZMjss8K9esEFF8BLL0FVVerUkmqgUeoAkiQVvRVXhMMOi9v06fGiuVGj4IEH4Jhj4j5t2sDGG0OPHnFbccU4tEJSUbH8SpJUG82axdkjevaECy+EyZPh0UfhkUdg9GgYMSLut+yys4tw9+6w5JJpc0sCLL+SJNXN0kvHVeT23DNOkzZhQizCjzwSl1weOjTut/LK8Je/xNkmNtgAllvOM8NSApZfSZIKJYRYcldeOc4h/NNPcTzwI4/EoRLDh8eZJCDOLzyrCG+wAayyClR4KY5U3yy/kiTVl8pKWHvtuB13XCzDr74KTzwRy/Cjj8LNN8d9W7eG9deP27rrQqdO8cI7SQVl+ZUkqaFUVsaFNNZcM148l2VxpbnHH4+F+Ikn4K674r6NGsFaa8E660C3bvHWi+ikOrP8SpKUSgiw0kpx22+/+Nwnn8CYMXHhjTFj4pjhSy+Nr7VuHUvwrEK89tpxyjVJNWb5lSSpmLRrB9tsEzeIQyVeey2W4VmF+L77Zu+/wgpxiESnTrEMd+oEbdumyS6VAMuvJEnFrLIS/vjHuPXuHZ+bOhWeew5eeGH2dscds9+zzDI/L8OdOsWp1hwyIVl+JUkqOYsuCptuGrdZpk6NM0vMKsPPPw/33BPHFUMcMrHGGrO31VePZ5UrK9P8DFIill9JksrBoovCRhvFbZZvv4Vx42IZfuWVuA0dGp8H2GefuGjHyJGzC/Eaa0CHDtC0aYIfQqp/ll9JkspVy5aw3npxm6WqCt5/PxbhF1+E776LM07cf388EwxxvuH27eN8xR07zt5WXjku4+zwCZWwOpXfEEIv4CKgEhiSZdnAgqSSJEn1o6IiXiS3wgrw5ZfxuXPPhR9+gDffjPMQjx8/exs9Or42S6tWs4vwyivHmSpWXNEhFCoZC1x+QwiVwGXApsAk4LkQwt1Zlr1eqHCSJKmBNG06+8K6Of30E3zwQSzCb745uxSPGjV76WaIQygaN4bBg2MZbt/+57dLLdWgP440L3U589sVmJhl2TsAIYRbgW2Aoiq/EydO5Ntvv2WjOcdASVKJWGuttQD8M0z1ok6/XyuvHEvt9Onw/fes1awZVFWx0YQJ8PLL8fk5hcBa++4LFRUsuvTSALz07LPxTHRFBRtdc83/nzku1t/7Ys1VrNZaay1atmyZOsav1KX8Lg18OMfjSUC3X+4UQugN9AZYdtll6/DtJElSUWnUKI4rbtly9gVya64Zb7MsFuDqcsz338dyW1U1+/3ffz/7/pNPxq/XtGm84K6iIo5NbtoUmjSZvTVu7Jhj1Uldyu/cfvOyXz2RZVcBVwF07tz5V6/Xt5VWWgmAQYMGNfS3lqQ6u+666wD/DFP9KOTvV02+1qx93n//fQDW6tYtjieePp1BAwbAhx/CBx9wXdOm8MMPDHrvvV9/kYoKWHzxOG9xu3bxdsklYYkl4uIebdrErXXreNu8eZ1/ttr8jJpt1n+vYlOX8jsJ+P0cj5cBPqpbHEmSlBvNmsVtkUXimOFZZpWmSy6BTz+Fjz+O2yefzL4/6/G4cXGfWTNV/FKLFrML8Zxb69ZxerhFFvn17SKLwO9+F89Eq+zU5ag+B3QIIawATAZ2Bf5WkFSSJEnNmsFyy8Xtt1RVwRdfxO3zz397e/vtePv11/P//i1b/rwUd+sWh27suScstND8txYt4rCNObcmTWbfL6UhHFkW/4ExYwb8+GPcZt2vPnv/q23KlKL8B8QCJ8qybGYI4TDgQeJUZ9dmWfZawZIVSLt27VJHkKQF5p9hqk+F/P2qydeatc+sYQ/zek+tc1VUxCEPbdvW/D0//gjffBNXxvv669m3c97/xXPtvvoqvu/pp+P8yLO2upizDM+636hR/JkqK2t3C7GkVlXV7vann2YX2l8W2znv11K7Xr3iUttFJmRZww3D7dy5czZ27NgG+36SJEn1KsvihXtzluE5txkz4pnRWducj+f12syZsZj+9NPPb+f23Jy3EEtwCLW/nXUxYePGNbs/5+NZw1d+uTVtGs+YJyjAIYTnsyzrPLfXiu9ctCRJUqkIIQ5vaNGidmeelUxF6gCSJElSQ7H8SpIkKTcsv5IkScoNy68kSZJyw/IrSZKk3LD8SpIkKTcsv5IkScoNy68kSZJyw/IrSZKk3LD8SpIkKTcsv5IkScoNy68kSZJyw/IrSZKk3LD8SpIkKTcsv5IkScoNy68kSZJyw/IrSZKk3LD8SpIkKTdClmUN981CmAK832DfcLY2wOcJvq8alsc5HzzO5c9jnA8e53xIdZyXy7Ks7dxeaNDym0oIYWyWZZ1T51D98jjng8e5/HmM88HjnA/FeJwd9iBJkqTcsPxKkiQpN/JSfq9KHUANwuOcDx7n8ucxzgePcz4U3XHOxZhfSZIkCfJz5leSJEkqz/IbQjgzhPByCOGlEMJDIYSl5rHf3iGEt6q3vRs6p+omhHBeCGF89bH+Vwhh0Xns914I4ZXq34exDZ1TdVOL49wrhPBmCGFiCKFfQ+fUggsh7BRCeC2EUBVCmOdV4X6WS1stjrOf5RIWQmgVQni4uls9HEJYbB77/VT9WX4phHB3g2Ysx2EPIYTfZVn2TfX9I4BVsyzr84t9WgFjgc5ABjwPrJ1l2VcNnVcLJoTQE3g0y7KZIYRzALIsO34u+70HdM6yzPkkS1BNjnMIoRKYAGwKTAKeA3bLsuz1hs6r2gshrAJUAVcCx2RZNtdi62e5tNXkOPtZLn0hhHOBL7MsG1j9j5fF5vF387dZlrVs+IRleuZ3VvGtthCx3P7SX4GHsyz7srrwPgz0aoh8Kowsyx7Ksmxm9cNngGVS5lH9qOFx7gpMzLLsnSzLZgC3Ats0VEbVTZZlb2RZ9mbqHKpfNTzOfpZL3zbA9dX3rwe2TZhlrsqy/AKEEM4OIXwI7A6cOpddlgY+nOPxpOrnVJr2A0bN47UMeCiE8HwIoXcDZlLhzes4+3nOBz/L5c/PculbIsuyjwGqbxefx37NQghjQwjPhBAatCA3ashvVkghhNFAu7m8dFKWZXdlWXYScFII4QTgMOC0X36Juby3/MaAlLj5HefqfU4CZgLD5vFl1suy7KMQwuLAwyGE8VmWPV4/ibUgCnCc/TwXuZoc4xrws1zkCnCc/SyXgN86zrX4MstWf57bA4+GEF7JsuztwiT8bSVbfrMs26SGu94M3Mevy+8kYKM5Hi8DPFbnYCqo+R3n6gsVtwR6ZPMYwJ5l2UfVt5+FEP5F/N9q/oVZRApwnCcBv5/j8TLAR4VLqLqqxZ/Zv/U1/CwXuQIcZz/LJeC3jnMI4dMQwpJZln0cQlgS+GweX2PW5/mdEMJjwJ+ABim/ZTnsIYTQYY6HWwPj57Lbg0DPEMJi1Vci9qx+TiUihNALOB7YOsuy/81jn4VCCAvPuk88zq82XErVVU2OM/GimA4hhBVCCE2AXYEGvXpY9cvPcm74WS59dwOzZtDaG/jVGf/q7tW0+n4bYD2gwS5qLMvyCwwMIbwaQniZ+AfkkQAhhM4hhCEAWZZ9CZxJ/KA9B5xR/ZxKx6XAwsT//flSCOEKgBDCUiGE+6v3WQJ4MoQwDngWuC/LsgfSxNUCmu9xrr4g7jDiP2DfAIZnWfZaqsCqnRDCdiGEScC6wH0hhAern/ezXEZqcpz9LJeFgcCmIYS3iLN2DISfdzBgFWBs9ef538DAhpzRoyynOpMkSZLmplzP/EqSJEm/YvmVJElSblh+JUmSlBuWX0mSJOWG5VeSJEm5YfmVJElSblh+JUmSlBuWX0mSJOXG/wEdBcwYslFPUQAAAABJRU5ErkJggg==\n",
"text/plain": [
"