{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 3 Part 2 - Non Linear Equations - Solutions\n",
    "\n",
    "\n",
    "## Exercise 1 - Relaxation method\n",
    "\n",
    "Solve the equations a) $x=2-e^{-x}$ and b) $x=e^{1-x^2}$. \n",
    "\n",
    "1. Write a short program that iterates equations a) and b).\n",
    "2. Start iterating a) from $x=1.0$ and b) from $x=0.5$.\n",
    "3. At every step, take the difference to the value at the previous step and plot this difference as a function of iteration number.\n",
    "\n",
    "\n",
    "## Talking points\n",
    "\n",
    "1. What do you observe?\n",
    "2. How many iterations do you need for 1E-3 accuracy?\n",
    "3. What is happening in case b)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x=  1.6321205588285577\n",
      "x=  1.8044854658474119\n",
      "x=  1.8354408939220457\n",
      "x=  1.8404568553435368\n",
      "x=  1.841255113911434\n",
      "x=  1.8413817828128696\n",
      "x=  1.8414018735357267\n",
      "x=  1.8414050598547234\n",
      "x=  1.8414055651879888\n",
      "x=  1.8414056453310121\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEUCAYAAADwYOuyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7DUlEQVR4nO3deXwUVbbA8d/pkEDCvoRFIAnIJCDoIKCAgGwuIKj4ZBGDCg4ThAcI4igYEVQi6LCKwBh4GIWMsjkiKosLm4oji3GUHTEJDKIimxAETO77ozoxS3fSId1dWc7386lPqHtvVZ1KQp/cqlu3xBiDUkop5Q8OuwNQSilVdmjSUUop5TeadJRSSvmNJh2llFJ+o0lHKaWU32jSUUop5Tfl7A6guKtVq5aJiIiwOwyllCoxdu7cecIYE+qqTpNOASIiItixY4fdYSilVIkhIinu6vTymlJKKb/RpKOUUspvNOkopZTyG006Siml/EaTjlKqRDl69CijRo2iffv2hISEICIkJyfnabdp0yZExO1y+vRprxxHFY4mHR9ITEwkIiICh8NBREQEiYmJdoekVKlx6NAhli9fTvXq1enUqVOB7V9++WW2bduWZ6lcubJXj6M8o0OmvSwxMZGYmBjS0tIASElJISYmBoDo6Gg7Q1OqVLj55pv58ccfAVi0aBEbNmzIt32zZs1o166dz4/jLRcvXqR8+fJ5yo0xXL58maCgIK/v25+0p+NlsbGxWQknU1paGiNGjGDFihUcPHiQjIwMm6JTquRzOPzzsVXU45w4cYLhw4dTv359ypcvT9OmTYmPj8/RJiEhARFhy5Yt9OvXj2rVqtG2bVvAekZw0KBBLF68mKZNmxIUFMT7778PwLp162jfvj3BwcFUrVqVPn36sH///hz77tKlCx07dmTNmjVcf/31lC9fnvnz5xfpnLxBezpelpqa6rL87Nmz9O/fH4BKlSrx5z//mZYtW3L99dfTsmVLmjdvToUKFfwZqlJlQkZGBr///nuOMhEhICDAZ8c8e/YsHTp04MKFC0yePJlGjRqxfv16hg8fzsWLFxk1alSO9tHR0QwcOJCVK1fmiHXjxo0kJSUxadIkateuTUREBOvWraNXr15069aNZcuWce7cOZ555hk6duxIUlIS9evXz9r+wIEDjB49mokTJ9K4cWNq1Kjhs3P2lCYdLwsLCyMlJe/DuA0bNuSdd94hKSmJpKQkvvrqK9544w3mzZsHQLly5WjWrBktW7bMsRSHXxKlSrLbb789T1nz5s359ttvfXbMOXPmkJKSwjfffMOf/vQnAG655RZOnz7Ns88+y/DhwylX7o+P3759+/LSSy/l2c+pU6fYuXMndevWzSobMGAAjRs3Zu3atVn7aN++PZGRkcyYMYOZM2dmtT1x4gQbNmygZcuWPjrTwtOk42VxcXE57ukAhISEMHXqVFq1akWrVq2yyjMyMvj+++/56quvspLRxx9/zJIlS7LahIWF5egRtWzZkvDwcETEr+elVEk1b948brzxxhxlwcHBWf/O3QsKCAgo8v+vdevW0bZtWxo1apRj/7fffjuLFi1iz549XHfddVnl99xzj8v9tGvXLkfCOX/+PLt27eKpp57KkbQaNWpEhw4d2Lx5c47tIyIiilXCAU06Xpc5WCA2NpbU1FTCwsKIi4tzOYjA4XBw9dVXc/XVV9O3b9+s8p9++omvv/46RzJ67733su4FVatWLSsBZSajZs2aERgY6J+TVKoEiYyMpE2bNm7rc/+/2bhxI126dCnSMX/66ScOHTrk9v/kL7/8kmO9Xr16LtvlLj916hTGGJft69atm+cqi7v92kmTjg9ER0cXaaRa7dq1ufXWW7n11luzytLS0vjmm2+yLs0lJSXx6quvcuHCBQCCgoJo0aJFjmR03XXXUaVKlSKfj1Kl2fbt23OsR0VFFXmfNWvWpHbt2syZM8dlfe5juOtZ5S6vXr06IsLx48fztD1+/Dg1a9b0aL920qRTQoSEhNC2bduskS0A6enpHDx4MEePaM2aNSxevDirzdVXX53j0tz1119PvXr1iuUvo1J2yK8XdKV69OjB3LlzCQsLo3bt2l7bb8WKFWndujUrVqxg8uTJWYMhUlJS+Pzzz/MMUCiONOmUYAEBATRt2pSmTZsycOBAwBrL/8MPP+ToESUlJbFy5cqs7UJDQ/PcJ4qMjPTpaB6lvCnz93nnzp0ArF27ltDQUEJDQ+ncuXOOtnv37qVSpUp59nHttddSsWJFrx0nu7Fjx7Js2TI6derE2LFjiYqK4vz58+zbt4+tW7eyevVqz082l+eff55evXrRu3dvRowYwblz55g0aRJVq1Zl3LhxV7xfvzHG6JLP0rp1a1ManDlzxmzdutXMnTvXPPzww6ZVq1YmKCjIAAYwwcHBpm3btmbYsGFmwYIF5osvvjDnz5+3O2ylXMr8vc29dO7cOavNxo0b3bYDzPbt271yHHdOnjxpxowZYyIiIkxgYKAJDQ01HTt2NLNmzcpq89prrxnAHDx4MM/24eHhJjo62uW+165da9q1a2cqVKhgqlSpYu666y6zb9++HG06d+5sOnToUGCcvgDsMG4+U8WqV+60adPGlNaXuF2+fJm9e/fmGMadlJSUNSeVw+EgMjIyz+W50FCXLwRUSikARGSnMcbldUtNOgUozUnHFWMMqampeS7PZR8Vc9VVV+W5PNe4cWO/PSmulCreNOkUQVlLOu6cPHmSr7/+Okcy2rNnD+np6QBUrlzZ5SwLds/zpJTyP006RaBJx73ffvuN3bt3Z/WGMpdz584B1iwL11xzTY5Lc3/+85+pXr26zZErpXxJk04RaNIpnIyMDA4fPpzn8tyxY8ey2oSHh+e5T9SwYUMdxq1UKaFJpwg06XjHjz/+mOfy3P79+8n8/atevXqeWRaaNm2qsywoVQJp0ikCTTq+c/78+axZFjKT0X/+8x9+++03AMqXL+9yloXsL99KTEz0aMohpZT/5Jd09OFQH3I1f1P//v0ZMWIEaWlp3HHHHXnqBw8ezODBgzlx4kSO+dgyDR8+nAEDBnDkyBEeeOCBPPXjxo3jzjvvZP/+/QwbNixP/dNPP80tt9xCUlISY8aMyVP/wgsvcNNNN/H555/z1FNP5amfPXs2LVu25KOPPmLKlCl56l999VWioqJYs2YNM2bMyFO/ZMkSGjZsyLJly1iwYEGOuuDgYL7//ntOnTrF7Nmz+eCDDzh8+DD/+c9/+L//+7+sdk2aNKFy5cocOXKEX375Jau3lP2FeT/88APvvfdenv2vXbsWsB6w+/jjj3PU16xZk1WrVgEwYcIEtm3blqO+QYMGLF26FIAxY8aQlJSUoz4yMjLrfSkxMTEcOHAgR33Lli2ZPXs2AIMGDeLo0aM56tu3b8/UqVMBuPfee/PMz9W9e3cmTpwIQM+ePbOmQMrUu3dvHn/8cUB/9wr7uwfWg6C1atUiISGBhISEPPUffPABISEhzJ8/n+XLl+ep37RpEwDTp08vFb97mefjbZp0VLGS+YqH9u3bZ72UyhjDpUuXOHfuHPfeey/ffvstmzdv5sSJE3m2T0tLIzY2lpEjR/o7dKWUB/TyWgH08lrx5XA4cPf7u2zZMvr06VOkV/sqpa5MfpfX9Gk+VWKFhYW5LA8ICGDAgAE0aNCAJ554goMHD/o5MqWUO5p0VIkVFxdHSEhIjrKQkBASEhJYt24dnTp1YubMmURGRtK9e3eWLVvGpUuXbIpWKQWadFQJFh0dTXx8fNabVMPDw4mPj2fQoEHcfvvtrFq1iiNHjhAXF8fhw4e57777tPejlM30nk4B9J5O6ZCRkcGHH35IfHw8q1evJj09na5duzJs2DD69Omj0/Uo5UV6T0eVeQ6HI0fv54UXXiA5OTmr9/O3v/1Nez9K+YEmHVXm1KtXjwkTJnDo0CHWr1/PzTffzOzZs4mMjKRbt2689dZbXLx40e4wlSqVNOmoMsvhcHDbbbexatUqUlNTs3o/AwcOzOr95H7ITilVNJp0lCJv76dz587Mnj2bqKgo7f0o5UWadJTKJrP3s3Llyhz3frT3o5R3aNJRyo26deu67f107dpVez9KXQFNOkoVwFXvJyUlhYEDB1K/fn0ef/xx7f0o5SFNOkoVQvbez4YNG+jatStz5szJ6v28+eab2vtRKh+adJS6Ag6Hg1tvvZUVK1Zw5MgRpk6dSmpqKvfff39W7ydzlmyl1B806ShVRHXr1mX8+PEcPHgwR++nadOmdOnSRXs/SmWjSUcpL8ne+zl69ChTp07lyJEjWb2fcePGae9HlXmadJTygTp16mT1fj788EO6du3Kyy+/nNX7+ec//5n1Wm6lyhJNOkr5kMPh4JZbbsnq/UybNo0jR44QHR1NgwYNGDduHPv27bM7TKX8RpOOUn5Sp04dnnzyyazeT7du3Xj55Zdp1qxZnt5PYmIiEREROBwOIiIiSExMtDl6pbyjUK82EJG6QBhQIXedMWaLF+MqNvTVBsqXfvzxRxISEoiPj+fw4cPUrFmTG2+8kY0bN+a4/BYSEkJ8fDzR0dE2RquUZ/J7tYFHSUdE6gNLgZtdVQPGGBNQpCiLKU06yh8yMjL45JNPiI+PZ8WKFS7bhIeHk5yc7N/AlLoC3kg67wLtgWnAN0Ce8Z/GmM1FjNNnRKQC8BYQhRX7j8BwY8zhgrbVpKP8zeFw4Or/pYiQkZFhQ0RKFU5+Saech/voBIw2xizxXlh+t8AYsx5AREYCi4Bu9oakVF5hYWGkpKTkKa9cuTJnz56lSpUqNkSllHd4OpDgAvCTtw4qIg1EZK6IbBORNBExIhLhpm1DEVkpImdE5KyIvC0iYYU5njHmt8yE4/QF0LgIp6CUz8TFxRESEpKjLCAggLNnzxIVFcXSpUtd9oSUKgk8TToLgQe8eNwmQH/gFLDVXSMRCQE+AZoCDzlj+BOwUUQqFuH4o4DVRdheKZ+Jjo4mPj6e8PBwRITw8HBef/11vvzyS8LCwnjggQfo1KkTSUlJdoeqVKF5ek8nBhgPJAMfACdztzHGLPb4oCIOY0yG899DsZJaI2NMcq52jwIzgShjzCFnWSPgIPCEMWams2wX1qg6V643xhzJts8JwF1Ad2NMWkGx6j0dVZxkZGSQkJDA+PHj+eWXXxg2bBhTpkyhRo0adoemVBZvDCQo6O7lFY9eKyDpfAxUMMZ0yFW+2XnQzoU81uPAfcAtxpjTnmyjSUcVR6dPn2bSpEnMmzePatWqERcXx9ChQwkIKJWDSFUJk1/S8fTyWqMCFl/dH2kOfOuifDdwTWF2JCKPAQOBWz1NOEoVV9WqVWPOnDl89dVXtGjRgkceeYQbb7yRbdu22R2aUvnyKOkYY1IKWnwUXw2s+z65nQSqe7oTEWkAzACqYd0PShIRt90XEYkRkR0isuPnn38uZMhK+c+1117Lxo0befPNNzl+/Dg33XQTgwcP5vjx43aHppRLhZoGR0RaiMj/ishEERkhIi18FVg2rq7/SaF2YMxRY4wYY642xrR0Li67fs728caYNsaYNqGhoYUOWCl/EhHuu+8+9u/fz/jx4/nnP/9JVFQUs2bN4vLly3aHp1QOHiUdESknIkuBr4G5wLPAK8DXIrJERHx1IfkUVm8nt+q47gEpVWZVqlSJqVOn8u2333LTTTfx2GOP0bJlSz755BO7Q1Mqi6c9nUlYQ5yfwbqHE+z8+gwwwPnVF3Zj3dfJ7Rpgj4+OqVSJFhkZyQcffMDq1au5cOEC3bt3p3///qSmptodmlIeJ51BwPPGmDjnPZyLzq9xwBTgQR/F9y7QTkSyBio4HyLt4KxTSrkgItx1113s3r2b5557jjVr1tCsWTPi4uL0PT7KVp4mnasAd8NiPnfWF4qI9BWRvkBrZ1FPZ1n2YdALsZ4NWi0id4vIXVgPdR4BXi3sMZUqa4KDg5k4cSL79u2jZ8+ePP3007Ro0YL333/f7tBUGeVp0jmG1btw5SZnfWGtcC6PONfnO9efzWxgjDmPNT/aAWAJkAh8D3Qzxpy7gmMqVSaFh4ezcuVKNmzYQGBgIL1796Z3794cOnTI7tBUGeNp0kkEYp2j1hqLSLCINHI+3R+LlRAKxTmazNXSJVe7VGPMvcaYKsaYysaYPrkfIlVKeebWW2/l66+/Zvr06WzZsoXmzZsTGxvL+fPn7Q5NlRGeJp3JwEqsXshB4BxwCIjLVq6UKgGCgoIYN24c+/fvZ8CAAbzwwgs0bdqU5cuXl4iJRI8ePcqoUaNo3749ISEhiIjL9wxt2rQJEXG7nD59Ot/jrF+/nm7dulG3bl3Kly9PgwYN6N+/P3v26BimovD04dDfjTH3A9cCI7FGq40EWhhjoo0x6T6MUSnlA/Xq1eONN97g008/pVatWgwYMIDu3buze/duu0PL16FDh1i+fDnVq1enU6dOBbZ/+eWX2bZtW56lcuXK+W538uRJWrduzSuvvMKGDRuYOnUqu3fvpl27di5fPaE8ZIzRJZ+ldevWRqnS7vfffzfz58831atXNwEBAWbMmDHm9OnTdoflUnp6eta/Fy5caADz/fff52m3ceNGA5gPP/zQa8fet2+fAcz06dO9ts/cfvvtN5flGRkZ5uLFiz7Zt7cBO4ybz1S3PR0RCRORwGz/znfxT4pUSvlCQEAAw4cP58CBAwwdOpQ5c+YQGRlJQkJCsXtbqcNRqIlUvKpmzZoABAYGFtj2xIkTDB8+nPr161O+fHmaNm1KfHx8jjYJCQmICFu2bKFfv35Uq1aNtm3bAhAREcGgQYNYvHgxTZs2JSgoKGvU4bp162jfvj3BwcFUrVqVPn36sH///hz77tKlCx07dmTNmjVcf/31lC9fnvnz53vj21Ak+b059HusV1R/iTVsuaCLvTq9rVIlXK1atfjHP/7BX//6V0aOHMmQIUN49dVXeeWVV2jdunXBOyiGMjIy+P3333OUiYjHM3Knp6eTnp5OSkoK48ePp27dutx33335bnP27Fk6dOjAhQsXmDx5Mo0aNWL9+vUMHz6cixcvMmrUqBzto6OjGThwICtXrswR68aNG0lKSmLSpEnUrl2biIgI1q1bR69evejWrRvLli3j3LlzPPPMM3Ts2JGkpCTq16+ftf2BAwcYPXo0EydOpHHjxsXjFRjuukBYL02r6fz3YOe628Xdfkr6opfXVFmVnp5uEhISTJ06dYyImJiYGLNgwQITHh5uRMSEh4ebpUuX2hqjJ5fXXC3Nmzf3+BitW7fO2q5JkyZmz549BW7z3HPPmfLly5sDBw7kKB86dKipWbOmuXz5sjHGmNdee80AZsyYMXn2ER4eboKDg80PP/yQJ54mTZpk7cMYYw4fPmzKlStnxo4dm1XWuXNnIyLmq6++8vhcvYV8Lq+57ekYY17P9u8EbyY6pVTx53A4eOihh+jTpw/PPvsss2fPzvyDFICUlBRiYmIA6y/14mrevHnceOONOcqCg4Oz/p27FxQQEIDIH3MKL1myhLNnz3L48GGmT5/OrbfeyqeffkpERITbY65bt462bdvSqFGjHPu//fbbWbRoEXv27OG6667LKr/nnntc7qddu3bUrVs3a/38+fPs2rWLp556inLl/vj4btSoER06dGDz5s05to+IiKBly5Zu47SDpxN+fiIiTd3URYqIziioVClVtWpVZs6cmePDL1NaWhqxsbE2ROW5yMhI2rRpk2Np3vyPKR0DAwNzLLk/uJs1a0bbtm0ZOHAgH3/8MefOnWPatGn5HvOnn35iy5Ytefbdr18/AH755Zcc7evVq+dyP7nLT506hTHGZfu6dety8mTOlzq726+d8runk10XoIqbuspAod7gqZQqedy9o6ekTyS6ffv2HOtRUVFu21arVo0mTZoUOJNDzZo1qV27NnPmzHFZn/sY2XtW+ZVXr14dEXH5szh+/HjWQIeC9msnT5MOuB9IcDXWw6JKqVIsLCzM5fMpVapUIT09vcS+KrtNG7ev1srjxx9/ZN++fQVeTuzRowdz584lLCyM2rVrFzXELBUrVqR169asWLGCyZMnZ33PU1JS+Pzzz/MMUCiO3CYdERkCDHGuGiBeRH7N1SwYaAF87JvwlFLFRVxcHDExMaSlpWWVBQQEcObMGXr27EliYiL+eunhypUrAdi5cycAa9euJTQ0lNDQUDp3znnhZe/evVSqVCnPPq699loqVqzo9hj33HMPrVq14rrrrqNKlSocOHCAWbNmUa5cOcaNG5dvfGPHjmXZsmV06tSJsWPHEhUVxfnz59m3bx9bt25l9erVhT3lLM8//zy9evWid+/ejBgxgnPnzjFp0iSqVq1aYFzFgrsRBlij0jY6lwxgZ7b1zGUdMBOo424/JX3R0WtK/WHp0qV5Rq8tXLjQlC9f3jRo0MBs27bNL3HgZlRa586ds9rkN3oNMNu3b8/3GNOmTTOtWrUyVatWNcHBwSYyMtLExMS4HCnnysmTJ82YMWNMRESECQwMNKGhoaZjx45m1qxZWW0yR68dPHgwz/bh4eEmOjra5b7Xrl1r2rVrZypUqGCqVKli7rrrLrNv374cbTp37mw6dOjgUazeRj6j18QYd1fN/iAiG4Hhxph93kl1JUebNm3Mjh077A5DqWJt165d9O3bl6NHjzJjxgxGjhxZLO8nKP8QkZ3GGJfXLT2de61rWUw4SinPtGrVip07d9KjRw9Gjx7NwIEDOXdOb/WqvDwaSCAiBb4Z1BjzRtHDUUqVVNWrV+edd97hxRdf5Omnn+brr79m1apVXHPNNXaHpooRTy+vuZt8KWtjY0zJHLpSAL28plThbdy4kfvuu4/z58+zcOFCBg4caHdIyo+KfHkNaORiacMf79dp64U4lVKlRNeuXdm1axctW7bk/vvvZ+TIkVy8eNHusFQx4Ok9nRQXyy5jzHPAm8Bjvg1TKVXS1K9fn40bN/LYY48xb948br755hL/IKkqOm/MEb4V6OWF/SilSpnAwEBmzJjBypUr2bt3L61atWL9+vV2h6VsVJgZCdxph85I4FKXLl3ylPXv358RI0aQlpbGHXfckad+8ODBDB48mBMnTtC3b9889cOHD2fAgAEcOXKEBx54IE/9uHHjuPPOO9m/fz/Dhg3LU//0009zyy23kJSUxJgxY/LUv/DCC9x00018/vnnPPXUU3nqZ8+eTcuWLfnoo4+YMmVKnvpXX32VqKgo1qxZw4wZM/LUL1myhIYNG7Js2TIWLFiQp37lypXUqlWLhIQEEhIS8tR/8MEHhISEMH/+fJYvX56nftOmTQBMnz6d9957L0ddcHAwa9euBawH7D7+OOczzTVr1mTVqlUATJgwgW3btuWob9CgAUuXLgVgzJgxJCUl5aiPjIzMel9KTEwMBw4cyFHfsmVLZs+eDcCgQYM4evRojvr27dszdepUAO69994883N1796diRMnAtCzZ08uXLiQo7537948/vjjQPH83Xv55ZeZMWMGPXv2JCwsjPDw8BzDqvV3r3j97mWej7d5OnrtGRfFQVizEfQCXvFmUEqp0qdBgwZ88cUXDBgwgPfff5+zZ8/SrFkzj16IpkqPooxeuwikAG8BU40xpfIuoY5eU8q7jDEsXLiQUaNGUadOHVauXJnn1QOqZPPGw6EOF0uwMaapMWZyaU04SinvExFiYmL47LPPcDgcdOzYkfnz5+PJH8Cq5LPvZeNKqTKtTZs27Nq1i9tuu43//d//ZdCgQTqLQRngNumISFhhFn8GrZQqHWrUqMG7775LXFwcb731Fm3btmXfPp1xqzTLr6eTDHxfiEUppQrN4XDw1FNPsWHDBn7++WduuOEGl6PDVOmQ3+i1h3H/4jallPKq7t2789VXX9G/f38GDBjAZ599xt///neCgoLsDk15kdukY4xJ8GMcSilF/fr12bRpE0888QSzZ89m+/btLF++nM2bNxMbG0tqaiphYWHExcUV+PZOVTx5NGQ6q7H1JNc1QA3gF2CvKeVDTnTItFL2WLFiBQ8//DAAly9fzjF3W0hICPHx8Zp4iilvTPiJiAwFfgD+A2wCvgGOichfvBGkUkpl169fP3bs2MHFixfzTBaalpZGbGysTZGpovB0RoJoIB74GFgKHAfqAtFAvIikGWPe9FmUSqkyKSoqit9//91lnU4eWjJ5OvfaE0CiMSb3hEuvi8gS4Ems2aaVUsqrwsLCSElJcVmuSh5PL69FYfVwXFnqrFdKKa+Li4sjJCQkR5nD4WDSpEk2RaSKwtOk8yvQwE1dA2e9Ukp5XXR0NPHx8VmzUtesWZOMjAwWL17MmTNn7A5PFZKnSWct8IKIdMpeKCLtgSnOeqWU8ono6GiSk5PJyMjgxIkTLF++nH//+99069aNn3/+2e7wVCF4mnSeAM4Am0QkVUT+LSIpwKfAWWe9Ukr5Rb9+/Vi9ejV79uyhc+fO/Pe//7U7JOUhT2eZPg60BB4FtmElmi+AUcD1xpgffRWgUkq50rNnT9atW8fRo0fp1KkThw8ftjsk5YFCPRxaFunDoUoVb9u3b6dHjx5UqFCBDz/8kGuuucbukMq8Ij8cKiK1cs8kLSLDRGSuiPT2RpBKKXUlbrjhBjZv3kxGRgY333wzO3futDsklQ9P7+ksBsZnrojIRGABcD+wWkQG+CA2pZTySIsWLdi6dSuVKlWiW7dubN261e6QlBueJp02WLMRZHoEeMEYUxOYBzzm7cCUUqowmjRpwqeffkq9evW4/fbbWb9+vd0hKRc8TTo1gB8BRKQF1hQ4rzvr3kEfDlVKFQMNGjRgy5YtREVFceedd7Jq1Sq7Q1K5eJp0fuGPh0O7AceMMQed64GF2I+tRGSIiBgR6WN3LEop36hduzYbN27khhtuoH///rz++usFb6T8xtO51z4CJotILWAcVu8mU1Mg78RIxYyIhAN/xRrqrZQqxapVq8aGDRvo06cPgwcP5tdff2XkyJF2h6Uo3MOhR4CpwHfAs9nqorEeEvWIiDRwjnrbJiJpzp5HhJu2DUVkpYicEZGzIvJ27lF0Hh7TAfwf1nNFFwtorpQqBSpWrMiaNWu4++67GTVqFC+88AL6iIj9POrpOB/+vNVN9S3Ab4U4ZhOgP7AT2Arc5qqRiIQAn2AliYewXp09BdgoItcZY84X4piPAZ8ZY3Za76FTSpUFFSpUYMWKFQwZMoTY2FjOnDnDtGnT0M8B+3h6ec0tY8zZQm6yxRhTB7JeDOcy6WBdCmsMRBljDjnb/wc4CAwDZjrLdgHuej/XA1WAvkAnN22UUqVYYGAgb7zxBpUrV+all17i7NmzzJs3D4ejRNyKLnWKnHQKyxiT4WHTu4AvMhOOc9vvReQz4G6cSccY0yq/nYjIcCAcOOj866Yu1ovn6hljFlzBKSilShiHw8H8+fOpWrUqL774Ir/++iuvvfYagYGBdodW5vg96RRCc2C1i/LdQD9Pd+JMLFnJRUQ2AbONMe8UMT6lVAkiIkybNo2qVavy1FNPce7cOd566y0qVKhgd2hlSnHuX9YATrkoPwlU9+WBRSRGRHaIyA6dNl2p0mXChAnMnTuX1atXc+edd3L+fGFuD6uicpt0ROQuEanqz2BccDXUpEh3AI0xXQrq5Rhj4o0xbYwxbUJDQ4tyOKVUMTRy5EgSEhL45JNPuO2221i4cCERERE4HA4iIiJITEy0O8RSK7/La/8C2gNfikg60N4Y86V/wgKsXk4NF+XVcd0DUkopjz300ENUqlSJ/v37s23btqzh1CkpKcTExADWy+OUd+V3ee0ckNnTsWN84W6s+zq5XQPs8XMsSqlS6N5776VWrVp5nt9JS0sjNjbWpqhKt/x6OjuBV0Vki3N9ooi4u8FhjDF/8W5ovAtMF5HGxpjDAM6HSDuQbcZrpZQqCnf3bVNTU/0cSdmQX9IZDswCbsa6t3IjcMlN20I95isifZ3/bO382tOZ0H42xmx2li0ERmK9OuFp5zGex5oZ4dXCHE8ppdwJCwsjJSXvTF5hYYWe/ER5wG3SMcbsB+4AEJEM4E4v3tNZkWt9vvPrZqCL8/jnRaQbVuJbgnWJ72NgjDHmnJfiUEqVcXFxccTExJCWlpZVFhAQwPPPP29jVKWXp0Omu+LF+yjGGHGzdMnVLtUYc68xpooxprIxpo8xJtlbcSilVHR0NPHx8YSHhyMiVKtWjfT0dLZv365ztfmAR0nHGLPZGHNORFqIyP+KyEQRGeF8t45SSpVo0dHRJCcnk5GRwalTpxg7dixz587lpZdesju0UsejGQlEpByQAAwk50g2IyL/BAYbY9K9H55SSvnf9OnT+eGHHxg/fjz16tXjwQcftDukUsPTy2uTsGaGfgZoBAQ7vz4DDHB+VUqpUsHhcJCQkEC3bt34y1/+wrp16+wOqdTwNOkMAp43xsQZY1KMMRedX+OwXjegfwYopUqV8uXL869//YvmzZvTt29fduzYYXdIpYKnSecqYJubus+d9UopVapUqVKFtWvXUqtWLXr16sV3331nd0glnqdJ5xjWQ5mu3OSsV0qpUqdevXqsX7+e9PR0br/9dn766Se7QyrRPE06iUCsc9RaYxEJFpFGIjIBiMV6jkYppUqlqKgo3nvvPY4dO0avXr04d04fFbxSniadycBK4FmsN3eeAw4BcdnKlVKq1GrXrh3Lly/nq6++om/fvly+fNnukEokT5/T+d0Ycz9wLdbUNM84v7YwxkTrcGmlVFnQu3dv/vGPf7B+/XqGDh2qD49egUK9OdQYsxtr9mellCqThg4dyrFjx5g0aRJXXXUVU6dOtTukEqU4v65aKaWKpYkTJ3Ls2DGmTZtG/fr1GTlypN0hlRiadJRSqpBEhHnz5nH8+HFGjx5N3bp16du3b8EbKo8HEiillMomICCAN998k/bt2xMdHc3mzZsL3khp0lFKqSsVHBzMmjVruPrqq7n77rv55ptv7A6p2Csw6YhIkIg8qjNKK6VUXjVq1GDdunVUrFiRnj176htHC1Bg0jHGXAKmATV8H45SSpU8YWFhrFu3jnPnztGjRw9Onjxpd0jFlqeX1/YCjX0ZiFJKlWTXXnst77zzDt999x3t2rUjLCwMh8NBREQEiYmJdodXbHiadJ4BJorItb4MRimlSrIuXbowbNgwDh48yJEjRzDGkJKSQkxMjCYeJ/HkiVoR2QpEAjWBZOAHIPuGxhjT2RcB2q1NmzZGpzRXSnkqIiKClJSUPOXh4eEkJyf7PyAbiMhOY0wbV3WePqeTDuzxXkhKKVU6uRtIoAMMLB4lHWNMFx/HoZRSpUJYWJjLnk5YWJgN0RQ/+pyOUkp5UVxcHCEhITnKRITY2FibIipePE46IlJPRKaLyHYR+U5EvhSRl0Skri8DVEqpkiQ6Opr4+HjCw8MREerUqQPA+++/T0ZGhs3R2c+jpCMikUASMBrrXTpfAueBR4EkEfmTrwJUSqmSJjo6muTkZDIyMjh+/DizZs1i9erVvPTSS3aHZjtPezovAmeBSGNMV2PMQGNMV6wRbWec9UoppVwYPXo0AwYMIDY2lk8++cTucGzladLpCkw0xiRnLzTGpGC9VbSrd8NSSqnSQ0RYtGgRUVFR3HfffRw9etTukGzjadIJAn51U/ers14ppZQblSpV4u233+bChQv079+fS5cu2R2SLTxNOknAKBHJ0V5EBBjhrFdKKZWPpk2bsnjxYrZt28bjjz9udzi28PTh0OeA94C9IrIMa0aCukA/4E9AL9+Ep5RSpUu/fv0YO3Yss2bNon379gwcONDukPzKo2lwAESkBzAFuB4QrGlwdmLd61nvswhtptPgKKW87fLly3Tr1o1du3bx5Zdf0rx5c7tD8qr8psHx5H06gSJyN7DfuZPKQEOgsjHmxtKccJRSyhcCAwNZvnw5lStX5n/+5384e/as3SH5jSfv07kMLAcinOtpxpj/GmPSfBybUkqVWvXq1WP58uV89913DBkyBE+vOpV0ng4kOAzU9mUgSilV1tx88828+OKLvP3228ycOdPucPzC06TzEhArIqG+DEYppcqaxx57jHvvvZcnn3ySLVu22B2Oz3k6eq0b1uuqvxeRL3D9Pp2HvB2cUkqVdiLC4sWL+fbbb+nfvz+7du3iqquusjssn/G0p9MJuAz8DFwNdHSWZV+UUkpdgSpVqrBq1Sp+/fVXBgwYwOXLl+0OyWc8SjrGmAhjTKN8lsa+DlQppUqz5s2bs2jRIj799FOefPJJu8PxGU+GTAeJyC4Ruc0fASmlVFk1cOBARo0axaxZs1ixYoXd4fiEJ0OmLwGNgN99H45SSpVt06dPp3379jz88MPs3bvX7nC8ztN7Oh8C2tNRSikfCwoKYvny5QQHB3PLLbcQFhaGw+EgIiKCxMREu8MrMk9Hr80FlopIOeAd8o5ewxhz2LuhKaVU2dSgQQOGDh3K1KlTs8pSUlKIiYkBrJfElVQezb0mItnfsepyA2NMgLeCKk507jWllB0iIiJISUnJUx4eHk5ycrL/AyqE/OZe87SnM8SL8fidiARhPeDaC7gEpBhj7rA3KqWUci81NbVQ5SWFR0nHGPO6rwPxsRewXjQXZYzJEJF6dgeklFL5CQsLc9nTCQsLsyEa7/F0IIFbIuIQkRqFaN9AROaKyDYRSRMRIyIRbto2FJGVInJGRM6KyNsiUqjvuIiEADHAeGNMBoAx5ofC7EMppfwtLi6OkJCQHGWBgYHExcXZFJF3uE06InJSRFplWxcReVdEcj8IegPWTAWeagL0B04BW/M5fgjwCdAUeAh4AOuFcRtFpGIhj3cKGC8i20Xkc+erGpRSqtiKjo4mPj6e8PBwRISQkBDS09O59tpr7Q6tSPLr6VQj5+U3B9DbWV4UW4wxdZz3VPJ7+umvQGOgjzHmHWPMauAuIBwYltnI+eDqCTdLQyAQCAO+M8bcADwMLBKRq4t4Hkop5VPR0dEkJyeTkZFBSkoKoaGhPPjgg1y6dMnu0K5YkS+vFVbmJS4P3AV8YYw5lG3b74HPgLuzlbUyxtRysxwBUrBG3C1xtt8HJGG9AVUppUqEWrVqsXDhQr7++mueffZZu8O5Yn5POoXQHPjWRflu4BpPd2KMOQGsB3oAOAcRtAC+cbeNiMSIyA4R2fHzz4W5cqiUUr5z5513MmTIEKZNm8YXX3xhdzhXpDgnnRpY92JyOwlUL+S+hgOPisg3wDrgcWPMfneNjTHxxpg2xpg2oaH6CiGlVPExe/ZsGjRowIMPPsj58+ftDqfQCko69UWksXPwQOPcZc7yBj6Mz9WDqFLonRiTbIzpboy51hjzZ2NMyZ9LQilVJlWpUoWEhAQOHjzI+PHj7Q6n0Ap6Tmeli7J3cq0LbmYpKKJTWL2d3KrjugeklFJlQteuXXn00UeZM2cOffr0oXv37naH5LH8ko7dsxDsxrqvk9s1wB4/x6KUUsXK1KlTWbduHUOGDOGbb76hatWqdofkEbdJpxjMQvAuMF1EGmdOJup8iLQDUPL6lEop5UXBwcG88cYb3HTTTTz66KMkJCTYHZJHbBlIICJ9RaQv0NpZ1NNZ1jlbs4VAMrBaRO4WkbuA1cAR4FW/BqyUUsXQjTfeyIQJE3j99dd555137A7HIx7NMu31g4q4O+hmY0yXbO3CgFnArVj3jj4Gxhhjkn0dYyadZVopVZxdunSJdu3acfToUb799ltq165td0j5zjJtS0/HGCNuli652qUaY+41xlQxxlQ2xvTxZ8JRSqniLigoiDfeeIMzZ87wyCOPYEdHojCK83M6SimlPNCiRQumTJnCv/71L5YuXWp3OPnSpKOUUqXAY489RseOHRk1ahRHjhyxOxy3NOkopVQpEBAQQEJCAr///jsPP/wwGRmeTnPpX5p0lFKqlLj66quZMWMGH330EQsWLLA7HJc06SilVCkSExNDjx49GDt2LPXr18fhcBAREUFiYvGY/UuTjlJKlSIiQq9evbh8+TLHjh3DGENKSgoxMTHFIvFo0lFKqVJm+vTpecrS0tKIjY21IZqcNOkopVQpk5qaWqhyf9Kko5RSpUxYWFihyv1Jk45SSpUycXFxhISE5CgLCgoiLi7Opoj+oElHKaVKmejoaOLj4wkPD0dECAoKokKFCvTu3dvu0DTpKKVUaRQdHU1ycjIZGRl89tln/Prrr0ycONHusDTpKKVUademTRtGjBjBvHnz2LVrl62xaNJRSqkyYMqUKYSGhvLII4+Qnp5uWxyadJRSqgyoVq0aM2bMYPv27cTHx9sWhyYdpZQqI+6//366devGhAkT+PHHH22JQZOOUkqVESLCvHnzSEtL429/+5stMWjSUUqpMqRp06Y88cQTLFmyhE2bNvn9+FLcX21qtzZt2pgdO3bYHYZSSnlNWloazZs3Jzg4mKSkJIKCgry6fxHZaYxp46pOezpKKVXGhISE8Morr7B3715mzpzp12Nr0lFKqTKoV69e3HPPPTz33HMkJyf77biadJRSqoyaM2cODoeD0aNH++2YmnSUUqqMatiwIZMmTWLNmjW8++67fjmmDiQogA4kUEqVZpcvX6ZVq1acPXuWPXv2ULFixSLvUwcSKKWUcikwMJAFCxaQmprK888/7/PjadJRSqkyrmPHjgwZMoQZM2awe/dunx5Lk45SSilefPFFKleuzIgRI/DlbRdNOkoppQgNDeXFF19ky5YthIaG4nA4iIiIIDEx0avHKefVvSmllCqxgoODcTgc/PLLLwCkpKQQExMDWC+F8wbt6SillALg6aefJiMjI0dZWloasbGxXjuGJh2llFIApKamFqr8SmjSUUopBUBYWFihyq+EJh2llFIAxMXFERISkqMsJCSEuLg4rx1Dk45SSinAGiwQHx9PeHg4IkJ4eDjx8fFeG0QAOg1OgXQaHKWUKhydBkcppVSxoElHKaWU32jSUUop5TeadJRSSvmNJh2llFJ+o6PXCiAiPwOngTNXsHkt4IRXA1LuVOXKfkbFXXE9L7vi8vVxvb1/b+2vKPu50m2L8vkVbowJdVWhSccDIhJvjIm5gu12uBs2qLzrSn9GxV1xPS+74vL1cb29f2/tryj7KW6fX3p5zTNr7A5AFai0/oyK63nZFZevj+vt/Xtrf0XZT7H6HdKejg9pT0cpVVJpT6dkirc7AKWUukI++fzSno5SSim/0Z6OUkopv9Gko5RSym806dhERK4WkU9F5ICIfCUiOuBAKVViiMhTIrJfRDJEpI+n22nSsc8/gARjTCTwBJAoImJzTEop5amPgTuALYXZSJOOh0SkgYjMFZFtIpImIkZEIty0bSgiK0XkjIicFZG3RSQsW30o0A54HcAY86GzqrWvz0MpVTZ58zMMwBjzb2PMd4WNQ5OO55oA/YFTwFZ3jUQkBPgEaAo8BDwA/AnYKCIVnc3CgGPGmMvZNk1xliullC948zPsipUr6g7KkC3GmDoAIjIUuM1Nu78CjYEoY8whZ/v/AAeBYcBMN9vppTWllC/5+jPMI9rT8ZAxJsPDpncBX2T+sJzbfg98BtztLEoFrhKRwGzbhTvLlVLK67z8GXbFNOl4X3PgWxflu4FrAIwxPwNfAoMBRORWrJ7OTv+EqJRSbhX4GVYUmnS8rwbWNdPcTgLVs60/AgwRkQPA34Foo9NDKKXs59FnmIg8LSJHgfbAIhE5KiJ1C9q53tPxDVfJI8c9G2PMQeAm/4SjlFKF4sln2BRgSmF3rD0d7zuF9ZdCbtVx/deDUkoVJz79DNOk4327sa6J5nYNsMfPsSilVGH59DNMk473vQu0E5HGmQXOB7A6OOuUUqo48+lnmL7aoBBEpK/zn92xBgKMAH4GfjbGbHa2qQh8DVwAnsa6Nvo8UBm4zhhzzt9xK6UUFI/PME06hSAi7r5Zm40xXbK1CwNmAZlDoT8Gxhhjkn0do1JKuVMcPsM06SillPIbvaejlFLKbzTpKKWU8htNOkoppfxGk45SSim/0aSjlFLKbzTpKKWU8htNOkoppfxGk47yOxEZ7Hw/+2kRqZ6rrpyzbrINcU12HrtYz74uIg4RmS0iP4hIhoi8k0/bHN9LEekjIo/5I053RGSMiPyPi/LJ+Ty8qEoJTTrKTlWBJ+0OogTqCzyK9R6mDsAT+bRtDyzKtt4HsDXpAGOAPEkHK872/g1F+Vux/otOlXobgFEiMtsYc9zuYPxBRMobYy4WcTfNnF9nF/QKYmPMF0U8VoG8dE4YY44CR70QkirGtKej7JT5AqjY/Bq5u+wiIgkikpxtPcJ5OekREZkqIsdF5FcRWSoiISLSRETWi8g5ETkkIg+5OWQzEdkoImnOS1jPiUiO/ysiUktEFojIf0XkoojsE5GYXG0yLyPeLCIrROQ08O8CzrWHiGwTkQsickZE3hGRqGz1ycBk52q6c/+D89lf1uU1EUkAHgLqO8tNru9fkc5JRG4QkZXON0heEJH9IvKCiATnij8ciM4WQ4KzLs/PWUSqiMgrInLMGdN+ERkrIpKtTRfnfu5ytj0hIj87f+7Vcu3vURHZ64zvlIjsEJF78vmRKC/Tno6y0w/AK8AYEZlujEnx0n4nAJuwPmCvAV4CMoDrgYXAdGA48JqI7DDG7M61/TvAYmAqcDsw0bn9ZLA+CIHPgGBn2ffOdgucf/XPzbW/ROBNrMtibv/PiUgP4H3gE2AAUAl4DvhURFoaY/4L3AOMBgbzx6Wo7zz4noA1U3AocANwl7PsohfPKQxIAhKAX7HeyfIM0Bi4z9nmHuADrFmMJzvLfnYVrDPRvw+0cu7nG6AXMNN5Hk/l2mQO8B5wPxCF9XNPx/o9QESigRlY39OtznO9DtcvLFO+YozRRRe/LlgfmAZogvUf/jSw2FlXzlk3OVv7ydavap79JADJ2dYjnNt+kqvd287yQdnKqgO/A5NyHwcYn2v7hVgfotWc6xOB34A/uWh3AiiX6zxnefh92QEczNzeWdYIuAzMzFY2xdX3w80+c38vE4CjLtp59ZywZiYuBwzCStg1s9UlA0tdbJPj5wz0dh5rcK52i7CSZS3nehdnu9dztXvFeU6SbX2X3b//ZX3Ry2vKVsaYk1h/fT6Y/TJSEa3Ntb7P+XV9tuOeAn4CGrrYfnmu9beweh0tnOs9sC4pfS/WaLtyzhFv64GaWL2r7P5VUMBivcOkFbDMGPN7tji/x+qBdC5oH0VU5HNyXgp7UUS+w0oKl4ElWAnoT1cQ081YCevNXOVLgSDyDjp4P9f6N0B5oI5zfTvQUkTmisgtIhJyBTGpItLLa6o4mAWMwrrsEe2F/eV+j/ulfMoruNj+Rzfr9Z1fa2P10i67OX7NXOs/uGmXXXWsD2dXbY9j3QfxJW+c02vALViXwpKA88CNwDxcf58LUgM4afIOUjierT67k7nWM7fLPPYbzn//BevlZZdF5APgMaPvuvIbTTrKdsaYcyIyFavH83cXTX4DEJEgY8ylbOW5Pwi9pQ5wONc6wH+dX3/B6iU96mb7/bnWPXn25JSzXV0XdXWdx/SlIp2TiFQA7sa6lDcnW/m1RYjpJFDDxc8983tUqO+Jsa6xvQq8KtbzYbdh/c4tA9oWIU5VCHp5TRUX87E+1Ke4qMscYJB5eQvnqKSbfBRL/1zr9wHngG+d6+uApkCqMWaHi+XXwh7QGHMe2An0E5GAzHIRCcc6z81XciIuXMS6gZ5bUc+pPBBA3p7S4ELEkNtmrM+ofrnKo7F6qVc8HNwYc8oYswzrUmqLgtor79GejioWjDEXReQ5IN5F9VrgDLBQRCZhfcA9gZUIfOGvzpFT27FGcA3F+gv+tLN+Ftbosq0iMgurF1AR60O7kzHm7is87kSs+xLvich8rPtIz2Kd+4wr3Gdue7B6D8OxBi78Zoz5hiKekzHmjIh8AYwTkR+wBh88zB+XJHPH0ElEemNdKjvh5vLWWuBT4B8iEgrsBu7A+nlMNcacKMyJi0g81oCQbVi9ukjgAaznxZSfaE9HFSevYY3eysH5Yd8b66bycqyhzHOBjT6K426sd8O/izX6agrWcOPMeM5g9T4+wJpRYT3WEOu7ixKTMWYd1pDgaljn+Q9gL9DRGHPsSvebyyKsgREvAF8Ca5zH9sY5DcTqrc3DGiV3HNeX6yZgJbXlWIl9squdGevB117A686Y3neuP0YBz3a58RnQGqtX/aFzH0txDqlW/pE5lFAppZTyOe3pKKWU8htNOkoppfxGk45SSim/0aSjlFLKbzTpKKWU8htNOkoppfxGk45SSim/0aSjlFLKbzTpKKWU8pv/B9SY1RgZgPHEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from math import exp\n",
    "from numpy import linspace\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# Set up lists for plotting\n",
    "delta = []\n",
    "it = []\n",
    "err1 = 0.1\n",
    "err2 = 1E-3\n",
    "\n",
    "# initialize iteration\n",
    "x = 1.0\n",
    "x2 = x \n",
    "N = 10\n",
    "for k in range(1,N+1):\n",
    "    x = 2 - exp(-x)         # iterate value\n",
    "    delta.append(abs((x-x2)/x)) # compute difference\n",
    "    x2=x                    # update x2\n",
    "    it.append(k)\n",
    "    print('x= ',x)\n",
    "    \n",
    "# Make the graph\n",
    "plt.rc('font',size=16) # set the font size\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.plot(it,delta,\"-ko\")\n",
    "plt.xlabel('Number of iterations')\n",
    "plt.ylabel('Error of solution')\n",
    "plt.hlines(err1,it[0],it[N-1],linestyles='dashed',colors='k')\n",
    "plt.hlines(err2,it[0],it[N-1],linestyles='dashed',colors='k')\n",
    "plt.text(5, err1+0.1, '1E-1 error')\n",
    "plt.text(5, err2+0.001, '1E-3 error')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x=  2.117000016612675\n",
      "x=  0.030755419069985038\n",
      "x=  2.715711832754083\n",
      "x=  0.0017034651847384463\n",
      "x=  2.71827394057758\n",
      "x=  0.001679913095081425\n",
      "x=  2.7182741571849562\n",
      "x=  0.0016799111168229455\n",
      "x=  2.7182741572030236\n",
      "x=  0.0016799111166579386\n",
      "x=  2.7182741572030253\n",
      "x=  0.0016799111166579221\n",
      "x=  2.7182741572030253\n",
      "x=  0.0016799111166579221\n",
      "x=  2.7182741572030253\n",
      "x=  0.0016799111166579221\n",
      "x=  2.7182741572030253\n",
      "x=  0.0016799111166579221\n",
      "x=  2.7182741572030253\n",
      "x=  0.0016799111166579221\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEUCAYAAADwYOuyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABSSElEQVR4nO2deXhV1dX/PyshZCIQCGEQCAGBxImioBJAoah11jrXAq2+bbHwE1+nt3VCca5WKr4qTrwKNdiC2FatBK0KCIKt2EKUGcI8EwKZIcP+/XHviTc3dzh3voT1eZ7z5J5z9tln3Zvc883ae621xRiDoiiKokSDhFgboCiKopw4qOgoiqIoUUNFR1EURYkaKjqKoihK1FDRURRFUaKGio6iKIoSNdrE2oB4pnPnziY3NzfWZiiKohxXfPPNNweNMdmezqno+CA3N5cVK1bE2gxFUZTjChHZ5u2cDq8piqIoUUNFR1EURYkaKjqKoihK1FDRURRFUaKGio6iKIoSNVR0FCXGzJ49m9zcXBISEsjNzWX27NmxNqkFsbIxUvcNZ7+h9BXstYFeF0j7iP+ujTG6edkGDx5sFCWSFBYWmrS0NAM0bWlpaaawsDDWpjURKxsjdd9w9htKX8FeG+h1gbQP12cDrDBenqsxf7DH86aio0Sa3r17N/uCW1vv3r1jbVoTsbIxUvcNZ7+h9BXstYFeF0j7cH02vkRHHOcVTwwZMsRocqgSSRISEvD0HRQRGhsbY2BRS2JlY6TuG85+Q+kr2GsDvS6Q9uH6bETkG2PMEI/22+5FUZSwk5OTE9DxWBArGyN133D2G0pfwV4b6HWBHI/K79qbC6SbDq8pkaewsNAkJSXpnI6X+7Zp0ybu53SSk5Nb1ZxOsO/HFXROR0VHiV9+8IMfND1ck5OT40pwLKZOndr0EEpJSYmajQMHDmw2rxCu+xYWFpqEhAQDmKysrJD6nTRpUtA2FhYWmvT0dAOYzMxM29f+8Y9/bLpnTk6O3+tmzZpl28Z7773XAEZEgv7MVXRUdJQ4paamxqSmpppJkyaZO+64w6Smpppjx47F2qwWvPfeewYwffv2NX379o3affPy8gxgunbtGtZ+q6urmx7CTz/9dEh9vfTSSwYw2dnZQV1/1VVXGcDcf//9tq85cOBAk/0VFRV+2x88eNAAJikpyTQ2Nvps++c//9kAZvXq1bbtcceX6OicjqLEkMWLF1NTU8Oll15KQUEBNTU1FBcXx9qsFhQXF5OQkMC1117L1q1bOXr0aMTv2djYyJYtW2jTpg379u2jtrY2bH3v3r276fXevXtD6mv//v0AlJaWBhXgUFpaCkBZWZntaw4fPtz0+siRI37bl5eXA1BXV+f3d1dRUQFARkaGbXsCQUVHUWLI/PnzSUlJYdSoURQUFACwfPnyGFvVklWrVtG/f3/OPPNMGhsb2bRpU8TvuXv3bo4dO8bZZ58NwM6dO8PW965du5peh0t0GhsbAxIOC0t0Dh06ZPsa1/u4CpA3XIXJEiBvqOgoSiumqKiIH/7wh6SmppKTk0P37t356quvYm1WC4qLi/nBD35AXl4eAOvWrYv4PUtKSgAYNWoUANu3bw9b35ank5WVFTbRATh48GDA10fT03F/7QlLdNq1a2fbnkBQ0VGUGLF582Y2btzIpZdeCjhyIQoKCuLO0ykvL6ekpISBAwcyYMAAANavXx/x+0ZSdCxPZ8iQIezZsyekvvbv309CguNRGqjoGGOaPJxoeTqWqHijoqKC1NRU2rSJzBqfKjqKEiOKiooAmkQHoKCggJKSkmb/Pcea7777DoCBAweSkZHBSSedFDXRSUhIYNiwYYhI2EUnPT2dAQMGhMXTOfnkkwE4cOBAQNceOXKEhoYGIL48nUgNrYGKjqLEjPnz59O/f3/69evXdCwe53WswIYf/OAHAOTn50dNdHr16kW7du3o1q1b2EXnpJNOonv37pSXl1NdXR10X/v27ePUU08FAvd0rKG1jIyMgDydQEUn0DkdFR1FaWXU1NSwcOHCZl4OwODBg0lKSoor0Vm1ahWZmZn06tULgLy8PNatW+fIuYggJSUl9O3bF3BkxIdbdHr06EG3bt0Ah3AEw9GjRzly5AinnXYaELinYwlNv379OHz4cJPX44+ysjJEBLA3vKaejqKc4CxevJja2toWopOSksKZZ54ZV6JTXFzMwIEDmx5yeXl5HDlyJOJDgNEUnWDndSyR6d27N+np6UF7Ov379wfseS3gEJqsrCySkpLU01EUxT/z588nNTWVkSNHtjhXUFDA119/TV1dXQwsa05jY2OT6FhYEWyRHGKrqqpi3759LUQnHN6VMYbdu3fTo0cPunfvDgQfNm0Jb5cuXejcuXPQomMNsdodYisrK6Njx4506NDBdiBBamoqYC+QQEVHUVoZrqHS7gwdOjRukkS3bt1KZWVlM9HJz88HIis6W7ZsAWgSnV69elFTU9P0kA6F0tJSjh071szTCZfoBDq85u7p2A0mOHz4MB07diQzM9N2IEH37t1JSEjw6+mUl5er6ChKa2LTpk1s2rSpxdCaRTwFE7gHEYDD60hJSYloro4VLu3q6UB4wqatcOmTTjqJ7OxsEhISgh5ecxWd7OzsoDwdEWl6n4F4OpmZmQF5Oh06dKB9+/a2htfat29vy45gUNFRlCjjKVTaFStJNB5EZ9WqVYhI00Q5ONZc6d+/f0Q9nWiITo8ePUhMTCQ7Ozumw2uZmZlkZ2cDkfV02rdvb1t01NNRlFbE/PnzGTBgQFNuhzvxlCRaXFxM//79SU9Pb3Y8Ly8v4qKTkZFBVlYWEDnRAejevXvQorNv3z6Sk5PJyMggOzs7qOG1Tp060bFjR8C+p3P48OEmT8duIIHl6fia02lsbKSqqkpFR1FaCzU1NSxatMirl2NRUFDAli1bgg7lDRerVq1qNp9jkZ+fz5YtWyJW+NOKXLMi5jp37kxKSkpYRGf37t2ISFMQQbdu3ULydLp27YqI0LlzZyorKwMqTFpaWkpWVlaT6NjxdIwxAQ+vWZ5ORkaGT0+nsrISiFzdNVDRUZSosmjRIo+h0u5Y8zqxrMNWWVnJ5s2bm83nWOTl5dHQ0MDmzZsjcm/XcGlweH/hCpvetWsXXbp0ISkpCXCITihzOl26dAEcwgiBJYgeOnSIrKwskpOTSUtLs+Xp1NTUUFdXF9Dwmt05nUgX+wQVHUWJKr5CpV2JhyTRb7/9FsCjpxPJsGlrSQNX0YHw5epYOToW3bt3Z9++fUEtS+AqOta8TCCiY3k6AJ06dbIlOpY3ZHk6lZWV1NfXe21vjLE9p6OioyitjKKiIkaPHk1KSorPdvGQJGpFrkVbdPbu3UttbW1EReekk05q2u/WrRv19fUBlaGxCNXTcRWdjh072hpes4bTLE8HfCd81tTUUF9fr56OopxobNy4kc2bN/sdWrOIdZJocXEx7du3p3fv3i3OtW/fnu7du0dEdNwj1yxycnLYs2dPyPNI7p5OsLk6xhiPomM3mODYsWNUVFSE7OmA71I41vBbhw4dyMjI8BlIoKKjKK0If6HS7sR6JVEriMCazHcnUhFsvkQHmi/AFihHjx7l4MGDHkUn0HmdI0eOcOzYsaCH1yyBcRWdYD0dX/M6lmdjDa9VVFR4HUpU0VGUVsT8+fPJy8tr8TD1RiyTRI0xTQu3eSNShT9LSkoQkRYeVjjCpi1hcZ/TgcA9HStHp2vXroBDBETEtuhY1Qhch9fseDqW6Lh6Or5Ex9XTsZI+rSg1d1R0FKWVUF1dbStU2pVevXpx0kknxUR0tm7dSkVFhcf5HIu8vDzKysqCWi3TFyUlJfTs2ZPk5ORmx8MhOu45OhD88JprYihAYmIiWVlZtofXLNHp1KlT0087nk6gw2vuno7rMXdUdBSllbBo0SKOHj0akOjEMknUVxCBRaRqsLmHS1v07NkTCL/otGvXjrS0tJBFBwioKoH78FrHjh2pqanxm+fj6unYGV7z5Ol4m9dR0VGUVkJRURFpaWmcf/75AV03dOjQmCSJFhcXIyKcfvrpXttEKoLNm+ikpqbSpUuXsIiOa/SaiASVq+NNdAL1dFzndMB/gmhZWRnp6ekkJSUF7OlYYuLL00lISCAtLc3WewgGFR1FiTDGGObPn28rVNqdWM3rrFq1ipNPPpl27dp5bdO7d2+Sk5PDWvizurqaPXv2eJ33CjVseteuXaSkpDRVALAIphSOJTpW1BoQUNFPT3M64L8UjlV3DQh6TseX6LRr185r8Eg4UNFRlAizceNGSkpKAhpas4hVkqi/IAJwzGH069cvrJ7O1q1bgZaRaxbhEJ0ePXq0eKgGUwpn//79dOzYkbZt2zYdC2R4rbS0lLZt2zbVtbM8HX+iY5XAAWjTpg3p6em2PJ2MjAxbohPJoTVQ0VGUiBNoqLQrsUgSraqqYtOmTT7ncyzy8/PDKjrewqUtQl3MzVq8zZ1gRGffvn3Nhtbge9GxY5+VGGoJoN36a66eDuC3FM6RI0dIT0+nTZs2tuZ0VHTCgIikiMjfRGStiKwUkY9FxF7cqqKEyPz588nPz6dPnz5BXV9QUMCKFSuiliT63XffYYyxJTp5eXmUlJSEzTY7olNVVWV7CQB33BNDLbp160ZZWVlAxTqtYp+uZGdnU19fb6semms1ArDv6VgVpi38VZo+cuRIk9jYmdNR0QkfrxhjTjHGDAI+BGbE2B7lBKC6uprFixcH5eVYRDtJ1NPCbd7Iy8ujvr4+bIU/S0pKSE9Pb0q0dCeUsGljTIsSOBZWrk4gARuu1QgsAimFYy1rYBFIIIG76PgbXrPmfvyJTqRXDYU4Fh0R6SkiL4rIchGpFhEjIrle2vYSkXkickREykXkLyKSY503xtQaYz52ueQrQD0dJeIsXLgw4FBpd6IdTLBq1SoyMjI8lr9xJ9wRbO5LGrgTiuhYnow3TwcCy9XxJTp2ItisCtMW7du3R0QCCiQAe8Nrlui0bduWlJQU9XS80A+4ESgDlnhrJCJpwOdAPvBzYBzQH1goIuleLpsEvB9WaxXFA8GGSrsS7STR4uJiBg4cSEKC/8dDpETHG6GIjqccHYtARaeuro5Dhw61EJ1ASuG4D68lJCT4LfrZ0NDAkSNHAvZ0XJef9lX080QXnS+MMV2NMZcB7/po9yscXsuPjTF/M8a8D1wF9AZuc28sIvcDA4D7I2CzojRhhUpfcMEFLbLrAyGaSaJW+Rs78zng+C+7a9euYREdY4xf0cnOziY5OTko0dm9ezfgW3Ts5upYohKsp2OMaSE64L8UjiUWwXo6gM+inxUVFc0EKhIEJDoi0k1EzhGR8923cBtmjLG7uMVVwFfGmE0u124BvgSudm0oIvcC1wGXGmOqw2Wronhiw4YNbNmyJaShNYtorSS6fft2jhw5Yms+x8KqwRYq+/bto6amxqfoJCQk0KtXr7B7Ol26dEFEbHs6nhJDwb6nU1lZSV1dXQvR8Vdp2rUEjoXl6XiLmDsuPR0R6SEiC4FdwHJgocu2yPkzVpwGfOfh+GrgVGtHRO4GbgYuMsYcjo5pyolMKKHS7kRrXmfVqlWA7/I37oSr2rS/yDWLYHN1LNGxggZcSUpKonPnzrZFxxJ/9+i1tLQ0UlJS/IqOe2Kohb/hNdcK0xaZmZnU1dV5jbxz93S8ic6xY8c4duxYfIgO8ApwOvAb4FJgtMv2Q+fPWNEJx7yPO4eAjuAISgCmApk45npWisgKT52JyHgRWSEiK+yWs1AUTxQVFZGfn09ubm7IfZ111llRSRK1Itd8lb9xJz8/n9LS0qYHabBEQ3Q6d+7sdagzkFI43jwdEbFVCseb6PjzdFzrrln4qkrQ0NBAZWWlLU8nGnXXANrYbHcecIcx5u1IGhMCnvzKpvAXY8xO132fHRnzOvA6wJAhQ8Jbs105YaiqqmLRokXcfvvtYekvWkmixcXFnHzyyQE9eFyDCYYNGxb0vS3R8SfSOTk57N69m7q6OpKSkmz37y1HxyKQBFFvogP2SuG4V5i28Fdp2tvwGjgEyZqbsrCExM6cTrREx66nUwPsj6QhIVCGw9txpyOePSAlRsyePZvc3FwSEhLIzc1l9uzZsTYpYixcuJBjx46FZWjNIhpJotbCbYFgiU6o8zolJSX06NHDb326nJycppybQPAnOoHUX9u/f3+zgpuu2CmF415h2sIaXvO2yJq34TXw7OlYx+LJ07ErOm/gCEWOR1bjmNdx51RgTZRtUbwwe/Zsxo8fz7Zt2zDGsG3bNsaPH99qhaeoqIj09HTOO++8sPVpJYla8y7hprq6mo0bNwYURAAOzyQpKSnkeR1/kWsWwYZNeyuBY2F5OnZK2Fg5Op7yiUIdXmtsbPQaXebP03HHEhc7czrxJjq7gGEi8rmI3Csi/+W+RdJIP3wADHUta+NMIh3uPKfEAQ8++CDV1c0DBqurq3nwwQdjZFHkCFeotDuRDiYIpPyNK23atKF///5xLTp1dXXs37/fr+gcO3bMVokdT4mhFqEMr/mrNH348GESEhKaCUMwns7Ro0c5duxYs7bxJjqvArnAKOBZHCVkXLc3ImAbInK9iFwPDHYeutR5bKRLszeArcD7InK1iFyFI/FzB/BaJOxSAsfbAyKUisHxyvr169m6dWtYh9Yg8kmidhZu80aoEWy1tbXs2rXLluj06tULCOxvZ8+ePRhj/IoO2EsQ9SU6nTt35siRIy0e6q6UlpbSvn37FnNS/krhWCVwXD2sYDwdaFn0M95Ep4+fLVIlZd51br927k937j9qNTDGVOGIntsAvA3MBrYAo40xnhcCV6LK0aNHSU1N9XjO+q+1NRHOUGlXrCTRr776Kqz9WhQXF9OuXbugCpPm5eWxadOmoOeb/C1p4EpaWhqdO3cOSHQ8Ld7mjhVKbUd09u3b1yJc2sLK1fEVzecpMRTseTruawHZ8XTcAwmgZf21uIpeM8Zsi6gV3u9rN+JsO46kTyXOOHz4MNdeey3V1dUkJSU1eyilpKTw5JNPxtC6yFBUVMQpp5xiq3ZZoBQUFPDee+/5fOgFy6pVqzjjjDNslb9xxyr8uWXLFgYMGBDw9XbDpS0CDZv2lRhqYdfTMcb49XTAkSDqKScIvIuOv0rT7hWmAdLT00lMTPQoOq6rhlp4W1Mn3jwdAETkdBH5fyIyWUQmioj9YH7lhGPHjh2MGDGCpUuXUlhYyFtvvUXv3r0RERISEujSpQvXX399rM0MK5WVlSxevJjLLrssIv1Hal7HKn8TaBCBRX5+PhB8DbZgRGfHjh22+/dVAsfCbimcyspKamtr/YqOr2ACf6Ljb3jNFRHxWn/Nk6dzXIiOiLQRkUJgFfAijuGtl4BVIvK2iCRG0EblOKS4uJihQ4eyY8cOFixYwJgxYxgzZgxbt26lsbGRv/71r2zfvp2HH3441qaGlUiESrsSqSTRHTt2cPjw4aDmcyD0wp8lJSWkpqba9t6C8XSSk5M9PugtOnToQEpKil9Px1eODtgrheO+rIFFMMNrlu3ePJ3ExETS0tKajvkSnbZt2zZbCTUS2PV0HsFR8flhHHM4qc6fDwM3OX8qCgCffvopI0aMICEhgaVLlzJ6dMuCFVdddRXjx4/n97//PYsWLYq+kRHCCpUeMWJERPpPSUnhrLPOCrvohBJEAI6HZXZ2dtC5Ov6WNHAnJyeH8vJyW4ulAU3r6PjqX0RsJYj6Ex07a+q4L2tgkZqaSkpKSkCeDnivNG0t4Ob6vi1PxlMgQaS9HLAvOmOBx40xTxpjthljjjp/Pgk8AfwsciYqxxN//OMfufTSS8nNzWX58uWcccYZXtv+4Q9/oF+/fowbNy7olSDjiUiFSrsTiSRRK/fH1+/LH6FEsNkNl7YINGza2+Jt7oRDdCwPxtvwWn19PYcPH/bqdfmqNO3N0/FWadp11VALX55OPInOSTgKfXpimfO8cgJjjOHJJ5/k5z//OSNHjmTJkiX07NnT5zXp6enMnj2bvXv3MnHixKDXvY8X1q1bx7Zt2yI2n2MRiSTR4uJi+vTpE1JZ+/z8/KBEx86SBu4EIzq+5nMs7NRf8yc6SUlJdOzY0aunY/2D5U10vJXCqa2tpba2NiBPx3XVUAtvohONVUPBvujsxpFs6YlhzvPKCUp9fT233XYbDz30EOPGjWP+/Pkey4N44uyzz2bKlCn8+c9/5p133omwpZElUqHS7kQimCCUIAKLvLw8Dhw44HflS3cOHDhAVVVVxETHKpljR3TslMKxKkx7Ex3wXZXAWzUCC2+ejqcSOBaBeDrp6emISNx7OrOBB51Ra31FJFVE+jgXRHsQR36McgJSWVnJ1VdfzRtvvMGDDz7IrFmzAp6IvO+++xg+fDgTJ05k27aYROeHhaKiIk499dSI5x717NkzrEmiNTU1bNiwIej5HItggwkCjVwDx5ICSUlJtkSnvLyc6upq255OaWmpz8TO/fv306FDB59DqL6qEvgTHW+ejqcK0xa+Agnc/wFMSEigXbt2cT+nMwWYhyNqbSNQCWwCnnQ5rpxg7Nu3j1GjRrFgwQJee+01nnjiCdsTwa4kJiby9ttvY4xh3LhxNDQ0RMDayFJZWckXX3wR8aE1CP9KoqtXr6axsTEsng5ER3QCWczNTo6OhRU2bQ2hecJXjo6Fr6KfwXo6nuquWWRmZlJeXt6iUKgnTwc811+LK9ExxtQbY34KnAHcjiNa7XbgdGPMGGPM8feUUEJi/fr1FBQUsHbtWt5//33Gjx8fUn99+vThpZdeYsmSJfz+978Pk5XR4/PPP49oqLQ7BQUFbN261XZVZF8Es3CbJ/r06RNU4U+7Sxq4YzdsOhjR8TWvY1d0/A2veQqZto4HOrzWoUMHjDEtvBdPng4cB6JjYYxZbYx5xRnF9ooxRqs4n4B8+eWXDBs2rGnNmCuuuCIs/Y4bN44bbriByZMn8+9//zssfUaLoqIi2rVrF7FQaXfCOa9TXFxMenp6QJ6GJ5KSkjj55JODEp3u3bs3yyWxQ6CiYyd6zU4pHDuiYw2veQqO8basgUWnTp2oqqpqMcTny9PxVn/NfdVQi7gUHRHJEZEkl9c+t4hbqsQF7733HhdccAGdO3dm+fLlnH322WHrW0R49dVX6dq1Kz/96U9bVKWOV4wxFBUVccEFF0Q8sc7CShINRx224uLioMvfuJOXlxdwrk6gkWsWOTk57Nq1i/r6ep/tAhEdO6Vw7Ho6x44do7KyZfnH0tJS2rRp4zVS0PJk3Od1/AUSQPP6a1YlaU/3cV/IzRjTYoXRSOHrr2wLcKbz9Vbnvq9NaeVMmzaNG264gcGDB7Ns2bKQ/zP2RKdOnZg1axbr16/nf/7nf8LefyRYu3ZtVEKlXQlXkqgxJqiF27xhFf70JwSuhCI6DQ0NfkOcd+3aRadOnbwWnXXFEhNvfTY0NHDw4EG/lRN8lcKxqhF4m//0VgrHXyCBaxvwXALHwt3Tqa6uprGxMSqejq+Cn/8FbHZ5fXwnUShB09jYyL333svzzz/PtddeS2Fhoa0vcLBccMEF3H333fzhD3/g8ssvj+rDPBiiFSrtTkFBAa+99lrAyza7smvXLsrKykIOIrDIz8+nrq6OrVu30q9fP7/tjx49ys6dO4MWHXCETVvLHXjC3+JtriQnJ9OpUyevno41ZGZneM1q7/7evNVds/BWCqesrIzU1FSPUXOePB1Pa+lYuItOtOqugQ9PxxgzyxhT6nw907nvdYu4pUpMqK2t5Sc/+QnPP/88d9xxB3Pnzo2o4Fg89dRTnHHGGdx6660+I4nigaKiIk477TSfD75IEI4k0XAFEVgEGsFmrSQbquj4wm6OjoWvXB1/iaEW/jwdX6Ljy9Px5OXA996Mq+h4WkvHIi5FxxXniqH5Xs4NEJHPw2uWEg8cOnSIiy66iHfffZepU6cybdo0EhOjU9s1OTmZ2bNnc+TIEX75y1/GbbWCiooKvvjii6h7ORCeYAKr5loo5W9cCVR0ggmXtrC7mFugouOrFI5d0fFV9DMUT8eb6FjHPQ2vefN0Kioqmr5XcSc6OFYM9TbDlAGM9HJOOU7ZsmULw4YN41//+hdz5szh7rvvDioHJxTOOOMMfve73/Hhhx/y+uuvR/Xedvn888+pq6uLyRBgr1696NGjR8iik5uba7uChD+ysrLIysqyHUwQiui0a9eOjh07+hSd+vp69u3bZyuIwMJXKZxAPR1vouMtXBp8ezqeggggcE8nIyODxsbGpmCdeBQd8D6nczKOZFGllfDNN99QUFDA/v37+fTTT7nxxhtjZssdd9zBRRddxF133RXScsiRwgqVHj7cW5WoyBJqkmg4gwgsAqnBVlJSQkpKSlPUWKD4C5veu3cvjY2NQXk6nrxru6KTkZFBUlKSx+E1bxWmLSyvJRBPp23btqSkpATk6cD3whQXoiMit4rIFyLyBQ7Bed3ad9m+BmYBSyJuqRIVioqKGDlyJCkpKXz55Zecd955MbUnISGBmTNnkpqaytixY8NaWTlUrFDpCy+8MGqh0u6EkiRaW1vL+vXrwxZEYBFItemSkhL69OkTdLi2P9Gxs3ibO927d6e2trZFHgs4RCcxMdGrx2EhIh5L4VRXV1NbW+tTdBITE+nQoUML0fHl6UDL+mv+5nRc28SF6ACNQINzE7d9aysFXgF+EVkzlWgwY8YMrrzySgYMGMDy5cs55ZRTYm0S4MiveP3111mxYgWPPho/FZfWrFnD9u3bYxpdF8q8jlX+JtyeTl5eHvv27fNY9didYMOlLfyJTiDVCCx85epYOTp2RNJTKRx/JXAsPNVf8xVIAC0rTR93no4zKu2HxpgfAouBMda+y3aJMeZuY8y+iFuqRAxjDA8//DC/+tWvuOiii1i8eLHXtd1jxXXXXcett97K008/zdKlS2NtDhC7UGlXzjzzzKBXEg114TZv2A0mCGZJA3dycnI4fPiwR68EQhMdT/M6+/bt8zu0ZuGpFI5d0XGvv9bY2BiUp5OSkuLRC3dfyC0uRMcVp8AEtySgEtfU1dVx66238vjjj/OLX/yCDz74ICp/eMHwwgsvkJuby7hx47w+ZKJJUVERp59+ut91gyJJKEmiq1atIi0tjZNPPjmsNuXnOwJd/YlOaWkpFRUVIYsOOJbb9sSuXbtISkpqmti3gx1Pxw6ehteC9XQqKytpbGz06+m45+l4CxDx5um0a9fOp13hwG7I9M/8bZE2VAk/5eXlXH755cyaNYtHH32UN954I+gkw2iQkZFBYWEh27dvZ9KkSTG1paKigiVLlsTUy7GwVhL1VY7fE8XFxZx++ulhD4Pv27cvbdq08Ss6oUSuWfjL1dm1axfdu3cPaM7IV/21QEQnnJ6Or7prFpmZmS2G17yVtfEkOunp6WEpheQPu3eY6WV7y2VTjiN27drF+eefz8KFC3nzzTd5+OGHox4SHQwFBQU89NBD/PGPf2Tu3Lkxs+Ozzz6LWai0OwUFBdTW1gaUJGqMCcvCbZ5ISkqib9++cSM6gQytgeOBn5SUFBZPp6ysrFlJoGA9HV911yzcPR1vFaahpehEa9VQsC86fTxsQ/h+fZ1zI2KdEhFWr15NQUEBmzdv5u9//zu33nprrE0KiIceeohzzjmHX//61+zcuTMmNhQVFZGRkRGzUGlXrGCCQIp/7t69m9LS0rDP51jYKfxpiU6fPn2Cvk/37t1JTEz0KjqBlMCxEBGPuTpVVVVUVVUF5OlA89Bn67WvPB3r/KFDh5rCtu14Op4CCbx5Op7mdOJKdIwx2zxs/zbGPAb8Cbg7smYq4WLRokUMHz6c+vp6vvjiCy6++OJYmxQwSUlJFBYWcuzYMX7+85+3WLgq0riGSsfDcGQwSaJWEEEkPB1wzOts2rTJ54J8JSUldO3alfT09KDvk5iYSM+ePcPq6YDnqgTWUJm/Yp8WnkrhlJaWkp6e7nPVUXB4NPX19VRVVQH2PJ3MzExqa2s5evQo4NvTSU5OJikpqdnwWlyJjh+WAJeHoR8lwvzpT3/i4osvblrq+Mwzz/R/UZzSv39/pk2bxueff860adOieu/Vq1ezY8eOuJjPsQg0SdQaigtX+Rt38vLyOHr0qM/lx0ONXLPwFjZdUVFBRUVFUKLjqf7avn2OIN1AhtegeVUCfyVwLCxPyPKMfFWYtnCvSuDL0xGRZvXXjjfRGYpWJIhrjDE8++yz/PSnP2Xo0KF8+eWX9O7dO9ZmhcwvfvELrr76au6///6Qil4GSjyESrsTaJJocXExOTk5Ph9ioWAnbDrSohPIOjruePJ07FYjsPBUCseu6LjXX7OG1/x5OvC96PjydID4Fh0RedjD9oSI/A3HvE7sZnQVnzQ0NHD77bfz29/+lptuuolPPvnEbzb18YKIMGPGDDp16sSYMWOora2Nyn2Lioo444wzYhoq7U6gSaKRCiKwsETH27zOsWPH2LFjR9hEZ+fOnS2G8oLJ0bHo1q0bBw4caBYEEKzouA+vBeLpWGJz+PDhJu/EG66eTmNjI+Xl5T7bZ2RkxK/oAFM8bPcA+cCTwH3hNkwJnerqaq677jqmT5/O//zP//DOO+/4HUs+3ujcuTNvvfUWq1ev5v7774/4/crLy1m6dGlceTngWEm0bdu2tkSntraWdevWRSyIABy/l06dOnn1dLZv305jY2PYRKe+vr6FZxJMCRyL7t27Y4xptqyG9doaNvNHuD2d9u3b+wxpdq00XVVVhTHGr6cTz4EECR62VGNMvjFmijHmaKQNVQLjwIEDjB49mg8++IAXX3yRZ599Niox+LHgkksuYdKkSUybNo1PPvkkoveKp1BpV5KTk20nia5du5aGhoaIejoi4rMGWzjCpS28hU2H6ulA81yd/fv3065dO9LS0mz10bZtW9q3bx82T8ffCIWrp+OrBI5FXA+vKccXmzZtYtiwYaxatYr33nuP22+/PdYmRZxnnnmGU045hVtuuaUpFyISFBUV0b59e4YNGxaxewSL3STRcC/c5o14EJ3MzEzbIuGKp1I4geToWLhWJWhsbKSsrMxvuDS0DCTwVWHawnXJal9LVVtYolNfX09NTY1PgQonvqpM5wSyRcVaxS///Oc/KSgooKysjM8//5xrrrkm1iZFhdTUVN555x0OHjzI+PHjI7LoW7yFSrtjN0m0uLiY1NRUW8tJh0JeXh579uzxWLKopKSEtm3bBjXJ744v0QnGywHvno7dcGkL16Kfhw8fxhhjy9NJS0sjKSkpIE/HNZDA+szteDqVlY44sHjwdLYCWwLYlBjz/vvv88Mf/pD27duzbNmypsnlE4VBgwbx5JNP8pe//IWZM2eGvf/vvvuOnTt3xt18jsXQoUMB/8EEkSp/446vGmyhLmngSvv27enQoYNH0QlW1DyJTiDFPi1cS+HYrUYAjuFJK0EU/FeYBodoiIhtTycjI6MprNzajwa+fuP/FeCmxJDp06dz7bXXcvrpp7N8+XIGDBgQa5Niwt13382oUaO444472Lx5c1j7jsdQaVfsJIkaYyKycJsnfIVNhytc2sJT2HQonk5KSgqZmZktPJ1QhtcCER1oXn+trKzMr6eTkJBA+/btm3k6/obXqqurm7ypaIlOG28njDEzo2KBEhKNjY088MADPPPMM1x55ZX86U9/CinD+3gnMTGRWbNmMXDgQMaNG8cXX3xBmzZe/8wDoqioiIEDBwb9IIsG/pJE9+7dy8GDByMaRGBx8sknk5iY2EJ0jDFs3rw5rPNi7qLT0NDA3r17Q/pduZbCaWxs5MCBA0F5OsGKjmv9NTueDnxff81uIAF8H+UXD55OC8TBaSJynoicKsdDhchWzNGjRxk7dizPPPMMv/71r/nLX/5yQguORU5ODq+++irLly/nqaeeCkuf8Roq7U5BQQHbtm3zuBYMRC+IABzRW3369GmRq1NWVkZ5eXlEPZ39+/fT0NAQsuhYns6hQ4dobGwMSnRqamqoqqoK2tOpq6ujqqrKluhYlabtejoQx6IjIr8E9gDFwCLgW2C3iOiqoTHg8OHDXHLJJfzpT3/i6aefZvr06WH7j7418JOf/IQxY8bw2GOP8c9//jPk/j799FPq6+uPC9EB78U/I7Vwmzfy8/NbeDrhjFyzyMnJ4dChQ02T4qGES1u4lsIJNDHUwrUUTrCejp26axauno6I+FwfxxIZ67OKK9ERkTHA6ziE5r+Ay5w/vwVeF5GbI2ah0oLt27czYsQIvvzySwoLC7nvvvuOi2UJos3LL79Mjx49GDt2bNPDKFjiOVTaFX9JoqtWraJXr15Rq0qRl5fHxo0bmxVljZTowPeLuYVDdFw9HUt0goleA0feXGlpKQkJCT69D1csT8dOhWkLq9K0tVSBr0ANy9OJS9EBfgPMNsZc5FzG+mPnzx8B7wC/jZyJiisrV66koKCAHTt2sGDBAsaMGRNrk+KWDh068Pbbb7N582buuuuuoPuxQqUvuuiiuAyVdsVfkmhxcXHUvBxwiE5tbW2zoa9wLGngjnvYdCh11yy6detGVVUVFRUVQXs6rlUJDh06RMeOHW1H7HXq1Iny8vKmOSE7/yhYS1b7KvZpEe+ikwcUejlX6DyvRJh//OMfnH/++SQkJLB06VJGjx4da5PinvPPP5/f/va3zJgxg7/+9a9B9fHdd9+xa9euuB9as/CWJHr06FHWrVsXlSACC0812EpKSsjOzg7rQ86T6CQmJgYsEq64hk0HWmHawn14ze7QGnyfILp161YgME/H11LVFvEuOhWAt+qGPZ3nlQgya9YsLrvsMnJzc1m+fHnEStK3Rh599FHOOussfvWrX3mdYPfF/PnzgfgNlXbHW5Lo2rVrqa+vj6qn4ylXJ9zh0uDwaBISEppEZ/fu3U0LvAWL67LV+/fvJyEhwVY1AVfch9cCER3Ls7E8Q7ueTnl5uS1Px3VOJzExkZSUFNu2hYJd0SkCnhKR81wPikgB8ITzvBIBjDE88cQT3HLLLYwcOZIlS5bEVXXj44G2bdsye/ZsqqurufXWWwOuVlBUVMQPfvCDsGTPRwNvFacjvXCbJ7Kzs8nMzIy46LRp04YePXo083RCDW139XT2799P586dAxaxzMxMEhMTQ/J0LNGx6+k0NDSwe/du257O/v37mxJLo0EgczpHgEUisl1E/iki24ClQLnzvBJm6uvrGT9+PJMnT2bcuHHMnz/f9iSk0pz8/Hyee+45Pv74Y1566SXb1x05coQvv/zyuPFyAHr27EnPnj1biM6qVatISUmJePkbV9wLf9bX17N9+/awiw40D5sOp+js2bMnqMRQcLx/K1cnVE/Hbsg0OIYZ7Xo67q8jjd0q03uBQcB/A8txCM1XwCTgTGPMvkgZeKJSWVnJVVddxYwZM3jwwQeZNWsWbdu2jbVZxzUTJkzgsssu4ze/+Q2rV6+2dc3xEirtjqck0eLiYk477bSoh9bn5eU1zens2LGDhoaGqIhOqJ5pVlYWbdq0afJ0gp0fskrhhOLptG3bltTUVL/XWP+UVldX+/0HNTExsSmvL+5EB8AYU22MeckYc5Mziu0mY8x0Y0x1JA08Edm7dy8jR47k448/5rXXXuOJJ57QkOgwICK8+eabZGRkMHbs2Ka15H1RVFREhw4djrs6dp6SRCO9cJs38vPz2b17NxUVFREJl7bIyclhx44dVFRUcOTIkZA9nYSEBLp27dokOoGGS1tkZ2ezc+dOqqurg/J0du7cSWZmpq1ngKvQ2BkVsbyhuBMdEensXklaRG4TkRdF5IrImHZism7dOgoKCli3bh3vv/8+48ePj7VJrYquXbvyf//3f6xcuZLJkyf7bGuMYcGCBcdFqLQ77sU/rQdnNIMILKwItg0bNkRcdOrq6vjPf/4DhJajY2Hl6oTq6VjDi4EEIliiY4yxnVflOgRnZ6kCS2ziTnSAN3FZHVREJgOvAD8F3heRmyJg2wnH0qVLGTZsGNXV1SxatIgrrlA9jwRXXnklt912G8899xwLFy702u7bb789rkKlXXFPEo1FEIGFa+HPkpISkpKSIlK/zgqbtt5zuERny5YtlJeXhyQ6VlmaQDydpKSkJjGwM58DrcjTAYYAn7ns/xp4yhiTBbwM3B1uw0405s2bx4UXXkh2djbLly/n7LPPjrVJrZqpU6fSv39/fvaznzVlfLtjhUpfcskl0TQtLFhJolY5HCt8Ohah9v369SMhIYF169ZRUlJCbm5uRJZViJTobNiwAQg8R8fCdXnrQEQHvvd2IuXpxLPodAL2AYjI6UA3YJbz3N/Q5NCQeP7557nxxhsZPHgwy5Yti8jQg9Kc9PR0Zs+ezd69e5k4caLHMOqioiIGDRp03IRKu+OaJFpcXEyPHj0CfuiFg+TkZPr06dPk6UTq79tddMLxe+vevXtTCZ9QPB2LQD9/azjuRPR0Svk+OXQ0sNsYs9G5nxRAPzFBRE4WkaUiskFE/iMiQ2JtEzjKpd91113cfffdXHPNNXz66acxeSicqAwZMoQpU6bw5z//mXfeeafZueMxVNod1yTRWAURWFhh05EUnQ4dOpCRkdGUdxKOB6kVNg2x9XTsik5KSkpTlOvxPqfzKTBFRG4H7sHh3VjkA9vCbFe4eRWYaYwZgLOOXKyXZaitreWmm25i2rRp3HHHHcydO9dWSKQSXu677z6GDx/OxIkTm8qNgKPkUENDw3EvOgCLFy9m7dq1MQkisMjLy2PNmjUcOnQoYqIjIk3eTrjmjFxFJ9jotXB4OnaH10SkycMJxNOxI1DhIpDk0B3A08Bm4FGXc2NwJImGDRHp6YyMWy4i1SJiRCTXS9teIjJPRI6ISLmI/MU10k5EsoGhOIcDjTH/cJ4aHE6bA6G0tJQLL7yQefPmMXXqVKZNmxbxpYMVzyQmJvL2229jjOFnP/sZDQ0NwPEbKu2KlST61ltvUVdXF3NPp66uDohM5JpFuEXHKoUDoQ+vpaamBvyPZaCejmvb43pOxxizz5mbk2GMGW2MOehy+kIcSaPhpB9wI1AGLPHWSETSgM9xeFs/B8YB/YGFImKtZpaDYziwzuXSbc7jUWfLli0MHz6cr7/+mjlz5nD33XdrDk6M6dOnDy+99BJLlizh97//fVOo9I9+9KPjfo2igoIC1qxZA0RvDR1PWDXY4PgSHcvTSUtLC3qBRGt4LdC6ba7XBLIURTCeTjRFB2NM3G1AgsvrXwIGyPXQ7r+BBqCfy7E+QD1wt3N/MLDe7bp/ANf6s2Pw4MEmnHz99dema9euJjMz0yxevDisfSuh0djYaG688UYjIiYrK8sAJisryxQWFsbatJAYM2aMcX5/TE5OTszez8svv9xkR69evSJmxw033NB0n969e4d8n//7v/8Lub8333wzqD4KCwtNZmamAUznzp1tXVdYWGhSUlJsfc6FhYWmY8eOAfVvF2CF8fZ893YiXjY/ovMZ8KWH44uBxc7X2TiqYCe5nN8ADPF373CKzkcffWTS0tJM7969zZo1a8LWrxI+XnvtNSMiTQ8IwKSlpR23wuP6AIrl+yksLDRpaWkRt6OwsNC0bds2bPcJh93B9hHMdYFcE+nfSWsWnb3Aax6OTwcOuOx/BvzK+foiYCMg/u4dLtF5/fXXTWJiojnzzDPN7t27w9KnEn569+7d7Evo+t/p8Ui8vJ9o2RHu+4Sjv2D7COa6QK6J9O+kNYvOMeB3Ho4/AdS77PcHljk9nJXAOT7uNx5YAazIyckJ+MMuLCw0vXv3NiJicnJyzI9//GMDmEsuucSUl5cH3J8SPdy9HGsTkVibFhTx8n6iZUe47xOO/oLtI5jrArkm0r8TX6LjNZBARK4SkQ7ezscRxsOxZjPzxpiNxphhxpgBxphBxph/ee3MmNeNMUOMMUNc4+vtMHv2bMaPH8+2bdswxrB9+3b+9re/MXLkSD744IPoTtYpAWNNQts9Hu/Ey/uJlh3hvk84+gu2j2CuC+SamP5teFMjHBP057i/jvaGb09nHzaG14LdAh1e8+ayBuMxKdEnWnMP0SJe3k8053TCeR+d0wkeghlew7Fo20XO141xKjqfA0s9HF+EM5AglC1Q0YmX4QwleFyHR8MR/RRr4uX9RMuOcN8nHP0F20cw1wVyTSR/J75ERxznWyIinwO5wBfAz4CPgAPeHSbzCy/nQkJEfgm8AfQxxmx1O3cn8BwwwBhT4jyWiyNQ4D5jzNRQ7j1kyBCzYsUK2+1zc3PZtq1lcYbevXs3y3ZXFEVpzYjIN8YYj+XGfGW+TQCeB87H8R/7OTgm7j3hWblCQESud760KgdcKiIHcAybLXYeewO4HcfyCg857XgcR/WE18Jtkz+efPJJxo8fT3X19+vapaWl8eSTT0bbFEVRlLjEq+gYY9YDlwGISCNwpfExAR8B3nXbn+78uRgYBWCMqRKR0TjE8W0cAQSfAXcaYyqjZGcTY8aMAeDBBx9k+/bt5OTk8OSTTzYdVxRFOdHxOrzWrJHISOCbWDzIY0mgw2uKoihK8MNrTVjDWc61dEbiWF+nFPjCGPNduAxVFEVRWje2REdE2gAzgZtpngNjROQd4BZjTEP4zVMURVFaE3aXNngER9Xnh3EU1Ex1/nwYuMn5U1GUOGbnzp1MmjSJgoIC0tLSEBGPUZWLFi1CRLxuhw8fDst9lBMTu3XbxwKPG2Ncw7C2AU+KSCJwKw5hUhQlTtm0aRNz585l8ODBnHfeeXzyySc+2//v//4vZ599dovj/iprBHof5cTCruicBCz3cm4Z8GB4zFEUJVKcf/757Nu3D4AZM2b4FYNTTjmFoUOHRvw+4eLo0aMkJye3OG6Moa6urmkZ53D2rQSO3eG13cBwL+eGOc8rihLHJCTY/brH9j4HDx5kwoQJ9OjRg+TkZPLz83n99debtZk5cyYiwhdffMENN9xAZmYm5557LuBI0h47dixvvvkm+fn5tG3blo8++giABQsWUFBQQGpqKh06dODHP/4x69evb9b3qFGjGDFiBB9++CFnnnkmycnJTJ8+HSU82PV0ZgMPOvN1ZgN7gG7AT3B4Oc9ExjxFUWJFY2Mj9fX1zY6JSESXVi8vL2f48OHU1NQwZcoU+vTpw8cff8yECRM4evQokyZNatZ+zJgx3HzzzcybN6+ZrQsXLmTlypU88sgjdOnShdzcXBYsWMDll1/O6NGjmTNnDpWVlTz88MOMGDGClStXNlttdMOGDdxxxx1MnjyZvn37BrXqp+IZu6IzBegLPOp8bSHAn5zHFUVpRVx88cUtjp122ml8913ksiReeOEFtm3bxrfffkv//v0BuPDCCzl8+DCPPvooEyZMaLaE+PXXX8+zzz7bop+ysjK++eabpuWmAW666Sb69u1LUVFRUx8FBQUMGDCAqVOn8oc//KGp7cGDB/nkk08YNGhQhN7piYvdPJ164Kci8iSOsjidgEM4imquiaB9iqLEiJdffplzzjmn2bHU1NSm1+5eUGJiIiLNVhUJmAULFnDuuefSp0+fZv1ffPHFzJgxgzVr1jBw4MCm49dcc43HfoYOHdpMcKqqqvj3v//NAw880Ey0+vTpw/Dhw1m8eHGz63Nzc1VwIoRdTwcAY8xqYHWEbFEUJY4YMGAAQ4Z4TCoHICkpqdn+woULGTVqVEj33L9/P5s2bWrRt0VpaWmz/e7du3ts5368rKwMY4zH9t26dWtRqNdbv0roBCQ6iqIoFl9//XWz/by8vJD7zMrKokuXLrzwwgsez7vfw5tn5X68Y8eOiAh79+5t0Xbv3r1kZWXZ6lcJHRUdRVGCwpcXFCyXXHIJL774Ijk5OXTp0iVs/aanpzN48GDeffddpkyZ0hQMsW3bNpYtW9YiQEGJHCo6inICMW/ePAC++eYbAIqKisjOziY7O5uRI0c2a7t27VratWvXoo8zzjiD9PT0sN3Hlbvuuos5c+Zw3nnncdddd5GXl0dVVRXr1q1jyZIlvP/++/bfrBuPP/44l19+OVdccQUTJ06ksrKSRx55hA4dOnDPPfcE3a8SIN5Wd9Mt8JVDFSXewcPKtoAZOXJkU5uFCxd6bQeYr7/+Oiz38cahQ4fMnXfeaXJzc01SUpLJzs42I0aMMM8//3xTm7feessAZuPGjS2u7927txkzZozHvouKiszQoUNNSkqKad++vbnqqqvMunXrmrUZOXKkGT58uF87Fe8QzMqhFiLSFseCbp+ZE6yitC5toCiKEji+ljbwmzpsjDkG/A5HmLSiKIqiBI3dehVrcSSHKoqiKErQ2BWdh4HJInJGJI1RFEVRWjd2o9d+C7QD/iMiW3HUXnOdDDLGGO8hKYqiKIqCfdFpALTcjaIoihISdmuvjYqwHYqiKMoJQHQW2FAURVEUAhAdEekuIs+JyNcisllE/iUiz4pIN/9XK4qiKIrN4TURGQAsAToCXwKbcCzi9t/Az0TkPGPMxohZeRziqdrujTfeyMSJE6muruayyy5rcf6WW27hlltu4eDBg1x//fUtzk+YMIGbbrqJHTt2MG7cuBbn77nnHq688krWr1/Pbbfd1uL8Qw89xIUXXsjKlSu58847W5x/6qmnGDZsGMuWLeOBBx5ocX7atGkMGjSITz/9lCeeeKLF+ddee428vDw+/PBDpk6d2uL822+/Ta9evZgzZw6vvPJKi/Pz5s2jc+fOzJw5k5kzZ7Y4P3/+fNLS0pg+fTpz585tcX7RokUAPPfcc/z9739vdi41NZWioiLAUQ7ls88+a3Y+KyuL9957D4D777+f5cubr87es2dPCgsLAbjzzjtZuXJls/MDBgxoWt1y/PjxbNiwodn5QYMGMW3aNADGjh3Lzp07m50vKCjg6aefBuC6665rUU35ggsuYPLkyQBceuml1NTUNDt/xRVXcO+99wL6t6d/e+H527PeU7ixG0jwDFAOnGuM2WodFJHewCfO89eG3TpFURSlVeG3DA6AiBwGfm2M+bOHczcD040xHcNvXmzRMjiKoiiBE1IZHCdtgQov5yqc5xVFURTFJ3ZFZyUwSUSatRfHSkcTnecVRVEUxSd253QeA/4OrBWROTgqEnQDbgD6A5dHxjxFURSlNWE3OXSBiFwBPAE8CAiOMjjfAFcYYz6JnImKoihKa8Gv6IhIEnAZUGyMGSIiaThCp8uMMdWRNlBRFEVpPdhZT6cOmAvkOverjTG7VHAURVGUQLEbSFACdImkIYqiKErrx67oPAs8KCLZkTRGURRFad3YjV4bjWO56i0i8hWe19P5ebiNUxRFUVoXdkXnPKAOOACc7Nxc8V/WQFEURTnhsRsynRthOxRFUZQTAL9zOiLSVkT+LSI/ioZBiqIoSuvFTsj0MaAPUB95cxRFUZTWjN3otX8A6ukoiqIoIWE3kOBFoFBE2gB/o2X0GsaYkvCapiiKorQ27IrOYufPu4G7vLRJDN0cRVEUpTVjV3RujagViqIoygmB3ZDpWZE2RFEURWn92A0k8IqIJIhIp3AYoyiKorRuvIqOiBwSkbNc9kVEPhCRvm5Nz8ZRqUBRFEVRfOLL08mk+fBbAnCF87iiKIqiBEzIw2uKoiiKYhcVHUVRFCVqtHrREZEUEfmbiKwVkZUi8rGHeSlFURQlCvgLme7h8oBOdDl22KVNz7BbFX5eMcZ8DCAitwMzcKwRpCiKokQRf6Izz8Oxv7ntCwGspyMiPYHfAkOAHwCpQB9jzFYPbXsBzwMXOe/zKXCnMWa73fsZY2qBj10OfQXca/d6RVEUJXz4Ep1IVSHoB9wIfAMswUshURFJAz4HjgI/xyFsTwALRWSgMaYqyPtPAt4P8lpFURQlBLyKTgSrEHxhjOkKICK/xHv16l8BfYE8Y8wmZ/tiYCNwG/AH57F/Azle+jjTGLPD2hGR+4EBwAVheB+KoihKgNitvRY2jDGNNpteBXxlCY7z2i0i8iVwNU7RMcac5eX6ZojIvcB1wIXGmOrArFYURVHCQTxHr50GfOfh+Grg1EA6EpG7gZuBi4wxh/20HS8iK0RkxYEDWmhBURQlnMSz6HQCyjwcPwR0tNuJM3BhKo5KCgudYdMrvLU3xrxujBlijBmSnZ0doMmKoiiKL6I+vBYgnqLiJKAOjNkZ6DWKoihKZIhnT6cMh7fjTkc8e0CKoihKnBPPorMax7yOO6cCa6Jsi6IoihIG4ll0PgCGupasEZFcYLjznKIoinKcEZM5HRG53vlysPPnpSJyADhgjFnsPPYGcDvwvog8hGN+53FgB/BaNO1VFEVRwkOsAgnedduf7vy5GBgFYIypEpHROMrgvI0jGOAzHGVwKqNkp6IoihJGYiI6xhhb0WTOGmvXRdgcRVEUJUrE85yOoiiK0spQ0VEURVGihoqOopwg7Ny5k0mTJlFQUEBaWhoiwtatW1u0W7RoESLidTt8+LDP+3z88ceMHj2abt26kZycTM+ePbnxxhtZs0YzHZT4r0igKEqY2LRpE3PnzmXw4MGcd955fPLJJz7b/+///i9nn312i+MZGRk+rzt06BCDBw9m4sSJZGdns337dn73u98xdOhQvv32W3r37h3S+1COb1R0FOUE4fzzz2ffvn0AzJgxw6/onHLKKQwdOjTg+9x8883cfPPNzY6dc8455OfnM2/ePO65556A+7TD0aNHSU5ObnHcGENdXR1t27YNe99K4OjwmqKcICQkxO7rnpWVBUBSUpLftgcPHmTChAn06NGD5ORk8vPzef3115u1mTlzJiLCF198wQ033EBmZibnnnsuALm5uYwdO5Y333yT/Px82rZty0cffQTAggULKCgoIDU1lQ4dOvDjH/+Y9evXN+t71KhRjBgxgg8//JAzzzyT5ORkpk+fjhIe1NNRFMUjjY2N1NfXNzsmIiQmJtq6vqGhgYaGBrZt28Z9991Ht27d+MlPfuLzmvLycoYPH05NTQ1TpkyhT58+fPzxx0yYMIGjR48yadKkZu3HjBnDzTffzLx585rZunDhQlauXMkjjzxCly5dyM3NZcGCBVx++eWMHj2aOXPmUFlZycMPP8yIESNYuXIlPXr0aLp+w4YN3HHHHUyePJm+ffvSqZOnMpBKUBhjdPOyDR482ChKa+SNN94wgNmyZUuLcwsXLjQ4KoC02E477TTb9xg8eHDTdf369TNr1qzxe81jjz1mkpOTzYYNG5od/+Uvf2mysrJMXV2dMcaYt956ywDmzjvvbNFH7969TWpqqtmzZ08Le/r169fUhzHGlJSUmDZt2pi77rqr6djIkSONiJj//Oc/tt+r0hxghfHyXFVPR1EUj7z88succ845zY6lpqY2vXb3ghITExH5Pu/77bffpry8nJKSEp577jkuuugili5dSm5urtd7LliwgHPPPZc+ffo06//iiy9mxowZrFmzhoEDBzYdv+aaazz2M3ToULp169a0X1VVxb///W8eeOAB2rT5/rHXp08fhg8fzuLFi5tdn5uby6BBg7zaqQSPio6iKB4ZMGAAQ4YM8XrefX5m4cKFjBo1qmn/lFNOAeDcc8/l0ksvJTc3l9/97ne8+uqrXvvcv38/mzZt8jr3U1pa2my/e/fuHtu5Hy8rK8MY47F9t27d2LZtm61+ldBR0VEUJSi+/vrrZvt5eXle22ZmZtKvXz82bdrks8+srCy6dOnCCy+84PG8+z1cPStfxzt27IiIsHfv3hZt9+7d2xTo4K9fJXRUdBRFCQpfXpA7+/btY926dYwZM8Znu0suuYQXX3yRnJwcunTpEqqJTaSnpzN48GDeffddpkyZ0hQMsW3bNpYtW9YiQEGJHCo6inICMW/ePAC++eYbAIqKisjOziY7O5uRI0c2a7t27VratWvXoo8zzjiD9PR0r/e45pprOOussxg4cCDt27dnw4YNPP/887Rp08Zvjs5dd93FnDlzOO+887jrrrvIy8ujqqqKdevWsWTJEt5///1A33ITjz/+OJdffjlXXHEFEydOpLKykkceeYQOHTpELHdIaYmKjqKcQNxwww3N9idOnAjAyJEjWbRoUbNzd9xxh8c+vv76a59eztChQ5k7dy5Tp07l2LFj9OrVi1GjRnH//ff7DCIA6NChA8uWLeOxxx7jmWeeYdeuXWRmZpKXl8d114VWcP6SSy7ho48+4tFHH+XGG2+kbdu2jBo1imeffZaTTjoppL4V+4gjuk3xxJAhQ8yKFStibYaiKMpxhYh8Y4zx+J+JViRQFEVRooaKjqIoihI1VHQURVGUqKGioyiKokQNFR1FURQlaqjoKIqiKFFDRUdRFEWJGpocGiFcCx9a3HjjjUycOJHq6mouu+yyFudvueUWbrnlFg4ePMj111/f4vyECRO46aab2LFjB+PGjWtx/p577uHKK69k/fr13HbbbS3OP/TQQ1x44YWsXLmSO++8s8X5p556imHDhrFs2TIeeOCBFuenTZvGoEGD+PTTT3niiSdanH/ttdfIy8vjww8/ZOrUqS3Ov/322/Tq1Ys5c+bwyiuvtDg/b948OnfuzMyZM5k5c2aL8/PnzyctLY3p06czd+7cFuet5MbnnnuOv//9783OpaamUlRUBDgy0z/77LNm57OysnjvvfcAuP/++1m+fHmz8z179qSwsBCAO++8k5UrVzY7P2DAgKaFxsaPH8+GDRuanR80aBDTpk0DYOzYsezcubPZ+YKCAp5++mkArrvuuhaFLS+44AImT54MwKWXXkpNTU2z81dccQX33nsvoH97+rcXnr8992ThcKGejqIoihI1tCKBD7QigaIoSuBoRQJFURQlLlDRURRFUaKGio6iKIoSNVR0FEVRlKihoqMoiqJEDRUdRVEUJWqo6CiKoihRQ0VHURRFiRqaHOoDETkAbHM73AE4YuPyzsDBsBvVOrD7GcaSWNgYyXuGq+9Q+wn2+kCv0+9p6ITyu+5tjMn2dEJFJ0BE5HVjzHgb7VZ4y8g90bH7GcaSWNgYyXuGq+9Q+wn2+kCv0+9p6ETq71GH1wLnw1gb0Ao4Hj7DWNgYyXuGq+9Q+wn2+kCvOx7+xuKdiHyG6ulECP0PSlHiH/2eRh/1dCLH67E2QFEUv+j3NMqop6MoiqJEDfV0FEVRlKihoqMoiqJEDRWdGCEiJ4vIUhHZICL/ERGdzFSUOENEHhCR9SLSKCI/jrU9rQEVndjxKjDTGDMA+A0wW0QkxjYpitKcz4DLgC9ibUhrQUXHJiLSU0ReFJHlIlItIkZEcr207SUi80TkiIiUi8hfRCTH5Xw2MBSYBWCM+Yfz1OBIvw9Fac2E83sKYIz5pzFmc1SMP0FQ0bFPP+BGoAxY4q2RiKQBnwP5wM+BcUB/YKGIpDub5QC7jTF1Lpducx5XFCV4wvk9VSJAm1gbcBzxhTGmK4CI/BL4kZd2vwL6AnnGmE3O9sXARuA24A9ertOhNUUJnUh/T5UQUU/HJsaYRptNrwK+sv6QndduAb4ErnYe2g6cJCJJLtf1dh5XFCVIwvw9VSKAik74OQ34zsPx1cCpAMaYA8C/gFsAROQiHJ7ON9ExUVFOePx+T5XIoKITfjrhGE925xDQ0WX/18CtIrIB+D0wxmh5CEWJFra+pyLykIjsBAqAGSKyU0S6RcnGVonO6UQGT+LRbM7GGLMRGBYdcxRF8YCd7+kTwBPRMefEQD2d8FOG478odzri+T8rRVGij35PY4SKTvhZjWO82J1TgTVRtkVRFM/o9zRGqOiEnw+AoSLS1zrgTE4b7jynKErs0e9pjNClDQJARK53vrwARyDAROAAcMAYs9jZJh1YBdQAD+EYN34cyAAGGmMqo223opxI6Pc0vlHRCQAR8fZhLTbGjHJplwM8D1ih0J8BdxpjtkbaRkU50dHvaXyjoqMoiqJEDZ3TURRFUaKGio6iKIoSNVR0FEVRlKihoqMoiqJEDRUdRVEUJWqo6CiKoihRQ0VHURRFiRoqOkpUEZFbnOvWHxaRjm7n2jjPTYmBXVOc947ryusikiAi00Rkj4g0isjffLRt9lmKyI9F5O5o2OkNEblTRK71cHyKj6ROpRWhoqPEig7Ab2NtxHHI9cB/41iDaTjwGx9tC4AZLvs/BmIqOsCdQAvRwWFnQXRNUWJBXP9Xp7RqPgEmicg0Y8zeWBsTDUQk2RhzNMRuTnH+nOZvaWZjzFch3ssvYXpPGGN2AjvDYJIS56ino8QKa2GsB3018jbsIiIzRWSry36uczjp1yLytIjsFZEKESkUkTQR6SciH4tIpYhsEpGfe7nlKSKyUESqnUNYj4lIs++JiHQWkVdEZJeIHBWRdSIy3q2NNYx4voi8KyKHgX/6ea+XiMhyEakRkSMi8jcRyXM5vxWY4txtcPZ/i4/+mobXRGQm8HOgh/O4cfv8QnpPInK2iMxzrqxZIyLrReQpEUl1s783MMbFhpnOcy1+zyLSXkReEpHdTpvWi8hdIiIubUY5+7nK2fagiBxw/t4z3fr7bxFZ67SvTERWiMg1Pn4lSgRQT0eJFXuAl4A7ReQ5Y8y2MPV7P7AIxwP2VOBZoBE4E3gDeA6YALwlIiuMMavdrv8b8CbwNHAxMNl5/RRwPAiBL4FU57EtznavOP/rf9Gtv9nAn3AMi3n9vonIJcBHwOfATUA74DFgqYgMMsbsAq4B7gBu4fuhqM02PhNwVFDOBs4GrnIeOxrG95QDrARmAhU41qp5GOgL/MTZ5hpgPo7qzlOcxw54MtYp9B8BZzn7+Ra4HPiD83084HbJC8DfgZ8CeTh+7w04/g4QkTHAVByf6RLnex2I54XclEhijNFNt6htOB6YBuiH4wt/GHjTea6N89wUl/ZTHH+mLfqZCWx12c91Xvu5W7u/OI+PdTnWEagHHnG/D3Cf2/Vv4HiIZjr3JwO1QH8P7Q4Cbdze5/M2P5cVwEbreuexPkAd8AeXY094+jy89On+Wc4EdnpoF9b3hKNicxtgLA7BznI5txUo9HBNs98zcIXzXre4tZuBQyw7O/dHOdvNcmv3kvM9icv+v2P996+b0eE1JXYYYw7h+O/zZ67DSCFS5La/zvnzY5f7lgH7gV4erp/rtv9nHF7H6c79S3AMKW0RR7RdG2fE28dAFg7vypW/+jNYHGu7nAXMMcbUu9i5BYcHMtJfHyES8ntyDoU9IyKbcYhCHfA2DgHqH4RN5+MQrD+5HS8E2tIy6OAjt/1vgWSgq3P/a2CQiLwoIheKSFoQNilhQIfXlFjzPDAJx7DHmDD0576+/TEfx1M8XL/Py34P588uOLy0Oi/3z3Lb3+OlnSsdcTycPbXdi2MeJJKE4z29BVyIYyhsJVAFnAO8jOfP2R+dgEOmZZDCXpfzrhxy27eus+79R+frX+BY1K1OROYDdxtdPyeqqOgoMcUYUykiT+PweH7voUktgIi0NcYccznu/iAMF12BErd9gF3On6U4vKT/9nL9erd9O7knZc523Tyc6+a8ZyQJ6T2JSApwNY6hvBdcjp8Rgk2HgE4efu/WZxTQZ2IcY2yvAa+JIz/sRzj+5uYA54ZgpxIgOrymxAPTcTzUn/BwzgowsIa3cEYlDYuQLTe67f8EqAS+c+4vAPKB7caYFR62ikBvaIypAr4BbhCRROu4iPTG8T4XB/NGPHAUxwS6O6G+p2QgkZae0i0B2ODOYhzPpxvcjo/B4aUGHQ5ujCkzxszBMZR6ur/2SnhRT0eJOcaYoyLyGPC6h9NFwBHgDRF5BMcD7jc4hCAS/MoZOfU1jgiuX+L4D/6w8/zzOKLLlojI8zi8gHQcD+3zjDFXB3nfyTjmJf4uItNxzCM9iuO9Tw2yT3fW4PAeJuAIXKg1xnxLiO/JGHNERL4C7hGRPTiCD/6L74ck3W04T0SuwDFUdtDL8FYRsBR4VUSygdXAZTh+H08bYw4G8sZF5HUcASHLcXh1A4BxOPLFlCiino4SL7yFI3qrGc6H/RU4JpXn4ghlfhFYGCE7rgYuAj7AEX31BI5wY8ueIzi8j/k4Kip8jCPE+upQbDLGLMAREpyJ432+CqwFRhhjdgfbrxszcARGPAX8C/jQee9wvKebcXhrL+OIktuL5+G6+3GI2lwcwj7FU2fGkfh6OTDLadNHzv278ZPb5YUvgcE4vOp/OPsoxBlSrUQPK5xQURRFUSKOejqKoihK1FDRURRFUaKGio6iKIoSNVR0FEVRlKihoqMoiqJEDRUdRVEUJWqo6CiKoihRQ0VHURRFiRoqOoqiKErU+P/2zEx/+2Mf/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from math import exp\n",
    "from numpy import linspace\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# Set up lists for plotting\n",
    "delta = []\n",
    "it = []\n",
    "err1 = 0.1\n",
    "err2 = 1E-3\n",
    "\n",
    "# initialize iteration\n",
    "x = 0.5\n",
    "x2 = x \n",
    "N = 20\n",
    "for k in range(1,N+1):\n",
    "    x = exp((1-x**2))         # iterate value\n",
    "    delta.append(abs((x-x2)/x)) # compute difference\n",
    "    x2=x                    # update x2\n",
    "    it.append(k)\n",
    "    print('x= ',x)\n",
    "    \n",
    "# Make the graph\n",
    "plt.rc('font',size=16) # set the font size\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.plot(it,delta,\"-ko\")\n",
    "plt.xlabel('Number of iterations')\n",
    "plt.ylabel('Error of solution')\n",
    "plt.hlines(err1,it[0],it[N-1],linestyles='dashed',colors='k')\n",
    "plt.hlines(err2,it[0],it[N-1],linestyles='dashed',colors='k')\n",
    "plt.text(5, err1+0.1, '1E-1 error')\n",
    "plt.text(5, err2+0.001, '1E-3 error')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2.1 - Newton's method\n",
    "\n",
    "1. We can now return to our earlier example of the root of the function  $f(x)=x - e^{1-x^2}$, for which the simple relaxation method failed. For the Newton method, we need to add the derivative, which we have already derived earlier $f'(x)=1+2xe^{1-x^2}$. Modifying your relaxation method code to incorporate the expression for the Newton method.\n",
    "\n",
    "## Talking points\n",
    "\n",
    "1. What do you observe?\n",
    "2. How quickly does Newton's method find the right solution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x=  1.0187680487630895  delta=  -0.5187680487630895  deriv= 3.117000016612675\n",
      "x=  0.999878200121373  delta=  0.018889848641716504  deriv= 2.961781414435417\n",
      "x=  0.9999999950561345  delta=  -0.00012179493476153045  deriv= 3.00024357008081\n",
      "x=  1.0  delta=  -4.943865522880825e-09  deriv= 3.000000009887731\n",
      "x=  1.0  delta=  0.0  deriv= 3.0\n",
      "x=  1.0  delta=  0.0  deriv= 3.0\n",
      "x=  1.0  delta=  0.0  deriv= 3.0\n",
      "x=  1.0  delta=  0.0  deriv= 3.0\n",
      "x=  1.0  delta=  0.0  deriv= 3.0\n",
      "x=  1.0  delta=  0.0  deriv= 3.0\n"
     ]
    }
   ],
   "source": [
    "from math import exp\n",
    "from numpy import linspace\n",
    "\n",
    "x = 0.5\n",
    "deriv = 1.0\n",
    "alpha = 1.0\n",
    "for k in range(10):\n",
    "    delta=(x - exp(1-x**2))/(1+2*x*exp(1-x**2))\n",
    "    deriv=(1+2*x*exp(1-x**2))\n",
    "    x -= alpha*delta\n",
    "    print(\"x= \",x,\" delta= \",delta,\" deriv=\",deriv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2.2 and 2.3. - Newton's method\n",
    "\n",
    "Use Newton's method to calculate the inverse (or arc) hyperbolic tangent of a number $u$. By definition, $\\tanh^{-1}{u}$ is the number $x$ such that $u=\\tanh{(x)}$. To put that another way, $x$ is a root of the equation $\\tanh{(x)}-u=0$. Recalling that the derivative of $\\tanh{(x)}$ is $1/\\cosh^2{(x)}$ the next Guess in Newton's method becomes $$x'=x-(\\tanh{(x)}-u)\\cosh^2{(x)}.$$ \n",
    "\n",
    "2. Write a function that calculates $\\tanh^{-1}{(u)}$.\n",
    "3. Plot $\\tanh^{-1}{(u)}$ from -1 to 1.\n",
    "\n",
    "## Talking points\n",
    "\n",
    "2. What do you observe?\n",
    "3. Does your function $\\tanh^{-1}{(u)}$ give the right solution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEQCAYAAACnaJNPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAk10lEQVR4nO3deXxU9b3/8deHQAgJW1b2kIQdBBTDJgqIW7UutbXaa29btS6t7dVa21+vtnXpYu+9va32Ya/eq72tSxet1datLaIIKDtcBQmyZWEnZGFLQtb5/v6YgcYYIIfMzJmZvJ+PxzySzMzhfM73cZj3nPM95/s15xwiIiId1c3vAkREJL4oOERExBMFh4iIeKLgEBERTxQcIiLiSXe/C4iGrKwsl5eX53cZIiJxZe3atZXOuey2z3eJ4MjLy2PNmjV+lyEiElfMbHt7z+tUlYiIeKLgEBERTxQcIiLiiYJDREQ8UXCIiIgnCg4REfFEwSEiIp4oOEREEtDqsmoefWsrdY3NYf+3FRwiIglo0eb9PPLWVnokhf9jXsEhIpKASitryc1IVXCIiEjHlFTUkp+VFpF/W8EhIpJgAgFHWZWCQ0REOmjf4XrqmwIKDhER6ZjSyloAChQcIiLSESWh4MjPVnCIiEgHlFbU0qtHEgP6pETk31dwiIgkmNLKGvKy0ujWzSLy7ys4REQSTGllbcT6N0DBISKSUBqbA+w8cDRiV1SBgkNEJKHsPFBHS8ApOEREpGNKKyJ7RRUoOEREEkqk7+EABYeISEIpqawlPbUH/VOTI7YOBYeISAIprayJaP8GKDhERBJKaWUt+Vm9I7oOBYeISIKobWim/HADBRHsGAcFh4hIwjjWMa5TVYCZXWNmL5rZdjM7amabzewnZtbH79pERGKFguOjvgW0APcCnwAeB74KLDCzeNkGEZGIKgndw5GXGdng6B7Rfz18rnDOVbT6e7GZVQNPA3OBhb5UJSISQ5Zuq2TswD70Sk6K6Hri4tt6m9A4ZnXo55Bo1iIiEosqaxpYvb2aSyYMjPi64iI4TmBO6OeHvlYhIhID3txYjnNw8YQBEV9XXAaHmQ0BfgC86Zxbc4L33Gpma8xsTUVFewcsIiKJY37RPoam92L8oL4RX1fcBYeZ9QZeBpqBG0/0PufcE865QudcYXZ2dtTqExGJtiP1TSzdVsUlEwZiFpnJm1qLl85xAMwsBXgFKADmOOd2+VySiIjvFm2uoLElEJX+DYij4DCzHsCLwDTgQufcBz6XJCISE+YX7SMzLZmzh6dHZX1xcaoqdK/G74ALgKuccyt8LklEJCY0NLewaHMFF40fQFKE5hhvK16OOP4L+CzwY6DWzGa0em2XTlmJSFe1rLiKmobmqJ2mgjg54gAuDf38LrC8zeNmv4oSEfHbq+v2kJacxDkjM6O2zrg44nDO5fldg4hIrNl98CivvL+Hf54xnJ7dI3u3eGvxcsQhIiJtPLmkBIBbZhdEdb0KDhGROFRV08Bzq3fwqbOGMKR/r6iuW8EhIhKHfrO0jIbmAF+ZMyLq61ZwiIjEmcP1TTy9vIxPTBjIyJzIThPbHgWHiEiceXb5do7UN3P73JG+rF/BISISR/YdqufxRcXMG5vDxKH9fKlBwSEiEkfuf2UDzYEAD1wxwbcaFBwiInHijaJ9zC8q584LRpObmepbHQoOEZE4UNPQzP2vFDF2YB9uPi/f11ri4s5xEZGu7j/+vol9h+t57PNT6JHk73d+HXGIiMS419fv5Znl27lpVj5n5UZn6PSTUXCIiMSwbfuP8O0/rePs4el85xNj/S4HUHCIiMSsmoZmbnt2LanJSfzX9VNI7h4bH9nq4xARiUEtAcfdf3yfsqo6fvvl6Qzsl+J3ScfFRnyJiMhxzjm+95cNzC8q57uXjWPmiOjNtdERCg4RkRjz0/mb+cOqHdw+dwQ3nevvpbftUXCIiMSQJ5YU89iiYq6fnsu3LxnjdzntUh+HiEgMcM7xy4Xb+NmCLXxy0iB+eNUZmJnfZbVLwSEi4jPnHD9+/UN+9W4pnz5rCP9xzSSSusVmaICCQ0TEV00tAe596QNeWLuLL80czv1XTKBbDIcGKDhERHxzoLaRr/5uLStKqrnjglHcdeGomD091ZqCQ0TEB1vKj3Dz02vYd7ieh6+bzNVnDfW7pA5TcIiIRNnL7+/m3pc+ILVnd56/dUZMjD/lhYJDRCRKjja28OCrRTy3eidnD0/nl9efxaB+vfwuyzMFh4hIFBTtOcRdz7/PlvIabp87grsuGu378OinKy6Cw8yGAt8BCoHJQC8g3zlX5mddIiKn0tQS4LG3i3l04VbS05J5+qZpzBmd7XdZnRIXwQGMBK4F1gLvABf7W46IyKlt2H2If31pPRt2H+bKyYN58MoJpKcl+11Wp8VLcCxxzg0AMLObUXCISAyraWjm4QVb+M3SUjLSknn881O4dOIgv8sKm7gIDudcwO8aREROxTnHa+v38tBfP2TvoXqun57Ldy4ZS7/UHn6XFlZxERwiIrHug12H+MFrRawuO8C4QX355fVTOHt4fF1m21EJGxxmditwK0Bubq7P1YhIotpZXcfDC7bw5/d3k5GazE8+PZFrC4fF9FhTnZWwweGcewJ4AqCwsND5XI6IJJiqmgYeW1TMs8u3Ywa3nlfA1+aNpG9KYp2Wak/CBoeISCRU1zbyxJISnlleRn1TC589exjfuGhUXN7Id7oUHCIiHVBxpIH/fbeUZ5eXUdfUwhWTBnPHBaMYmdPb79KiTsEhInISO6vrePKdEp5fvZOmlgCfnDSYO+aNZNSAPn6X5pu4CQ4zuyb069mhn5eaWQVQ4Zxb7FNZIpKg1u08yBPvlPC3D/aS1M34zJSh3DZnBPlZaX6X5ru4CQ7ghTZ/Pxb6uRiYG91SRCQRNbcEmF9UzlPLSllddoA+Kd25ZXYBN5yT16X6ME4lboLDOZe417aJiK8qaxp4fvVOfrtiO3sP1TMsoxffv3w8100dRu+ecfMxGTVqERHpkpxzrCyt5ncrd/D3DXtpanGcMyKTH1x1BvPG5iT0fRidpeAQkS6lqqaBF/9vF8+t3klJRS19UrrzzzOG8/npw7vkFVKnw3NwmNlAIBdIafuac25JOIoSEQmn5pYAizZX8MLanbz14X6aA47C4el89ZoRfHLSIFKT9R3aiw63lpkNAX4LzG7vZcABSWGqS0SkU5xzFO05zJ/f283L7++hsqaBrN7J3Dgrj2sLh3Xpy2k7y0vMPg6cAfw/4AOgISIViYh0wq4Ddbyybg8vv7eHzeVH6JFkzBubw2emDOX8sTlxO+teLPESHOcBdzjnno1UMSIip6OypoG/fbCXV9btYXXZAQCm5Pbnh586g8snDkqIyZNiiZfgOArsj1QhIiJeVNc2Mr9oH6+v38uy4koCDkbl9OZbF4/mqjOHMCwj1e8SE5aX4HgS+AIwP0K1iIicVMWRBuYX7ePvG/axvKSKloBjeGYqt88dyRWTBzNmoPotouGkwWFmN7X6cxfwBTNbCPwVqG77fufcr8Nbnoh0dTur63hjYznzi/axuqwa5yA/K43bZhfwyUmDGD+oL2a65yKaTnXE8at2nsuj/SE+HKDgEJFOcc6xce9hFmwsZ8HGcor2HAZgzIA+3DFvFJdOHMiYAX0UFj46VXDkR6UKEenSGppbWFlSzVsflvPmh/vZffAoZjAlN517LxvLxeMHkqfBBWPGSYPDObc9WoWISNdScaSBtzfv5+1N+1mypYLaxhZSenTj3JHZ3HnBKOaNyyGrd0+/y5R26HZJEYmKQMCxfvch3t60n0Wb97Nu1yEABvZN4cozh3DhuBxmjcwipYfuI451Xu4cTwbuAf6J4JAjbb8KOOecgkhEjquqaWDJ1goWb65gydZKqmsb6WZwVm46d180mnnjctS5HYe8fND/FPga8DfgJXTnuIi00dQS4L0dB1mypYIlWyv4YPchnIPMtGTmjM5m7phsZo/K1g15cc5LcFwD3O+c+3GkihGR+OKco6yqjne3Bo8olhdXUdPQTFI346xh/fnmhaOZPTqbiUP60U3DlCcML8HRG1geqUJEJD4cqG1kWXEV726r4J2tlew6cBSAoem9uPLMwcwelcXMEVn069XD50olUrwEx6sER8ZdGKFaRCQGHW1sYXVZNUuLK1m6rZKiPYdxDvr07M7MEZncNmcE543MYnhmqvoquggvwfEo8IyZBTjxneMl4SpMRPzR1BJg3c6DLCuuYum2St7bcZDGlgA9koyzctO568LRnDsqi0lD+tFdI812SV6C49hpqgeA+0/wHl1HJxJnWgKOjXsOs7ykkmXFVawqraausQUzGD+oLzfMymPWyCym5qVrwiMBvAXHTQSHFRGROBYIODbtO8LykiqWF1exqrSKw/XNAIzITuMzU4Yyc0QmMwsydfWTtKvDweGceyqCdYhIhAQCjq37a1hRUnX8caCuCYDcjFQuPWMQ54zMZEZBJgP6fmxGaJGP0XGnSIIJBByby4+wsqSKFSXVrCz9R1AMTe/FBeMGMLMgkxkjMhnSv5fP1Uo88hQcZpZD8M7xMUDbrybOOfflcBUmIh3TEnB8uPcwK0urWVlSxaqyag6GgmJI/17MGzuA6QUZzCzI1ORGEhZehhwZA6wg2AGeBlQCGaG/DwCHIlFgq/UPAx4GLgIMeBP4hnNuRyTXKxJrmlsCFO05zKrS4NHEqtLq430UwzJ6cdG4AUwvyGR6foaCQiLC65Ajq4BPAbXApcB64IvAg8DV4S7uGDNLJXj/SAPwJYKd9D8C3jazSc652kitW8Rvjc0BPth9iJWlVawsqWbt9gPUNASDIj8rjcsmDmJ6QQbT8zMZrFNPEgVegmMq8BX+MUZVN+dcM/BrM8sCHgHOD295x90CFABjnHPbAMxsPbAVuA34eYTWKxJ19U0tvLfjIKtKq1lVVsXa7QeobwoAwTm1P3XWYKblB48o1JktfvA65Ei1cy5gZoeArFavrQHuC2tlH3UlsOJYaAA450rNbClwFQoOiWM1Dc2s3X6AVaHTTut2HqKxJYAZjBvYl89NzWVGQQZT8zLI1PwUEgO8BEcZMDD0+2bgs8DfQ39fDhwMW1UfNwF4uZ3ni0J1iMSNg3WNrC4LBsXK0mo27D5EwEH3bsaEIf24cVYe0/IzKMzL0HhPEpO8BMcCgh3TLxD8hv+cmZ0LNANjgUiOmptBsAO+rWogvb0FzOxW4FaA3NzcyFUmcgr7j9SzuvQfQbG5/AjOQXL3bpw5rD9fO38k0/MzOSu3P2k9dYW8xD4ve+k9hCZvcs790cyOAtcBqcAvgCfDX95HtHfX+glHVHPOPQE8AVBYWKg73iVqdh88GgyJkmpWlVZTUhm8diM1OYmzh6dz+aRBTM3LYPKw/prtTuKSlzvHG2g1eZNz7lWCI+ZGwwGCRx1tpdP+kYhIVDjn2FFdx8qSalaEwmL3weAw431TujM1L4PPTRvGtPxMJgzuSw8NCigJwMt9HCXA1c65de28dgbwinOuIJzFtVJEsJ+jrfHAxgitU+RjnHMUV9QevzR2VWk1+w7XA5CRlsy0vAxuPi+f6fmZjBnYhyRNXiQJyMupqjw+Ps/4MSnA8E5Xc2KvAP9pZgXHhm43szxgFvCvEVyvdHHOObYdG+eptJqVJdVU1gQPvLP79GR6fsbxm+1G5fTWfBTSJXjtiTtRX0Ehkb2q6kng68DLZva9UB0/BHYC/xPB9UoX49xHBwRcWVJNVW0jAAP7pnDuyMzjQZGflaagkC7ppMFhZncBd4X+dMCrZtbY5m29CPY/PBf+8kIrdq7WzOYRHHLkWYKd4m8RHHKkJlLrlcTX+ohieZugGNwvhTmjs5lekMGMgkxyMzTDnQic+oijhOAHNASH+lgDVLR5TwPBfoZfhbe0jwqNSfWZSK5DEp9zjtLK2uNzUaxodeppcL8U5ozJZkZBJjPyMxmW0UtBIdKOkwaHc+5lQjfehf4D/cA5VxqFukTCZvfBoyzbVsny4iqWFVcd78we0Lcn547MDE1alKWgEOkgL30ctwHt3sZqZmlAo3OuKSxViXRCdW0jy4urWFpcybJtlZRV1QHBq55mFgSD4pwRmeqjEDlNXoLjSYLBcX07r/0P0EhwelmRqDra2MKqsmqWbqvk3a2VbNx7GIDePbszoyCDL8zMY9bITEbn9KGbLo8V6TQvwXE+8O0TvPYKwWHXRSKuJeAo2nOId7YGg2Lt9gM0tgRITurGlOH9+dbFozlnZBaThvSju264Ewk7L8GRA+w/wWsVwIDOlyPSvj0Hj/Lu1kqWbK1g6bbK41OhjhvUlxtm5TFrZBbT8jLolawhPEQizUtw7AcmAm+389pEoCosFYkQnJNiVWk1i7dUsGRLBVv3B6+6zunTk3ljBzB7dBbnjMgiu4+GGReJNi/B8RrwfTNb5Jxbf+xJM5sIfBf4c7iLk65le1Utb2/az6ItFawoqaK+KUBy925Mz8/g2sJhzB6dzegBujtbxG9eguM+gsOqrzWz1cAuYAgwDSgFvhf+8iSRNTYHWFVazcJN+1m0ef/xUWTzMlP53NRc5owO3lOh008iscXL6LiVZjYV+CbBADkTqCQ4D8fDzrlDEalQEkrFkQbe3rSftzaV8+7WSmobW+jZvRszCjL54szhzB2TQ15Wmt9lishJdCg4zCwZ+Hfg9865+4jsNLGSQI6N/bRgYzlvfljO+zsP4hwM6pfCVWcN4YKxOZwzIktHFSJxpEPB4ZxrNLPbUD+GdEBLwLF2+wHeKNrHgg/L2R66AW/S0H7cdeFoLhiXw/hBfdVXIRKnvPRxvEfw6qklEapF4lhjc4BlxZXML9rHgo3lVNY0kpzUjXNGZnLr7AIuHDeAAX1T/C5TRMLAS3DcDfzBzLYDrzvnNB1rF9fQ3MI7Wyr564a9vLmxnMP1zaQlJ3H+2BwumTCQuWOy6ZPS7ig1IhLHvATHC0A/goMeNpvZfj46P4dzzkVyMieJAY3NAZZuq+TVdXtYsLGcIw3N9E3pzkXjB3LpGQM5d1SW5tEWSXBeguMtTjyRkySwQMCxorSKV9ft4W8b9nGwrom+Kd35xBkDuWzSIGaNyCK5u4b2EOkqvFyOe0ME65AYtHHPYf783i5eXbeXfYfrSU1O4uLxA7hi8mDOG5WtsBDporxOHSsJbv+Rel5+bw8v/t8uNu07Qo8kY87oHL77yXFcOG6ALpsVEe/BYWaTgTHAxy6Rcc49E46iJLqaWgIs3LSfF9bs5O3NFbQEHGcO688Pr5rA5ZMGk56W7HeJIhJDOhwcZtYfeB2Yceyp0M/W/R4KjjhSWlnLc6t38OLaXVTWNJLdpye3nFfANWcPZWROb7/LE5EY5eWI4yEgE5gNvANcDRwiOHnTTOBzYa9Owq6pJcCCjeX8dsV2lhVXkdTNmDc2h89NHcac0dmav0JETslLcFwCPAisCP29yzm3FlhkZo8DdwJfDHN9EiYVRxr43crt/H7lDvYfaWBI/17cfdForp06TDfmiYgnXoJjEFDinGsxs3qgT6vXXgKeC2tlEhYbdh/i1++W8tr6vTS2BJgzOpuHrh7O+WNzSNI0qiJyGrwExz6gf+j37QRPTy0K/T0yfCVJZwUCjsVbKnhiSQnLS6pIS07i+um5fHHmcAqy1XchIp3jJTjeJRgWrwHPAvebWR7QDHyJ4Lzj4qPmlgCvrd/LY4u2saW8hoF9U7j3srF8bloufTX0h4iEiZfgeBAYHPr9pwQ7yq8DUgmGxr+Et7R/MLNvAucDhcBA4EHn3AORWl+8aWwO8Ke1u3h88TZ2Vh9l9IDePHzdZC6fNJge6uwWkTDzcud4MVAc+r2J4KCHd0eorrZuAQ4DfwG+EqV1xrymlgAvrt3Fowu3sfvgUSYP6899l0/ggrE5dFP/hYhESLzcOT7BORcws+4oOAgEHK+u38PP3tjCjuo6Jg/rz4+vPoM5o7M1x4WIRFxcBIdzLuB3DbHina0V/NvfNlG05zDjBvXl1zcUcv6YHAWGiERNXASHQElFDT96/UMWbtrP0PRePHLdmVw5ebBOSYlI1CVscJjZrcCtALm5uT5Xc/pqGpr5xZtb+M3SMlJ6JHHPpWO5YVYePbtrsEER8UfUg8PMLgQWdOCti51zc093Pc65J4AnAAoLC+NuHhHnHPOLynnw1SL2Hqrn2sKhfPuSsWT36el3aSLSxflxxLEMGNeB99VFupBYVX64nu/+eQNvfljO2IF9+OX1Uzh7eLrfZYmIAD4Eh3OuDtgU7fXGA+ccf3l/N/e/XERjS4B7Lh3LTefm614MEYkpCdvHEW+qaxu556X1zC8qZ0puf/7zs5M1PIiIxKS4CA4zKwTygGNfvceb2TWh3/8aOoqJW6tKq7njD+8Fw+PSsdx8XoEGIBSRmBUXwQF8neB4WMd8NvQAyAfKol1QOAQCjscXF/PzBVsYlt6Ll24/hzOG9PO7LBGRk4qL4HDO3QDc4HMZYVXb0Mxdz7/PGxvLuWLyYB66+gz6aCBCEYkDcREciWZndR23PLOGLeVHuO/y8dw4K093fotI3FBwRNn7Ow/y5adW09gS4KkbpzF7dLbfJYmIeKLgiKKl2yq55Zk1ZPXuyR9vnMoIXTUlInFIwREl84v28S+/f4/8rDSe/fI0cjTPt4jEKQVHFLy2fg93Pvc+E4f046kbp9I/NdnvkkRETpuCI8IWb6ngruff5+zcdH5z41TSeqrJRSS+aSyLCFq7vZqvPLuWUTl9+NUNhQoNEUkICo4I2VJ+hBt/s5oBfXvy9E3T6Kt7NEQkQSg4IuBwfRO3PrOGnj2SePbL0zUUuogkFAVHmAUCjrv/uI5dB47y2OenMCwj1e+SRETCSsERZv+9pJgFG8u597JxTM3L8LscEZGwU3CE0bJtlfzn/M1cMXkwN87K87scEZGIUHCESV1jM//vxfXkZabxb5+eqLGnRCRh6frQMPnFm1vZdeAoz986Q5fdikhC0xFHGBTtOcSv3i3lusJhTC/I9LscEZGIUnB0UkvAce9LH5Ce2oN7LhvrdzkiIhGn4Oik36/czrpdh/j+5eM1BpWIdAkKjk5oaG7h0YXbmJ6fwZWTB/tdjohIVCg4OuHFtbvZf6SBOy4YpauoRKTLUHCcpuaWAP+9uJjJQ/txzgh1iItI16HgOE2vf7CXHdV13H7+SB1tiEiXouA4Dc45Hl9UzKic3lw0boDf5YiIRJWC4zQs3LSfTfuO8NW5I+jWTUcbItK1KDhOw1PLyhjSvxdX6EoqEemCFBweHaprYnlxFVdMHkyPJDWfiHQ9Mf/JZ2ajzewXZrbezGrMbK+ZvWJmk/2oZ+HmcpoDjksmqG9DRLqmmA8O4GLgfOBp4ArgdiAbWGlmZ0e7mPkbysnp05PJQ/tHe9UiIjEhHoZxfQ74L+ecO/aEmS0EyoA7gS9Gq5D6phYWb6ngM2cPUae4iHRZMR8czrnKdp47ZGZbgCHRrGXJlgqONrVwyYSB0VytiEhMiYdTVR9jZhnAGcCH0Vzv/KJy+qZ0Z4aGTheRLiwugwN4FDDgkRO9wcxuNbM1ZramoqKi0ytsbgnw1qZyLhg3QFdTiUiXFvVPQDO70MxcBx6LTrD8PcD1wNedc9tOtB7n3BPOuULnXGF2dnan615VWs3BuiZdTSUiXZ4ffRzLgHEdeF9d2yfM7CvAQ8D3nHO/DndhJ/PGxnJ6du/G7NGdDyERkXgW9eBwztUBm7wuZ2ZfAB4Dfuac+3HYCzuFlaXVTMvPIDU55q8nEBGJqLg4WW9mVwO/AX7lnPtWtNcfCDjKKmsZldMn2qsWEYk5Mf/12cxmA38A1gNPmdmMVi83OOfei3QN5UfqOdrUQn52WqRXJSIS82I+OIB5QE/gLGBpm9e2A3mRLqC0ohaAgiwFh4hIzJ+qcs494JyzEzzyolFDSWUwOPIVHCIisR8csaC0spaUHt0Y2DfF71JERHyn4OiA0spa8jLTND6ViAgKjg4praylQB3jIiKAguOUmloC7KiuU/+GiEiIguMUdlbX0RJw5Gf19rsUEZGYoOA4hVJdUSUi8hEKjlM4Fhy6h0NEJEjBcQollbX0T+1Belqy36WIiMQEBccplFbU6jSViEgrCo5TKK1UcIiItKbgOInahmb2Ha5X/4aISCsKjpMoqzp2RZUuxRUROUbBcRK6FFdE5OMUHCdxbDj1vKxUnysREYkdCo6TKK2sZVC/FE0XKyLSioLjJEp0RZWIyMfoq/RJnD08nUH9NAeHiEhrCo6T+P7l4/0uQUQk5uhUlYiIeKLgEBERTxQcIiLiiYJDREQ8UXCIiIgnCg4REfFEwSEiIp4oOERExBNzzvldQ8SZWQWw3eNiWUBlBMpJZGoz79Rm3qnNvDvdNhvunMtu+2SXCI7TYWZrnHOFftcRT9Rm3qnNvFObeRfuNtOpKhER8UTBISIinig4TuwJvwuIQ2oz79Rm3qnNvAtrm6mPQ0REPNERh4iIeKLgEBERTxQcIWb2TTN71cz2mpkzswc8Lv8pM3vPzOrNbLuZfc/MkiJUbkwws25mdo+ZlYW2e52ZfaaDyz4Vaue2j0ciXHbEmdkwM/uTmR0ys8Nm9pKZ5XZw2RQz+2loPzxqZsvNbHaka/ZbJ9usvf3ImdmZES7bV2Y21MweDe0jdaFtzuvgsp3azxQc/3ALkAP8xeuCZnYJ8CKwGrgU+AXwPeChMNYXi34IPAD8kuB2rwBeMLPLOrh8BTCzzePh8JcZPWaWCiwExgJfAr4AjALeNrOOTGD/vwT3xfuAy4G9wPxE/hAMQ5sBPMXH96UtYS82towErgUOAO94XLZz+5lzTo/gBQLdQj+7Aw54wMOy7wGL2zx3H9AIDPR72yLUXjlAA/Bgm+ffAtZ3YPmngF1+b0cE2uVOoAUY2eq5fKAZ+OYplp0c2vdubPVcd2Az8Irf2xaLbRZ6rwN+5Pd2+NBu3Vr9fnOoHfI6sFyn9zMdcYQ45wKns5yZDQPOBH7b5qVngR4Ev4knokuAZD6+3b8FJppZfvRLiglXAiucc9uOPeGcKwWWAld1YNkm4PlWyzYDzwGXmFnP8JcbEzrTZl3W6X5mEYb9TMHReRNCPze0fjK049cB46NeUXRMIHjEsa3N80Whnx3Z7hwzqzSzZjPbYmbfSYB+oQm02RdCijh1m0wASp1zde0sm0zw1EQi6kybHfNVM2sInetfaGbnha+8hNPp/ax7JKrqYjJCPw+089qBVq8nmgzgoAsd57ZS3er1k3kfWEtwZ00BrgZ+QvDc9s3hKzPqMmh/X6gG0jux7LHXE1Fn2gyCR7mvAXuA4cC3gYVmdpFzblG4ikwgnd7PEjI4zOxCYEEH3rrYOTe3s6sL/WzvTkpr57mYdBptZnRim51zj7R56q9mVgN8w8z+3Tm3tSP/Tow63XbpVJvGuc7sS19o9ec7ZvYywSOYHwHnhqG2RNPp/SwhgwNYBozrwPvaHqqdjpOldP9Wr8c6r21WDaSbmbU56khv9bpXfwC+ARQC8RocJzrKTKf9b3mtVQPtXYLamTaNB51ps49xzh0xs9eBL3e2sATV6f0sIYMjdO5uU5RWd+yc/gRg+bEnQ9dTpwIbo1RHp5xGmxUBPYERfLSf49g56dPZ7pMdvcWLIv7R79XaeE7dJkXA1WaW2ub883iCV+i17U9KFJ1psxM50bdqCcN+ps7xTnLO7QDWAZ9v89I/E7xy4W9RLyo6/k5wJ2tvuzeELg7w6nqC/9lXd7I2P70CzDCzgmNPhL5EzAq9dqplewCfbbVsd+A64A3nXEPYq40NnWmzjzGzvsAngZXhKjDBdH4/8/ta5Fh5EDw9cg3BG2oc8MfQ39cAqa3e9xawrc2ylwEB4H+AucBdQD3wU7+3K8Jt9m+h7fxmaLsfD7XDFW3e95E2I9iBuQS4HbgYuAL4dWjZx/3erk62SRrBb2wfELyU9EqCXyxKgN5t2qAZuK/N8s8RPD1zM3AB8KdQG0/xe9tisc2AbwFPEvzSMZfgDYQfEPxSc57f2xaFtjv2GfV46HPrq6G/50RyP/N9w2PlQfCGNHeCR16r9y0CytpZ/tOhnb0B2EHwBsAkv7crwm2WRPAO+e2h7V4PXNPO+z7SZgTPZ/8ltFw9cBT4P+DrtLqpKV4fBM8fvwgcBo6EtjWvzXvyaOdGU6AX8HNgX6htVgJz/d6mWG0zgl86lhKcFrUJqCL4jXqa39sUpXY70WfWokjuZxpWXUREPFEfh4iIeKLgEBERTxQcIiLiiYJDREQ8UXCIiIgnCg4REfFEwSEiIp4oOERExBMFh4iIeKLgEIkiM0szs01mtsrMerR6/mIzC5jZ1/ysT6QjNOSISJSZ2VnACuBh59y/mlkOwXG+VjnnrvS3OpFTU3CI+MDM7gJ+RnB04G8BE4HJzrlKXwsT6QAFh4gPzMyA14F5QDJwkXPuLX+rEukY9XGI+MAFv7E9S3AWxXUKDYknCg4RH5jZQOARgvOQTDazO/2tSKTjFBwiURY6TfU0wVnqLiIYIP9uZpP8rEuko9THIRJlZnY38B/APOfcYjNLJniVVU+g0Dl31NcCRU5BRxwiURS6FPch4CfOucUAzrlG4J8ITvP5c/+qE+kYHXGIiIgnOuIQERFPFBwiIuKJgkNERDxRcIiIiCcKDhER8UTBISIinig4RETEEwWHiIh48v8B+18YQ/Rv2NgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from math import tanh,cosh\n",
    "from numpy import linspace\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "accuracy = 1e-12\n",
    "\n",
    "def arctanh(u):\n",
    "    x = 0.0\n",
    "    delta = 1.0\n",
    "    while abs(delta)>accuracy:\n",
    "        delta = (tanh(x)-u)*cosh(x)**2\n",
    "        x -= delta\n",
    "    return x\n",
    "\n",
    "# collect data for the graph\n",
    "upoints = linspace(-0.99,0.99,100)\n",
    "xpoints = []\n",
    "for u in upoints:\n",
    "    xpoints.append(arctanh(u))\n",
    "\n",
    "# Make the graph\n",
    "plt.rc('font',size=16) # set the font size\n",
    "plt.plot(upoints,xpoints)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('arctanh')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3.1 - Golden ratio search\n",
    "\n",
    "The Buckingham potential is an approximate representation of the potential energy between atoms in a solid or gas as a function of the distance $r$ between them: $$V(r)=V_0\\left[\\left(\\frac{\\sigma}{r}\\right)^6-e^{-r/\\sigma} \\right] .$$ \n",
    "\n",
    "1. Plot the Buckingham potential for $\\sigma$=1.\n",
    "\n",
    "## Talking points\n",
    "\n",
    "1. What do you observe?\n",
    "2. What can you say about the Buckingham potential?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEQCAYAAADlK+DYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtfUlEQVR4nO3deXxU9b3/8ddnJssESCACYV8EBAUVF9Ta2latLdblqnWrV2s3tZv39tpbW1Fva2ttbeu9Wn92EWrFularVWutrWvVWhFwBRXZBWQPS0L2mc/vjzOBMCQhkEnOmcn7+XjkMZMz58z5TIC8+X7P93y/5u6IiIhETSzsAkRERFqjgBIRkUhSQImISCQpoEREJJIUUCIiEkkFYReQLwYMGOCjR48OuwwRkZwyd+7cDe4+sLXXFFBZMnr0aObMmRN2GSIiOcXMlrf1mrr4REQkkhRQIiISSQooERGJJAWUiIhEkgJKREQiSQElIiKRpIASEZFIUkCFbPaySm742wKakqmwSxERiRQFVMhee38Ttzy7iLomBZSISEsKqJCVFMYBqG1IhlyJiEi0KKBCVpwOqLpGBZSISEsKqJAlFFAiIq1SQIWsZHtA6RqUiEhLCqiQJQqDP4K6JrWgRERaUkCFLKFBEiIirVJAhaxE16BERFqlgArZji4+XYMSEWlJARWy4oJ0C0pdfCIiO1FAhaykKB1QGiQhIrITBVTIdB+UiEjrFFAhSxQEfwS1DboGJSLSkgIqZAXxGIVxUxefiEgGBVQEJAri6uITEcmggIqA4kIFlIhIJgVUBJQUxTQXn4hIBgVUBKiLT0RkVwqoCEgUxqlVQImI7EQBFQElugYlIrILBVQEFBfqGpSISCYFVAQk1IISEdmFAioC1MUnIrIrBVQEJNTFJyKyCwVUBGgUn4jIrhRQEaAuPhGRXSmgIqC4ME59U4pUysMuRUQkMhRQEVCSXhOqXsu+i4hsp4CKgERh8Megbj4RkR0UUBGwfVVdrQklIrKdAioCmrv4ahsUUCIizRRQEbCji0/XoEREmimgIqBYXXwiIrvImYAysxFm9kcz22JmW83sITMb2cFjE2b2czNbbWa1ZvYvM/tYK/vFzGyamS0zszoze8PMzsz+p9lZcxdfnbr4RES2y4mAMrNewDPA/sDngc8B+wHPmlnvDrzFbcDFwPeAU4DVwN/M7JCM/a4FrgFuAT4NvAw8YGYndf5TtE2DJEREdlUQdgEddDEwBpjg7osAzOxNYCHwFeD/2jrQzCYD/w58yd1vT2/7BzAf+CHwb+ltFcC3gevd/Yb04c+a2TjgeuDxLvhcwI5rULUNugYlItIsJ1pQBCHycnM4Abj7UuCfwGkdOLYR+EOLY5uA+4CpZlac3jwVKALuyjj+LuAgM9u3U5+gHdu7+HQflIjIdrnSgpoEPNLK9vnA2R04dqm717RybBEwLv18ElAPLGplP4CJwNI9qLnD1MUnItni7iRTTjL92JRykskd3+/05U6qeZ+U4872/Vq+jzskU04q43nKIeXOgUP7MrJ/r6x/llwJqH2ATa1srwTKO3Fs8+vNj5vdPXNCvMz9tjOzS4BLAEaO7NB4jVYlCnQflEiucHcakinqGlPUNyapbUxS35SiruVjY4r6phQNyeB5QzK147EpRUNTisZki8dkisak05j+vjEVPG9KpWhIOk3JFE1JpzEVPCZTTmMytdNjc8g0hTCn549OP5AL+o/K+vvmSkABtPZTtw4cZx08tqP77SjIfTowHWDKlCl7/bciURT0tGouPpHsq29KUlXXRHVdE9X1TWyta9z+fFt9E9saksFjfZKahuD7mvomahqS1DQmqW1oorYxSW1DavvzzmZAUTxGYdwoKohRVBCjMN78ZRTGYxTEYxSln5cUxSiMGQVxS2+PEY8ZBc3bYrH08+Axln4t3uKx+XnzazHbsT0eM+IWvBZsZ/vrMbMWz9m+j8GO12MwuCyRlT+rTLkSUJtopQVD0HpqrXXUUiXQWvOmvMXrzY/lZmYZrajM/bKuKB7DTNegRNpT09DExuoGNtU0sKmmkU3bdjzfUtPA5tpGtqS/ttY2srWuia21jR36j58Z9C4qoFdRnN7FwWOvojj9SgoZUpagV1GcRFGcksLgK1EYI1EYb/EVI1EQpzi9vTgdPMUFO54H38fS/9478n9ryZWAar5GlGki8HYHjj3DzHplXIeaCDSw45rTfKAYGMvO16Emph93d569ZmYkCuLq4pMep64xyfqqetZurWNdVT3rttaxvrqe9VX1bKhuYEN1PRurG9i4rb7dmVZKEwWU9yqib0khfUsKGdq3hLKSAkoThZQlgsfS9GOf4oLgK1FA7+I4fYoLKCmMKzQiKFcC6lHgBjMb4+5LAMxsNPAR4IoOHPsDgsEUd6SPLQDOBf7u7vXp/Z4gCKzz0/s3uwCYlx412GVKiuIaJCF5JZly1m6tY+WmWlZtruGDzXV8sLmW1VvqWL2ljjVbatlU07jLcfGYMaBPEQP6FNO/TzHjBvahf58i9uldzD69C7c/9utVRHmvIsoSBRTEc2VAsuyJXAmoGcClwCNmdjXBtaJrgRXArc07mdkoYDHwQ3f/IYC7v25mfwBuMrNCgpF4XwP2JQgj0vutM7MbgWlmVgW8ShBix7P7oeydliiIaS4+yTl1jUmWb6xh6YZtLN+4jeWVNayorGH5xho+2Fy7ywX78l6FDO1XwrB+CQ4b2Y8hfRNUlCUYVJagorSYitJiynsVEYupNSM5ElDuvs3MjgduBO4kGLjwNPBf7l7dYlcD4ux6f9cXgeuAHwH9gDeAE9391Yz9rgKqgW8Cg4EFwDnu/uesfqBWJArj1OoalETUpm0NLFxXzcJ1VSxcW82SDdtYvK6aD7bU0vKKbXmvQkb2783Bw/tyysFDGF7ei2HlJQwvL2Fo3xJKiuLhfQjJOTkRUADu/j7Q7rx47r6MVkbduXst8K30V3vHJwlC7Ed7XeheShTGqVdASciakikWr9/G26u38M7qKt5ZvZUFa6pYV1W/fZ9eRXHGDuzDlNHljBkwgn0H9mZ0/16MHtCbskRhiNVLvsmZgMp3iUJ18Un3SqacReuqeWPlZt5auYU3V23hndVbaUiPeisqiLFfRR+O2W8A+w8uZb9BpYwfVMqQsoS64KRbKKAiQl180tW21jXy6vJNzFm2iVff38SbK7dQXd8EQJ/iAg4cVsaFHxrFpGFlTBralzEDemvwgYRKARURJYVxttTuOqJJZG9trmlg1tJKXl6ykVlLKnlnzVbcg1FyBwwp5YxDh3HIiH5MHtGPMQN6q1UkkaOAiohEYVw36kqn1DclmbtsEy8s2sCLCzcw74MtuAfdx4ePKuebn9iPI0bvwyEj+tG7WP/0Jfr0tzQiinUNSvbCmi11PP3uWp59dz0vLd5ATUOSgphx2Mhy/usT4/nwuP4cPLwvxQUaPSe5RwEVESVqQUkHLVxbxV/nreGpd9by5sotAAwvL+Ezhw3j2PEVHD22v1pIkhf0tzgi1MUn7Vm4topH3/iAx99azeL12wA4dGQ/Lp86gU9OHMR+FX00VY/kHQVURCQKY9Q2JnF3/aIRIOi+e/j1VTzy+ge8s3orMYOj9u3P5z88mqmTBjOoi2aQFokKBVRElBTGSTk0Jp2iAgVUT1XXmORv89fwx7kr+eeiDaQ8aCldc+pETj54KANLi3f/JiJ5QgEVES1X1S0q0L0nPc2idVXc+8oKHnx1JZtrGhnWr4RLjxvHZw4bzugBvcMuTyQUCqiIKG4OqIakpovpIZIp5+l31jLzpWW8tHgjBTFj6qTBnHfkSD48tr/uS5IeTwEVESXNAaWh5nlvW30T981ewe3/XMrKTbUM6Zvg8qkTOGfKCHXhibSggIqIRGHQrac1ofLX+qp6bv/nUu56eTlb65o4YnQ5V510AJ+cOEhTCom0QgEVEYn0jZRaVTf/rN1ax63/WMI9ryynvinFiZMGc8nHxnDoyPKwSxOJNAVURDSvk6N7ofLH+qp6fvnsIu555X2SKef0Q4bxjePGMmZgn7BLE8kJCqiI2NHFp2tQuW5LbSPTn1/M715cRkMyxZmHDePS4/ZjZP9eYZcmklMUUBFRrC6+nNeYTHHPrPe56an32FTTyKmTh3LZCfupxSSylxRQEdHcxVevQRI56dl313HtX95myfptfHhsf6486QAOHNY37LJEcpoCKiK236ira1A55f2NNfzwsfk89c46xgzszW8vnMInDqjQdFUiWaCAiohEevYIdfHlhvqmJL95bgm/fG4RBTFj2qf354sf2VezgIhkkQIqIraP4tMgicibu3wTVzz4JgvXVXPywUO4+uQDGNK3JOyyRPKOAioimu+DUhdfdNU2JPnpE+9yx7+WMaQswe++MIXj9x8UdlkieUsBFRGxmFFUECy5IdHz2vub+O/732DJhm1cePQovnPi/vTRooAiXUr/wiIkURCjXnPxRUpjMsXNTy/kl88uYnBZgnsuOooPjxsQdlkiPYICKkIShXENkoiQlZtq+M97X+PV9zfzmcOGcc2/TdJM8yLdSAEVISVFcU0WGxFPzFvNd/74JimHX3z2EE47ZFjYJYn0OAqoCEkUxDVIImRNyRQ/feJdZrywlMnD+3LzeYcyqr8WDBQJgwIqQhKFMWp1DSo066vqufSeV5m1tJILjx7F1SdP1H1NIiFSQEVIolAtqLC8tXILF/9+DptqGvi/cybzmcOGh12SSI+ngIqQRGGczTUNYZfR4zz+1mq+df/r9O9dzENf/zCThmoOPZEoUEBFSNDFpxZUd3F3bnlmEf/75HscNrIft35uipZcF4kQBVSElBTGqdM1qG7RlExx1Z/m8Yc5Kzjj0GH85DMHbZ+wV0SiQQEVIboG1T1qGpr4xt2v8uyC9fzn8eO47JPjNfu4SAQpoCIkURhXF18Xq9zWwBdnzuatlZu57owDOf+oUWGXJCJtUEBFSKIwrqmOutC6rXVccNsslm2s4TcXHM6nJg0OuyQRaYcCKkIShTEakimSKSceU5dTNq3aXMv5M15mXVU9M79whObTE8kBCqgIabmqbm/NlJ01yzdu499nzGJrXSN3fvkoDh9VHnZJItIBuk0+Qkq07HvWrais4bzpL7OtoYl7L/6Qwkkkh3Tov+lmNgqYCOyT3lQJvO3uy7uqsJ6oOaBqGpL0D7mWfLBqcy3nzXiZ6vom7rn4Qxw4TDfgiuSSdgPKzE4ErgcOAjIviriZzQOmufvjXVRfj9InEfxxVNc3hVxJ7luzpY7zpr/MltpG7r7oKIWTSA5qs4vPzD4D/AXYAHwZOArYDxiffv5lYD3wqJmd0fWl5r/SdEBV1SmgOmPTtgY+d9ssKrc18PsvHcnBw/uFXZKI7IX2WlDfA6a7+9faeH02MNPMfg18H/hTtovraUrTi+FV1TWGXEnu2lbfxBdnzmZ5ZQ13fPFIDh2pa04iuaq9QRITgPs68B73pfeVTlILqnMamlJ89a65vLlyM//vvEM5eqyu5InksvYCag1wWAfe47D0vl3GzGJmNs3MlplZnZm9YWZn7sHxp5vZa+ljl5vZ1WYWz9jnGjPzVr4ezvoHasOOgFILak+5O9998E1eWLiB6888mKm6CVck57XXxTcd+LGZlQJ3u/vili+a2RjgfGAa8KOuKxGAa4FvA1cBc4HPAg+Y2Sm7G6BhZlOBB4HbgG8BhwI/BkqB77ZyyDFAy3HelZ2uvoPK0l18W9WC2mM3PrWQP722im9/ajznTBkRdjkikgXtBdT1QB/gSuD7ZlYPbAYcKAeKgUbgRuAnXVWgmVUQhNP17n5DevOzZjYuXePuRhBeD7zo7pe0OLYPcLWZ3ejuma2/We4eSkIkCuMUxWNsVQtqj/xx7kpufnoh50wZzjeOGxd2OSKSJW0GlLs7cJWZ3QScyI77oIygVTEfeMLd13dxjVOBIuCujO13Ab8zs33dfWlrB5rZCOAQ4JKMl+4EfgB8Grg9q9V2UmmiQNeg9sC/Fm9k2kNv8pFx/bnujIM0K7lIHmkzoMysxN1r0wF0ZzfWlGkSUA8sytg+P/04EWg1oNLHAsxrudHdl5pZTfrYTCvSrbaVBANArnH32r0pfG8ooDpuRWUNX797LqP69+ZX5x9OYVwTo4jkk/b+Ra81s9+Z2ce7rZrW7QNsTrfoWqps8Xp7xwJsauW1TRnHLgKuAD5P0GK8H7gMeLStNzezS8xsjpnNWb8+Ow3J0kShBkl0wLb6Ji7+/RySKee3F06hb0lh2CWJSJa1dw3qEeBs4PNmtgL4PXCnuy/szAnN7ATgyQ7s+g93P5agSzEznGDXmS1aPV36cbfHu3tmF+KTZrYSuMnMTnD3pzLfwN2nEwwmYcqUKa2dY4+pBbV77s7lf3yD99ZWMfOLRzJ6QO+wSxKRLtDeNajPmVlv4CzgQoLBEleZ2SxgJnC/u2/ei3O+BBzQgf1q0o+VQLmZWUYrqrzF621pr5XVbzfHAtwL3AQcAewSUF2hNFHAhg3buuNUOetXzy3m8bfWcNVJB/Cx8QPDLkdEuki7c/G5+zbgDuAOMxsGfA64APgNQcvisfTrf3X3Dq205+41wLt7UON8ghGDY9n5OlTz9aO3d3MsBNei/tW80cxGA712c2xLWWkddURZolAtqHa8tHgD//v3Bfzb5KFc9NF9wy5HRLpQh68qu/sqd7/e3Q8EjgRmAB8nuEazqovqA3gCaCC456qlC4B5bY3gA3D394E32ji2Efjrbs7dfNysDlfbSaWJQrbW6hpUa9ZtreM/732dMQP78JPPaMSeSL7bq1Xx3H1OuvuvnOCXeEVWq9r5XOvM7EZgmplVAa8C5wLHA6e13NfMngZGuXvLm2GuBB4zs1sJuuwOBa4GftHyHigze43gOtsCghbTJ4H/IBhK/2xXfb5MpYkCtjUktapuhqZkikvvfY1t9U3ce/FRWtBRpAfYo3/lZjaBHd18I4CtBDM03JH90nZyFVANfBMYTBAi57j7nzP2i5Pxmdz9cTM7i2BC2y8Aawlmkrgu49gFwKXAkPT7LAZ+CPwsmx9kd5qnO6qua6JvL41Ma3bjU+/xytJKbjx3MvsNKg27HBHpBrsNKDPrD5xHEExTgBTBgIErgIfdva5LKwTcPUkwnVK7UyqlR/21tv0h4KHdHPvZva0vm3ZMd9SogEp7afEGfvXcYj57xAjOOHR42OWISDdp70bdswhC6USgkGDAwRXAXe6+unvK63k0o/nONm1r4Ft/eIN9B/Tme6e2dl+1iOSr9lpQ9wMbgVuBO9x9bveU1LNpTagd3J1pD73Fxm31/PbzH6FXka47ifQk7f2Lvxr4ubvrN2U3KitRC6rZfbNX8MT8NVx50v5asl2kB2pvmPm1wPtmdquZnWhmuiDSDUpbXIPqyd7fWMO1j73NR8b156JjxoRdjoiEoL2AGkYwim0U8DCwwcz+YGafNbOy7iiuJ9I1KEilgqmM4mb8/KzJxDTcXqRHajOg3H2Nu//a3U8EBgJfJbg/6DfAOjP7m5l91cyGdlOtPYJW1YWZLy1j1tJK/ueUiQztVxJ2OSISkg7NJOHuVe5+b3oo9kDgDIIlLv6HYHmKWWZ2RRfW2WMUF8QpKoj12BbUkvXV/Oxv73LchIGcPUVDykV6sjYDKr3Uxscyt7t7o7v/1d2/6u7DCJZIf45gQlnJgrJEQY9c9j3o2nuToniM6888WFMZifRw7bWgziVYHn2pmf0gvcT6Ltz9X+7+XXfXTSpZ0lPXhLp71nLmLt/E906dxKCyRNjliEjI2guoQcBFwDKCIecLzOxFM7vYzDTmtwuV9cA1oVZvqeWnTyzgmHEDOPOwYWGXIyIR0N4giWp3v93djwNGE1xv2ofgxt3VZnavmX3azLTOdpaVJgp71DBzd+d7j8ynKZXix2dolnIRCXR0kMQKd/9xuhvvQ8DvgE8AjwGrzOyGLqyxx+lpq+o+MW8NT769lstOGM/I/r3CLkdEImKPWz/u/oq7X0pwn9SNBEttXJbtwnqyIKB6Rguqqq6R7z86n0lDy/jyMVqAUER22OPJzdKDJS4kWHJjFFAFPJDlunq00h60qu4vnlrI+up6pl84hYK4eotFZIcOBZSZ7UMwqu9CgtV0nWDJjauAP3XHkhs9SWmigJqGJE3JVF7/0l6wporbX1rGZ48YySEj+oVdjohETHvLbRQCpxCE0qeBIuBtYBrBkhsfdEuFPVDzfHzV9U3061UUcjVdw935n0fmUZoo4DtTJ4RdjohEUHstqDVAP6ASmEGw5Mac7iiqpytrMR9fvgbUI69/wCtLK/nxGQdR3js/P6OIdE57AfUCwVLuj2nJje6V7zOaV9c3cd3j7zB5eF/OPWJE2OWISES1GVDufno31iEtNLegttbm50CJXz+3iPVV9Uz/3OHENVO5iLQhf6/A57B8XlV3RWUNM15YyumHDOXQkeVhlyMiEaaAiqB8XhPq+ifeJWbwnRP3D7sUEYk4BVQE5euaUHOWVfKXN1fzlY+N1TpPIrJbCqgI2tHFlz8tqFTK+eFjbzO4LMFXPq4l3EVk9xRQEVRUECNRGKOqPn8C6rG3VvPmyi1cPnUCvYr2eAITEemBFFARlU9rQjU0pbjhbws4YEgZZxyqpTREpGMUUBFVmijIm2Hm98xazvuVNXz3xAnENKxcRDpIARVR+bImVFVdIzc/s4ijx/Tn4+MHhl2OiOQQBVRE5cuqujOeX0Lltgau+PT+WohQRPaIAiqi8mFNqPVV9cx4YSknHzSEyZqtXET2kAIqokqLc39NqF8/t5iGZIr//tT4sEsRkRykgIqospLc7uJbvaWWu2Yt58zDhjFmYJ+wyxGRHKSAiqjSRCG1jUkak6mwS9krtzyzCHfnP47fL+xSRCRHKaAiqnm6o+ocbEWtqKzhD7NXcO4RIxixT6+wyxGRHKWAiqhcXhPq5qcXEosZlx6n1pOI7D0FVEQ1rwm1pTa3Amrphm08+OpKPvehUQzumwi7HBHJYQqoiBpYWgzAhur6kCvZM7c8s4iighhf/fjYsEsRkRyngIqoirKg9bF2a+4E1PKN23j49VWcf9So7QErIrK3FFARNbBP8At+XQ4F1K+eXUw8ZnzlY1pOQ0Q6TwEVUUUFMcp7FbKuqi7sUjpkRWUND766kn8/cuT21p+ISGcooCKsojTBuqrcaEH9+h+LiZlpMUIRyRoFVIRVlBXnREB9sLmWB+as4JwjhjOkr5ZyF5HsUEBF2MDSYtZvjX4X3/Tnl+CORu6JSFYpoCKsojTB+up63D3sUtq0sbqe+2a/z2mHDGN4uWaNEJHsyYmAMrOYmU0zs2VmVmdmb5jZmR089lQzu8fM3jOzlJk9186+x5jZS2ZWa2ZrzOz/zCy0PquK0mIak86mmujerDvzpWXUN6X42rG69iQi2ZUTAQVcC1wD3AJ8GngZeMDMTurAsacDh6SPWdnWTmZ2MPAksA44Bbga+CIwc6+r7qRB6dFwUR3JV1XXyMyXljF14mDGVZSGXY6I5JmCsAvYHTOrAL4NXO/uN6Q3P2tm44Drgcd38xYXu3sq/V4vtrPfDwgC7Gx3b0zv3wDcYWY/dfdXO/M59kZFWXAv1Nqt9ew/uLvPvnt3z3qfqromvn6crj2JSPblQgtqKlAE3JWx/S7gIDPbt72Dm8OpPWZWCJwI3N8cTmn3Aw3AaXtUcZZUlDbfrBu9FlRdY5LbXlzKMeMGcPDwfmGXIyJ5KBcCahJQDyzK2D4//TgxC+cYCySAeS03unsdsDhL59hjFaXNXXzRG2r+4KsrWV9Vz9ePVetJRLpGLgTUPsBm33UoW2WL17NxDoBNrbxW2dY5zOwSM5tjZnPWr1+fhTJ2VlIUp7S4gPURC6hkypnx/BImD+/L0WP7h12OiOSpbg8oMzvBzLwDX881HwK0Ns7asllW+nGPzuPu0919irtPGThwYBbL2WFgWXHkBkk8+fYalm2s4ZKPjcUsm38MIiI7hDFI4iXggA7sV5N+rATKzcwyWlHlLV7vrPZaY+Xs6E7sdhWlxZGaMNbd+fU/ljByn16ceGAER26ISN7o9oBy9xrg3T04ZD5QTHCdqOV1qObrQm9noazFBNe5JrXcaGYJYAzwQBbOsVcqShO8vmJzWKffxStLK3ljxWauPf1A4jG1nkSk6+TCNagnCEbSnZ+x/QJgnrsv7ewJ3L0hfZ5zzKxlaJ9FEI6PdvYce6uiNOjii8psEtOfX0L/3kWcffjwsEsRkTwX+fug3H2dmd0ITDOzKuBV4FzgeDKGf5vZ08Aodx/XYtso4Ij0t/2BlJmdlf5+trsvTz+/BvgXcL+Z/RIYDfwc+KO7z+2Kz9YRg8oS1DWmqKpvoixRGFYZACxcW8XT767jshPGkyiMh1qLiOS/yAdU2lVANfBNYDCwADjH3f+csV+cXT/TccDtGduau+y2zxTh7q+b2VTgp8BfgC3A74Ers/MR9k7zzbrrttaFHlAzXlhCojDG544eFWodItIz5ERAuXsS+FH6q739jm1l20w6OF2Ruz8PHL3HBXahgaU7VtYNczqhdVV1PPzaB5x7xAj26V0UWh0i0nPkwjWoHi0qN+v+/qXlNKZSfPmYdifuEBHJGgVUxG3v4gvxXqiahibumrWcTx4wiNEDeodWh4j0LAqoiCstLiBRGAv1XqgH565kc00jF39MS2qISPdRQEWcmVFRmgitiy+Zcm57cSmTR/Rjyqjy3R8gIpIlCqgc0HwvVBiefHstyzbWcPFH99W0RiLSrRRQOaCirDi0FtRtLy5hWL8STpykaY1EpHspoHJARWmC9SFcg3p9xWZmL9vEl47Zl4K4/qqISPfSb50cUFFWTFV9EzUNTd163t++sITS4gLOmaJpjUSk+ymgcsD2e6G6sRW1clMNf523hvOOGklpyDNYiEjPpIDKAc1Lv6/txqXfZ/5zGQCf//DobjuniEhLCqgcMLp/cHPs0g3buuV8W+sauW/2Ck4+aAjD+pV0yzlFRDIpoHLA8PISSgrjLFhb1S3nu3/2Cqrrm7joo5rWSETCo4DKAbGYMX5QH97rhoBqTKb43YtLOXL0Phw8vF+Xn09EpC0KqBwxflApC9ZUd/l5Hn9rNR9sqeMSTWskIiFTQOWICYNL2VBdz8bqrhvJ5+7MeGEJYwb25vj9K7rsPCIiHaGAyhHjBwVrQb23tutaUf9aspF5q7Zy8UfHEItpWiMRCZcCKkdMGNwcUF13HWrG80sY0KeIMw4d1mXnEBHpKAVUjqgoLaZvSWGXjeR7b20Vzy5Yz4VHjyZRGO+Sc4iI7AkFVI4wMyYMKuW9NV0TUNOfX0KiMMYFHxrVJe8vIrKnFFA5ZPzgPixYW4W7Z/V9V22u5eHXVnHulBHs07soq+8tIrK3FFA5ZMKgUqrqmliT5SmPZjy/BIBLPj42q+8rItIZCqgc0jySb0EWu/k2Vtdz3+z3Oe2QYZrWSEQiRQGVQ3YMNc9eQM18aRn1TSm+dqxuzBWRaFFA5ZDy3kVUlBZnbUaJqrpGZr60jKkTBzOuojQr7ykiki0KqBwzYXBp1lpQd896n6q6Jr5+nK49iUj0KKByzPhBpSxcV0Uy1bmRfNX1TUx/fgkf3W+AJoUVkUhSQOWYCYNKqWtMsaKyplPvc/uLS6nc1sB/f2pClioTEckuBVSOGZ+e8mj+B1v3+j021zQw/fklfGriIA4Z0S9LlYmIZJcCKsdMGlpGea9C/jZ/zV6/x63PL6G6oUmtJxGJNAVUjimMxzjxwME89c5aahuSe3z8uqo6bv/nUk6bPHT7BLQiIlGkgMpBpx48lJqGJM8uWLfHx/7ymUU0Jp3/OmF8F1QmIpI9CqgcdNSY/gzoU8xjb36wR8fNW7WFO19eznlHjmD0gN5dVJ2ISHYooHJQPGacdNBgnnl3HdX1TR06pimZYtpDb9G/TzGXT92/iysUEek8BVSOOnXyUOoaUzz9ztoO7T/zpWW8tWoL3z91In1LCru4OhGRzlNA5ajDR5YzuCzBn99Yvdt9V26q4f+efI/j96/g5IOGdEN1IiKdp4DKUbGYcfLBQ3j+vfVsqW1sc7/GdNeeO/zwtEmYWTdWKSKy9xRQOezUyUNpSKa47YUlrb6eSjmXP/AGLyzcwP+cMpHh5b26uUIRkb2ngMphk4f35fRDhnLzM4v4bUZIuTtXPzKPh1//gMunTuDfjxoZUpUiInunIOwCZO+ZGTecPZn6phQ/+ss7FBfE+Mxhw3ljxWYeem0Vf5y7kq8dO5ZvHDcu7FJFRPaYuXduVmwJTJkyxefMmRPKuRuaUnz97rk89c464jHbPtP5Rcfsy1UnH6DrTiISWWY2192ntPaaWlB5oKggxi/PP4xbnlkEwOGjyjl0ZLmGk4tITlNA5YnigrgmfxWRvKJBEiIiEkk5EVBmFjOzaWa2zMzqzOwNMzuzg8eeamb3mNl7ZpYys+fa2O8aM/NWvh7O5mcREZGOyZUuvmuBbwNXAXOBzwIPmNkp7v74bo49HTgEeBlIdOBcxwAt17Go3NNiRUSk8yIfUGZWQRBO17v7DenNz5rZOOB6YHcBdbG7p9Lv9WIHTjnL3Ts2A6uIiHSZXOjimwoUAXdlbL8LOMjM9m3v4OZwEhGR3JILATUJqAcWZWyfn36cmOXzrTCzpJktN7OfmllJlt9fREQ6IPJdfMA+wGbf9Y7iyhavZ8Mi4ArgNcCBTwGXAYcBn2ztADO7BLgEYORITSUkIpJN3R5QZnYC8GQHdv2Hux8LGEFg7PJW2azL3TO7EJ80s5XATWZ2grs/1cox04HpEMwkkc16RER6ujBaUC8BB3Rgv5r0YyVQbmaW0Yoqb/F6V7kXuAk4AtgloFqaO3fuBjNb3oW1dKcBwIawi+gB9HPuevoZd4/O/JxHtfVCtweUu9cA7+7BIfOBYmAsO1+Har729HaWSmvPbltH7j6wG+roFmY2p625sSR79HPuevoZd4+u+jnnwiCJJ4AG4PyM7RcA89x9aReeu/mcs7rwHCIi0orID5Jw93VmdiMwzcyqgFeBc4HjgdNa7mtmTwOj3H1ci22jCLroAPoDKTM7K/39bHdfnt7vNeD3wAKCFtMngf8AnnD3Z7vq84mISOsiH1BpVwHVwDeBwQQhco67/zljvzi7fqbjgNsztj2QfvwiMDP9fAFwKTAk/T6LgR8CP+t8+TlnetgF9BD6OXc9/Yy7R5f8nLUelIiIRFIuXIMSEZEeSAElIiKRpIASzOwsM3swPb1TrZktMLOfmFlp2LXlMzN7Ir2ky4/CriXfmNlJZva8mVWb2VYzm2Nmx4ddVz4xs4+Y2d/NbF36Z/yqmX0pm+dQQAkEs8UngSuBE4FfA18jmE1Df0e6gJmdB0wOu458ZGZfAR4hWJrnDOBsgoFRvcKsK5+Y2cEEkxcUAhcDZwKzgdvM7GtZO48GSYiZDXT39RnbLgTuAD7h7s+EU1l+MrN+BDerXwbcA1zn7leHWlSeMLPRwDvANHe/Kdxq8peZ/ZjgP7b7uHt1i+0vA+7uR2fjPPrfsZAZTmmz04/DurOWHuJnwHx3vzfsQvLQl4AU8JuwC8lzRUAjUJuxfTNZzBUFlLTl4+nHd0KtIs+Y2THAhcDXw64lTx1D0Dr9rJktNrMmM1tkZt8Iu7A8MzP9eLOZDTWzfmZ2MfAJ4MZsnSRXbtSVbmRmwwhuUn7K3eeEXU++MLNC4FbgBndfEHY9eWpo+uvnBNdUFxNcg7rFzArc/RdhFpcv3H2emR0L/Ikd/9lqBL7q7vdl6zwKKNmJmfUhuMDcRDDThmTPd4ES4LqwC8ljMaAU+IK7P5Te9kz62tQ0M7u5lbXlZA+Z2X7AgwSTeX+VoKvvNOA3Zlbn7ndn4zwKKNnOzBLAo8AY4OPuvjLkkvKGmY0kmLLrIqDYzIpbvFycHjhR5e7JMOrLIxuB/dh1zbm/E4xQHQJ80N1F5aEfE7SYTnH3xvS2p82sP/ALM7vX3VOdPYmuQQmwvfvpQeBI4CR3fyvkkvLNGCAB3AVsavEFwWioTcBB4ZSWV+a3sb15gdNO/9IUIPi7+kaLcGr2CsGk3BXZOIkCSkjf63Q3wQXO09z95ZBLykevE0xcnPkFQWgdx87rncne+VP6cWrG9qnASndf08315Ks1wCFmVpSx/SigjiwtJKsuPgH4JcGF5OuAbWb2oRavrVRXX+e5+2bgucztZgaw3N13eU32yuPAs8CtZjYAWAKcBXwKXVPNplsIbn7+s5n9iuAa1L8B5wE3untDNk6iG3UFM1tG28su/8Ddr+m+anoWM3N0o25WmVkZ8BOCYConGHZ+vbvfE2phecbMPk0w8GcSQff1YoJlN27N1rVUBZSIiESSrkGJiEgkKaBERCSSFFAiIhJJCigREYkkBZSIiESSAkpERCJJASUSMWZ2nplVZczX19a+o9PLxn+hG0oT6VaaSUIkek4HnnD3+g7suxo4muAmSZG8oht1RSLAzIrdvT49t9l64OvZWrJAJFepi0+km5nZNeluuQPN7G9mVg3cn375EwRrRv0lve9gM7vDzD4ws3ozW21mj5lZRfr1Vrv4zOybZrbMzOrM7BUz+3D6+5kt9vlC+tgPm9n96W7FtWY2Lf36iWb2mpltM7PZZnZ4xjk+ZWaPp2uqMbN5ZvbfZhbvqp+d9Czq4hMJzyPAbcBP2bEMxOnAP9KTywLcSTBP4uXACmAQQYj1autNzewi4Kb0ez8AjAXuAfq1ccgdwO8J5lE7G/hxen2qkwgmEK4GfgY8bGZjW0wEOgZ4Gvh/BDNYTwGuAQYCV+z+44u0TwElEp6bWy5BbsHU5qey84q7RwNXZnT3PdDWG6aXTvk+8Fd3v6jF9jUE63215k53vza933PAGcC3gPHuvrTF+z6SrucfAO7+m4zaXwCKgG+b2ZXZWLBOejYFlEh4/pTx/YcIVnx9pMW22cDl6QB4Bpi3myXLh6e/vpex/RGgqY1j/tr8xN2bzGwR0Lc5nNLeTT+OaN5gZkMIWkwnAkPZ+fdJBcGaQSJ7TdegRMKzOuP704E5GetvnQs8CnwHeBNYZWbfS7doWjMk/biu5cb08gcb2jhmU8b3DW1sg2BZheYW1aPAKcCPgOOBI9jR+ku0cS6RDlNAiYQnsyV0GvDwTju4r3P3b7j7MGB/YCbwA+Arbbxnc+jttOR2euDCgE7W29JYgmtO33X3Ge7+grvPAbKyDpAIKKBEIsHM9gcmkBFQLbn7Ane/kqB1c2Abu61Mf52dsf10stul3zxIo7F5g5kVAudn8RzSw+kalEg0nAEscvf5zRvMrC/wFHA3wTWgRoJWVjnw99bexN1TZvYDYIaZ/ZZgQMUYglF1W9gxWrCz3gGWA9eZWTJd22VZem8RQAElEhWns2vrqQ54FbiYYKh5ClgAnO/uj9AGd/+tmfUhCIwLgHkELZs/E4RUp7l7g5mdDtxCMES9Evgd8D4wIxvnENFMEiIhS4+GWwV81N3/2UXnOAJ4BbjQ3e/sinOIZJsCSiTPmNm+wDcI7kvaChwAXEkwEu9Ad68JsTyRDlMXn0j+qSUYRHEhwfWqTQTXsq5QOEkuUQtKREQiScPMRUQkkhRQIiISSQooERGJJAWUiIhEkgJKREQi6f8D3tI45XfXEb8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from math import exp\n",
    "from numpy import linspace\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "xpoints = linspace(1.2,8.0,num=100)\n",
    "ypoints = []\n",
    "for x in xpoints:\n",
    "    ypoints.append((1/x)**6-exp(-x))\n",
    "\n",
    "# Make the graph\n",
    "plt.rc('font',size=16) # set the font size\n",
    "plt.plot(xpoints,ypoints)\n",
    "plt.xlabel('r/sigma')\n",
    "plt.ylabel('V/V0')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3.2 and 3.3. - Golden ratio search\n",
    "\n",
    "The Buckingham potential is an approximate representation of the potential energy between atoms in a solid or gas as a function of the distance $r$ between them: $$V(r)=V_0\\left[\\left(\\frac{\\sigma}{r}\\right)^6-e^{-r/\\sigma} \\right] .$$ \n",
    "\n",
    "2. Complete the golden ratio example program below to find the minimum of the Buckingham potential.\n",
    "3. Check your computational against the analytic solution.\n",
    "\n",
    "## Talking points\n",
    "\n",
    "3. How does the number of iterations depend on the specified accuracy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The minimum falls at 1.630516067174875 nm\n"
     ]
    }
   ],
   "source": [
    "from math import exp,sqrt\n",
    "\n",
    "# Constants\n",
    "sigma = 1.0             # Value of sigma in nm\n",
    "accuracy = 1e-6         # Required accuracy in nm\n",
    "z = (1+sqrt(5))/2       # Golden ratio\n",
    "\n",
    "# Function to calculate the Buckingham potential\n",
    "def f(r):\n",
    "    return (sigma/r)**6 - exp(-r/sigma)\n",
    "\n",
    "# Initial positions of the four points\n",
    "x1 = sigma/10\n",
    "x4 = sigma*10\n",
    "x2 = x4 - (x4-x1)/z\n",
    "x3 = x1 + (x4-x1)/z\n",
    "\n",
    "# Initial values of the function at the four points\n",
    "f1 = f(x1)\n",
    "f2 = f(x2)\n",
    "f3 = f(x3)\n",
    "f4 = f(x4)\n",
    "\n",
    "# Main loop of the search process\n",
    "while x4-x1>accuracy:\n",
    "    if f2<f3:\n",
    "        x4,f4 = x3,f3\n",
    "        x3,f3 = x2,f2\n",
    "        x2 = x4 - (x4-x1)/z\n",
    "        f2 = f(x2)\n",
    "    else:\n",
    "        x1,f1 = x2,f2\n",
    "        x2,f2 = x3,f3\n",
    "        x3 = x1 + (x4-x1)/z\n",
    "        f3 = f(x3)\n",
    "\n",
    "# Print the result\n",
    "print(\"The minimum falls at\",0.5*(x1+x4),\"nm\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4.1 - Gauss-Newton method\n",
    "\n",
    "Return to the Buckingham potential from the previous exercise: $$V(r)=V_0\\left[\\left(\\frac{\\sigma}{r}\\right)^6-e^{-r/\\sigma} \\right] .$$ Find the minimum of the Buckingham potential for $\\sigma$=1:\n",
    "1. For the Gauss-Newton method. Start from $r=\\sigma$.\n",
    "\n",
    "## Talking points\n",
    "\n",
    "1. What do you observe?\n",
    "2. What happens when you start from $r=4\\sigma$ and why?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The minimum falls at 1.6305159287146183 nm\n"
     ]
    }
   ],
   "source": [
    "from math import exp,sqrt\n",
    "\n",
    "# Constants (V0 is set to 1)\n",
    "sigma = 1.0             # Value of sigma in nm\n",
    "accuracy = 1e-6         # Required accuracy in nm\n",
    "\n",
    "# Function to calculate the Buckingham potential\n",
    "def f(r):\n",
    "    return (sigma/r)**6 - exp(-r/sigma)\n",
    "\n",
    "# Function to calculate the first derivative of Buckingham potential\n",
    "def fp(r):\n",
    "    return -6*(sigma**6/(r**7)) + exp(-r/sigma)/sigma\n",
    "\n",
    "# Function to calculate the second derivative of Buckingham potential\n",
    "def fpp(r):\n",
    "    return 42*(sigma**6/(r**8)) - exp(-r/sigma)/(sigma**2)\n",
    "\n",
    "# Initial position\n",
    "x = sigma\n",
    "err = 1.0\n",
    "\n",
    "# Main loop of the search process\n",
    "while err>accuracy:\n",
    "    xp=x-fp(x)/fpp(x)\n",
    "    err=abs(x-xp)\n",
    "    x=xp\n",
    "\n",
    "# Print the result\n",
    "print(\"The minimum falls at\",xp,\"nm\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot of first derivative of the Buckingham potential."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEQCAYAAABr8amkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAk1klEQVR4nO3deZhcZZn38e9d1Xtn34EkdEJQIYA6BATHmRFcEMdXcGHUkVHcGGeY11FnceNlUVDHmXlxH0FBBESRQUUEGUEQV2QJgiBbyAKBJJ2k00m6q7pru+ePc6q7utJLJV2nqvvU73NdfVXVWaruFFz96+c8z3kec3dERERqIVHvAkREpHEodEREpGYUOiIiUjMKHRERqRmFjoiI1ExTvQuY6hYsWOBdXV31LkNEZFq5//77d7j7wvLtCp0JdHV1cd9999W7DBGRacXMNo22XZfXRESkZhQ6IiJSMwodERGpGYWOiIjUjEJHRERqRqEjIiI1o9AREZGa0X06MdebyrB+Rz/dewbYunuAvQM58u4UHJoSxozWJma0NTGnvZlFs9pYPKuVBTNaaU7q7xERqT6FTswUCs5dT2znzse7uWdDD49t3bvf75EwOGh2O8vmtbN8XgcrFsxgxYJOVi7s5ND5HbQ2JSOoXEQagUInJlKZHDesfZZv/moD63f009GS5NhD5/KXRx/EkQfPYsnsNpbMamN2ezPJhGFm5PIF+gZz7B3I0ZvK0r13gG17BtmyO80zPSme7klxx2Pd7OjbPPQ5yYSxfF4Hhy3sZNWimRy+aAaHL57BYQtn0Nmq/51EZHz6LREDv31qJx+87gG27RnkmKWz+eLbXsxrVi+hpWn8S2RNyQRzOlqY09HCsnkAs0c9bs9Alo07+lm/vZ/12/tYt72Pdd193PXEdrL54ZVnD5nTzvMWz2DVohkcvmgmq8Lns9qaq/ivFZHpTKEzjeXyBb74syf50p3rWLGgk+vOfjHHr5iHmVX1c2a1NXPM0jkcs3TOiO3ZfIGne1I8uW0vT27r48nuPp7YtpdfP7WTTK4wdNyima2sWhS0hlYu7Bx6PHh2O4lEdWsVkalNoTNNpTN53nXlPdy9voc3H7uUC1+/uuaXt5qTCQ5bGITJa44a3p4vOM/0pHiyu4+nwlbRuu4+fvj7Z9k7kBs6rrUpQdf8TroWdNC1oJOu+UGf0aHzO1kyq42kAkkkdhQ601C+4Hzwugf43YYePvfmY/irNcvqXdIIyYQFIbKgk1exeGi7u7O9b5CnuvvZuLOfDTvCy3Xdfdz52HYy+eHWUUsywdK57Syd18Gyue0sndvB0rntHDK3naVz2lkwo1WtJJFpSKEzDX36lkf5n0e2cd7rjpxygTMeM2PRzDYWzWzjxMPmj9iXLzjP9abZtDPFpp5+nu5J8UxPimd60jy0uZfeVHbE8c1JY8nsNg6e3c5Bs9tYEj4untXK4lltLJ7VxoIZrRP2a4lIbU2b0DGzZcAlwKsAA24HPujuT1dwbhvwKeBMYA7we+Aj7v6LqOqNylW/3cjlv9rAWS/t4t0vW1HvcqommTCWzetg2bwOXsaCffb3DeZ4dleazbtSPNub5rneAZ7rTbNld5r7Nu1i254tIwY1FM3taGbRzDYWzGxh4YzgHqT5M1pZMKOFBTNamdvZwvzOFuZ1ttDRkqx6f5iIjDQtQsfMOoA7gEHgnYADFwF3mtkx7t4/wVtcDvwl8C/AeuAc4H/M7ER3/31khVfZ41v3cuFNf+SVRyzi/73uyHqXU1MzWpt4/pKZPH/JzFH3FwrOzv4M2/YM0L13gK27B9m+d5DuvQN07x1kR98g9z+9i+17BxnIFkZ9j5ZkgjkdzcztaGFOR3Pw097C7I5mZrc3M6s9fGxrYlb4OLOtmZltTbQ3K7BEKjEtQgd4H7ASeL67rwMws4eAJ4G/Bf7/WCea2QuBvwbe7e7fDLfdBTwCfBJ4fbSlV8/FtzxKZ0uS/zjjhepkL5NIGAtntrJwZitjDf0uSmVy7OzLsKNvkF2pDDv7MvT0Z9iVyrKrP8OuVIbedJaNO1LsSvWyO51lMDd6UBUlE0ZnS5KZbc10tibpbG1iRmsTnS1NdLY20dmapKOlic6WJO0twfOOliRtzUk6wm3tzcHr9pYkbU0J2sLX+m8tcTJdQuf1wN3FwAFw9w1m9mvgNMYJnfDcLHBdybk5M/su8FEza3X3wYjqrpqfP97NL57Yzrl/eQRzOlrqXc601tHSRMe8JpbN66j4nIFsnj3pLHsGcuwZyA493zuQZW/42D+YZ+9Ajr7B4edbdw+QyuTpz+RIDeZHDJaoVHPSaG1K0hoGUWtTgpamBK1NCVqbkrSUvG5pStCcDB5bkgmakza0rTl83ZxM0JRM0JI0mhIJmorbEsV9RjJ8nkwYTYnh44qvk+G2RAKaEgmSZiSTRtKMRILgdXgTskip6RI6q4EbR9n+CHBGBeducPfUKOe2AKvC51NWLl/g4psfpWt+B+84save5TSkYqtj0azJvU8uXyCVzZMazJPK5Ehn86Qz+RGPg9kC6WyegWyewVyBgWyegWyBwdzw60yuwGAu2JbK5NiVKpDJFcjmg8dM3snk8mTzTjZfIFfYt7+rFhIWtAITYQgFoWThNkiYhT9BazUxFFZBcCUseJ4Iwyx4XXouGMPHlD6aGRbWUHxuxX0w8ljGPgeA0mPCc4fPCQ8Y2j5yf3FveQAXjyk9b+h5ybHlNZSeX7q/9O1HO27Ee5XtsNE3c/qLDqFrQSfVNF1CZx6wa5TtPcDcSZxb3D+CmZ0NnA2wfPnyyquMyHfvfYYnu/v42pnHajTWNNeUTDArmaj5LA2FgpMtFIIQyhWGnufy4WOhQK4koLL5AoUCZMPt+ULwkysUwkcPtrtTCF8PPboPHVMoDB+TLwQTzRb3Fzz8KTB0TMGdfHiMh8e5M3Seu5c8Z+g93Bl6T6e4HSh57c7QOTD8PDgs3B8+LzgEexg6j+JxxW0Mn1M8Dka+T7Bl5PnFbcUTxzo2fLfhGkr/g/q++8t2jfi8fbeP8T9KmRctm9OwoQNl33mokra77e+57n4ZcBnAmjVr6vMnYqh/MMcltz3B8SvmccrqxROfIDKKRMJoTSRpbQJa612NTFVeaRpNwnQJnV2M0iIhaOWM1oop1QOM1lyZW7J/yrrlD1vY2Z/hv179fF0fF5FI1eJ3zHS5VvMIQd9MuSOBP1Zw7opw2HX5uRlg3b6nTB03rN1M1/wOjuua6CqiiMjUN11C50fACWa2srjBzLqAPw33TXRuMyUDDsysCXgL8NOpPHLtmZ4Ud6/v4Y1/slStHBGJhekSOl8HNgI3mtlpZvZ6gtFszwCXFg8ys0PNLGdm5xW3hTd/Xgd83szea2avAL4LrADOr90/Yf/94IFnAXjDiw+pcyUiItUxLUInnHHgZOAJ4Grg28AG4GR37ys51IAk+/673gV8k2AWg5uBZcBr3H1txKUfMHfn+2s3c8LKeft1P4mIyFQ2XQYSEM6x9qYJjtnIKKPS3D0NfDj8mRbu37SLjTtTnHPSqnqXIiJSNdOipdOIbli7mfbmJKcefVC9SxERqRqFzhQ0kM3z4we3cOpRS5hR44XZRESipNCZgn67fid7B3OcpgEEIhIzCp0p6Hfre2hOGsd3jXY/rIjI9KXQmYJ+t2EnxyydQ3tLst6liIhUlUJnikllcvxh825eskKtHBGJH4XOFLN2Uy+5gnO8QkdEYkihM8X8bsNOkgljjfpzRCSGFDpTzO/W93DUwbM0VFpEYkmhM4UMZPP8/pleXrJyfr1LERGJhEJnCnng6V4y+YKGSotIbCl0ppB7NvRgBsdpEIGIxJRCZwr53YadHLFkFrPbm+tdiohIJBQ6U0QmV2Dt07s0VFpEYk2hM0X84dleBrIFTlip0BGR+FLoTBGPbtkLwNFL59S3EBGRCCl0poh13X10tCQ5eHZbvUsREYmMQmeKWNfdx6pFMzDbZ+FTEZHYUOhMEeu6+1i1cEa9yxARiZRCZwrYM5Bl654BVi1W6IhIvCl0poCnuvsA1NIRkdhT6EwBT4ahc/jimXWuREQkWgqdKeCp7j5akgmWzW2vdykiIpHar/nzzewY4M+B+cCl7r7VzFYB29x9bxQFNoJ13X2sXNhJU1J/A4hIvFUUOmbWClwDvBEwwIGbgK3A54AngI9GVGPsPdndx9FLZ9e7DBGRyFX6p/XFwCuBvwEWEwRP0U+AU6pcV8MYyOZ5ZldKgwhEpCFUenntbcC57n6tmSXL9m0AuqpaVQN5ansf7nC4hkuLSAOotKUzH3h0nPdorU45jWddcbj0IoWOiMRfpaGzAThxjH3HA49Xp5zGs667j4TBigWd9S5FRCRylYbOVcBHzeztQEu4zc3sJOBDwBVRFNcI1nX3cej8Tlqbyq9aiojET6Wh8zngZuBqoCfc9ivgduBWd/9SBLU1hCfDiT5FRBpBRQMJ3D0PvNXMvkIwUm0RsJMgcO6KsL5Yy+YLbNzRz6uOXFzvUkREamK/bg51918Cv4yoloazaWeKXME5XC0dEWkQFV1eM7O1ZvZBM9Of5FWkkWsi0mgq7dPZRtCv84yZ3WJmbzUzLXE5Sc/2pgFYNrejzpWIiNRGRaHj7qcCS4F/JejPuRbYZmaXhyPY5ABs2zNAS1OCOR3N9S5FRKQmKp5h0t273f3z7r4GWA18BXgFcLuZbYqqwDjbsnuAJbPatES1iDSMA5rW2N0fBT4JfAJ4jqAVJPtp2+4BlszWVUoRaRz7HTpmdrKZfZOgn+cqYDPwf6tdWCPYuido6YiINIpKlzY4CjgT+GvgEGAT8AXgand/Mrry4svdg9BRS0dEGkil9+k8BOwGricIGt2rM0m9qSyZXIHFaumISAOpNHTeAvzI3QejLKaRbN0zAKDLayLSUCqdBuf6qAtpNEOho8trItJAxgwdMzsP+Ia7Pxc+H4+7+6eqW1q8bd2t0BGRxjNeS+cC4FaCIdEXTPA+Dih09sPW3QOYwaKZWv9ORBrHmKHj7onRnkt1bNszwPzOVpqT+mpFpHFUOuHncjMbda4WM2sys+XVLSv+guHSauWISGPZn+WqXzzGvheG+2U/bN2tG0NFpPFUGjrjTQ7WDBSqUEtD2aYbQ0WkAY03em0OMK9k0yFmtrLssHbgncDW6pcWXwPZPLtSWbV0RKThjDd67R+B8wlGpjnw32McZ+FxUqFt4T06mo1ARBrNeKHzQ2AjQahcAVwEPFV2zCDwR3d/KIri4kr36IhIoxpvyPSDwIMAZubAj919Z60KizNNgSMijarSaXC+FXUhjWSbpsARkQZV6YSfxeUN3gM8Hyj/benu/opqFhZnW3cP0tmSZGablqkWkcZS6c2hLwHuA04FTgHmAiuBlwOrGH9I9aSZWcLMPmZmG81swMweNLM3VXjulWbmo/x8Psqax7NtzwCL1coRkQZU6X06nwa+D6wmCJj3uHsX8EogSTDIIEqfIpj/7csEwXc3cL2ZvbbC87cDJ5b9XFL9MiuzZXda/Tki0pAqvbx2DMH9OB6+TgK4+x1mdhHwGeAl1S8PzGwR8M/AZ939P8LNd5rZKuCzwC0VvE3G3e+Oor4DsW3PIC9ZMW/iA0VEYqbSlk4z0O/uBaAHOKhk3+PAUdUurMQpQAtwTdn2a4CjzWxFhJ9ddYWC6/KaiDSsSkPnKeCQ8PlDwLvDfpYE8C6inZFgNcH9QOvKtj8SPh5ZwXssMrMdZpYzsyfM7CNmlqxqlRXa2Z8hV3AOUuiISAOq9PLaTQSDBq4l6N+5GdgD5IEZwAeiKC40D+h1dy/b3lOyfzy/B+4nCKk24A0ElwMPB9472glmdjZwNsDy5dWdQFuzEYhII6v0Pp0LSp7fbmYnAG8COoBb3f2nlX6gmb0SuK2CQ+9y95cTDFwoDxyocMScu3++bNMtZtYHfNDM/s3dnxzlnMuAywDWrFkz2mcfsKHZCBQ6ItKAKr5Pp5S7PwA8cICf+RvgiAqOS4WPPcBcM7Oy1s7ckv376zvAB4E1wD6hE6WtujFURBpYRaFjZt8HrgJudvfsZD7Q3VPAY/txyiNAK3AYI/t1in05fzyAMoqtpKq2Yiqxqz8DwLzOllp/tIhI3VU6kOAFBPfpbDGzr4SX12rlViADvL1s+5nAw+5+IAvI/TVB4Nw7ydr2W286S2dLUstUi0hDqrRP50gzOxb4G+AtwPvNbD1B6+fb7r4+qgLdvdvMLgE+ZmZ7gbVhDScDp5Uea2Y/Aw5191Xh60OBq4HvErSSWgkGEpwFXOru5bNmR253OsucDrVyRKQxVdyn4+73A/eb2T8R3DtzJvAR4AIz+427/1lENQJ8AugjWONnCcG9QX/l7jeVHZdk5L9pL0Gfz0eAxQStm0cJRtt9NcJ6x9SbyjKrXXOuiUhj2u+BBO6eJ5gF4BYzezVwOfDSahc2ymdexATT7YSj3Upf9wCnR1bYAdidzjBHoSMiDWq/OxbM7DAzO9/MngB+QtAp/59VryymdqezzFboiEiDqnT02lyCfpS/AU4gGM78A+Ac4PZRbtyUMfSmsszpUOiISGOq9PLaVoL+kjsIOuFvCIc+y35SS0dEGlmloXMucI27b4mymLgbyOYZzBWYrZaOiDSoSodM/3vUhTSC3engvlq1dESkUY0ZOmb2DoIZCHaGz8fl7ldVtbIY6k0FoTOnXffpiEhjGq+lcyXBoIGd4fPxOMGNojKO3lQwBY4GEohIoxovdFYAW0qeyyTp8pqINLoxQ8fdNwGYWTPwIuChA5znTEK9Ch0RaXAT3hwazir9PaAr8mpibk8xdHR5TUQaVKUzEqwHFkVZSCPoTWVJJoyZrQe0jJGIyLRXaeh8DviEmS2Mspi4601nmNXWhFlFi56KiMROpX9ynwzMAzaY2d0EAwxKp75xd39ntYuLm93pnJY1EJGGVmnovAzIAtsJVvA8rGy/5l6rQG8qo0EEItLQKp2RQEOmq2CPFnATkQanNZNrqDetGaZFpLFVHDpm1mlmHzCz/zazO83s8HD7W83sBdGVGB+aYVpEGl2l6+ksA34OLAUeA44CZoa7TwJeCbw3gvpio1BwdqezWjVURBpapS2d/wQGgcOBYwlWCy26C/jzKtcVO3sHcrjDLIWOiDSwSkevvQo4292fNrNk2b5ngUOqW1b8FOdd00ACEWlklbZ0WoC9Y+ybTTCcWsbRmw5mmFafjog0skpD5yHgTWPsOxW4vzrlxNdwS0ehIyKNq9LLa/8O/Hc4fcu14bYjzew04D3A6yOoLVaGF3BT6IhI46r05tDvm9nfA58F3h1uvorgkts/uPutEdUXG1rWQESk8pYO7v41M7saOJFgxumdwG/cfay+HilRXNZAo9dEpJHt1xz77t4P3B5RLbHWm8rQ1pygrbl88J+ISOMYM3TMbL/uvXH3X0y+nPgKbgzVcGkRaWzjtXR+zvDs0cbEM0nrT/hx9KY0BY6IyHihc1LJ8znAl4CHge8C24DFwNuA1cA5EdUXG73prJapFpGGN2bouPtdxedmdiXwU3cvn1/tKjO7HHgjcFMkFcbEnnSWZfM66l2GiEhdVXpz6GnAdWPsuy7cL+PoTWmyTxGRSkMnAawaY9/hqD9nQru1lo6ISMWhczPwGTM7ozjhp5klzeyvgIuAH0dVYBwM5vKks3kNJBCRhlfpfTofAJYRXErLmdkuYG54/q/C/TKG4rxrszXDtIg0uEqnwdkB/JmZvQo4ATgI2AL81t11s+gEdqc0BY6ICOz/jAS3AbdFVEtsFedd00ACEWl0lfbpyCSopSMiElDo1ECv1tIREQEUOjWxW8saiIgACp2aKIbOzDaFjog0NoVODfQP5uhoSZJMWL1LERGpq4pCx8zyZnb8GPuONbN8dcuKl1QmT0eLJm0QEam0pTPen+hJJl72oKGlMznaFToiIuPfp2NmCYYDJxG+LtUOnArsiKC22Ehl8nQ079ctUSIisTTeyqHnA+eFLx349Tjv89VqFhU36Wyejla1dEREJlo5FIKWznnA5cDmsmMGgT+iCT/HVRxIICLS6CZaxO0uADNz4Bvu/mytCouTVCbPvM7WepchIlJ3lU74eWH5NjM7EjiCYNLP56pdWJyksxq9JiIClQ+Z/rKZfa3k9RuBB4HrgT+a2XER1RcLGjItIhKodMj0qcBvSl5fSNCP80LgHuD8KtcVK+lMno4WjV4TEak0dJYAGwHMbCmwGviMu/8B+CKgls4Y3J1URgMJRESg8tBJAzPC538B7AHuC1/3ATOrXFdsDOYKFBzdHCoiQuWLuK0FzjGzp4FzgNvcvRDuW0GwiqiMIpUJZghSS0dEpPLQ+QRwK8HggV7g/SX7Tifo15FRpDI5QKEjIgIVXl5z93uB5cDxwAp3f6hk92VEPJDAzD5sZjeZ2RYzczO7YD/PP93MHjCzATPbZGbnmllNUiA91NLRQAIRkQlDx8xazGwt8Kfufr+77ynd7+43u/sTkVUYeB+wCPjh/p5oZqcANwD3EozC+wJwLvDpKtY3Jl1eExEZNuGf3+6eMbMVQK4G9YxltbsXzKyJkZf2KvFZ4Ffufnb4+k4zmwGca2aXuPvWqlZapj+8vKaBBCIilY9euw14dZSFjKdk0MJ+MbNlwIuAa8p2XQ00E7R8IqXLayIiwyr9Tfgl4JqwpfFDgtFqI9bQcff11S2tKlaHjw+XbnT3DWaWAo6MugBdXhMRGVZp6NwVPn4Y+NAYx0zF36rzwsddo+zbVbJ/BDM7GzgbYPny5ZMqIK3QEREZUmnovKtaH2hmryS4XDeRu9z95ZP9uPBxtJVNx1wN1d0vIxiVx5o1aya1KurwkGldXhMRqXSW6W9V8TN/QzA79URSVfisnvBxtBbNnJL9kUll1dIRESmq+Z/f7p4CHqvRxz0SPq4GflvcaGZdQAfBAnSRSg3mMYPWpkrHbIiIxNd4y1VfAXwq7HS/YoL3cXd/T3VLmzx3f9rMHgTeDnyjZNeZQBb4SdQ1pDJ5OpqTmI15NU9EpGGM19I5ieBGSoCTGb1fpGhS/R4TMbM1QBfDQ7yPNLM3h89vCVtPmNnPgEPdfVXJ6R8HfmxmlwLfAV5McHPoF6K+Rwcgnc3R0ar+HBERGH+56hUlz7tqUs3Y/gF4Z8nrM8IfCCYc3Rg+T1L2b3L3W8KAOh84C9hGMBvBxdGVO0wLuImIDBuzo8HMeszsT8LnV4SzEtSFu5/l7jbGz8aS414+WkC6+/fd/YXu3uruy939k+6er0XtqUye9maFjogIjD8jQSfQGj4/C1gYeTUxpAXcRESGjdfZsAl4n5kVg+fFZtY21sHu/ouqVhYTqUyeTt2jIyICjB86nwUuJehLceCrYxxn4X79OT+KdCbPwhmtEx8oItIAxhtIcIWZ/QR4HnAn8AHg0VoVFhcaSCAiMmzc6z7uvgXYYmbfAm529w21KSs+Upk87bq8JiICVD4NTtXmXms0aQ0kEBEZorlZIuTupLK6vCYiUqTQidBAtoC7ZpgWESlS6ERoeFkDtXREREChE6niqqHtCh0REUChE6m01tIRERlBoROh/kFdXhMRKaXQiVA6U2zpaCCBiAgodCKVyujymohIKYVOhFLq0xERGUGhE6F0OGRa0+CIiAQUOhEaurymRdxERACFTqSGQqdVoSMiAgqdSKUyOZIJoyWpr1lEBBQ6kUpl8nQ0JzGzepciIjIlKHQilM7kNQWOiEgJhU6EtGqoiMhICp0IpTI5zUYgIlJCoRMhtXREREZS6EQopT4dEZERFDoRSqulIyIygkInQqms+nREREopdCKklo6IyEgKnQj1Dyp0RERKKXQiUig46WxeM0yLiJRQ6ERkIKe1dEREyil0IqJVQ0VE9qXQiUh6KHR0eU1EpEihE5H+cNVQtXRERIYpdCJSvLymGQlERIYpdCKS1lLVIiL7UOhEJKU+HRGRfSh0IpIq9um0qqUjIlKk0ImIhkyLiOxLoRORodBp1uU1EZEihU5E0uHlNY1eExEZptCJSCqTpylhtDTpKxYRKdJvxIhoqWoRkX0pdCJyxEEzOfWog+pdhojIlKJe7oi85bjlvOW45fUuQ0RkSlFLR0REakahIyIiNaPQERGRmlHoiIhIzSh0RESkZhQ6IiJSMwodERGpGYWOiIjUjLl7vWuY0sxsO7Cp3nVUyQJgR72LiDl9x7Wh77k2JvM9H+ruC8s3KnQaiJnd5+5r6l1HnOk7rg19z7URxfesy2siIlIzCh0REakZhU5juazeBTQAfce1oe+5Nqr+PatPR0REakYtHRERqRmFjoiI1IxCJ8bM7M1mdoOZbTKztJk9bmafMbOZ9a4tzszsVjNzM7uo3rXEjZm91sx+YWZ9ZrbHzO4zs5PrXVecmNmfmtlPzaw7/I7Xmtm7q/X+Cp14+2cgD3wceA3wX8DfAbeZmf7bR8DM3ga8sN51xJGZ/S1wI3A/8AbgDOB6oKOedcWJmR0D3A40A+8D3gTcC1xuZn9Xlc/QQIL4MrOF7r69bNs7gG8Br3D3O+pTWTyZ2RzgMeBDwLXAxe5+bl2Ligkz6wIeBT7m7p+vbzXxZWafJvhjdZ6795Vsvxtwdz9xsp+hv3ZjrDxwQveGj4fUspYG8TngEXf/Tr0LiaF3AwXga/UuJOZagCyQLtveS5XyQqHTeP4ifHy0rlXEjJm9DHgH8Pf1riWmXkbQinyrmT1lZjkzW2dm59S7sJi5Mnz8opkdbGZzzOx9wCuAS6rxAU3VeBOZHszsEOCTwO3ufl+964kLM2sGLgX+w90fr3c9MXVw+PPvBH2UTxH06XzZzJrc/Qv1LC4u3P1hM3s58AOG/4DKAu939+9W4zMUOg3CzGYQdMLmgHfVuZy4+QjQDlxc70JiLAHMBM5y9++H2+4I+3o+ZmZfdHVQT5qZHQ7cADwCvJ/gMttpwNfMbMDdvz3Zz1DoNAAzawN+BKwE/sLdN9e5pNgws+XAJ4D3Aq1m1lqyuzUcXLDX3fP1qC9GdgKHA7eVbf8pwcjMg4Dnal1UDH2aoGXzOnfPhtt+ZmbzgS+Y2XfcvTCZD1CfTsyFl35uAI4HXuvuf6hzSXGzEmgDrgF2lfxAMApoF3B0fUqLlUfG2G7h46R+EcqQo4EHSwKn6B5gPrBosh+g0Imx8F6cbxN0Ap7m7nfXuaQ4+j1w0ig/EATRScC6ulQWLz8IH08p234KsNndt9a4nrjaCrzIzFrKtr8EGAB6JvsBurwWb18h6Gy9GOg3sxNK9m3WZbbJc/de4Ofl280MYJO777NPDsgtwJ3ApWa2AFgPvBl4NeqjrKYvE9xwe5OZfZWgT+f1wNuAS9w9M9kP0M2hMWZmG4FDx9h9obtfULtqGouZObo5tKrMbBbwGYKwmUswhPqz7n5tXQuLGTM7lWBwzGqCS8dPESxxcGk1+iYVOiIiUjPq0xERkZpR6IiISM0odEREpGYUOiIiUjMKHRERqRmFjoiI1IxCR6RGzOxtZra3bH62sY7tCpe8PqsGpYnUjGYkEKmd04Fb3X2wgmO3ACcS3JgnEhu6OVQkQmbW6u6D4VxW24G/r8b08CLTlS6viVSJmV0QXhI7ysz+x8z6gO+Fu19BsObOzeGxS8zsW2b2nJkNmtkWM/uxmS0K9496ec3M/tHMNprZgJndY2YvDV9fWXLMWeG5LzWz74WX9LaZ2cfC/a8xswfMrN/M7jWzY8s+49VmdktYU8rMHjazfzKzZFTfnTQOXV4Tqb4bgcuBf2N4yv3TgbvCCUIBriaYF+9fgGeAxQTB1DHWm5rZe4HPh+99PXAYcC0wZ4xTvgVcRTBv1hnAp8P1fV5LMAlsH/A54IdmdljJZI4rgZ8BXyKYWXgNcAGwEPjoxP98kbEpdESq74ulyydbMOX0/2HkyqInAh8vu9R2/VhvGC5TcT7wE3d/b8n2rQTrJY3manf/VHjcz4E3AB8GnufuG0re98awnrsA3P1rZbX/EmgB/tnMPj7ZRbyksSl0RKrvB2WvTyBY2fLGkm33Av8S/lK/A3h4guWWl4Y/55VtLy5BPpqfFJ+4e87M1gGzi4ETeix8XFbcYGYHEbRsXgMczMjfE4sI1lwROSDq0xGpvi1lr08H7itbv+gtBEuI/yvwEPCsmZ0XtjxGc1D42F26MZxqfscY5+wqe50ZYxsEU9gXWz4/Al4HXAScDBzHcCutbYzPEqmIQkek+spbLKcBPxxxgHu3u5/j7ocALwCuBC4E/naM9ywG2YjlgsPO/QWTrLfUYQR9OB9x96+7+y/d/T5g0uuoiIBCRyRSZvYC4PmUhU4pd3/c3T9O0Ao5aozDNoc/Z5RtP53qXiYvDmTIFjeYWTPw9ip+hjQw9emIROsNwDp3f6S4wcxmA7cD3yboU8kStIbmAj8d7U3cvWBmFwJfN7NvEAw6WEkwmmw3w6PkJutRYBNwsZnlw9o+VKX3FlHoiETsdPZt5QwAa4H3EQybLgCPA2939xsZg7t/w8xmEITAmcDDBC2QmwiCZ9LcPWNmpwNfJhhu3QNcATwNfL0anyGNTTMSiEQkHAX2LPBn7v7riD7jOOAe4B3ufnUUnyFSTQodkWnCzFYA5xDcN7MHOAL4OMEItKPcPVXH8kQqostrItNHmmCgwTsI+n92EfQNfVSBI9OFWjoiIlIzGjItIiI1o9AREZGaUeiIiEjNKHRERKRmFDoiIlIz/wtJw/VmnIsRyAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from math import exp\n",
    "from numpy import linspace\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# Function to calculate the first derivative of Buckingham potential\n",
    "def fp(r):\n",
    "    return -6*(sigma**6/(r**7)) + exp(-r/sigma)/sigma\n",
    "\n",
    "\n",
    "xpoints = linspace(1.2,8.0,num=100)\n",
    "ypoints = []\n",
    "for x in xpoints:\n",
    "    ypoints.append(fp(x))\n",
    "\n",
    "# Make the graph\n",
    "plt.rc('font',size=16) # set the font size\n",
    "plt.plot(xpoints,ypoints)\n",
    "plt.xlabel('r/sigma')\n",
    "plt.ylabel('first derivative')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4.2 - Gradient descent\n",
    "\n",
    "Return to the Buckingham potential from the previous exercise: $$V(r)=V_0\\left[\\left(\\frac{\\sigma}{r}\\right)^6-e^{-r/\\sigma} \\right] .$$ Find the minimum of the Buckingham potential for $\\sigma$=1:\n",
    "\n",
    "2. For gradient descent.\n",
    "\n",
    "## Talking points\n",
    "\n",
    "2. What do you observe?\n",
    "3. What is a good value for $\\gamma$ in gradient descent?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x:  1.0  err:  1.1264241117657114\n",
      "x:  2.1264241117657114  err:  0.01774835466002811\n",
      "x:  2.1086757571056833  err:  0.017806619876254537\n",
      "x:  2.090869137229429  err:  0.017846940009032863\n",
      "x:  2.073022197220396  err:  0.017867204371291745\n",
      "x:  2.055154992849104  err:  0.017865177022855505\n",
      "x:  2.0372898158262487  err:  0.017838512376462834\n",
      "x:  2.019451303449786  err:  0.017784778197627915\n",
      "x:  2.001666525252158  err:  0.017701487333570798\n",
      "x:  1.9839650379185871  err:  0.0175861394658805\n",
      "x:  1.9663788984527066  err:  0.017436274015973474\n",
      "x:  1.9489426244367332  err:  0.017249534998721616\n",
      "x:  1.9316930894380115  err:  0.017023748074423528\n",
      "x:  1.914669341363588  err:  0.016757009256169964\n",
      "x:  1.897912332107418  err:  0.016447783669723215\n",
      "x:  1.8814645484376948  err:  0.016095011448709284\n",
      "x:  1.8653695369889856  err:  0.015698216338227766\n",
      "x:  1.8496713206507578  err:  0.015257610994461324\n",
      "x:  1.8344137096562965  err:  0.01477419149419612\n",
      "x:  1.8196395181621003  err:  0.014249812457408417\n",
      "x:  1.805389705704692  err:  0.013687233729292592\n",
      "x:  1.7917024719753993  err:  0.013090130050958049\n",
      "x:  1.7786123419244413  err:  0.012463056785137594\n",
      "x:  1.7661492851393037  err:  0.011811367623238533\n",
      "x:  1.7543379175160652  err:  0.011141084139807855\n",
      "x:  1.7431968333762573  err:  0.010458721694200035\n",
      "x:  1.7327381116820573  err:  0.009771080903696605\n",
      "x:  1.7229670307783607  err:  0.009085017998583522\n",
      "x:  1.7138820127797771  err:  0.008407210108757113\n",
      "x:  1.70547480267102  err:  0.007743932399506104\n",
      "x:  1.697730870271514  err:  0.007100862768632865\n",
      "x:  1.690630007502881  err:  0.0064829267178345695\n",
      "x:  1.6841470807850465  err:  0.005894190543876743\n",
      "x:  1.6782528902411697  err:  0.005337805905794912\n",
      "x:  1.6729150843353748  err:  0.004816003906306321\n",
      "x:  1.6680990804290685  err:  0.004330132749509552\n",
      "x:  1.663768947679559  err:  0.003880730233141172\n",
      "x:  1.6598882174464178  err:  0.0034676209493804677\n",
      "x:  1.6564205964970373  err:  0.003090027996776712\n",
      "x:  1.6533305685002606  err:  0.002746689964558957\n",
      "x:  1.6505838785357017  err:  0.0024359755752207946\n",
      "x:  1.6481479029604809  err:  0.002155990301008215\n",
      "x:  1.6459919126594726  err:  0.0019046712055486026\n",
      "x:  1.644087241453924  err:  0.0016798679937304328\n",
      "x:  1.6424073734601936  err:  0.0014794096634189735\n",
      "x:  1.6409279637967746  err:  0.001301157201559544\n",
      "x:  1.639626806595215  err:  0.0011430434702242032\n",
      "x:  1.6384837631249909  err:  0.0010031018339919928\n",
      "x:  1.637480661290999  err:  0.00087948525202175\n",
      "x:  1.6366011760389771  err:  0.000770477560826377\n",
      "x:  1.6358306984781508  err:  0.0006744985646891521\n",
      "x:  1.6351561999134616  err:  0.0005901043767673197\n",
      "x:  1.6345660955366943  err:  0.0005159842505380041\n",
      "x:  1.6340501112861563  err:  0.00045095493301183787\n",
      "x:  1.6335991563531445  err:  0.00039395337388148377\n",
      "x:  1.633205202979263  err:  0.00034402844750047734\n",
      "x:  1.6328611745317625  err:  0.0003003321913770396\n",
      "x:  1.6325608423403855  err:  0.0002621109365237295\n",
      "x:  1.6322987314038617  err:  0.00022869660026936245\n",
      "x:  1.6320700348035924  err:  0.00019949832863752626\n",
      "x:  1.6318705364749548  err:  0.00017399461025613405\n",
      "x:  1.6316965418646987  err:  0.0001517259340337862\n",
      "x:  1.631544815930665  err:  0.00013228802572173848\n",
      "x:  1.6314125279049432  err:  0.00011532567147942707\n",
      "x:  1.6312972022334638  err:  0.00010052711752472909\n",
      "x:  1.631196675115939  err:  8.761902205778327e-05\n",
      "x:  1.6311090560938812  err:  7.636192743887271e-05\n",
      "x:  1.6310326941664424  err:  6.654621587331988e-05\n",
      "x:  1.630966147950569  err:  5.798850967408953e-05\n",
      "x:  1.630908159440895  err:  5.052847681374928e-05\n",
      "x:  1.6308576309640812  err:  4.402600336783635e-05\n",
      "x:  1.6308136049607134  err:  3.835869617696552e-05\n",
      "x:  1.6307752462645364  err:  3.341968129588757e-05\n",
      "x:  1.6307418265832405  err:  2.911566632635143e-05\n",
      "x:  1.6307127109169142  err:  2.5365237364960436e-05\n",
      "x:  1.6306873456795492  err:  2.2097363958195615e-05\n",
      "x:  1.630665248315591  err:  1.925008801251238e-05\n",
      "x:  1.6306459982275785  err:  1.6769375053238278e-05\n",
      "x:  1.6306292288525253  err:  1.4608108513058937e-05\n",
      "x:  1.6306146207440122  err:  1.2725209826980333e-05\n",
      "x:  1.6306018955341852  err:  1.1084869047994772e-05\n",
      "x:  1.6305908106651372  err:  9.655872439395807e-06\n",
      "x:  1.6305811547926978  err:  8.41101507820241e-06\n",
      "x:  1.6305727437776196  err:  7.3265879170225645e-06\n",
      "x:  1.6305654171897026  err:  6.3819300151202185e-06\n",
      "x:  1.6305590352596875  err:  5.559037768332331e-06\n",
      "x:  1.6305534762219192  err:  4.8422239711243265e-06\n",
      "x:  1.630548633997948  err:  4.217820418039864e-06\n",
      "x:  1.63054441617753  err:  3.6739185331757795e-06\n",
      "x:  1.6305407422589968  err:  3.200143204207251e-06\n",
      "x:  1.6305375421157926  err:  2.7874555967866144e-06\n",
      "x:  1.6305347546601958  err:  2.427981257602241e-06\n",
      "x:  1.6305323266789382  err:  2.1148602786791315e-06\n",
      "x:  1.6305302118186595  err:  1.8421167049531562e-06\n",
      "x:  1.6305283697019546  err:  1.604544722200174e-06\n",
      "x:  1.6305267651572324  err:  1.3976094759282631e-06\n",
      "x:  1.6305253675477565  err:  1.2173606478427246e-06\n",
      "x:  1.6305241501871086  err:  1.060357148752189e-06\n",
      "x:  1.6305230898299599  err:  9.236015059421732e-07\n",
      "The minimum falls at 1.630522166228454 nm\n"
     ]
    }
   ],
   "source": [
    "from math import exp,sqrt\n",
    "\n",
    "# Constants (V0 is set to 1)\n",
    "sigma = 1.0             # Value of sigma in nm\n",
    "accuracy = 1e-6         # Required accuracy in nm\n",
    "gamma = 0.2            # approximate of inverse 2nd derivative\n",
    "\n",
    "# Function to calculate the Buckingham potential\n",
    "def f(r):\n",
    "    return (sigma/r)**6 - exp(-r/sigma)\n",
    "\n",
    "# Function to calculate the first derivative of Buckingham potential\n",
    "def fp(r):\n",
    "    return -6*(sigma**6/(r**7)) + exp(-r/sigma)/sigma\n",
    "\n",
    "# Function to calculate the second derivative of Buckingham potential\n",
    "def fpp(r):\n",
    "    return 42*(sigma**6/(r**8)) - exp(-r/sigma)/(sigma**2)\n",
    "\n",
    "# Initial position\n",
    "x = sigma\n",
    "err = 1.0\n",
    "\n",
    "# Main loop of the search process\n",
    "while err>accuracy:\n",
    "    xp=x-gamma*fp(x)\n",
    "    err=abs(x-xp)\n",
    "    print('x: ',x,' err: ',err)\n",
    "    x=xp\n",
    "\n",
    "# Print the result\n",
    "print(\"The minimum falls at\",xp,\"nm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4.3 - Gradient descent\n",
    "\n",
    "Return to the Buckingham potential from the previous exercise: $$V(r)=V_0\\left[\\left(\\frac{\\sigma}{r}\\right)^6-e^{-r/\\sigma} \\right] .$$ Find the minimum of the Buckingham potential for $\\sigma$=1:\n",
    "\n",
    "3. For gradient descent with numeric 1st derivative.\n",
    "\n",
    "## Talking points\n",
    "\n",
    "1. What do you observe?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x3:  1.842986204589765  err:  0.8429862045897649\n",
      "x3:  1.9337285205513042  err:  0.09074231596153925\n",
      "x3:  1.9256366474164555  err:  0.008091873134848715\n",
      "x3:  1.9171395568522098  err:  0.008497090564245635\n",
      "x3:  1.9087059192629927  err:  0.008433637589217113\n",
      "x3:  1.900342659676813  err:  0.008363259586179783\n",
      "x3:  1.892055113507541  err:  0.008287546169271964\n",
      "x3:  1.8838486810373356  err:  0.00820643247020536\n",
      "x3:  1.875728817146367  err:  0.008119863890968526\n",
      "x3:  1.8677010127639664  err:  0.008027804382400694\n",
      "x3:  1.8597707746191081  err:  0.00793023814485827\n",
      "x3:  1.8519436033763719  err:  0.007827171242736242\n",
      "x3:  1.8442249702467572  err:  0.007718633129614627\n",
      "x3:  1.8366202921960026  err:  0.007604678050754599\n",
      "x3:  1.8291349059081174  err:  0.0074853862878851984\n",
      "x3:  1.82177404069769  err:  0.0073608652104275585\n",
      "x3:  1.8145427906004077  err:  0.007231250097282205\n",
      "x3:  1.807446085905991  err:  0.007096704694416589\n",
      "x3:  1.8004886644303817  err:  0.006957421475609404\n",
      "x3:  1.7936750428531982  err:  0.00681362157718346\n",
      "x3:  1.7870094884710803  err:  0.006665554382117911\n",
      "x3:  1.7804959917362289  err:  0.006513496734851465\n",
      "x3:  1.7741382399611736  err:  0.006357751775055265\n",
      "x3:  1.7679395925744794  err:  0.006198647386694178\n",
      "x3:  1.7619030583069786  err:  0.006036534267500837\n",
      "x3:  1.7560312746736126  err:  0.0058717836333659346\n",
      "x3:  1.7503264900918905  err:  0.005704784581722189\n",
      "x3:  1.7447905489444127  err:  0.005535941147477796\n",
      "x3:  1.7394248798504008  err:  0.005365669094011816\n",
      "x3:  1.7342304873605752  err:  0.005194392489825672\n",
      "x3:  1.7292079472323443  err:  0.00502254012823089\n",
      "x3:  1.7243574053795498  err:  0.004850541852794432\n",
      "x3:  1.7196785805249886  err:  0.004678824854561281\n",
      "x3:  1.715170770516399  err:  0.0045078100085895745\n",
      "x3:  1.7108328621998425  err:  0.004337908316556449\n",
      "x3:  1.7066633446804085  err:  0.004169517519434063\n",
      "x3:  1.702660325741001  err:  0.004003018939407577\n",
      "x3:  1.6988215511372982  err:  0.0038387746037027437\n",
      "x3:  1.6951444264424433  err:  0.0036771246948548164\n",
      "x3:  1.6916260410796329  err:  0.003518385362810461\n",
      "x3:  1.6882631941553894  err:  0.00336284692424349\n",
      "x3:  1.6850524216911718  err:  0.003210772464217637\n",
      "x3:  1.6819900248460788  err:  0.0030623968450929073\n",
      "x3:  1.6790720987283123  err:  0.002917926117766534\n",
      "x3:  1.6762945614069005  err:  0.002777537321411838\n",
      "x3:  1.6736531827569985  err:  0.002641378649901993\n",
      "x3:  1.6711436128005708  err:  0.002509569956427704\n",
      "x3:  1.6687614092379768  err:  0.0023822035625939986\n",
      "x3:  1.6665020639036026  err:  0.0022593453343742187\n",
      "x3:  1.664361027918721  err:  0.0021410359848814675\n",
      "x3:  1.6623337353557126  err:  0.002027292563008487\n",
      "x3:  1.660415625268787  err:  0.0019181100869256618\n",
      "x3:  1.6586021619858606  err:  0.0018134632829263087\n",
      "x3:  1.6568888535936805  err:  0.0017133083921800907\n",
      "x3:  1.6552712685829227  err:  0.001617585010757816\n",
      "x3:  1.653745050651037  err:  0.0015262179318857871\n",
      "x3:  1.6523059316881035  err:  0.0014391189629334722\n",
      "x3:  1.6509497429944753  err:  0.0013561886936281908\n",
      "x3:  1.6496724247984258  err:  0.0012773181960494817\n",
      "x3:  1.6484700341579113  err:  0.001202390640514528\n",
      "x3:  1.6473387513424609  err:  0.0011312828154503851\n",
      "x3:  1.6462748848000022  err:  0.0010638665424587135\n",
      "x3:  1.6452748748189672  err:  0.0010000099810349994\n",
      "x3:  1.644335295998996  err:  0.0009395788199710875\n",
      "x3:  1.6434528586440795  err:  0.0008824373549165543\n",
      "x3:  1.642624409190566  err:  0.0008284494535135689\n",
      "x3:  1.641846929779443  err:  0.0007774794111230499\n",
      "x3:  1.6411175370779354  err:  0.0007293927015075141\n",
      "x3:  1.6404334804500915  err:  0.0006840566278438764\n",
      "x3:  1.6397921395699884  err:  0.0006413408801031473\n",
      "x3:  1.6391910215647552  err:  0.0006011180052332143\n",
      "x3:  1.6386277577675459  err:  0.0005632637972092791\n",
      "x3:  1.6381001001539695  err:  0.0005276576135764266\n",
      "x3:  1.637605917528315  err:  0.0004941826256543713\n",
      "x3:  1.6371431915194479  err:  0.00046272600886720916\n",
      "x3:  1.6367100124396685  err:  0.00043317907977935555\n",
      "x3:  1.6363045750535399  err:  0.00040543738612863756\n",
      "x3:  1.635925174298235  err:  0.0003794007553048484\n",
      "x3:  1.6355702009911381  err:  0.00035497330709688413\n",
      "x3:  1.6352381375557834  err:  0.0003320634353547458\n",
      "x3:  1.634927553792511  err:  0.00031058376327242065\n",
      "x3:  1.634637102715945  err:  0.00029045107656600067\n",
      "x3:  1.6343655164779316  err:  0.0002715862380133771\n",
      "x3:  1.6341116023908298  err:  0.0002539140871018475\n",
      "x3:  1.6338742390635737  err:  0.00023736332725610154\n",
      "x3:  1.6336523726593704  err:  0.00022186640420329695\n",
      "x3:  1.6334450132828635  err:  0.00020735937650684377\n",
      "x3:  1.6332512315012806  err:  0.00019378178158291348\n",
      "x3:  1.6330701550029463  err:  0.0001810764983343205\n",
      "x3:  1.6329009653951263  err:  0.00016918960781997505\n",
      "x3:  1.6327428951415286  err:  0.00015807025359770854\n",
      "x3:  1.6325952246387652  err:  0.00014767050276343063\n",
      "x3:  1.632457279430349  err:  0.0001379452084162569\n",
      "x3:  1.6323284275558598  err:  0.00012885187448907764\n",
      "x3:  1.6322080770322465  err:  0.00012035052361336618\n",
      "x3:  1.6320956734638914  err:  0.00011240356835506837\n",
      "x3:  1.631990697777742  err:  0.00010497568614931119\n",
      "x3:  1.6318926640791225  err:  9.803369861960398e-05\n",
      "x3:  1.6318011176238276  err:  9.154645529485528e-05\n",
      "x3:  1.6317156329025844  err:  8.548472124325635e-05\n",
      "x3:  1.6316358118321612  err:  7.982107042314723e-05\n",
      "x3:  1.631561282050166  err:  7.45297819952917e-05\n",
      "x3:  1.6314916953070249  err:  6.958674314105906e-05\n",
      "x3:  1.6314267259516473  err:  6.496935537758652e-05\n",
      "x3:  1.6313660695066496  err:  6.065644499764211e-05\n",
      "x3:  1.631309441327713  err:  5.66281789367018e-05\n",
      "x3:  1.6312565753423998  err:  5.286598531317743e-05\n",
      "x3:  1.6312072228669403  err:  4.9352475459452094e-05\n",
      "x3:  1.6311611514930522  err:  4.6071373888123546e-05\n",
      "x3:  1.6311181440440095  err:  4.300744904273124e-05\n",
      "x3:  1.6310779975941538  err:  4.0146449855615884e-05\n",
      "x3:  1.6310405225505122  err:  3.747504364159937e-05\n",
      "x3:  1.6310055417895828  err:  3.498076092944835e-05\n",
      "x3:  1.6309728898504647  err:  3.265193911805575e-05\n",
      "x3:  1.6309424121762222  err:  3.047767424257941e-05\n",
      "x3:  1.6309139644057273  err:  2.844777049482161e-05\n",
      "x3:  1.6308874117084844  err:  2.655269724294307e-05\n",
      "x3:  1.6308626281619845  err:  2.4783546499929088e-05\n",
      "x3:  1.6308394961702095  err:  2.3131991774949512e-05\n",
      "x3:  1.6308179059177381  err:  2.15902524713929e-05\n",
      "x3:  1.6307977548597234  err:  2.015105801467243e-05\n",
      "x3:  1.6307789472446599  err:  1.880761506356521e-05\n",
      "x3:  1.630761393667314  err:  1.7553577345896798e-05\n",
      "x3:  1.6307450106521817  err:  1.638301513229301e-05\n",
      "x3:  1.630729720261839  err:  1.5290390342759252e-05\n",
      "x3:  1.630715449732808  err:  1.4270529030913082e-05\n",
      "x3:  1.6307021311332117  err:  1.331859959630144e-05\n",
      "x3:  1.630689701044225  err:  1.2430088986770116e-05\n",
      "x3:  1.6306781002621957  err:  1.1600782029219658e-05\n",
      "x3:  1.6306672735176415  err:  1.0826744554215395e-05\n",
      "x3:  1.6306571692189185  err:  1.010429872305707e-05\n",
      "x3:  1.630647739204943  err:  9.43001397546439e-06\n",
      "x3:  1.6306389385189215  err:  8.800686021492865e-06\n",
      "x3:  1.6306307251968688  err:  8.21332205269698e-06\n",
      "x3:  1.630623060068815  err:  7.665128053835346e-06\n",
      "x3:  1.630615906572724  err:  7.153496091039102e-06\n",
      "x3:  1.6306092305804385  err:  6.675992285432031e-06\n",
      "x3:  1.6306030002393515  err:  6.230341087043456e-06\n",
      "x3:  1.6305971858167199  err:  5.814422631589267e-06\n",
      "x3:  1.6305917595640926  err:  5.4262526272275124e-06\n",
      "x3:  1.6305866955792285  err:  5.06398486410653e-06\n",
      "x3:  1.6305819696900012  err:  4.7258892272861175e-06\n",
      "x3:  1.6305775593302  err:  4.410359801143571e-06\n",
      "x3:  1.6305734434449097  err:  4.115885290367771e-06\n",
      "x3:  1.630569602380238  err:  3.841064671661343e-06\n",
      "x3:  1.6305660177924655  err:  3.5845877726092823e-06\n",
      "x3:  1.6305626725605515  err:  3.3452319139204434e-06\n",
      "x3:  1.6305595507070443  err:  3.1218535072596865e-06\n",
      "x3:  1.6305566373212976  err:  2.913385746650121e-06\n",
      "x3:  1.630553918486741  err:  2.7188345566031558e-06\n",
      "x3:  1.6305513812144814  err:  2.537272259628409e-06\n",
      "x3:  1.6305490133819114  err:  2.3678325700160485e-06\n",
      "x3:  1.6305468036767354  err:  2.2097051759484287e-06\n",
      "x3:  1.6305447415429095  err:  2.0621338259196875e-06\n",
      "x3:  1.630542817124191  err:  1.924418718601828e-06\n",
      "x3:  1.6305410212272642  err:  1.7958969267262148e-06\n",
      "x3:  1.6305393452707513  err:  1.6759565129031984e-06\n",
      "x3:  1.630537781241573  err:  1.5640291783913796e-06\n",
      "x3:  1.630536321671583  err:  1.459569989981091e-06\n",
      "x3:  1.6305349595818706  err:  1.3620897123534093e-06\n",
      "x3:  1.6305336884622714  err:  1.2711195991332858e-06\n",
      "x3:  1.630532502242503  err:  1.1862197684120446e-06\n",
      "x3:  1.630531395245893  err:  1.1069966099341855e-06\n",
      "x3:  1.6305303621837153  err:  1.0330621778020799e-06\n",
      "x3:  1.6305293981237465  err:  9.640599687443796e-07\n",
      "The minimum falls at 1.6305293981237465 nm\n"
     ]
    }
   ],
   "source": [
    "from math import exp,sqrt\n",
    "\n",
    "# Constants (V0 is set to 1)\n",
    "sigma = 1.0             # Value of sigma in nm\n",
    "accuracy = 1e-6         # Required accuracy in nm\n",
    "gamma = 0.1            # approximate of inverse 2nd derivative\n",
    "\n",
    "# Function to calculate the Buckingham potential\n",
    "def f(r):\n",
    "    return (sigma/r)**6 - exp(-r/sigma)\n",
    "\n",
    "# Function to calculate the first derivative of Buckingham potential\n",
    "def fp(r):\n",
    "    return -6*(sigma**6/(r**7)) + exp(-r/sigma)/sigma\n",
    "\n",
    "# Function to calculate the second derivative of Buckingham potential\n",
    "def fpp(r):\n",
    "    return 42*(sigma**6/(r**8)) - exp(-r/sigma)/(sigma**2)\n",
    "\n",
    "# Initial position\n",
    "x1 = sigma-0.1\n",
    "x2 = sigma\n",
    "err = 1.0\n",
    "\n",
    "# Main loop of the search process\n",
    "while err>accuracy:\n",
    "    x3=x2-gamma*(f(x2)-f(x1))/(x2-x1)\n",
    "    err=abs(x3-x2)\n",
    "    print('x3: ',x3,' err: ',err)\n",
    "    x1=x2\n",
    "    x2=x3\n",
    "\n",
    "# Print the result\n",
    "print(\"The minimum falls at\",x3,\"nm\")"
   ]
  },
  {
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}