HW (Derivatives).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl4VOXZP/DvTZDFhVVATEIA2YrggriL4FahWrGuYK1tXVAsVV9rrVb76s9e4qvVVhQV0KJWBcStLKJIiCwCKqGIgBRFZAlLQgIJSNgCz++Pe57rnBlmMjOZOTmzfD/XNdeZOZmEGwh88+xijAEREVGsGvhdABERpRcGBxERxYXBQUREcWFwEBFRXBgcREQUFwYHERHFhcFBRERxYXAQEVFcGBxERBQXBgcREcWlod8FeOHYY481HTt29LsMIqK0smTJknJjTJto78vI4OjYsSOKi4v9LoOIKK2IyPpY3seuKiIiiguDg4iI4sLgICKiuDA4iIgoLgwOIiKKC4ODiIjiwuAI2L8feO01YMQIIPQ03UmTgC+/9KUsIqKUI5l45njfvn1NvOs4Dh0C2rYFKiqA1auBbt30/oYNQEEB0L49sGkTIOJBwUREKUBElhhj+kZ7H1scAQ0aABddpM8LC537ZWV63bIFWL68/usiIko1DA6XSy7R66xZzr2dO53n7vtERNmKweFy8cV6/fRToKZGnzM4iIiCMThcOnYEunQBqqoAO0TiDo5584C9e30pjYgoZTA4QthWhx3ncAfHnj3AwoX1XxMRUSphcISw4xyhwWFnU4V2V1VUAM88A2zbVj/1ERH5jcER4oILNCQWLgQOHHCC49xz9bpggfPeQ4eA668H7rsPGDOm/mslIvJDygeHiHQWkX+KyLv18eu1bAm0aqWhsX07sGuX3u/VS68VFc57R48GZs/W56Wl9VEdEZH/PA0OERkvImUisiLk/kARWS0ia0Tkgdq+hjFmrTHmFi/rDNWqlV537HBaHB066LWyUq/ffw/86U/O57jHQoiIMpnXLY7XAAx03xCRHAAvABgEoCeAoSLSU0R6i8j0kEdbj+sLq2VLvW7fHjk4Pv5YZ1g1b66vGRxElC08PTrWGDNPRDqG3D4DwBpjzFoAEJFJAAYbY54AcLmX9cQqXIujfXsgJweortZ9rbZv1/s/+Qnw+ecMDiLKHn6MceQC2Oh6XRK4F5aItBaRMQBOFZEHa3nfMBEpFpHibQlOcbItDndwNG8OtGihz6uqnODo2NG5R0TkhYMH9QfWVOFpiyOCcNsERtxp0RhTAeCOaF/UGDMOwDhANzmsc3VwWhzurqpjjtHgqKjQ7qodO/S+DQ62OIjIC/v3AzfeqM8nTtSeD7/50eIoAZDvep0HYLMPdUQUrsXRrJnT4qisdFocBQV6DQ2OVauAsWN1yi4RUV3s2QNcdRXwzjvAzJnAd9/5XZHyIzgWA+gqIp1EpBGAIQCm+lBHROFaHJGCI1JX1T33AHfcASxa5Hm5RJSBdu0CLrsM+PBDoHVroKgI6NHD76qU19NxJwJYBKC7iJSIyC3GmBoAIwDMBLAKwGRjzEov64iXbXGUlurMqZwcoGnT4OCwXVW5uUDDhsC+ffqw1q51vgYRUTx27NBdLD79VCfmzJ0LnHaa31U5vJ5VNTTC/RkAZnj5ayfCtjjWr9drs2a6mjxci6NVK/24XSzYuLGeILhpk/NeIqJYlZdraHz1lXaFz54NnHCC31UFS/mV436wLQ53cABOcOzY4bQ4bHAATndVZaX2TbrvERFFU16uB8p99RXQtSvw2WepFxoAgyMs2+LYskWvocGxZYt2SzVurF1YoYsAN7uG+hkcRBSLbduACy8Evv4a6N4dmDMHyMvzu6rwGBxh2BaHFRocP/ygVxsw9uM2OGw3FcDgIKLoyso0NJYv1wHwTz8Fjj/e76oi82MdR8qLFhx24Nu+j8FBRHVlQ2PlSic0jjvO76pqxxZHGE2bAk2aOK+jtThsV5UNCQYHEcWitFSPcli5UrcvmjMn9UMDYHBEZEMB0FXjgBMcP/4Y/J7aWhycVUVE4ZSWakvjm2+Anj21pdGund9VxYbBEYG7uyq0xRH6ntDg4OA4EdVm61ZtaXzzDXDiiekVGgCDIyJ3iyNScLCriojitWWLhsaqVXpA3KefAm19OUCi7hgcESTS4ogWHCZkC8bKSqBfP2DcuLrXS0Spz4bGf/8L9O6t24i0aeN3VfFjcEQQrsVx1FHBO1OGG+M4cCB4m5GqKicojAEGDNCHOzxmzdKFPi+/nOzfBRGlirIyXdy3ejVw0knpGxoAgyOicC0O97Yj7ve4u6pKSzUUjjsOaNRIt0Teu1c/vnq17jkzb17w3vorAgfrbt3qvK9XLz1lkIjSX0UFcPHFTvfU7NnAscf6XVXdMTgiCNfiAIKDI1yLw3ZT5eYePvYxZ47zudXVzvPly/W6datuw/7eezo9b9SohH8bROQzu2GhXdxXWJjeoQEwOCIK1+IAwrc4Yg2OuXOdz92923luWxw1NbpZov0aCxboPSJKT1VVwKWXAkuXAl26aEsjnWZPRcLgiCCeFod7ryr7n/7xxwcHhzHBLQ4bHNXVwJo1zv0tW5yvsWuXbnZGROln1y7gZz8DFi8GOnXSMY1U3kYkHgyOCNwtDrsAEKi9q6qqylnDkZsbfEb5t986YxiA01W1alXwQPnWrUBJifN63rzEfh9EVP+qq4Gf/xxYuBDIz9fQyM+P/nnpgsERQSwtDvvc3VW1YYM+D+2qcndTAU6Lw3ZTWe4WB8DgIEo3e/YAgwfrv/njj9d1Gvak0EzB4IggWoujeXNnam7jxsARR+hMKRsQJ53kBEdlZXA3FeC0OOzAeMPAdpMbNwZP550/n+eWE6WLffv0jPDCQh3LKCpKzfM0EsXgiMDOemjWLHjthg0Od7CIOCGxaZMGzcknB7c45s/X51276jW0xXHWWXpdulS7rtq316bt9u06w4qIUtv+/cB11+k0+mOP1YHw7t39rsobDI4IWrYEnnoKePbZ4Ps2ONxdWUBwd9a552oLwgbHxo06btGkiXNucGhwXHKJXouL9ZqbC5x/vj6fPTvx3w8ReaemBrjhBmDqVP2/o7BQ96DKVAyOWvzxj8Bvfxt8L1yLAwgODvsfvn3vF1/otXt3p9urulrnd2/apNu4n3223rfH1ebl6eAaoFuRhG5TQkSp4eBB4Fe/0vVXzZvrThAnn+x3Vd5icMTppJOABg2APn2C79vWBeAEh723dKlee/QAjjxSn+/e7YRE587awnDLzdW+0txcnXn1ySfJ/X0QUeIOHQJuuQWYNEl/KJw50+lVyGQMjjiddJKeDfzkk8H3bYujSROgb199boNj3z699uih+10B2uKwmyK2aHH44S25uTrgPmKEvv7HP5L7+yCixBw6BNx+O/D66/rvesYM4Mwz/a6qfjA46qBVKx0Qd7PBcdZZOssKCG6FAHrClw2O3bt1gZD93JYtnc8DnBbIsGHaSpk5U/fuJyL/GQP8/vfAK6/oD4vTpgHnned3VfWHwZEkrVvr1XZTAYcHR2hXlW1xHHOMBpG71ZGXp9dWrYAbb9Tnb72V/LqJKD7GAH/4A/Dii/rD3pQpulV6NmFwJMmIEfoTyD33OPfciwVFgG7dgruq3C0OIDg43GMeV12l1w8/TH7dRBQ7Y4AHH9Su4yOOAN5/H/jpT/2uqv4xOJKka1fgueeCZ1u5WxwdO+rsKXdXlbvFAejaDcsdHP37a0tl2bLg7UiIqH49+qiObzZsCLzzju5FlY0YHB5yB0ePHnoN11UV2uJo1gw4+mjnc5s00b38AR2Ai6S4mJsiEnll5Ejgscd0VuWECbqtSLZicHioUSP9Tx/QgXGg9q4q2+IInZoLAJddptdI3VUHDujpYueeC5SXJ147ETmefhp46CHtcn7jDeDaa/2uyF8MDo/ZVkdtLY7Qrio7MO5mm8SFhc6Jgm47dujXq67mEbREyfTcc7oYGADGj9cV4tmOweExO0BugyPcOg7b4ujXT2dR2daFW16erkatrtbzyUPt2OE8f/FFbYEQUWLGjAHuvlufjx0L/OY3vpaTMlI+OETkJyIyRkTeFZHhftcTr3vvBYYMcTYxjLSOA9BwKS93vlFD2Xniy5Yd/jF3cJSUAP/+d+K1E2Wz8eOB4YH/cZ5/XtdUkfI0OERkvIiUiciKkPsDRWS1iKwRkQdq+xrGmFXGmDsAXAegr5f1emHYMGDiRJ26B9TeVQUcvrDQzY6TrFp1+MdscDQI/I2OHl33momy3ZtvArfeqs+fecbZwYGU1y2O1wAMdN8QkRwALwAYBKAngKEi0lNEeovI9JBH28DnXAHgMwBpv09sbYPj0fTsqddwK8htcAwapFMFP/vM+fpEFLvJk4Ff/1rXbIwcqb0GFMzT4DDGzAOwPeT2GQDWGGPWGmP2A5gEYLAxZrkx5vKQR1ng60w1xpwD4Jde1lsfapuOG427xRG6W+72wJ9yhw7AKafoPjpffpl4vUTZ5IMPdPD70CHgkUd0sR8dzo8xjlwAG12vSwL3whKRASLynIiMBRBxFYOIDBORYhEp3rZtW/KqTbJGjbRFUFMDVFToPXdXVW3atdMFhpWVweeXA06Lo2VLnZIL6HnHRBSb6dOB66/XbdIffFCDg8LzIzjC9eJHPG3CGDPHGHOXMeZ2Y8wLtbxvnDGmrzGmb5s2bZJSqFfcA+RA7MEhEnmcwx0c55yjzxkcRLGZORO4+mqdjXjvvcDjj9c+3pjt/AiOEgD5rtd5ADb7UIdvbHcVoJukNWoU++fGExyLFvG8cqJoioqAK6/Uo19HjNDFfgyN2vkRHIsBdBWRTiLSCMAQAFN9qMM3tsUBxD6+YUUaIHcHR16enldeVcWt2IlqM3++nrS5d6/OgBw1iqERC6+n404EsAhAdxEpEZFbjDE1AEYAmAlgFYDJxpiVXtaRatwtjniDI5YWB8BxDqJoFi3SHRmqq3Vh30svOdPZqXZez6oaaoxpb4w5whiTZ4z5Z+D+DGNMN2PMCcaYx72sIRW5Wxyxjm9Y7hbHvn3O9iOhwcFxDqLIiouBgQOBH3/UWVSvvMLQiAf/qHyQSFdVfr62WEpLdTfd/HydZWWDo1UrvdoWx6xZOkskkupqba6HTu8lylRffaVnaOzcCVxzjR79mpPjd1XphcHhg0S6qho0AE48UZ9XVuoWJV9/fXiL49RT9QyQzZuBefMif72HHtJTC7lFCWWDFSuASy7Rfy+DB+v26A0b+l1V+mFw+CCRrioAePZZ4IEHnGNqv/tOp/bm5DjneIg4u3jWduTs9Ol6nT8//jqI0sl//6tHD5SX69jG2287WwFRfBgcPkikxQHo+MUTTzjdUUuX6rVly+AZIb8MrLN/910dDwm1eTOwZo0+//rr+OsgShdr1gAXXgiUlWmL4733dCo81Q2DwweJjHG4deyo1//8R6/uY2sBHUg/5RSdlhvu5MC5c53ny5ZxnIMy0w8/aGhs2QIMGKDdsvaANaobBocPEu2qsmxw2G3WQ4MDcFodEyYc/rE5c5zn5eU64E6USTZs0NDYuFFb6NOmBbf4qW4YHD5ItKvKssFRXa3XcMFx/fV6nTHDeZ9lWxz2lEJ2V1Em2bxZQ2PdOuDMM/XfgB0DpMQwOHyQrBZHhw7Br8MFR34+cMYZGhoff+zc37oVWL1aa7nuOr23fHndayFKJaWlGhrffw/06aPf+4n8kEbBGBw+SFaLo0kT55xywFnDEerqq/X63nvOPTtF95xzgNNO0+dscVAm2LZNZ0+tXg2cdBLwySfOEc6UHAwOHyRrcBxwuquA8C0OALjqKr1On+7MrvrXv/R6wQX6jwtgcFD6275dZ02tXKmTQ2bNAlq39ruqzMPg8EGyuqqA2IKjSxcNh507gdmzdY+eDz/U/t5bbwV69dL3ffONbivtdvCg7hpKlOoqK3VF+LJlQLdu+r3etq3fVWUmBocPktVVBcQWHIDTXXXvvcBdd+nze+4B2rTR8OrUSQPi22+DP69fPw2WmprE6iTy0q5demzykiVA5866Vfpxx/ldVeZicPigvruqAGD4cN1Zd/Vq3eCtRQvgD39wPt63r17fftu5t3Ontk6++06nMxKlot27dSX4558DBQUaGrkRzxSlZGBw+MDd4qiPripAWxbFxcDvfqd784wcGTxg+Pvf63XUKGffq++/dz6+fn1idRJ5obpaz9P47DM9h6aoSMODvMXg8EF9j3FYRx4JjB6t/9iGDw/+WL9+OhNl507dCwtwtiMBGByUevbuBX7xC+DTT7VbqqhIu6nIewwOH9jgOOqoxLdzdq/liDQdN1Skjd0eeUSvzz6r25SEC45583SzOCI/7d+vW6J/8om2pouKgK5d/a4qezA4fNCmjf6EdMopiX+tJk10gV/79kC7dol9rX79gLPP1lbHnDmHB8fGjbqoyk7vJfLDgQPAkCE6M7B1a509ZU/GpPrBneh90LixDjg3apScrzd/vv4Eloyvd/75OiC+eLHWaK1bp5spHjyox9ZWV3PPH6p/NTXAjTcCH3ygY3SzZgG9e/tdVfZhi8MnRx+dvOBo1Ch5e/CcfrpeFy8+vMWxYoXzOnTaLpHXDh4EfvtbYPJknY34ySd6YBnVPwYHBbHBsWiRbkNtx0M2bgxeWb56df3XRtnr0CFg2DDgzTd1bPCjj5zvVap/DA4Kkp+vq2137dLXXbromMz+/ToAaTE4qL4Yo9PIx48HmjbVXW7POcfvqrIbg4OCiDiLAQENDjsvvrzcuc+ZVVQfjNEdDsaM0Ykg06Y5RyaTfxgcdBh3F4A7OADnaFq2OMhrxgD33w8895yO433wga41Iv8xOOgwtQVHv356Xb06+lGzEybo3Hr37CyiWP3lL8DTT+tOB+++Cwwc6HdFZDE46DChweFenX7++TrmsXs3sGlT7V9n4kSdmTVpkidlUgb761+Bxx/XBbKTJum2IpQ6GBx0mLZtgR499Ce9E08MbnH07g10767Po3VXlZTo9csvvamTMtOTTwL/+79AgwY6i8ru7Eypg8FBYU2dqqvHc3ODg6NXLw0VIPwAuT0oCnB21P3yy+jdWkQA8I9/AA88oGNpr76qK8Qp9TA4KKyuXYFzz9XnHTtql0HTpno/UovjpZd008a33gL27AEqKvR+WRmwYcPhv8aUKdqdRQQAL76o58UAwLhxwE03+VsPRcYtRyiq5s21n/nII3VBYLgWx+zZujX7wYPAzJm6f5bbl18Gt1xqavSnyb17dRyFi7my2yuv6FoNAHjhBT2ZklJXyrc4RGSAiMwXkTEiMsDverLVNdfoYTmAc9Ts0qXaBbVxI3DttRoaQPiDn0LHOTZt0tAAgD//2bu6KfW9/rquCge0q+rOO/2th6LzNDhEZLyIlInIipD7A0VktYisEZEHonwZA+BHAE0AlHhVK8WuoEBnVpWX6+aHb7yhhz/ZhYPffecMjNut3kODw32+R2Fh8Kp0yh6TJgE336w/gDz5pC72o9QXNThEJEdE/lbHr/8agKDZ1yKSA+AFAIMA9AQwVER6ikhvEZke8mgLYL4xZhCAPwH4f3Wsg5JIBDjzTH3+xRfA3Ln6/N57tTurosLZ12rwYL0WFwefW26Dw270+PDD3tdNqeW993Sn20OHgMce08V+lB6iBocx5iCA00TsmuHYGWPmAdgecvsMAGuMMWuNMfsBTAIw2Biz3BhzecijzBhzKPB5OwA0jrcG8oYdw1iwQB8AMGCAjlcAeioboLuXduqk27BfeSUwdqz+dGmD47bbdGffRYvCD6BTZpo2Tce4Dh4EHnpIF/tR+oi1q2opgCki8isRuco+6vhr5gJw94CXBO6FFfi1xgJ4A8DoWt43TESKRaR427ZtdSyNYmWD4403dDFgt256mJQ9hW3pUr3m5zsHP334IXDHHXquhw2Onj2BSy5xPk6Z7+OPdcyspga47z5d7EfpJdbgaAWgAsCFAH4eeFxex18zXMsl4ix/Y8z7xpjbjTHXG2Pm1PK+ccaYvsaYvm3atKljaRQrOwuqqkqv/fvr1QaHXbeRlwf87W96+NPFF+u9pUud4CgoAC67TJ8zODLf7Nl6Tvj+/cBddwFPPeXsf0bpI6bpuMaY3ybx1ywBkO96nQdgcxK/PtWDVq2C96EaMECvoec+5+frfww9emhwFBbq+Me6dfrxggKgTx99XlSk6z+aNq2P3wHVt/nzgSuu0Nl0t9+uZ9szNNJTTC0OEckTkQ8CM6RKReQ9Ecmr46+5GEBXEekkIo0ADAEwtY5fi3zkXqsR2uIAdM2Hu/F30kl6XbbMGc8oKNAurj59NDTs2Ahlls8/1+nc1dXAb36ji/0YGukr1q6qV6H/uR8PHY+YFrhXKxGZCGARgO4iUiIitxhjagCMADATwCoAk40xK+tSPPnLzqw64QTdmgQIDo68PN1vyLLB8fnnujVJq1a60hxgd1UmW7JEd7b98Ufghht0sV+DlF9BRrWJ9a+vjTHmVWNMTeDxGoCoAwnGmKHGmPbGmCOMMXnGmH8G7s8wxnQzxpxgjHk8gfrJR1deqaFhV/wCQLt2zvnn+fnB7z/+eA2L/fv1tXsl+eWBEbP33tOWB2WGZct08kNVlQ6Iv/66bl9D6S3W4CgXkRsDazpyRORG6GA5ZbH8fN02/X/+x7kn4rQ68kI6M0WcVgcQHBynn67dVaWlOmWX0t/KlTqutWOHjm1MmKA7LlP6izU4bgZwHYCtALYAuCZwj+gwNjhCWxxA5OAQAR59VJ//3/9pXzilr2+/1dP6ysuBQYOAyZN1zIsyQ0wrxwFcbYy5whjTxhjT1hhzpTFmfbTPpexkB8rPOuvwj0UKDkC7q/r21VbHmDHe1Ufe+v574MIL9e/xoou0+7Exl+5mlFhXjg+uh1ooQ9x5J7Btm46BhHIHh/tkQSC41fH448627Nbu3cmskrywfr2GxqZNeszwlCmcXp2JYu2qWiAio0Wkn4j0sQ9PK6O0duyx4e+feKIzDTO0xQHolM0LLwS2b9dT4KyJE3XQ/e23k18rJUdJif7dbdgAnH22zpA76ii/qyIvxBoc5wA4EcBjAJ4JPJ72qijKXEceCfz0pzp91x4I5SYCjBqlM2/GjNFZOYCuPgecDRUptWzdqt1Sa9dqd+NHHzlTrSnzxDLG0QDAS8aYC0IeF9ZDfZSBpk7VfvBIP4326qVTfA8dAu6+W8PD7n21aVP91Umx2bZNQ+Pbb4GTT9aDvJo397sq8lIsYxyHoAv2iJKiUaPog6WPPqrdXXPn6tbbFoMjtWzfrus0vvlGuyFnzXLOYKHMFWtX1SwRuU9E8kWklX14WhlltZYtgZEj9fkK1zFgDI7UUVmp3Y7LlunuyIWFwVvMUOaKZx3H7wDMA7Ak8Cj2qigiQE+Gsxsg9uih4x+lpcCBA/7WRcCuXbo+Y8kS3T2gqAg47ji/q6L6ElNwGGM6hXl09ro4ym45OcDLLwO9e+v03HbtdLv2rVu1FfLww8655VR/du/WvcU+/1xnxhUVOXuVUXaoNThE5H7X82tDPjbSq6KIrD59dBv2q65y/nPatElD4/HHgVejbrVJybRnj24fMn++/n0UFQEdOvhdFdW3aC2OIa7nD4Z8bCCI6pE7OOy4R2Ghf/Vkm337NMBtt1RREdCZ/Q5ZKVpwSITn4V4TecoGx5o1ul4A0P+8Dh70r6ZssX8/cO21euxrmzZ6kl+3bn5XRX6JFhwmwvNwr4k8ZYOjqMg5mray0lnjQd6oqdFzNKZN06m2hYV6Vjxlr2ibHJ8sIjuhrYumgecIvG7iaWVEIWxwzJ8ffL+wUFcrU/IdPAjcdJNuVNi8OfDJJ8H7jVF2qrXFYYzJMcY0M8YcY4xpGHhuX3OTZKpXNjjsQU+2q2T2bH/qyXSHDgG33ursE/bxx8Bpp/ldFaUCHuBIaSN0yqc9efCzzzgtN9mMAYYPB157TfcX++ij8NvkU3ZicFDaCA2O/v11jcfevboQjZLDGN0jbNw4oEkTHds47zy/q6JUwuCgtNGsmf70C+gq8m7dnDM9tm3zrayMYgxw//3A88/rnmL//rdulU7kxuCgtCHitDo6d9YDgpo109c7d0b+PIrdX/4CPP20ng3+7rvApZf6XRGlIgYHpRUbHHY6qN2+u6rKn3oyyV//qqvxc3L0wKyf/9zviihVMTgorUQLjueeA+66i4sC4/XUU3riYoMGwJtv6gpxokgYHJRWfvpT7aIaNEhfhwbHI49o/3xRkT/1paNRo4A//Um7Al99FRgyJPrnUHZjcFBauekm3dK7f3997Q6Ogwd1JTkAvP66P/WlmzFjgHvu0edjx+qfL1E0DA5KOzk5znN3cNjQAID339eAAXSm0MiRuoCNHOPH61oNABg9GrjtNn/rofTB4KC0ZoNj505gxw7n/p49OisIABYsAB56yPlPEtAzPQ4dqr86U81bb+mqcAB45hlnMSVRLBgclNbcLQ53cABOd9WCBXpdt07ft2gR0L69hkk2eucd7ZKyLbF77/W7Iko3DA5Ka3YdR1UVsH27Pj/9dF3xPHcusGULsHCh8/7ly50uq6+/rt9aU8GUKbrT7aFDOovqwdBTdohiwOCgtBauxVFQAFx0kT6fPj04OL7+GvjyS30e2kLJdDNm6JkaNTU6i+rRR/2uiNJVygeHiPQTkTEi8oqILIz+GZRNwgVHq1bO4rVRo4Dycuf9y5YBixfrc/dgeqYrLNS1GQcO6CyqJ57Q6bdEdeFpcIjIeBEpE5EVIfcHishqEVkjIg/U9jWMMfONMXcAmA6AkywpyDHH6H+AP/7oBETLlsDll+vzlSv12qaNXj/8EKio0OfZ0uKYM0fPCd+3TycI/P3vDA1KjNctjtcQcja5iOQAeAHAIAA9AQwVkZ4i0ltEpoc82ro+9QYAEz2ul9JMgwYaHgCwfr1eW7bUFeZ9+jjvu/lmvW7a5NzLhuBYsEBDdM8enUU1ejRDgxLnaXAYY+YB2B5y+wwAa4wxa40x+wFMAjDYGLPcGHN5yKMMAESkA4CqFcs0AAAO00lEQVQqYwy3sqPD2O4qGxytWunVvdfS4MHAcccFf96+fZl9jscXX+gK+927dRbV2LEatESJ8uPbKBfARtfrksC92twC4NXa3iAiw0SkWESKt3GP7axig+OHH/TasqVebXA0bqytj3BHnrpbHa+/rq2XTDjbY8kS3dl21y5g6FBd7MfQoGTx41spXEPZ1PYJxphHjDG1DowbY8YZY/oaY/q2sR3alBVscGzYoFcbHH366N5Vo0dreLiDo3VrvbqDY/RoHSuZOdP7mr20bJnu6VVVBVx9NfCvfwWvtidKVEMffs0SAPmu13kANvtQB2UIGxwHDujVBodI8JRTGxwFBdptVVHhzKwqKwOKi/X5Rnd7OM2sXAlcfLGuabniCmDCBD1bgyiZ/GhxLAbQVUQ6iUgjAEMATPWhDsoQdhGgZcc4Ql16KdCjB3D77U642BaHu5VhWy5//ztwxhm6PUk6WL1a16+Ul+vYxuTJeoofUbJ5PR13IoBFALqLSImI3GKMqQEwAsBMAKsATDbGrPSyDspstsVh2VAI1bYtsGqVrpa277Etjo8+ct5nWxwvv6xrPv7xj+TW64UfftDQKC3VFsf772v3HJEXPG3EGmOGRrg/A8AML39tyh7u4HBPz61NixZ63bFDt2P/5BPnYxs36r21a/X12LHAww/H9nX9sHmzhsWmTUC/frqtSJMmfldFmYzzLCjtuYOjRYvYZg+5u6qKi3W8o1Mn/Q+3slJbJvv363uqqvSAo1S0bZuGxtq1ukfX9OnAkUf6XRVlOgYHpT13cEQa3wjl7qqaPVufDxwI5AembdgTBJs21euzz6becbSVlTpus2oV0KuXdreFjvcQeYHBQWnPHRyRxjdCubuq1qzR56ecAnTooM9tcFx3nc7A+uEH532pYPdu4LLLgKVLgS5dgFmznCnGRF5jcFDaq0twuFsc69bp84ICp8Uxd65eu3fXQXUAqK5OuNSk2LtXV8IvXKj1FhYeviqeyEsMDkp7ibY47FYlHTs6LQ4726pLF6e7ym5P8v77wFSfJpCXlwOXXKLda+3a6bWgwJ9aKHtxaRClPXe/frxjHBUVzrqNDh2cFofVpYszQ2nvXh0wHzpUD0IqLdWB6N/9DhgwAPjVrxL6bUS1bp2Gxpo1uonjxx8DXbt6+2sShcPgoLSXSFfVt9/qwUZt22rLwrY4rBNOcIJjzx7trrKzrexJguPHa3eRl8GxebOu01i7Fjj1VGDaNA0PIj8wOCjtJdJVtW+fXjt21Ku7xdG2rbZm3F1Ve/Y4H582zRn32LAB2LnTm1lN27drS2PtWuC007R7KnTRI1F94hgHpb26dFU1bx58LoUdJ3AHR5cuenV3Vbm3YZ8xw2l1AMA33wDffacD6hMmxF5/bQ4cAK65Rr/2iSfq1igMDfIbg4PSXsOGwFFH6fNYWxwNGgQHjg2Oo492vkZocOzZE9zi2LnT6bYCgBUrgLfe0u6vSMGxbBkwcqTT0omkulrPRr/5ZuDTT3XW1Mcfc8otpQZ2VVFGaN5c1zbEGhyAvreqSp/bripAxzl27HAGniN1VVldu2pLY8UK56haO+Ae6o9/1DUXDRsC998f/j1r1wLnnQds2aKvmzTRbUTy8mL/vRF5iS0Oygj2CJZ27WL/HHfIuKe0du6s1+7d9eruqrLBYe81bKj7WAHAV18BixbpczvFN9RXX+n1738PH0J21taWLboFyvXXa0vjjDNi/30ReY0tDsoIo0YB//mPbpseKztADgQHx2OP6SryK67Q1+6uKjvGcdZZOuZwwglA//56b/58naYLaDdWVVXweERpqe4tZZ+/+ipw553BNT38sHZRdeigp/jF04Iiqi8MDsoI/fs7/4HHKlKLo1cvfVjhuqqOOkpPDAQAY3Rs5Mcfg7/+hg1A797O6xUrnK+3Zw/w1FPAbbcBRxyh95ctA55+Wk/rmzCBoUGpi11VlLVsi6Nly9qn0YbrqrJhAujsLHfQWKHdVcuX6/WXv9RusPXrg7dzv/9+DaERI4Bzz43v90JUnxgclLXsT/TugfFwws2qcgcHEBwc9j/90AFy2+I4+WTghhv0+ZQpep01S0OkWTNnzIQoVTE4KGvZ4Ii211O4rqpIwdGihW7PDkRucfTuDVx5pT6fMkVXrtsZVg8+CBx7bHy/D6L6xuCgrGXXaZxySu3vC7cAMDQ4+vXTLquf/cxpwWzYoIPl69bp1U7V7dVLw6NTJ6CsDPj973W2VV4ecPfdyfidEXmLwUFZ69prgc8/Bx54oPb3heuqCj2atU8fbVGMGePsd7VhA/D88xoQt96q60zat9dFfCJOq2PMGL3+9a+HBxJRKmJwUNZq0AA480ygcePa3xdLVxWg03OPOcYJjvXrgTfe0Of26Fn3WIgNDkBbIF7vrkuULJyOSxRFtFlVoXJzNZQ2bwY2bQr+mHt67rnn6kaKZWXAk0/qNFyidMDgIIoilllVbkccARx/PFBSoq8vu0zXeMydC5x9tvO+nBwdHF+/3hlQJ0oHDA6iKNxdVZEGx0MVFDjBceWVOv120SLggguC33fWWfogSicc4yCKora9qiJxHwh12WV6UuBFF2kXFlG647cxURS2dRFrVxXgrA05/XSdSUWUSRgcRFHEOzgO6Il9ADB8uHd1EfmFYxxEUcSyADDUhRfqYU2NGnlbG5Ef2OIgiiJcV1W0MQ6AoUGZi8FBFIVdILh/v67+BrjCm7JbygeHiPQUkcki8pKIXON3PZR9RJzwqKzUK4ODspmnwSEi40WkTERWhNwfKCKrRWSNiETZKQiDADxvjBkO4CbPiiWqhQ2KHTuCXxNlI68Hx18DMBrAv+wNEckB8AKASwCUAFgsIlMB5AB4IuTzbwbwBoBHROQKAK09rpcoLDumwa4qIo+DwxgzT0Q6htw+A8AaY8xaABCRSQAGG2OeAHB5hC/1u0DgvO9VrUS1CR0Mj2VwnChT+TEdNxfARtfrEgBnRnpzIHj+DOAoAH+r5X3DAAwDgA7uZbtESRDawmCLg7KZH8EhYe6ZSG82xqxDIBBqY4wZB2AcAPTt2zfi1yOqC3cLo2FDfRBlKz9mVZUAyHe9zgOw2Yc6iGLmDg62Nijb+REciwF0FZFOItIIwBAAU32ogyhm7rBgcFC283o67kQAiwB0F5ESEbnFGFMDYASAmQBWAZhsjFnpZR1EiXK3ODgwTtnO61lVQyPcnwFghpe/NlEysauKyJHyK8eJUgG7qogcDA6iGLDFQeRgcBDFgGMcRA4GB1EM2FVF5GBwEMWAXVVEDgYHUQwYHEQOBgdRDNhVReRgcBDFgIPjRA4GB1EM2FVF5GBwEMWAXVVEDgYHUQzY4iByMDiIYsDgIHIwOIhi4A4LDo5TtmNwEMWALQ4iB4ODKAYcHCdyMDiIYsAWB5GDwUEUAy4AJHIwOIhiwK4qIgeDgygG7KoicjA4iGLA4CByMDiIYnDEEUCDwL8WBgdlOwYHUQxEnFYHB8cp2zX0uwCidHH66cCGDUDr1n5XQuQvBgdRjGbPBmpqgIb8V0NZjv8EiGKUk6MPomzHMQ4iIooLg4OIiOLC4CAiorgwOIiIKC4MDiIiiguDg4iI4sLgICKiuIgxxu8akk5EtgFY73cdAI4FUO53ERGkam2sKz6pWheQurWxrsgKjDFtor0pI4MjVYhIsTGmr991hJOqtbGu+KRqXUDq1sa6EseuKiIiiguDg4iI4sLg8NY4vwuoRarWxrrik6p1AalbG+tKEMc4iIgoLmxxEBFRXBgcREQUFwYHERHFhcFRz0Sks4j8U0Tedd0bICLzRWSMiAxIobp+EqjpXREZnkJ1HXbPTyLSU0Qmi8hLInKN3/VYItIv8Pf3iogs9LseKxW+3yNJhe/5cFLte57BEQcRGS8iZSKyIuT+QBFZLSJrROSB2r6GMWatMeaW0NsAfgTQBEBJqtRljFlljLkDwHUA4l6Y5GFd4f4M6yQZNQIYBOB5Y8xwADelSl3GmPmBv7/pAF5PlbqQ4Pe7l7Ul+j3vYV1J+55PCmMMHzE+AJwPoA+AFa57OQC+B9AZQCMAywD0BNAb+g/W/Wjr+rx3Xc8bBK7tALyVKnUFXl8BYCGAG1Kprkj3/Kgx8HgBwN8ALEjB77XJAJqlSl2Jfr97/WeWyPe8H/8O/Hj4XkC6PQB0DPkGOBvATNfrBwE8GMPXCfcfYaO6fmN4WVfg/oepVley/hElscYcAFNS6XsNQAcALyerpiT/edX5+93r2gLvrdP3vMd/ZikRHOyqSlwugI2u1yWBe2GJSGsRGQPgVBF5MHDvKhEZC+ANAKNTqK4BIvJcoLYZKVTXYfeSLN4aO4rIOAD/grY6vBJXXQG3AHjVs4pUvH9eXny/J6s2L77nk1GX19/zcWnodwEZQMLci7iq0hhTAeCOkHvvA3g/BeuaA2BOUqtKTl2H3UuyeGtcB2CYZ9U44qoLAIwxj3hUi1u8f15efL9HEm9tc5D87/lwEv534Ce2OBJXAiDf9ToPwGafanFjXXWXqjWyrvilam2pWldMGByJWwygq4h0EpFGAIYAmOpzTQDrSkSq1si64peqtaVqXbHxe5AlnR4AJgLYAuAA9CeGWwL3fwbgW+gsiYdYV2rXlQ41sq7MqS1V60rkwU0OiYgoLuyqIiKiuDA4iIgoLgwOIiKKC4ODiIjiwuAgIqK4MDiIiCguDA6ieiIij4rIfX7XQZQoBgcREcWFwUHkIRF5KHBYTyGA7oF7i+3JdyLyhIg87meNRPHi7rhEHhGR06B7EJ0K/bf2HwBLAPwGwLsicheAgQDO9KtGorpgcBB5px+AD4wx1QAgIlMBwBizUkTeADANwNnGmP0+1kgUN3ZVEXkr0mZwvQFUQo9PJUorDA4i78wD8AsRaSoixwD4OaAn4AFoDT2L+jkRaeFjjURx4+64RB4SkYcA3ARgPXRL7W+gpwVeZIzZGBjnOM0Y82sfyySKC4ODiIjiwq4qIiKKC4ODiIjiwuAgIqK4MDiIiCguDA4iIooLg4OIiOLC4CAiorgwOIiIKC7/H1aEBfP4N5QBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1892b8cea58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def derifw(y,x,dx): #forward difference method for a derivative\n",
    "    return (y(x+dx)-y(x))/dx\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x=(np.pi/4) #set the x value to evaluate cos'(x) at.\n",
    "dxarr=[]\n",
    "yarr=[]\n",
    "for dx in 10**np.arange(-16., -0.5, 0.1): #create a loop to store 10^-16,10^-15,...,10^-1 in an array\n",
    "    y=derifw(np.cos,x,dx) \n",
    "    dxarr.append(dx) #stores the dx value used\n",
    "    yarr.append(y) #stores the numerically evaluated y value using the specific dx\n",
    "    \n",
    "plt.plot(dxarr, abs(yarr+np.sin(x)), 'b-', lw=2) #plots the error (defined as true value-numerically calculated value) vs the corresponding dx value    \n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('dx')\n",
    "plt.ylabel('Error')\n",
    "\n",
    "plt.plot(x,10^-5)\n",
    "\n",
    "plt.show()\n",
    "#the error as Error~dx^n shows that n is 1. It is 1st order. The optimal value for dx is around 10^-7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4VGX6PvD7SQhFQARpSkcQQUGRiFhQFAsqiq4NlNWvsotY1rYW1LWtsqiILgoWVARZG2JDRBFsqKBSVBQQ6RJRSkB6z/v745n3d85MZpJJ5sw5Zyb357pynZmTSfIECHfeLsYYEBERlVdO0AUQEVFmY5AQEVFKGCRERJQSBgkREaWEQUJERClhkBARUUoYJERElBIGCRERpYRBQkREKWGQEBFRSioFXYAf6tata5o3bx50GUREGWX27NnrjDH1SntdhQiS5s2bY9asWUGXQUSUUURkRTKvY9cWERGlhEFCREQpYZAQEVFKGCRERJQSBgkREaWEQUJERClhkCTh55+BRx4Btm8PuhIiovCpEOtIUrFpE3DGGcDy5cCKFcCIEUFXREQULmyRlOKGGzREAOCpp4BPPgm0HCKi0GGQlOCtt4DRo4GqVYErrtB7/foBW7YEWhYRUagwSBIwBhg6VB8PGQI8+yzQsaO2Tux9IiJikCQkAnz0EfDkk8C11wJ5ecCwYfq+oUOBdeuCrY+IKCwYJCWoXh247joNFQDo2hXo0QPYvBl4+OFgayMiCgsGSRk9+KBehw8HVq0KthYiojBgkJRRp05Ar17Ajh3Ayy8HXQ0RUfAYJOVw+eV6ff31YOsgIgoDBkk59OgB1KgBzJ4NLFkSdDVERMFikJRDtWravQUA48YFWwsRUdAYJOV08cV6ZfcWEVV0DJJyOu00oFYt4IcfgEWLgq6GiCg4GRkkItJSRF4QkfFB1VClCnDccfp4/vygqiAiCp7vQSIio0RkjYj8FHO/h4gsFJHFIjKwpM9hjFlqjOmX3kpLV7++Xu0q9zfeAK68UqcGExFVFEG0SEYD6OG+ISK5AEYAOANAOwB9RKSdiLQXkYkxb/X9Lzm+unX1aoNk8GDgxReBiRODq4mIyG++n0dijJkmIs1jbncGsNgYsxQAROQ1AL2MMYMB9PS3wuTVq6dXGySrV+v1s8+ACy4IpCQiIt+FZYykEYCVrucFkXtxicj+IvIMgI4ickeC1/QXkVkiMmvt2rXeVhvhbpEYA6xZo88//TQtX46IKJTCEiQS555J9GJjTKExZoAx5qBIqyXea0YaY/KNMfn1bNPBY+4g+fNPYM8efT5/vtM6cfv+e+C999JSChFRYMISJAUAmrieNwYQ+i0R3UFiWyPW558Xf/3FF+tCxt9/T39tRER+CUuQzATQWkRaiEhlAL0BTAi4plKVFCSx3Vu7dul6E2OA337zpz4iIj8EMf33VQAzALQRkQIR6WeM2QPgOgCTASwAMM4YM8/v2srKBsnatU5X1gEH6DU2SFas0BABgMLC+J9v40btIiMiyiRBzNrqk+D+JACTfC4nJfvtB+TkaADYVsYZZ+i2KQsXaheWDZZly5yPixcke/cCRx8NbN8OLF6sJzISEWWCsHRtZaScHGD//fXxggV6bdRIz3YHgJ9/dl67dKnzOF6QzJ6t4fPrr8DKlcXfT0QUVgySFNkJYXablPr1gSaRaQMFBc7rSmuRfPhh/NcSEYUdgyRFdpzEHSSNG+tjd5CU1iJhkBBRpmKQpMgGiQ2H8gTJ+vXAN984zxkkRJRJGCQpskFiJQqSkrq2pk4FiooAkeKvtYqKgHPPBa66KvWaiYi8xCBJUTJBsmGDvlmxQWK7tc48U6/xguTXX4F33wVGjoz+XEREQWOQpMgdJDk5QJ06xYPEBkONGnp1B0lRkRMkV18d/Xo390yub79NvW4iIq8wSFLkDpJ69TRMGjQAcnN1tfvOnU4wHHmkXt1BMmuWrjdp0gQ4/XRdP7J6ta4ncXMHiXs8hYgoaAySFLn3g7QHXeXmAgceqI9XrXIG2jt21Pdt3qxbpgDAO+/o9dxzgUqVgKZN9fny5dFfxx0kX3/t6bdARJQSBkmK3C2S+q4jt9zdWzZIDjpIu74AnakFRAcJADRvrtdly4BHHgGGD9fnsS0Sk3BvZCIif/m+RUq2KS1IVq50urZatNCV8GvXavfWxo26In6//YCuXZ3XAMDbbwPPP68zua64IjpI1q/XbVRat07f90VElCy2SFKUTJDMi2w/2bKls6VKYaHOwgKAnj2dvbVskLz4ol6N0a1TbJA0ihz35R4nefddYNw4YNMmb74nIqKyYJCkqHp1oEoVfRwvSCZN0u6tBg2AQw6JDpLYbi3ACZK9e5178+c7QXL++Xq14ySzZ+vHX3yxjtc88IB33xsRUTIYJCkScVol8YJk2jS9nn129CaPBQU6jTcnBzjtNOfjbJC4zZmjZ57k5QHnnKP3pk/X66hRem3QQAfwn33Wm++LiChZDBIPlBQk1tln69UGyaefaqujbVugZk3nde4g6dlTrx99pNdGjYAuXXQ9ynff6fqTV1/V99kjfP/4I7o1Q0SUbgwSD+Tna2uhQwfnnjtIqlYFTjlFH9sg+fhjvXbqFP256tfXcDnoIODuu/WeHWNp0kS70m6/XZ9fdJGucu/YETjqKA20vXt1MJ+IyC8MEg+MHKktATt1F9ADrezeWaeeCuyzjz62QWIHxmODRERbGz/8oMGU4/obstvT33yzrlPZvFmfX3GFXu1API/yJSI/MUg8YLdGccvLc05HtOMagBMkll3t7lalirY8qlbVmV6WDZJ99gEefFAfV64MXHKJPraLIBkkROQnriNJo3POAT74IHpWljtIRIAjjij5c7Rtq2tGACdIAOCyy3RacKtWzue0LZJVq1KvnYgoWQySNHr6aV0HYru4gOggOeQQZyPHRNq1cwbS3UGSmws89FD0a9kiIaIgsGsrzdwhAkQHSez4SDxt2zqP3UESD1skROS2YoU/X4dB4jN3kMQbH4nVrp3zONkgYYuEiBYtAtq0AS69NP1LAjIySESkrYg8IyLjReTqoOspi8qVne6sZFskVaroYH7sQH2sRF1b69YBY8fq6Yqffqr3jAH+97/iuwwTUeYzBrjuOj3GIi9Pu8LTyfcxEhEZBaAngDXGmMNc93sAGAYgF8DzxpiHEnwKGGMWABggIjkAnktzyZ478URg7lxdf1KaGjV0wL5ateLdZLHidW0VFGgYbdmiz7/+WqcWT5gA/PWvwAUXAG+8Ub7vg4jCafx4XchcuzYwZEj6v14Qg+2jAQwH8JK9ISK5AEYAOBVAAYCZIjIBGiqDYz7+SmPMGhE5B8DAyOfKKBMmaFPTbtRYmpNOSu51devq51y/Xg/GqlZNQ2PLFu0WW70a+PFH4M8/gU8+0Y9ZuLB83wMRhdPmzcCNN+rjwYOjz0xKF9+7towx0wCsj7ndGcBiY8xSY8wuAK8B6GWM+dEY0zPmbU3k80wwxhwL4FJ/v4PU5eQkHyJl/bx27crvv+v1jz/02r27dqUZo62SL77Q++7t6Yko8913n/ZKdO4M/O1v/nzNsIyRNALg/i+tIHIvLhHpJiJPiMizACYleE1/EZklIrPWVqA9Q2LHSWyQNGwIHHecPv7gA22pANo6sd1eRJTZ5s4Fhg3TXyqffjr9YyNWWIIkXu9/wjMAjTGfGWOuN8ZcZYwZkeA1I40x+caY/Hp+tO1CInacxB0kxx+vj194ASgqcj6moMC/+ogoPYqKgKuv1m7za69NblaoV8ISJAUA3JNbGwPgaohyiJ0CbLu4GjYEjj1WH2/dGv0x7N4iynyjR+vxEg0a+H8uUViCZCaA1iLSQkQqA+gNYELANWWkRF1bBxygg25t2jivtdOJGSREma2wELjtNn382GNArVr+fn3fg0REXgUwA0AbESkQkX7GmD0ArgMwGcACAOOMMfP8ri0blNS1BTjjJIBO/QXYtUWU6QYO1DA5+WSgTx//v77v03+NMXG/TWPMJCQYOKfklTTYDmiQjBqlK+YPP1zvJWqRFBbqPPScsLRbiaiYGTOA55/XmaAjRpS+3iwd+F9ElmnaVK/Llul88q1bdT2JPYXxgguA887TKYJ2y5V4QTJlinaFPfqoL2UTUTns2aMD7ABw6626EWwQGCRZpnlz3YaloMDZfr5hQ+e3lH33Bd56C7jwQucUx3hBMmqUrjl5801fyiaichgxQqfyN2sG3HVXcHUwSLJMpUpA69b6+PPP9WoXKcayLZLYMZJdu4BJkU7G777TVfJEFC6rVjnHcT/5pHMKaxAYJFnINm/tBo12fCRWnTra7bVpk3P0LwB89pnzfPduYPbsxF9r3Tqdarg+dq8CIkqrm2/W7utevYCzzw62FgZJFrJnmEybptdEQSISv3vr3Xf1ardxmT498df65z+Be+7xZ2M4IlJTpgCvv66/CA4bFnQ1DJKsZFskf/6p10RBAhQfcDfGCZIBA/Q6Y0b8jy0s1H/MAPDll+Wvl4iSt3OnrlwH9Je4Zs2CrQdgkGSl2JkbicZIgOLjJDNn6tThRo2Af/xD702frgETa/Ro/UdtP27nTt1heOhQYOPGlL4FIkrgkUf00Kq2bbV7KwwYJFnIvXodKLlF4u7aKioCbrlFn194IdCqlW5Nv2YNsHRp9McVFQHPPquP8/I0RL77TlfX3nKLnnUSL3yIqPyWLAEGDdLHTz2lMzTDgEGShWrUcAICSK5r67PPgIcf1u3lGzTQ2SAizv5cseMkn36qvxU1bgz07av3pk7VqcUA8N57wDPPePLtEBH0F7Prr9df2vr2Bbp1C7oiB4MkS9kBd6DkIOnWTX+rmTYNuPNOvTdihM7oApwtVV55JbqFYacHX345cMIJ+vjRR3VLeruH1803a9gQUereeUd/7mrVCt/kFgZJlnKPkzRokPh1bdoA336rB18BwEUXAeef77z/ssv0H+6HH0YvTpw7V6+dOzutFjsucs89uoJ+xw6nhUJE5bdli7ZGAO3aKumXwyAwSLKUDRJ7/G5JDj9cpxP+8gvw8svR72vYEHjoIX18/fVOWPz4o17bt9cFkHXr6vOcHA0j2+y2q+uJqPz+/W+dENOpkzObMkwYJFnKBklJM7bcRDQQKsXZxrN/f+CYY/Rsk8ce08H31at1/65mzaLHUrp31/Bp1UqfM0iIUvPTT8Djj+vPmZ+nHpYFgyRLde0K9OsH/OtfqX+unBxnK4YPPnBaI4cd5uwMfNllQNWqukARYJAQecEY4JprdHPGAQOAo44KuqL4fN9GnvyRl6dbS3vlhBO0tTJ7ts7sAoAOHZz3n39+9NhKs2b6+oIC3aurWjXvaiGqKF56SX/e6td3pv2GEVsklJTq1YEuXaLXj7Rvn/j1lSrpTsSAbmlPRGWzfr1uDQ/ojMjatYOtpyQMEkqandllD8sqKUgAdm8RpeKuu4C1a7U3wK7VCisGCSXt5JOjnzNIiNLj22+15V+pkq5gD+LUw7JgkFDSunRxxjoaNy69qV1SkEyZAlx8MbBhg7c1EmW6vXv11ENjdFHvoYcGXVHpGCSUtMqVdTYYED3QnkhJQTJkCDBuHDBmjHf1EWWDp58G5szR7YvuuSfoapLDIKEysQfoHH986a8tKUjsvalTvamLKBv88YdzZO4TT+gkl0zAIKEyueYa/c8/me2rmzfXdSYrVgCzZulvV1u36lG+K1boaz77TE9htMaMAc46C9i2LR3VE4XbLbfo6aRnnaUnH2YKBgmVSU6Ozt6qUqX011apAjRtqlOGjztOj+R9+WVg+XK9B2iwfPONPt62DbjhBt2Yzu42PHmy/oa2d29avh2i0PjkE/35qFpVWyNhH2B3y8ggEZFuIvKFiDwjIt2CrocSs91bu3bp9fvvi3d1ffyxXl97zdnLy7ZI7r4b+M9/dEU9UbbatUtb+4DuRtGyZbD1lJXvQSIio0RkjYj8FHO/h4gsFJHFIjKwlE9jAGwBUBVAQbpqpdR16qTXjh31+sMPTpDYM1PsOIn7/BIbJJs2Rb+GKBsNHQosXAgcfLBzuFwmCaJFMhpAD/cNEckFMALAGQDaAegjIu1EpL2ITIx5qw/gC2PMGQBuB3C/z/VTGdx7r3ZT2fNLfvxRdxkGgCuu0K6yr7/WGVwzZzofZ4PEXhkklK2WL9duX0DPAkqm2zhsfN9ryxgzTUSax9zuDGCxMWYpAIjIawB6GWMGA+hZwqfbACAD/9grjmrVdOdgQHcF/uMPJxSOPFLPM/n6a11TAmi/sDHFg2TePN19ONndjIkyxfXX6350vXsDp5wSdDXlE5YxkkYAVrqeF0TuxSUifxGRZwGMBTA8wWv6i8gsEZm1du1aT4ul8rFrTxYu1GurVrqH0Fln6Qyvpk2Bv/xF37d9u17ds7fYKqFsM2GCHktds6Ye0ZCpwhIk8eYnmDj39B3GvGWMucoYc7Ex5rMErxlpjMk3xuTXq1fPqzopBbGLGFu21NlcEyfqxo4rVgDt2un7tm3TmV02UAAGCWWXrVudUw8ffDCzW9thCZICAE1czxsDWBVQLZQmhx/uPG7UCNhnn+Kvsfe2bdOjet2mTIk+N54okw0apL88HXGEM2MrU4UlSGYCaC0iLUSkMoDeACYEXBN5zN0isdOCY9m9vLZtc7q1atfWc+d//91ZX0KUyRYs0G5de+phvJNJM0kQ039fBTADQBsRKRCRfsaYPQCuAzAZwAIA44wx8/yujdLrkEOc8+MTBYltkWzf7gRJ9ep62iMA3H67dnndeitw+ukMFso89tTD3buBv/9dN0PNdEHM2uqT4P4kAJN8Lod8VLky0LYtMHdu6UHibpHssw9w223AyJHAV1/pzJZPP9X3ffSR7pQ6YoT+djd9OrB5s4YMURi98opuDVS3ri62zQZh6dqiCsIejmWnBMdKFCS1aumaFEBDJC8P+NvfdM7900/r+Q07dgA9eugssN9+S+/3QVQef/4J/POf+viRR4D99w+2Hq+UGiQikisiQ/wohrLfoEHATz8BJ54Y//2JggQArrpKz2YQAcaOBZ57TsMEAN5/X/cq2rxZ9+XiDC8Ko3/9C1i9WmcrXn550NV4p9QgMcbsBdBJJJO2EKOwqlat5IN64g222yDJywNmzACWLnUWMJ51ll7ff1/n41tTpnhbN1GqZs/W0w5zc7UVnZNF/UHJjpF8B+BdEXkDwFZ70xjzVlqqogor3mC7e5pwzZr6ZnXrpuEzZ45uNWFNnaqDmvz1h8LAferhjTeWfkx1pkk2E+sAKARwMoCzI28lbV1CVC4ldW3FU62aM+6yfj1w4IH6tnq17utFFAbPPad7yTVq5Iz1ZZOkWiTGmCvSXQgRUPYgAYCePXV1PKAnOG7fDrz0knZvJXMkMFE6rV4N3HGHPv7vf6Nb1NkiqRaJiDQWkbcj27+vFpE3RaRxuoujiqekMZJEzjzTeXzOOcCpp+pjDrhTGNx2m87W6tEDOP/8oKtJj2S7tl6ErjQ/ELqZ4nuRe0SeKm2MJJ4mTXTn1KOOAk4+2dlB9fPPo/fqIvLb559r67hKFeDJJ7N3zC7ZIKlnjHnRGLMn8jYaAHdCJM+5WyRbI9M6SgsSAHj1VV1LUrWqblffqZOGyOTJ6auVqCS7dzt7aN1xR+JFuNkg2SBZJyJ9I2tKckWkL3TwnchTlSrpCnhjgA0b9F4yQRLrggv0On68d7URlcXjjwPz52uA3H570NWkV7JBciWAiwD8AeB3ABdE7hF5zgbHunXRz8vCBsmECcV3EXabPds5sZHIK7/+CtwfObt1xAhtKWezpFa2AzjfGHOOMaaeMaa+MeZcY8wKH+qjCsh2bxVG2rzlCZJWrXR77s2bEy9OLCwEjj9eN82zrR8iL9x4o3bPXnghcNppQVeTfsmubO/lQy1EALxpkQD6QwwAb7wR//2zZ2trZcMGYAg3ASKPvP8+8PbbQI0a2r1VESTbtfWViAwXka4icqR9S2tlVGF5FSTu7q2iouLv/+475/GwYXqePFEqtm0D/vEPfXz//boAsSJIdouUYyPXf7vuGehKdyJPeRUkBx8M7LsvsHGjzuOvUyf6/XPm6HXffYFNm/S40+HDy/e1iABg8GA9NrpDB+cY3YogmTGSHABPG2NOinljiFBa2DESO0he3iABnG26168v/j7bInniCb2+8gqP8qXyW7hQt4YHdHPGTD/1sCySGSMpgp5eSOSL2OBIJUhsK8QO3M+Zo3twbdoELFqkU4379NH+7A0bnNcBGiovv+y0XIgSMQa47jpg1y7gyit1m/iKJNkxkikicouINBGROvYtrZVRheVlkNgWSWEhsHOn7hZ87LHAxx/r/cMO0zBp3VqfL1rkfOyECUDfvvofA1FJXn9dt+SpUwd4+OGgq/FfWdaRXAtgGoDZkbdZ6SqKKrZ0BMn69cCaNTodeMsW4IYb9H7Hjno9+GC92iApKnJ2aV24kF1elNjGjcBNN+njhx/WI3QrmmR3/22R7kKIrHS1SNasce6vXKnXIyNzD22LxC5OfPtt4Icf9PGOHTqj64ADyl8HZa9779V/H126VNzWa4ktEhG5zfX4wpj3Zcmx9RQ2drDd8mqMxB0klm2RuLu2ioqA++7T53aTvWXLyl8DZa/vv9fNGHNysu/Uw7Io7dvu7Xp8R8z7enhcCxGA6ODIydExjPKK7doCgAYNnM9tzytxd23NnKnnyh94INArshSXQUKxior01MOiIl07csQRQVcUnNK6tiTB43jPfSMiXQFcCq2/nTHm2FI+hDKIO0j22Se1rbfdXVtr1+rj3r11Z9aGDYHq1fWeu0VizzE55xygVi19vHRp+Wug7PTCC8DXX2uX57//Xfrrs1lpQWISPI73PCkiMgp6TO8aY8xhrvs9AAwDkAvgeWPMQwmLMuYLAF+IyLkAZpanDgqv2CBJhe3acrdIGjYEBg6Mft3++wO1a+sU4Jdf1nvduzvTgdkiIbe1a50dfR9/XBe1VmSldW0dLiKbRGQzgA6Rx/Z5eY+vH42YbrHIxpAjAJwBoB2APiLSTkTai8jEmLf6rg+9BMCr5ayDQsrLIIk32F4vwUk6tntrwQJtBZ10EtAiMs2EQUJuAwfqLx2nnAJcdFHQ1QSvxBaJMSbX6y9ojJkmIs1jbncGsNgYsxQAROQ1AL2MMYOhrZdiRKQpgI3GmE1e10jBcg+2exkktmurfv34r23dGvjmG318xBH6sQwSivXVV8CoUTp2N2JE9p56WBZhmWPQCMBK1/OCyL2S9EMJx/2KSH8RmSUis9ba/0EoI6S7a6ukILG6d9dr06b6H8XKlboGpX9/4M03ndctX87jfCuS3buBAQP08e23O63Yii4sQRIv00scgzHG3GuMmV7C+0caY/KNMfn1EvVlUCh5GSS1aunsrE2bgFWr9F5pXVuAEyRVqugOrkVFuo/Sc885CxXnzwcOOgi44orUaqTM8eSTOqOvRQs9PpdUWIKkAEAT1/PGAFYFVAsFzMsgyclxWiU2SEprkeTlAV27Ovdt95bd3HHxYg2W6dP1On68021G2augwPklYvjw4uudKrKwBMlMAK1FpIWIVIauX5kQcE0UEC/HSIDo7eOrVXOm/Mbq0AE4+2zg1lujX2ODZFNkNG7nTu3qWrhQn+/dq3stUXa76SbdXue884Azzwy6mnDxPUhE5FUAMwC0EZECEelnjNkD3WF4MoAFAMYZY+b5XRuFg5ctEsAZcAe0WyvR4Ghenm7UOGhQ9P0WcTYIWrQo+qx3O2XYKizk/lzZ5MMPteVZvboegkbRfA8SY0wfY8wBxpg8Y0xjY8wLkfuTjDEHG2MOMsYMKu3zUPZKZ5Ak6tYqiTtI7H5bixY5LRJAF6YtWaKPZ83Sr8M+9OywY4duEQ9o11aTJiW/viIKS9cW0f+X7hZJWbVsqdf99tNZW4CuNbHBcf75erWtkvnzdexk4sTy1Uvh8vDD+nfdrh1w441BVxNODBIKHa+DxD1GUp4WybHHAtdco5vyHRbZi2HKFGDPHv3t9JJL9N70yBzCbdv0On++M65i7d5d9q9PwVm8WI/PBfTvPy8v2HrCikFCoeP1YHuqXVu5ubrwrHdvZ2bXzz/rtU0b5/yJrVv1aoPEGGD2bOfzjBmjfeyPPVb2Gsh/9tTDnTuByy8HTjgh6IrCi0FCoZOX55x3HYauLbdWraKfH3ywU6MNEHsFnJXy770H9OunLZL7749/hjyFy5tvApMna5emPYud4mOQUCjZ/5zD0LXlVr26bi9vtWnjTBWObZEAGiRz5+p+THv3ai2bNgFDh+rrJ03S33gpXDZvdsZDBg9O/d9NtmOQUCh5GSSpdm3Fcm+lEq9FYgMF0CC5806d+XPZZc4A/LBhQNu2wFlnAX//e+o1kbfuuw/47Tegc2f+/SSDQUKhZMdJwta1BUQHSWktkt9/B95/X7+PRx8FjjlGF7Nt3eoc9zt2LDBtWup1kTfmztWgz8kBnnpKx8ioZEmd2U7kt7B2bQFOkFSpops62q6peGMk1jXXOCE2bJiOAZ1xhm7b8sADwLXXAnPmcFZQ0Oyph3v36kB7p05BV5QZGCQUSjZAvNjPKF0tklat9LdVW+OOHfofkA2SQw8F5s3T9996q/PxrVoB776rj7dvB/73P90I8MordYppjRqp10jlM3q0TuNu0AB48MGgq8kc7NqiUOrbF+jSxZvfCKtXBy69VMcovAimU04BTj5Zz+kGdMsVG3zbtztB8te/AlWrAvfck7glVK2a7ihctaoGSseOwIoVqddIZVdYCNx2mz5+7DHnmGUqHYOEQun664EZM4CaNb35fP/7n67j8ELNmsDHHwNXXeXcs+Mk27Y5QXLssToWEnusb6zu3XVblfbtdQGcXWfy88/A//2f7jpL6TdwoIbJyScDffoEXU1mYZAQecC2SLZudQbd99lHB2yTceihwH//q4+//FKvDzyg4We3Lqf0mTEDeP55HaPiqYdlxyAh8kC8FklZJwocfbQOwn//va41mTpV77/2GrBxo3e1UrQ9e3SAHdCxrEMOCbaeTMTmsiG2AAASd0lEQVQgIfKAu0VS3iCpXh048kidOfTcc87RwNu2Aa+84l2tFG34cOCHH4DmzYG77gq6mszEICHygBctEgA4/ni92i05GjTQ68iROlX4m294zomXfvsNuPtuffzkk95MN6+IGCREHvCiRQI4R/za1sgDD+j05e+/17Pju3QB3n479XpJ3XyznnrYqxfQs2fQ1WQuBgmRB2xobNmiU4CB8k01Pu646OdnnQUMGBB9j6vgvfHRR8C4cfp3x1MPU8MgIfKA7doqLNRrtWrJz9hyq1dPt10B9CClAw/UfZ9++EF3owV0Cw9KzY4dupsAoOt8mjULtp5MxyAh8oBtkaxdG/28PGz31imn6LVSJaBDB91AENAgcY+TjB6toWP37qLSPfKIrtlp1w646aagq8l8DBIiD9gWiRdBcuedenbJ7bdH32/UCKhdW1s9v//u3H/6aT3694MPyv81K5IlS4D//EcfP/UUULlysPVkAwYJkQdscKxbF/28PFq00MVx7nNPAF0k16GDPrbdW7t3a7cXACxfXv6vWVG4Tz3861+BE08MuqLswCAh8oCXLZKSxAbJvHnO7sMMktK99Rbw4Yd66uGQIUFXkz0YJEQeiG2R2GDxWmyQuM+EX7YsPV8zW2zeDNxwgz7+z3+cNTqUuowLEhFpJyLjRORpEbkg6HqIgOBaJLNmOe9ji6Rk99+vCxDz84H+/YOuJrv4GiQiMkpE1ojITzH3e4jIQhFZLCKl7JWKMwA8aYy5GsBlaSuWqAy8HCMpyaGH6ljJggXArl3RLZI//nDWsFC0H3/UTTFFdHICTz30lt8tktEAerhviEgugBHQgGgHoE+k1dFeRCbGvNUHMBZAbxEZAmB/EIWAbZHs3q3XdAVJ9ep6MNaePdoqsQPt9rwTnmVSnPvUw2uu0RYJecvXIDHGTAOwPuZ2ZwCLjTFLjTG7ALwGoJcx5kdjTM+YtzWRt2sBDASwzs/6iRKJDY507tlku7fuvFNbJa1bA4cdpvfidW/t2QNMnAhcfDFw7rm6GK8iGTMG+OornnqYTmE4arcRAPdSqgIARyd6sYg0B3AngOoAEs67EJH+APoDQNOmTT0okyix2MH1dAbJDTcAEyYAU6bo8/x85+vFC5J+/YCXXnKeT5sGnHZa+uoLk8JC55jjoUN1thZ5LwyD7fGOkEm4v6kxZrkxpr8x5lJjzJclvG6kMSbfGJNfz4uDuolKEBsc6Zq1BejK97FjncOXOnXSLdCB4kGydy/wzjv6uF07vS5Zkr7awubOOzVMunUDLrkk6GqyVxiCpABAE9fzxgBWBVQLUbn42bUFaDfV888DxxwDXHSREySxU4DtIVktW+qZ9QCwaFH0a4qKsnNr+q+/1u338/J0BTtPPUyfMATJTACtRaSFiFQG0BvAhIBrIioTP7u2rCuvBKZPB5o00dXwQPEWyeef6/XEE3WQHtA9ptwuvhho2lR3Ls4W7lMPb7kFaNs22Hqynd/Tf18FMANAGxEpEJF+xpg9AK4DMBnAAgDjjDHz/KyLKFV+t0hiJeraKi1Ili0Dxo8HCgqAX35Jd5X+eeopbY01awb8619BV5P9fB1sN8b0SXB/EoBJftZC5KWgg+SAA3TzwTVr9HCt6tW1y+qLL/T9J54I1K2rj5cs0bGT3Fw9j8Oyiykz3apVTng88QRPPfRDGLq2iDJeXp6+WX7/55WT45ypYdeSzJ0LbNig95s3B2rUABo21CnDv/2mr3ntNedzZEuQ/POfuh3KOefoG6Ufg4TII+5xknTO2kokdsDd3a1lubu3fv5Zu3+sbAiSqVM1HKtV46mHfmKQEHnE3QoJojvFLpeyB1zZbq1u3ZzXuIPk9df1sT3JcV2GL+/dudM59fDuu51gpfRjkBB5xN0KCUOQ/PyzXo84wnmNDZJfftG1KADQq5deM71FMmSIfl+HHKLdW+QfBgmRR4JukTSJrMb69VddF2LHSty/mdsgGTNGB92bNdPpv0BmB8nSpcCgQfqYpx76j0FC5JGgWyQ2SFauBNav13Uh++4bvS2IDRLbjXXddToADzhBsmdPZq0pMQb4xz90D7FLLwVOOinoiioeBgmRR4Jukbi7tmxrpFmz6BXdBx3kPK5eXffhsjsI2SA59VQNnFWl7C/x/vvOfl9BeucdYNIkoFYt4NFHg66mYmKQEHkk6FlbjRvrdeVKZ+aWnRJs7befs57k8suB2rWdIFm3TlsjX3wBrF6t+1QlsmqVTq097zz9mKBs2eKcejhokNO6In8xSIg84m6FVKsWzNfff389E+Xbb/VevJlL3boBNWsCN96oz+vU0VbL+vW6Mn7vXr0/Zgwwc2b8r/X227rgcevW0lsu6fTAAxqcnToBAwYEV0dFxyAh8ohthVSpEtwJfLZ7y079jW2RAMArr2jXV+vW+jw3VwMIiD66F9Cwibeh45tvOo+DOiv+p5+Axx7jqYdhwCAh8ohtkQS5JYcdcLeBEK9FkpenXVputnvLtmR69dKWyvTpxTd5XLvWWewIBBMkxuhph3v2aEvkqKP8r4EcDBIij4QpSOyRv/FaJPHYcRMbJIcf7ixk/Oqr6Ne++652a1nxDtNKt7FjtdVVv74z7ZeCwyAh8ojt2goySGIPA012dbdtkcyZo9cWLYBjj9XH06dHv9Z2a3Xpole/WyQbNujW8IDO0optXZH/GCREHglTiwTQAX/b0iiNDZLt2/XqDhJ3i2THDuDjj3Vcwg7W+x0kd96p3WsnnAD07evv16b4wnBmO1FWsC2SIKb+Wu4gad48+VMBY0+jbtECaNBAJw7Mn6+tgNq1dYB79249KMqOS/jZtTVrFvDss0ClSjz1MEzYIiHySBhaJO6urWTHR4DoIKlUCWjUSEMkP1/vzZihV7tbcMeOGlo5OXoo1q5dqdWdrIEDdaD9xhuBQw/152tS6RgkRB6pU0ev7i1J/Hbggc5uvmXZ/dYdJE2bOlNpjztOr3ac5Lvv9HrEETr7q3Fj/Y/dbhS5axdw+unOMbdemjpVu9X226/kxZLkPwYJkUdOOQW4/37gnnuCq6FSJQ0ToGwtEvdYij3/HSg+TmKDpGPH6NfacZLJk4GPPgKeeQb48cey1V4SY5zwuO02DrCHDYOEyCNVqmiItG8fbB12nKS8LRJ3kBxzjF6/+QbYtk1PXQScreljg8R94uKIEcl//ZLs3g3cdJOusm/QALj+em8+L3mHQUKUZfr21fED94FWpUkUJPXrA0ceqbO5hgzRLVEaN3ZaMDasli/XoJkwwfnYsWOBP/8s5zcR8eefuonksGHa2ho+PNjJDBQfg4Qoy1xzjc6uKssGhom6tgDgb3/T60MP6dV2a7lfu2yZ7sC7ZYvO5ureXYNlzJiy12/t3avbwn/+OXDAAXq94ILyfz5KHwYJEaFKFT27BCgeJJdcojPRduzQ5/GCxH3iYu/ezpG3jz+ugVIe996r4VSnjo7R2PEaCp/QB4mItBSRF0RkfEn3iCg1HTtqmLRtG32/Vi3nFEX7Ost2bc2e7XRrXXghcPbZQIcOujlkeSYffPmlbn2Sk6Nny8eGG4VLWoNEREaJyBoR+Snmfg8RWSgii0VkYEmfwxiz1BjTr7R7RJSaDz/UI2tr1Sr+vv79ncfuM+AbNQLOPFN3D27RArj1Vh3sr1QJeOEFDYLHH0+8HX0ijzyi19tv19lwFG7pXtk+GsBwAC/ZGyKSC2AEgFMBFACYKSITAOQCGBzz8VcaY9akuUYiAlC1qr7Fc/TRehDW7t3R04pzcvSkxHjy83W21dChuoAwdvPHRBYuBN57T7vb7DYsFG5pDRJjzDQRaR5zuzOAxcaYpQAgIq8B6GWMGQygZzrrIaLyEQFGjy77x913n860mjEDKCx0zj0pyeOP6/Wyy3TWGIVfEGMkjQCsdD0viNyLS0T2F5FnAHQUkTsS3Yvzcf1FZJaIzFprD6MmIl/VqKGD5MYAn36q288PHRp9nonbypXOTK+bbvKvTkpNEJs2xttmLc4ZbJF3GFMIYEBp9+J83EgAIwEgPz8/4ecnovTq3l1D5JNPdOzkllt00aZd3GgtXarjITt26GB97KA/hVcQLZICAK49StEYQICnPhNROnXvrtePPwaef14fx+4YvHYt0LWrrkc56ijgxRd9LZFSFESLZCaA1iLSAsBvAHoDuCSAOojIB/n5QM2autZk0SK9t3kzsHGjM0Ns7Fhg1Sqgc2fdnLFmzeDqpbJL9/TfVwHMANBGRApEpJ8xZg+A6wBMBrAAwDhjzLx01kFEwalUCTjxRH1sXJ3MBQXOY7tH1623MkQyUbpnbfVJcH8SgEnp/NpEFB7duwMTJ+rjffcFNm3SIDn0UGDJEl1nUqOGrkmhzBP6le1ElPlOP12nELdsCfTqpffsGSbjxum1V69gDwWj8mOQEFHatW2rs7Y++shZ0Gi7tmy3lnsbFsosPLOdiHxht7Vv3FivK1fqAPzcuXrq4WmnBVYapYgtEiLylT14q6DAOQu+e3fdEoUyE4OEiHzlbpHMmaOPO3UKrh5KHYOEiHzlbpHEngFPmYlBQkS+2m8/nZ21eTPw7bd6j0GS2RgkROQrEadVsnMncOCBQIMGwdZEqWGQEJHv7DgJABx5ZHB1kDcYJETkuyaubVvZrZX5GCRE5Dt3kLBFkvkYJETkO3fXFlskmY9BQkS+sy2S2rWBpk2DrYVSxyAhIt/l52uIXHSRzuKizMa9tojId/XqAYWFQVdBXmGQEFEg2BLJHuzaIiKilDBIiIgoJQwSIiJKCYOEiIhSwiAhIqKUMEiIiCglDBIiIkqJGGOCriHtRGQtgBVB1xFRF8C6oIsoQZjrC3NtQLjrC3NtQLjrC3NtQHrra2aMqVfaiypEkISJiMwyxuQHXUciYa4vzLUB4a4vzLUB4a4vzLUB4aiPXVtERJQSBgkREaWEQeK/kUEXUIow1xfm2oBw1xfm2oBw1xfm2oAQ1McxEiIiSglbJERElBIGCRERpYRBQkREKWGQhICItBSRF0RkvOteNxH5QkSeEZFuIautbaSu8SJydVC1lVBfsXtBEpF2IjJORJ4WkQuCrieWiHSN/H0+LyLTg64nVlh+FuIJ089CPH79LDBIUiQio0RkjYj8FHO/h4gsFJHFIjKwpM9hjFlqjOkXexvAFgBVARSEqTZjzAJjzAAAFwEo90KoNNYX788zsBoBnAHgSWPM1QAu86IuL+szxnwR+fucCGBM2OqDBz8L6arNq5+FNNbn2c9CiYwxfEvhDcAJAI4E8JPrXi6AJQBaAqgM4AcA7QC0h/6wut/quz5uvOtxTuTaAMDLYaot8vwcANMBXBK2P7uS7gVRY+RtBIAhAL4K8b+/cQD2DVt9XvwspPPPzoufhSB+PjytNZ2fvKK8AWge85d9DIDJrud3ALgjic8T7z/Dyqn8I0hnbZH774f4z86THx4Pa8wF8G4Y//0BaArgOa9r8/jPL6WfhXTWFnltSj8Laf6zS2uQsGsrPRoBWOl6XhC5F5eI7C8izwDoKCJ3RO79RUSeBTAWwPCQ1dZNRJ6I1DfJw9q8qq/YvYBrbC4iIwG8BG2VpFuZ6ovoB+DFtFUUrax/fun6WfCitnT+LMST8s9HOlRK1yeu4CTOvYQrP40xhQAGxNx7C8BbHtcFeFPbZwA+87Qqhxf1FbvnsbLWuBxA/7RVU1yZ6gMAY8y9aaolnrL++aXrZyGestb2GdL3sxBPyj8f6cAWSXoUAGjiet4YwKqAaokV5tqA8NcHhL9G1ld+Ya4NCGl9DJL0mAmgtYi0EJHKAHoDmBBwTVaYawPCXx8Q/hpZX/mFuTYgrPWlcwCmIrwBeBXA7wB2Q39b6Be5fyaAX6AzLO5ibZlXXybUyPqys7ZMqM/9xk0biYgoJezaIiKilDBIiIgoJQwSIiJKCYOEiIhSwiAhIqKUMEiIiCglDBKigIjIfSJyS9B1EKWKQUJERClhkBD5SETuihxKNBVAm8i9mfbkPxEZLCKDgqyRqKy4+y+RT0SkE3RvpI7Qn705AGYD+D8A40XkegA9ABwdVI1E5cEgIfJPVwBvG2O2AYCITAAAY8w8ERkL4D0AxxhjdgVYI1GZsWuLyF+JNrdrD+BP6HGyRBmFQULkn2kAzhORaiJSE8DZgJ4ACGB/6BndT4jIfgHWSFRm3P2XyEcicheAywCsgG4NPh96emJ3Y8zKyDhJJ2PM5QGWSVQmDBIiIkoJu7aIiCglDBIiIkoJg4SIiFLCICEiopQwSIiIKCUMEiIiSgmDhIiIUsIgISKilPw/Nx3qczvgL1cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1892b84f3c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def dericd(y,x,dx): #central difference method for a derivative\n",
    "    return (y(x+dx/2)-y(x-dx/2))/dx\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x=(np.pi/4) #set the x value to evaluate cos'(x) at.\n",
    "dxarr=[]\n",
    "yarr=[]\n",
    "for dx in 10**np.arange(-16., -0.5, 0.1): #create a loop to store 10^-16,10^-15,...,10^-1 in an array\n",
    "    y=dericd(np.cos,x,dx) \n",
    "    dxarr.append(dx) #stores the dx value used\n",
    "    yarr.append(y) #stores the numerically evaluated y value using the specific dx\n",
    "    \n",
    "plt.plot(dxarr, abs(yarr+np.sin(x)), 'b-', lw=2) #plots the error (defined as true value-numerically calculated value) vs the corresponding dx value\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('dx')\n",
    "plt.ylabel('Error')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "\n",
    "#This scheme is much better because not only is the optimal value for dx larger (thus making it more\n",
    "#manageable), you are also able to achieve much greater accuracy (the error value is smaller here)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\kevin\\Miniconda3\\lib\\site-packages\\matplotlib\\ticker.py:2206: UserWarning: Data has no positive values, and therefore cannot be log-scaled.\n",
      "  \"Data has no positive values, and therefore cannot be \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAADB1JREFUeJzt3VGIpfddxvHn1w2hRXSlTRRMoptACAaKxC6VCpWgN1shKba1JFKqNRhSDF4VSeiFIqgXvauNqatto6WkxpjSRANFL0JL24tsKkJiiMRiyJJqNimuWC9i6t+LOTHTyW92Zrrn3fc9k88HDjvn3Zndh5l998uZM+ecGmMEAHZ6w9wDAFgmgQCgJRAAtAQCgJZAANASCABaAgFASyAAaAkEAC2BAKB10dwDzscll1wyjh07NvcMgI3y2GOPvTDGuHSv99voQBw7diynTp2aewbARqmqZ/bzfr7FBEBLIABoCQQArcUEoqquqqpPVdX9c28BYOJAVNWnq+r5qnp8x/ETVfVUVT1dVXckyRjjm2OMW6bcA8D+TX0L4p4kJ7YfqKojSe5K8q4k1ya5uaqunXgHAAc0aSDGGF9O8u0dh9+e5OnVLYaXknw+ybun3AHAwc1xH8RlSZ7ddv10ksuq6i1V9ckk11XVnbt9cFXdWlWnqurUmTNnpt4K8Lo1xwPlqjk2xhgvJrltrw8eY5xMcjJJjh8/Pta8DYCVOW5BnE5yxbbrlyd5boYdAJzDHIF4NMnVVXVlVV2c5KYkD86wA4BzmPrHXO9N8vUk11TV6aq6ZYzxcpLbk3wpyZNJ7htjPDHlDgAObtL7IMYYN+9y/OEkD0/5dwNwfhbzSGoAlkUgAGgJBACtjQxEVd1QVSfPnj079xSAQ2sjAzHGeGiMcevRo0fnngJwaG1kIACYnkAA0BIIAFoCAUBLIABoCQQALYEAoCUQALQEAoDWRgbCU20ATG8jA+GpNgCmt5GBAGB6AgFASyAAaAkEAC2BAKAlEAC0BAKAlkAA0BIIAFobGQhPtQEwvY0MhKfaAJjeRgYCgOkJBAAtgQCgJRAAtAQCgJZAANASCABaAgFASyAAaAkEAC2BAKC1kYHwZH0A09vIQHiyPoDpbWQgAJieQADQEggAWgIBQEsgAGgJBAAtgQCgJRAAtAQCgJZAANASCABaAgFASyAAaAkEAK2NDITXgwCY3kYGwutBAExvIwMBwPQEAoCWQADQEggAWgIBQEsgAGgJBAAtgQCgJRAAtAQCgJZAANASCABaAgFASyAAaAkEAC2BAKAlEAC0NjIQXnIUYHobGQgvOQowvY0MBADTEwgAWgIBQEsgAGgJBAAtgQCgJRAAtAQCgNaegaiqI1X1sQsxBoDl2DMQY4zvJnlbVdUF2APAQly0z/f7hyRfrKq/SvKdVw6OMR6YZBUAs9tvIN6c5MUkP7/t2EgiEACH1L4CMcb40NRDAFiWff0UU1VdXlVfqKrnq+rfq+qvq+ryqccBMJ/9/pjrZ5I8mOTHklyW5KHVMQAOqf0G4tIxxmfGGC+vLvckuXTCXQDMbL+BeKGqPrB6TMSRqvpAtu60BuCQ2m8gfj3J+5P8W5JvJXnf6hgAh9SeP8VUVUeSvHeMceMF2APAQuz3kdTvvgBbAFiQ/T5Q7qtV9Ykkf5nvfST1NyZZBcDs9huIn139+nvbjo187yOrAThE9nMfxBuS3D3GuO8C7AFgIfZzH8T/Jrn9AmzZt6q6oapOnj17du4pAIfWfn/M9e+q6iNVdUVVvfmVy6TLzmGM8dAY49ajR4/ONQHg0NvvfRCvPObhN7cdG0muWu8cAJZiv8/meuXUQwBYlnN+i6mqfnvb27+84/f+YKpRAMxvr/sgbtr29p07fu/EmrcAsCB7BaJ2ebu7DsAhslcgxi5vd9cBOET2upP6p6rqP7N1a+FNq7ezuv7GSZcBMKtzBmKMceRCDQFgWfb7QDkAXmcEAoCWQADQEggAWgIBQEsgAGgJBAAtgQCgJRAAtAQCgJZAANASCABaAgFASyAAaAkEAC2BAKAlEAC0BAKAlkAA0BIIAFoCAUBLIABoCQQALYEAoCUQALQ2MhBVdUNVnTx79uzcUwAOrY0MxBjjoTHGrUePHp17CsChtZGBAGB6AgFASyAAaAkEAC2BAKAlEAC0BAKAlkAA0BIIAFoCAUBLIABoCQQALYEAoCUQALQEAoCWQADQEggAWgIBQEsgAGgJBAAtgQCgJRAAtAQCgJZAANASCABaAgFASyAAaAkEAC2BAKAlEAC0BAKAlkAA0BIIAFoCAUBLIABoCQQALYEAoCUQALQEAoCWQADQEggAWgIBQEsgAGgJBAAtgQCgJRAAtAQCgJZAANASCABaAgFASyAAaAkEAC2BAKB10dwDXlFVP5Dkj5O8lOSRMcbnZp4E8Lo26S2Iqvp0VT1fVY/vOH6iqp6qqqer6o7V4fckuX+M8RtJbpxyFwB7m/pbTPckObH9QFUdSXJXkncluTbJzVV1bZLLkzy7erfvTrwLgD1MGogxxpeTfHvH4bcneXqM8c0xxktJPp/k3UlOZysSk+8CYG9z/Ed8WV69pZBsheGyJA8keW9V3Z3kod0+uKpurapTVXXqzJkz0y4FeB2b407qao6NMcZ3knxorw8eY5xMcjJJjh8/Pta8DYCVOW5BnE5yxbbrlyd5boYdAJzDHIF4NMnVVXVlVV2c5KYkD86wA4BzmPrHXO9N8vUk11TV6aq6ZYzxcpLbk3wpyZNJ7htjPDHlDgAObtL7IMYYN+9y/OEkD0/5dwNwfvw4KQAtgQCgJRAAtDYyEFV1Q1WdPHv27NxTAA6tGmNzH2tWVWeSPDP3jpVLkrww94hzWPK+JW9Llr1vyduSZe9b8rZk2n0/Mca4dK932uhALElVnRpjHJ97x26WvG/J25Jl71vytmTZ+5a8LVnGvo38FhMA0xMIAFoCsT4n5x6whyXvW/K2ZNn7lrwtWfa+JW9LFrDPfRAAtNyCAKAlEAC0BAKAlkBMqKquqqpPVdX9245dX1VfqapPVtX1C9v2k6td91fVh+fado59rzk2p6q6tqruq6q7q+p9c+/Zqareufp6/llVfW3uPdst5TzYzZLOhZ0u5HkgELuoqk9X1fNV9fiO4yeq6qmqerqq7jjXnzHG+OYY45adh5P8V5I3ZuvV9RazbYzx5BjjtiTvT/J9P0Bnwn3d53O2jUneleSPxhgfTvLBdexa574xxldWX8+/SfLnS9qWNZwHU+5b17kw0ba1nQd7GmO4NJckP5fkp5M8vu3YkST/kuSqJBcn+cck1yZ5a7ZOwu2XH9n2cfdve/sNq19/NMnnlrRtdf3GJF9L8itL+9yd69gcG1eXu5J8LMlXF/zv774kP7Skbes4D6b+3K3jXLjQ58a6L5P+4Zt+SXJsxxfyHUm+tO36nUnu3Mef0/0nd/H5fIGn3LY6/rcL/tyt5cRY48YjSb64xH9/SX48yZ8ucdvq/c7rPJh63+p9z+tcmPBzN3kgfIvpYC5L8uy266dXx1pV9Zaq+mSS66rqztWx91TVnyT5bJJPLGzb9VX18dW+db/i3zr2vebYzBuPVdXJJH+RrVsRUzvQvpVbknxmskWvOujnbqrzYDcH3TfluXC+26Y+D/7fpC85eghVc2zXRxqOMV5MctuOYw8keWDNu5L1bHskySNrXfWqdex7zbE1O+jGf01y62RrXutA+5JkjPE7E23Z6aCfu6nOg90cdN8jme5c2Om8z42puAVxMKeTXLHt+uVJnptpy05L3pYsf1+y/I1L3rfkbcmy9y12m0AczKNJrq6qK6vq4iQ3JXlw5k2vWPK2ZPn7kuVvXPK+JW9Llr1vudumvpNjUy9J7k3yrST/k63C37I6/otJ/jlbP3XwUds2b98mbFzyviVvW/q+JW/rLp6sD4CWbzEB0BIIAFoCAUBLIABoCQQALYEAoCUQsGZV9btV9ZG5d8D5EggAWgIBa1BVH1294MvfJ7lmdezRV14trar+sKp+f86NcFCezRXOU1W9LVvPn3Ndts6pbyR5LMmvJbm/qn4ryYkkPzPXRvh+CAScv3cm+cIY47+TpKoeTJIxxhNV9dkkDyV5xxjjpRk3woH5FhOsx25PavbWJP+RrZfWhI0iEHD+vpzkl6rqTVX1g0luSLZeNS3JW7L1OsQfr6ofnnEjHJhnc4U1qKqPJvlgkmey9TTO/5StV5v7hTHGs6v7Id42xvjVGWfCgQgEAC3fYgKgJRAAtAQCgJZAANASCABaAgFASyAAaAkEAK3/A3fn04F9zKDqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x242b4fdc400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def dericd(y,x,dx): #central difference method for a derivative\n",
    "    return (y(x+dx/2)-y(x-dx/2))/dx\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x=(0.0) #set the x value to evaluate cos'(x) at.\n",
    "dxarr=[]\n",
    "yarr=[]\n",
    "for dx in 10**np.arange(-16., -0.5, 0.1): #create a loop to store 10^-16,10^-15,...,10^-1 in an array\n",
    "    y=dericd(np.cos,x,dx) \n",
    "    dxarr.append(dx) #stores the dx value used\n",
    "    yarr.append(y) #stores the numerically evaluated y value using the specific dx\n",
    "    \n",
    "plt.plot(dxarr, abs(yarr+np.sin(x)), 'b-', lw=2) #plots the error (defined as true value-numerically calculated value) vs the corresponding dx value\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('dx')\n",
    "plt.ylabel('Error')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "#The results differ because there is no round off or truncation error here. The derivative formula and\n",
    "#the true value of cos(0) will return 1 every single time, no matter what value of dx there is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}