P03-TaylorSeries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taylor Series\n",
    "=============\n",
    "\n",
    "Suppose you have some function that may be expensive or difficult to evaluate and so you’d like to find an easy approximation for that function in some limited domain. One particularly nice way to handle that is with a polynomial approximation since they are easy to compute. The question is: How can we find such a polynomial? One easy answer is called the \"Taylor Series\" set up like so:\n",
    "\n",
    "$$f(x) = A + B(x-x_0) + C(x-x_0)^2 + D(x-x_0)^3 + \\cdots$$\n",
    "\n",
    "So.. $x_0$ is the point around which we are interested in finding an approximation, but what are the values coefficients $A$, $B$, $C$, etc. of the binomials $(x-x_0)^n$ ? There are many ways to answer that question but *Taylor* *Series* answer is that the value of the function, and all of its derivatives *must* *match* the value of the *series* and all of *its* derivatives at the point x0. Let’s try it and see what happens: Just substitute $x \\rightarrow x_0$ in the expression.\n",
    "\n",
    "What happens? Well.. $x_0 - x_0$ is nothing but $0$ so all but the first term go to zeor and we just get:\n",
    "\n",
    "$$f(x_0) = A$$\n",
    "\n",
    "We have our first coefficient, $A=f(x_0)$, so $A$ is simply equal to the *value* of the function at the point of interest.\n",
    "\n",
    "What about $B$, $C$, and all the rest? Let’s demand that the first derivative of $f(x)$ match the first derivative of the series and see what we discover. What is the first derivative of the series?\n",
    "\n",
    "$$f'(x) = 0 + B + 2C(x-x_0) + 3D(x-x_0)^2 + \\cdots $$\n",
    "\n",
    "Now, put $x \\rightarrow x_0$ and see what we get:\n",
    "\n",
    "$$f'(x_0) = 0 + B + 0 + \\cdots $$\n",
    "\n",
    "so $B = f'(x_0)$, see how easy this is!\n",
    "\n",
    "Now try the second derivative:\n",
    "\n",
    "$$f''(x) = 0 + 0 + 2C + 6D(x-x_0) + \\cdots $$\n",
    "\n",
    "putting $x \\rightarrow x_0$ and we see:\n",
    "\n",
    "$$f''(x_0) = 2C$$\n",
    "\n",
    "so $C=f''(x_0)/2$. If you carry this on for a bit you can see you get:\n",
    "\n",
    "$$f(x) \\approx f(x_0) + f'(x_0)(x-x_0) + \\frac{1}{2}f''(x_0)(x-x_0)^2 + \\cdots$$\n",
    "\n",
    "or more generally:\n",
    "\n",
    "$$f(x) = \\sum_{n=0}^{n=\\infty} \\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n$$\n",
    "\n",
    "Morse Potential\n",
    "===============\n",
    "\n",
    "Let's use this technique to analyze the motion of nitrogen atoms in a nitrogen molecule. The atoms live in an attractive potential of the form:\n",
    "\n",
    "$$U(x) = U_m \\left( \\left(1-e^{-\\alpha(x-x_0)}\\right)^2 - 1\\right)$$\n",
    "\n",
    "Let's load up sympy and define all the symbols we need: $U_m$, $x$, $x_0$ and $\\alpha$:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Um,x,x0,alpha=sp.symbols('Um x x_0 alpha', real=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let's define the Morse Potential function $U(x)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAAyCAYAAACztwV0AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAMWElEQVR4Ae2di3XVOBCGSQ4FhFDBkg6ArQDogEcFsB3AoYI90AFLBTw6IFRASAeECgjpgP0/oTGybF9Zzs29lq90jmNbz3lpZjSSb/Z+/fp1bQlpb2/vlfD4JHyOl4BPxWFeFPDydSCobuk60/VCsnYxLyh/QyNYgfOt4HuUgm8/VaGEcmNOnfwlcKs8GL18vZF8/aPrgTA41PV5rpgIRhTTG8H9NQVj8QrAM+cWzEkhW8srBSZS4GHU7l+935bs4Q3MMmk+4AknlUDRCkAMgDHPdD2dJRcqUEuiQDjZzfUP82aHq5TAfwLq3BvJXviu9+YWkCmk3DpHoD7yLk8BUFcQS6SA5Osogvu2fz+J8uf4Shzgu+bLF+HxMQZwT5lxXhHvQuiTAD0T/NX1L4JjywFSsvdN2BATeF0CVoIXT5mg4I0Y3iIVgEfog5C5IaTMHYtxW/muPlg6XFN73KSadoQC4rvzHMX3ZIS8jyRqz27TgdoXZXgEN0rro+B+EeJVqgIAmeOpTBAxnqv9E7W/ExKjPu8GBcT/+8L0VS7/1Q6jcWeq3G2TuoIdLwCjeST42cZ0qbggoGcCwRc0cXbyzH+phveyG9cGi6CAJgAR8neSBSbEqOTlhsnjLD/vuiwWMKqPbVYS3Kz/3fmFEI7iPAARHevP2p/92Kyktrh/33XdU/vTrMa18uIoIHlgn5y1/MploJ/ob1WX7T9LKIKiAtDCAw/mja5m6VyUAvCMgGkPxLTsE39q7zS+2k5a/2ncnU+iIWcuGhdyjgRBTgRjUsGDi+DHoDQTog8f1fupfIxHK2mMvVbGzF+EBziAC6cYXQCzNAWA9nos4DvRzBTtEQrVQXm01kCpdrX8DwVEQ9bNrSDSn9L5PAlO1rsXgjVpJFR3p4yC8HWnA0UbF/8qLQbwWIxNMnVAFHHhiILO2noNwL71bD9RUMCzT+Ix691HghkLn0q49Q9Vt2PhUw0LLX8nuJtTjMUoADEICw6T2P/PSl4QaD9765WF2IYqe9oTdylJecLrZKBYOLFU4HLbwhsi6TaHscNA7IRcK0YBCFYHsO5TPACEoTQBFsizSeyaYDmKSZrYFwL2llf+KbjxbMBx8ckrcWjjguglHQUGYNZ1U6wQS4cwgrt4Ro9BUJMDC/lD101dzeEW5fOMkFjqBNWCthyTZQKhoP/W9a/aJgNwqteb1twvSgu4Vkb5VY5ROdDY9wX7FAOj5kWlE0HbbGGyTiIQQjSUc8FECXGzn4kYnJRzl8+jnItAwnMr28Rd4zm4csdSO/ADZoS4wWfXn0UP+I3AG3+RAcdT3Qn2WT70+2bv3KmnCyXBM+2++udWn2GbMc/qJ9mv6rAMRHHhsnN3cPT1rzJ3+KWvLM5TXeSrwTsuX9I7eOoSSuKfIWaZujdCYWV2VxkM2ujEd0D+ZjrSxp5tA/OYZ7VBUFR1fDvVR3AQiEHhyulvbnWFFxabH09paKl36IRix2I2MuDf47qNMvVtWsYi7DfnWX0l+1UdjJMpH5RBC7ZwPJWBy2B5VLdRZGH+nJ+F2yQ5VTs3J6DPdf2xBLEQiFUuEO7eNgJpWCESFiY3sd2RXDbg9qoeOwXnug518T77JLgRAg6lpFJ4aIUlUfwJNTSCztSL+7sIO1d56OLjSnZkZgpcY/rVWCgnBw93jXMX3lleCKeejZdRdu/rF+VCy1kncBWAl5VTlgCk26EFwMKu1JYq/yBCN2029axxYQzwPcwdU21QGivxivtUfTwdxluqBwBuHOhpeKl3JnIHZ5/fSz9rY/3wbs+XuQ/1q3wU1M+wb73D395xlY9RGyWzqmcy1qJLONbcngXzJDlVO5QIvH6zrz/XpFWc9dcj7lVv8ponaUl7G18+E4tMmjI+QjOlnRtwoX+w4NDFJc/bJ7xIyLGqCIgl8wzcu8qIrCM8JNqEtHV9uJLMP1fUbw7vDQ8U4aITPPYIHtoSwM7Vd1y5gBK4jR0FIcahPGgPsW15YK4UUWF+iOB1UI+oM/UJuKwaT1WaZAJpgDcFIx+mLB1Gdn211UQ3cGeLCrqRWIZBOxNYl5n5h6PQL9Q3wg5N8XR4JxJOgAh32O0XKx+FgJEwNxu3+qPesT5OHvwzOwmX2WlJ9gvOwKFxwoRxGKIFyz/W9mOS9WHGZkyb0uscOBdQWDBBRN8/LmH8rPLeAJzyXeRUd4SBqwkS6hlmYS0QljCwhIJouXLxeOG76rqope5Z7pnq2/hZQSq1m+RahTCv41lwoFzhTYM3z7o2Gq3WeNAje/m1DhrEfQgOZMwtzbjzHtexd5W53Ql7T91V38lqqt5cygXvZDlVW4LcX697jYpQpaxxRzOqLRYEa0GijwsRJ/yVFGtzU/lh/+THmpw+hpLVxUrkJBs/p80s6oq20BPr9VS0c9bJ8woviyj+xhI81dhMPDtFtrGxewbCe3kleNixwBuKA5auicoxMlPohCezK+mAJQBWhjToKomYCKNNdFfZ/+F03anKmaDUsaWE1UFBkGLXkHqhQnCVVvxxE1ljXayo01c0VXH09TWY5/H/rAo23mBdXwAeTGzW4kMJr4d0qP7R9CTabev36PmFWTyp1KEaB+hV/fEy0DvpbUzPD74YXVnP6kf3sTyMmnVfr0guugNNy8GYuhgA63SSbQ38fmv/hZAdbeqZQU2nRPQeT2omOkoCwQ0T9WOlEJav+/lKPQGPH+vNdSYsGB8vbXXCGUKCg3U/H80QcV+luKzJNu/sc1s8amtwXJFcrBUfPACXEkxlDWpBEmsS3pnofULBRG+5jRIg8tCyLh8t6QmlrMGEtgoDUYMVo4JY8UTFs39dRfONAy8+tXi5cQBGDig4w2XoqFbIoa9YusyMwtdXOt/Xg3PtRQBc+E5SPkuDlDZ9rDot6+8JSp+sHcPEGu4UheLHRCGkkjEl15I7xaHOjbmpceZUbjh3YBLdbEnQKasZkylgslXsjlEm5uB7se+1Ota7NcklZHwzzOTn445BS+QnMRMsnugoBbYWWopBWXd1WR7/ammMVbEtsKyJrL5tEpUY2GGJ1FKOojXbdCzFxtBM1WrKoIDJlhmNjKbFVj2/DuiaKHckWERWmcRMGojAByBY61TCyrPOt0lt9SFon6Ai2A80Fuu0d1Y5cbeJbFo6Ub1VTNte76ZVSy9+cjGGTbzPykP5sdW00bW4xiPybjBhlRwNlE88pqb1U8Bky2Rt/SOsqcc1yemBwDkr4ifBhDABMbyR7N8CVFuU2qEmzrqDdOq2pqVQQHKCQcK7Wvn7gEvAV7gy+TkH8Hq/EIRsCXJ7Arwsb0Z5ABP6rk2WQwHOFHCOZfYewBpIbvPhRxEKQExhEpNgUm4iyMnaGa1XU6XAEAXwMuNl7FDd0vNNAZwVoQA8tfECCCDmJmOqC0rmNq71l08BbxyYFHEge6nI29mf05IUABM5ewngXTraxqcUl8rcilc+BSzoa8Yiv4eyWjgPQHOjKA+As9/uMNAEWhNANCZPaF6bLJwC7iAbE2LheBp6zAWn7EryAE489FMm8nu1JQ7AOq+mSoGYAiwPO0fd40pLeNccwPoTD3PLnWIUgA8EEqF9kssIvwzgoNPL3La1/rIpoAnhTlVKRjZ6zmOLVDUj6M7oFHEOwIglZqGl+SJt0v9kU3v2Pus/BjWC1jtLSmSCLyx3QgEIXyw/3/a4HbViPAAvq6zlYdqUZQBN8QLe8lBTpYDkCOvPKdZdmfy4/sydZrlTlAIQowhcEKgZc0RZ1drJGO0Z3y6sbztFAckAa2GWhJNkqVBicdqR1Ci8opYAQC7GXerIptqjBdlR4Oev7YARXde0IxQIZADXv+97lUVSQniz3DkWzo3SK04BwBmPyHshMunDGLVH+7MW4tuCXdn6gXQ1iQLiP7zn58IbS7h0wghngn8soY9CmS9VARgykz/c8EqAHyOpXsDSpT/AT3zHA7wrvu/KoR+HvfDG6z0R3i2jWaQCAKMhhBy29U+lQKVAQwHNFQwmwe+/pADYSm9SyQoAN57v5O9UK97wsz5UCrQo4D2e78ok5tXxeoraBQgxEzK2G+C2BsOy+lwpUCnQUID5QbysM/mpUawCAHghRQT3WFquKgEIUlOlQEABzQu37ad50lr3B1XKVgAgYsgJWXekM0SuPlcK7CoFNB848MNvbq78CrZoD8CYKyTZ1zzySFt2vVcK7CQFNA/Y6cDq30sR4H+38vVAx71ioQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle Um \\left(\\left(1 - e^{- \\alpha \\left(x - x_{0}\\right)}\\right)^{2} - 1\\right)$"
      ],
      "text/plain": [
       "   ⎛                  2    ⎞\n",
       "   ⎜⎛     -α⋅(x - x₀)⎞     ⎟\n",
       "Um⋅⎝⎝1 - ℯ           ⎠  - 1⎠"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = Um*((1-sp.exp(-alpha*(x-x0)))**2-1)\n",
    "U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, Let's take the derivatite of $U(x)$ with respect to $x$: $U'(x)$. \n",
    "\n",
    "Remember the force is $-U'(x)$ so the point in space where $U'(x)=0$ is the equilibrium position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python expression for Force:\n",
      "Force is: -2*Um*alpha*(1 - exp(-alpha*(x - x_0)))*exp(-alpha*(x - x_0))\n",
      "Latex: 2 Um \\alpha \\left(1 - e^{- \\alpha \\left(x - x_{0}\\right)}\\right) e^{- \\alpha \\left(x - x_{0}\\right)}\n",
      "Pretty math\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASkAAAAmCAYAAACYuG3jAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAKe0lEQVR4Ae2di3EUORCGzy4HYEwEZzIAEwGQAY8IgAygiICCDHxE4IMMgAg4yAAyAJwB93+yNJY00uyMvTuP3VbVeKRWq9XqbrWes9778+fPXxbKEtjb23tGjmT0TxnDoNsmAen8UG16J50/2ra2LbU9B0tlfNN8y1hfqI4nMtY7m67L6M9HAtL3uXR/quer6X4eetmTIubByYy4kIHeFzvv9fyN0c6INWNlJAn4Qequ9G8zqpFkXqtmv5axq3AZJ9N9HNQ9c1DDrUDyOx5eavwS4vN2V63S/VvlHwvPLfm7cJeSt1TdmJNqW9g7gT7JSL+1swzSJQF1gjeS248unBnl4YCYMXcFZlEs/Ri4Fh2WrBtzUpHpSZGMrg/1vIzAFu0hAcmO2edpD9RZoMiZfhAjj8R3debnHS54DFyLDUvXjTmp1PQwxg8Lmg2k3E+U8s79xwLlxmD0ZoXYXiv/odq4yNnUNujGnJS3UCmTEZWZlM2ivEwGvF4J92wA/ixQ5VQ5FGHZ1zWbYtnPs9S9qcXrxk73fHeRobJUuS/DveVBO/+STJhl/NRzU8+hZPMcoQhOnA7ugtLfS3KLyiNTJ1+97+p5Lfwr7/mtk65ocdXkXPxU78IJBwfFftsNvWcRIhlsvW64qGjPxTWM37K+FyaPC3uQLL7rwWk7+1CcPScnH73psAHOLAQnldgRuHpwZnq5stw7Ip7QzcutSqv8SrrCYWmGg3XOhXSNrvLYg3xfyweuQBuJNPLowt90nvhIZKj0VuvGXeaUV0YJYZkTpr4vJexmtPM4KJ4lETiMpP/pQcHNKCS8j4KFUxPKnymf49zZBt82DPvTbJkckTHJg1kPe0yxPNDrc+WhU+IhYAulEz1OSMNsCxxo4siuO1PtQ5dO+4j6xS96Jf1ATynAIzjVIDo/RAc8aMQyqZbZVMYu6ubAK5GRsbm0JhjO6KveDwR3SkFRgnEaQh6jGUbQUphgD4RDPgY5a+cEjz44pyp+G6ccMkpvtY/Rl0322V/29Ly6ZVqpLRHMdWqffqz30yiPKPrH2YCX06MDJyGTJQNby1auwlsfuqqLGY/jibfqOdGTLFEjZn8pfhSla1H4D4NvDWcQ/CrtVwU7qRtGmdZ0WDCWP7+l5Hwa/1VwgVN4nFY+o2aLZowzp7jnt7VkiXmkPXqQFW3jzfR/MW2M27Iq7tt2HOMJhqNptdnDP8a4cTyUCzDSIX6dd42u4DjSxG6VZnlUrFdwHE/ncg8+FRh4Fa3b/Rh58KBnp3SzrwajJBxPHhg5DuXtMc44FEfFGEHxIymsNbpmOHNKIoPSkqXhkfboCbOIL03GdkaYUdLZXWAWosgTEsjBp12e/oQZVkizsc6JmevZAlIulq2j0yAPiGyILu2M+atx5HDgoYYwEnzndHPgFXRbwmdW0OlYhENnJsR7EhcQ/xc6ivZRelJu4gSG11qOTMzTlNWz9H8pXTIgYRPYBmkGLZb7OGkuOQanlZ/4sYT6IFxmH85WfJyTKO4dXTWspCs+2T/CBuNwpETNJvmAnJnxqhDKI5MQX1VmE/k7p5vqFQQpmikynfdGcF6Chf2oO4Lh0VtBOJyo5JuuDk95ODk2H6EbNurZ3yFwNP1FdN9GeBx/g8+eWdGJeJ7AI7Apy9S9iOswCn9Eg1Gfg4Jee2jCp/Mhi0Y2BbKjgMQLHfKVnlgGyGu0juTlgc6d4xql4R2ViB8co9tj8/LBJrCjVlB+r1878HTYAnkuWs1BUYtgBPBlTDepTIbrRgJnNEwe0cSZ0HEx9iZPaRyXQJewPK780xwW0oGe3jDK0xz5K05no04cQHPUqzhOLNljgJ4Cji4/ioUGhuTW7HqTfhjqL709DgSflfJLMOG6/QnKlvLHgql+9OQGk1Cn0sgl0VvI2+RbdVb3pTZZb4m2eEHv7B2GKwjJHk4oo3xsa4jenX2G8l1v0TXdFPyE5DJYN40DigUuQhh+spnoiaOkTmPMywW6gjNNdg5Db+i7ezNRPp0L+kkHUxpDE1rqGAXDGSW44Cjg/Jyj1Btn0ulIlB/qHWKskzspzzcyaJwwbdVD5yxuEucyXGcaPvT0luE6674KLfHqZDWkrMq07LNUXnjYlOkm67MlWZVguW4OBEiCpqgYOXdR8mNmRgZCdf2usiintqnMcuCbnwKDl0+/cWKEfM8CvGT55nnEyHJcyrPMOSGicFN1nl9Eq3+hQ2C/YyPBt/mziIe6VtUDz0+RVwciy03CkejjNAmUY9nKe9SgOtmD4hs3HGQX36Py1VEZg1/YcuhAa2X10aHppiW2QYBEN4mTCsYuI8sdFDWwZ0TgAmctUA4n1wpRx3HOTunE8agAzghHlncw8HNn9Fiw+FKfkk2gPKdLGEperkEqRNhc3UjwbWKDdp2BmQsfQ/faH1lnxTVa4mUWe1I1/mK4eO21/xiXGRA33QwQVo6a66ZxUoyCQr4lhMZBCcaMh7VWswmreNcoydq/wc0r92mcUYkGzigxctUPjJHLwZUOoxjvEg2BXSCfTfhzn+569cHpKj9l3ipZT8nbVtUd2V5fezHdrMkC9qEjBbDU4nZ546A8fRxXCG4ZJ1znuAIwvAVnGdhn+uxmQaEcb28A0GU/KQ4ct37D8fl6m5vEgrOvVQ0qkzi8KuLlMg/HtqRQ7SySVVj+Lak9c+c1zLQ77c43wnSzRm3u+86Pg+GS3hv/8GuEwDhudSOC7/TMXhJHJDzuWIHLl+2do4fwcEQ4g9wZ4biYseVLwBOBAwxeguNhiZPvaeHsWMs2zsbz1qSV1wqiGQyKOzxLCixlmWk2Qe3lHhPL7SCnJs8i15ZAsKNfPSiZbnoIqS8Kyz0cDM6jNPomTkcd+o46AY4MJ0PnRmF8TsKMp0+gntIdKgyg1LFQNt8C4nzO9Lig+vjQ1TlUAZqRTXC3P6M8HBZ7Uj8F67P3QFvgrTOIJg6AETU4h8+CISNOPEfdG6Jdvp3whAxoA44+nw0DtnB9CYSZlJNzFznTTZd0hudVL3MOJ7XcEursOF0+5Vn35vZyhWKcJxKQjTBQMiBMfoE3YWwHEvs70MY+TWQZu3Im1YeQ4WytBPia4VwD2cqZ1NZKYKKGmZO6EPwXvdjPOZxID1bt/CXAIVLYH50/t1vEoTmpC2UG43Mb+FukX2vKGiTgBy9m2vmBzxqoG4lVEjAnJQn5KTyOqnViuEqAlr8TEggHJWEw24lGz6WR5qQuNcEpZzDGS6jFTAIXg5e7r2fCGF8C5qQuZf6vouxLxRdYL3MttssSYBuAkz0LE0jAnJQXul/ycVH11QR6sCpnKgENWu7+oOxj1HtwMxXHJGzZPalM7DJKfmLjnoySawkWdlwC3h74ZQlzUhPZgs2k2oJnNvWuDTbIrknAz6L4QsIc1ITKt5lUQfgyTv4xxZmMs88nNQUKBlq6BGQDx2oDdsBPZSefhy29bUvj32ZSZY3dE5jvA/l1CAs7JgHpnUu93InihwfNQU2sf5tJVRTgR1IMlZ+wMUOtyGkbwdI9eufns22ZNwMFm5PqUIJ3VPxGum2id8hpm7L8LOpEOreLmzNR7P/t6IpdbP8tlAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle 2 Um \\alpha \\left(1 - e^{- \\alpha \\left(x - x_{0}\\right)}\\right) e^{- \\alpha \\left(x - x_{0}\\right)}$"
      ],
      "text/plain": [
       "       ⎛     -α⋅(x - x₀)⎞  -α⋅(x - x₀)\n",
       "2⋅Um⋅α⋅⎝1 - ℯ           ⎠⋅ℯ           "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Up=U.diff(x)\n",
    "print(\"Python expression for Force:\")\n",
    "print(\"Force is:\",-Up)\n",
    "print(\"Latex:\", sp.latex(Up))\n",
    "print(\"Pretty math\")\n",
    "Up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can solve for the value of x where $\\frac{dU}{dx} = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAB0AAAAUCAYAAABxnDbHAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABtElEQVRIDe2U/U0DMQzFU9QBKtigbHASG8AGZQQ6QlEnQGWDY4V2gzICdIOyQVE3OH4vlwQrSVUhHfyFpZc4ieNnOx/OObcAe/AOmq7r3NDA7zz4F88tcCspQxPV/EWuC7H+tfyT/mrFx6e8j0ajCWu6dZIb8AimoAFX4MBleab/sVRJA+ESpyJyjEW+Ba2IGOvqH0EiZU5PzwejPu5FL6RKitUSPBnra3Rl+RLmFMwu6DGoREQAM7CG+D7a5H3xTjGYsCF9Eoz1cWztnNVZ848+m+sqfjxX9cmwWaWzonNUeQshI529qvCZLcqH/32yeVcltUY4FaHkte+KVoSSPFAFcelXsqZKGs5EGUh8tGRvz3ARMuwt+jbPVLPRh7UrM8WZSNYglubO7mBdmR0rR5BnlY+Tm3HSvpU31A2YQqBnoBs4R2/pdaF0weIt1vBDDZJnpXGqjrcITUEaMsivenqPdrN02ROQiGuZKYFCqmdaWJ2fUBVSoAQxY7ypHIH3VGR63n9pgXP9UqtwBEcsVNqH0rKfGYRUriD2X+YpIjuv8h5AS5R70NjFoXT86onp11LZ3RfAVCvDyzKnkwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left[ x_{0}\\right]$"
      ],
      "text/plain": [
       "[x₀]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.solve(Up,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aha! So the function has a minimum at $x=x_0$. Now take the *second* derivative of $U(x)$ with respect to $x$: $U''(x)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAAmCAYAAADKt5LtAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAN+UlEQVR4Ae2di7XVthKGLywKOPekgks64FEBoQNIKiB0EBYVsJIOuFSQCx0QKjiBDkIHJKcD7v8JyUf2lm35bW2P1vK2JY1Go3/GGj1s71tfv379l4XtELh169bP1C49/Hc7KazmtRGQ3i9U5xvp/enadVt9hoAhYAjcMQi2Q0AO4BfV/pMcwP3tpLCat0BAOr+W/l/r+Gj630IDVqchcGwEbqnjOTYCG7Venf4Pqvqtjv/gCDYSw6rdGAE/AHwoG7AVgI11YdUbAkdC4PaRGruXtqrDZ8kXx//IHP9wrQi/u8NLrV9Cct7rq1X6/000d0Xrtn/66EvIL0E/ObopAWuTcXsESrB3UGravDn/gbYjAH/VwXLte3/GkQ8Nb1TgD3X8n4YWPDo9+Au3z4XggFNnhacvMOvHpsbYUh/vVfML0k+ublbF75wrwzZ0TO07R0O0RP0F2Tu41WzenP8AU0LRIn8t5/Ncx2NdX+r4MIBFGH09UZkXQ8oZrRu5slryuhQsZCPvJOtT2U3nSoUfzEDLoLDYoHYWo59c3RSrjJ0JPkffOaVJS9Rfkr2DXdPmzfkPsyicdhxeKXJPRtDZuccFdE0H/853+I0si7YhIIxZQv9cIG4M8hg09gVs6YnaWeTsv1D95OqmT3eW34/AHH1nfy3tFLPWX6i9g05l8+b8242lLSd29NeeKE5rK8esHzqcmM36W1FqzXipnN9bc3eaocEKNsJyW6eNiI4tII5S9/6L00+ubnZqWiWKFd8Dg/rOmRo7Z/3F2TsYxjZvr/oNsCoB932DPDzQ9WcjvS2K0y9x9trWnsnpcorMir/o+E7HhTB+DlOlcx06CJLuKX7yjERUHt2wJcAe+0Mdr1L0Ss8KM/Nl0IJcfd9yQH7w4CHAXYQIh3PVT65udqGPUoXQvTi175zU9Nz6D2Dv4PjN5gUKowE7RmAgEP/S8UsufqL9Zwh9Lt9S6Tx+PwT5FWfP2OGpMw/2ObvUNSP2v0I8nKHVwSBBJ1eWd+a5Ri8V30Cfe1bZXr6iudCBo2amztnJkapDeSw5vk3lxWmioZ00YLTsMb+p15KjhqPiZ6cftSlLN1OxtPJ1H+NtK7vvnBu/VP0+rbr3FD87ewdHBWfzszh9MaPTYtbC8d4fzNQq/kqDBjDpUJAARwjtzw060sjn+KhjMwOJ5WpeSy46fB7+q9rYdS3a0LHXcOkqc855YKfjfdxGxXGk6JxZcnwTEq/RUk6hwlLXlKvZUsx7yHUOX9Fgp2HgwUDgRL5Qp/KS8of8+Cxa7otq4BPnrXktGQ6hnyG6WQJ/1W99Z2YfOhf+3F/Yd8yPuI7aPaz4WfdHd9TASUHLJG4GJCCrj5T4pZOPOj9W+h9UoDOvZ/HkM8Azs3oa8nRdBaU9Fg35lNnN8mcloC4kH0ZBx++WqOO8jmscAG06WbpuK6N6GKG90bHrDwF5OXOwQOdhKf9HtetZo+3YCJ0hdE1+oVxVpIElWzDO1ioCXYyRLYevWDM4cTJxVj0PdGATJ3KK9m8dlzpyAm1wtpJDnEMzBgPxPYp+hugmB+5sGuxFxAz0rO/MRm3cPR3YC/O2vvso9g4U32xehpc1c22jEyNm8ydLnkpjBvNPs5zSmKEpub1e5TMKO+HZVWatPMlFx1zNzHy8moG2ySE62nSydN2kFw0dAphCz5kp7i6xaMo+JO7bdTcuozQc+El7fXptVJ4qF9KgD9dTzkGewCPw1ZkBSs22FWdFK1mv0rGZ3mV/6lFg4KvL9vtjjTxk0HH2+lEbs3XThrt4MEgfvEKpMtZ3ttwzbVhPSfe6TvbdyjuEvYOfx+Htba4mBm4eHHozMIO50EgrPBQX8pMztJDpz5cSMjWDapCtG/VtYeXiSte8lsVN7x7iy5AEnHo/TkO7dYSZ71UG31JJWAHBibogLBn0/EQEDHzc5eknrAiEODN6nqB3lqxEysXYOj4V8YCLhfjSzli+LokcHXJ0Ea2QdxT9DNHN3LBb35l/X0zCXvcTfqer7z6KvYOjs/k7kxD9VpjOinfdmZ12OmzRYOwE9kuTAT7KyO0okzwWTPwg3sjHiL0Kfe32hAB+sixdMTneBUudL/xNid1gP8QZMLqbVGl8+CYMBnD4sY2xdPVOacyUnT35a55Kf6VjbOjlKzk/I0ujgkvF2+yWP26q2UyjbBwNPOiswnWcv9b1UfQzRDdzY49+re+cG9U0v76++yj2DjrO5hf7Yx91jiyD4vD+HZyj0sJ+/32lMdI6CaJhT4bX4U4cpfIYPPBlPfgy4yYw+yY81HGlcr9FdF+UBj1LPSf8lI5DQSboCLyOwvJsktZRjPxRPcxSXyBfLguVwbEhX4Vhbtk56SQHju6ljhgnMF3NOXkssAs3IJizfWN4SR4GHO4ZBo8PdoOtnQTlZ/9zn+fFlhlfkex7NdDV5cuYfiLkhUmWfoboJmJfuxQP+iC2SLLv7RqDRkT8Dt13qv2b9zcNlZxEJSN9c9n9kQyWmdWsh0DBSePsqv0V6lDAqHXZXp/ya09hxrTKc/x05sbmqPbZdI3BUAlKiZ8U58as7c96WRhAIE9MCw86XrfXqTPxJ7EMY649H2Qb9DS66GkL5ZjxzqqjXH6qG126ziiUURzsaroNeUueVWfrvv+S9aZ4SxZsg+cyGKwyQKvtj4cySsf+hurd2XHg0XUWb9NP4t4QLr36GaOblC48n6ovStHkpnl9ov/a/aX4IfpO3/5d9Dd9OpOsRfdHizgUb6i1B5yUxs2IUXcCpvxauaAApbMM6hyxzhiHe6c7yschpW4aOmeR1Z2n0pKvVCmdQYUbgOiM853seMUjyDbUCWzq/L3c4FQNgMBDB04v+XBbE+c548ihYxCGc9Y/lJdkdViNKHdixyke4o9dmX4a93YKq2baWN00+RBXwC7ncv70bbU+UPFD9J1q567sOaXrOE3yFt0f3fFLLB+8genUG65F8UwgtC3b4xj4x7rmq1rMUAite5+SBeW3PeTGEssnLy90zSVWBgeE5n4vdLVlfPFARm6oJq2S3D7rAy4UvlOdtHdqoC4C+8mLhbl1KUGZ0RIuxZuBCAE82L6YAxfHMPdHdbLHz4OWDDyS9pfLayU6Bp5he2polcFmusqZfrrQ6c4brBvZHXiHfizmjq64R1IPmtJvsZ/cG1T+6H3nruy5T2HF90dqwGyzf4GFg0gu2ysdxTJMbp0xeprk8mmQUzSMthSty600bpyTV+mUxsyoNir3ackVCOUhp1sV0HnyrB85FRiYcDFo1ir6rWf+yFybhTRxt3jdDufAw9tK8j6K+Xs600+jL4gxWuNaepg88xePw/edZs/z9yVd9n9HgM8SNGrlBvhelVUzfqUxQ8dTVw+G6bprxobjr2hbBGMmn+LBiLz2MJjqJ41RuUtXPMymOKd4KNkF8q8ky7WPTz3NxWeqHGPK9+ljDE8r04JAZKO5NmP6acGylGTp3PrOG2WZPd9gsejV7Tm4y3iZ2fI1v8rxe74YdQhuOV+0bkAQEsNZ6WwH5CyR/ii65jI+zhq+7NfHgeW2TwwofL3VV9mUzt5aa1CZ2kCilTAv429PhpwlhVYHJDyZqViYH4FLz7LTPj2N6Wd+/FflqPvI+s4bxM2eb7BY/Gqy8/dOFcfNh1l+9cdrnUnjdSU3kvPOlNl2zcGLjvdcoeVf2DpHfaLDweNAm06eAQErDLVBgZIe6AhpyBIcOq9QNZ8Z4LU/9gHh74KXrYqH9KFn1RuMmnfQSwo8E1Hb4xQmvIfPFkvAsqT2lCBrsLcwYOyS2fTThc7O83Qf0Z9Z33mjJ7PnGywWv5r8nr8MmK/7MXpNBR52qf2Vo+jZU4ceh0gHxz591vuxKosjYi+0yZNZ6EOl1x6sET0rDzh5ZPxT+Qw+XPBycF3NsJTv3qtWHs4N+b4oLUs2GHUF8eQ5Ah6ErMmYKuPrZwYYti2Qm4ERzylkvfud4jsmTbKALXiDkxvErC3DGLlLLSO80TmD2+p/MbraYvrpQmedPOmAfoYty0F9hcpZ39lQkdlzA5AFo5Od/4KynRVrGTUdOp8t5utKFgyBJAKyE1afGHxu+mGnpHCWmERAOhvl/JPMLNEQWAmBycv+K8l5DtUwe2eZz4Ih0IUAqyzXGiS6VZYuQsvbDQJsLdo22G7UYYLkIGDOPweleWiuxIb98ot52BmXM0WAWWR4TuVMm3hezWKgpqPzeaXzarG15hwQMOe/nhZDh+4eTlyvWqupFAT8wJDVoeYDraU0weQ0BAyBQhAw57+SovwyLgOAk7cMVhLBqtk/AuHNijBQ3L/EJqEhYAgUiYA5/3XVxms9oYNft2arrQQE3AesbAm5BFWZjIZA2QiY819Xf/9Tdez7s69rwRBoIsCWEE/6WzAEDAFDYFEEzPkvCm+duV/65yNHL+s5Fjs6AhoQ8j0FPlS16nccjo67td8QOCoC9p7/BppXR88Hfx6po68+OrSBGFbljhDwNsG/JZrz35FeTBRD4FwRsJn/Nppl9v9mm6qt1r0h4Gf9fA3THP/elGPyGAJnioDN/DdSrDp8Pu35uzr8QZ8E3Uhcq3YhBGQHd8UaW7gvW7B3xRfC2dgaAoZAHQGb+dfxWDP2SJU9V+ff9r8Ia8pidW2AgHTPB594p/+ZOf4NFGBVGgIHRsBm/hsq38/66Pz5Exeb9W2oiy2qlv7R/Vvp3pb7t1CA1WkIHBgBc/4bK98PAC7kAOzhv411sWb1ftb/QHq3D/qsCbzVZQgYAg6B/wPdFfSsns2GIgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle - 2 Um \\alpha^{2} \\left(1 - e^{- \\alpha \\left(x - x_{0}\\right)}\\right) e^{- \\alpha \\left(x - x_{0}\\right)} + 2 Um \\alpha^{2} e^{- 2 \\alpha \\left(x - x_{0}\\right)}$"
      ],
      "text/plain": [
       "        2 ⎛     -α⋅(x - x₀)⎞  -α⋅(x - x₀)         2  -2⋅α⋅(x - x₀)\n",
       "- 2⋅Um⋅α ⋅⎝1 - ℯ           ⎠⋅ℯ            + 2⋅Um⋅α ⋅ℯ             "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Upp=Up.diff(x)\n",
    "Upp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, we can use the `subs` method to evalute $U''(x)$ where $x=x_0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEEAAAAUCAYAAADStFABAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADbElEQVRYCeWX8VFbMQyHCccA9LpBskEoG8AGtCPQDWCEHmxAV6AbQCcAsgFsQGED+n2OZewXp++ufxSu0Z0jW/5ZlmVJfpm8vLxsbQpNJpMzzrpLm9IeaKec/3myKU7IDrjgvB5+i/ElbMp4b5OccM+hDysnzBnf0Wbb/GwSmQZBz7kz3QnJ/86JgNngjEaCdLtpkbA89vLXIpkK4xYeUmSYXOR2BbfNnYvGWIzFxNxy0RNN3HFgsi5lztvMuZN6/j30sUkHWCTj/OnJuAyBPIPsHNTysbnAsu7EFuP3xLHrmFYckM+Ubnd3aChAb/qpI/d2Eb9GybDvJrQVnUPcvx5j0wHtLPbN41QbPOx9TARHZugbDcO0UHYVuB53bU/+ljLPQfMCj6pm6u76OvjxMOfjwZuLZwPRKoHRk5KLu6QeJtIHSRfwdsKfbK1tXm4hz7zDz16RtJ14QuoDHWbIdQttRl8YdZ2UnagOi+xpXuXNSPu0G+w5r3CPyMQbwt09wVrkxEk+g0bhChbZh4To/TBZXoDog/PGDfuSP85B6WUIXI+DaYpOjQl9cJ1kK8WTvrfkJhbVUpDp66RebdKR2lNj1WF6+zkMSzqPaht6/RUH5MUqb/I6K1TzX9UD1hlZySC4+u9qgxh7KPUPHW81B9peFjIP22Cz7To3XQRch44WaGtCQ4SXlf0ahV+biWV0KGpyqsaw1oPc1LKq/4DOBRhvS1ykVkAi/b6FIHNxTXhnG9UzxLrE9P1kB/rIns/L7h9+9V40YHquG87IzT1vqnktYq08Y1Io1vK6D8bwRvS6r31I5/deKW+8pE3GKutGJHLtTFECH40C9ZVIwLsaN0NYIiDfrBaX4kh/AW4d6YCCXQPyZns6rEM/6jXsr8wbT3LG9iV5T0eazPMW2fEoALztKpQbiv7NLA5QDumYoBTm4ZgQBkdumkTFD3GP+3oMw9tDmSLDV+UzsoWOzftaBONg1pW1BK5x6FpgnnBzFRpG0QxND9WEJ2M/Npp0YawDxa5NEwxKKQfGvYz9UtGdg9YVP/dLxS94xif7Qm/wrMc5REmvto2mxASQG0VRotuQxaz5C8qN6Cjx3sgvmo46h48Saw1vnTjUaS3aR+7NFwJvJJo+2njLfEmBbIfYEhHMf1fAnI7Qvkdko7b9BlwxXIMxsvAoAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 2 Um \\alpha^{2}$"
      ],
      "text/plain": [
       "      2\n",
       "2⋅Um⋅α "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Upp.subs(x,x0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we see now that:\n",
    "\n",
    "$$\\frac{d^2U}{dx^2}\\bigg|_{x_0} = 2 U_m \\alpha^2$$\n",
    "\n",
    "so we can connect the parameters of the Morse Potential to the Taylor Series coefficients near the minimum of the Potential.\n",
    "\n",
    "Project 3\n",
    "===========\n",
    "\n",
    "See the podcast! But basically the project is to use the experimental data, the analytical approximation (Taylor Series) to the Morse Potential and the Heun Method to:\n",
    "\n",
    "* Verify/show that the experimental data is consistent with the Morse Potential model.\n",
    "\n",
    "* Use the Heun Method to compute the period of small oscillations about the equilibrium position and compare that to the experimental data *and* the Taylor Series result.\n",
    "\n",
    "$U_m$ = 7.37 eV\n",
    "\n",
    "$x_0$ = 1.2 A\n",
    "\n",
    "$\\alpha = 2.287\\ A^{-1}$\n",
    "\n",
    "$f$ = $5.19\\times 10^{13}$ Hz\n",
    "\n",
    "IMPORTANT NOTES\n",
    "---------------\n",
    "\n",
    "1) Without careful attention to keeping the namespaces separate it's not really practical to use sympy and pylab at the same time in the same notebook. If you're interested in this, I have a notebook that illustrates how this can be done, but I don't want overload you with too many new concepts in this one lesson. If you're already confused just take my advice and don't mix them. Doing this is easy: Create *one* notebook for sympy. Take the *result* of that notebook and use it to produce the python expression you need for the numerial algorithm. Then create a *second* notebook and simply *paste* python expression produced by the first notebook in the evaluation of the force for the second notebook. Simple.\n",
    "\n",
    "2) This project is meant to add two tools to your toolbox 1: sympy, 2: Talor Series expansion. Both of these are really just analytical tools that work with algebra and calculus. We're going to use these tools to understand and approximate the low amplitude behavior system interacting under the Morse Potential. Then we're going to *check* our understanding and approximation by calculating the (arbitrary amplitude) behavior of a system interacting under the Morse Potential numericaly and evaluating the frequency at various amplitudes. So:\n",
    "\n",
    "    a) Use sympy/Taylor series to estimate the low-amplitude frequency of the system\n",
    "    b) Use sympy to get a python expression for the Force of interaction\n",
    "    c) Use a numerical method (e.g., Heun Algorithm) to compute the time evolution (frequency) of the system\n",
    "    d) Compare the results of (c) and (a)\n",
    "    e) Once you've validated the code in step \"d\" experiment with increasing the amplitude of motion. \n",
    "    f) Describe how the period and trajectory of the particle changes.\n",
    "\n",
    "Please ask if you have questions!\n"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}