{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Particle in  multistable potentials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "\n",
    "## Phase portrait for a one-dimensional system\n",
    "\n",
    "The Newton equation  for a point particle in a one-dimensional potential $U(x)$ can be written as a set of two first-order differential equations:\n",
    "\n",
    "$$ \\begin{cases} \\dot x = v \\\\ \\dot v = - \\frac{1}{m}U'(x) \\end{cases} $$\n",
    "\n",
    "We can draw a phase portrait, i.e. parametric solutions $(x(t),v(t))$ and a vector field defined by the right hand sides of equations in the $(x,v)$-phase space.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Example: motion in the $U(x) = x^3-x^2$ potential\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('v')\n",
    "m = 1\n",
    "U(x) = x^3-x^2\n",
    "xmax,xmin = sorted([s_.rhs() for s_ in solve(U.diff(x)==0,x)])\n",
    "Emin = U(xmin)\n",
    "Etot = 1/2*m*v^2 + U(x)\n",
    "\n",
    "plot(U(x),(x,-0.4,1.1),figsize=4) +\\\n",
    " point([xmin,U(xmin)],color='red',size=40)+\\\n",
    " point([xmax,U(xmax)],color='red',size=40)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The potential $U(x) =x^3-x^2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pkt = point((xmin,0),size=20,color='black')\n",
    "plt = sum([ implicit_plot(Etot==E0,(x,-1/2,1.2),(v,-1,1),color='blue')\\\n",
    "          for E0 in srange(Emin,0.0,0.02)])\n",
    "plt +=implicit_plot(Etot==0,(x,-1/2,1.2),(v,-1,1),color='black') +pkt\n",
    "plt +=sum([ implicit_plot(Etot==E0,(x,-1/2,1.2),(v,-1,1),color='green')\\\n",
    "           for E0 in srange(0.02,-2*Emin,0.04)])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The phase curves $(x, v)$ depends on initial conditions $(x_0, v_0)$. There are two types of phase curves: closed (periodic motion of the  particle in a bounded interval) and open (the motion is unbounded: the particle can escape to -infinity or can return from -infinity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vector_field = vector([v,-U.diff(x)])\n",
    "plt  + plot_vector_field(vector_field.normalized(),(x,-1/2,1.2),(v,-1,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vector field shows direction of motion of the particle on the $(x, v)$-plane. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "\n",
    "## Harmonic  oscillations  limit  for  one-dimensional systems\n",
    "\n",
    "Consider a conservative one-dimensional system. In this case, the force $f(x)$ can always be represented as gradient of the potential $U(x)$, namely, \n",
    "$$f(x) = -\\frac{\\partial U(x)}{dx}.$$\n",
    "Consider a certain potential that has a minium at some point $x_0$. The necessary condition for minimum of the function is its zero first derivative  at this point.  Let's expand the potential in the Taylor series around the minimum. We obtain:\n",
    "$$ U(x) = U(x_0) + \\underbrace{U'(x_0)( x-x_0)}_{=0}+\\frac{1}{2} U''(x_0)(x-x_0)^2+...$$\n",
    "For small deviation from the minimum this series can be approximated by the function \n",
    "$$U(x) = \\frac{1}{2} k (x-x_0)^2,$$\n",
    "\n",
    "The Newton equation for such  motion is as follows:\n",
    "\n",
    "$$m \\ddot x = m  a  = F = -U'(x)  =  -k (x-x_0)$$\n",
    "\n",
    "For the new variable $y=x-x_0$ it takes the form \n",
    "\n",
    "$$m \\ddot y =-ky$$\n",
    "\n",
    "This is the already known equation for the harmonic oscillator with the shifted equilibrium point $x_0$. \n",
    "\n",
    "Now,  let us return to the system with the potential $U(x) = x^3-x^2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('x v')\n",
    "Etot = 1/2*v^2 + U(x)\n",
    "Elin = 1/2*v^2 + U(xmin)+1/2*U.diff(x,2).subs(x==xmin)*(x-xmin)^2\n",
    "show(Etot)\n",
    "show(Elin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Emin = Etot(x=xmin,v=0)\n",
    "Emin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Let's have a look at  the trajectories for the  exact system with $U(x) = x^3-x^2$ and the linearized system with $U(x)=(1/2) x^2$. The blue line below is a separatrix - i.e. a solution with $E=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt = sum([ implicit_plot(Etot==E0,(x,.4,.91),(v,-.3,.3),color='red') \\\n",
    "          for E0 in srange(Emin+1e-3,-0.1,0.005)])\n",
    "plt += implicit_plot(Etot==0.00,(x,0,1.1), (v,-.6,.6),color='blue') \n",
    "\n",
    "plt_lin =sum([ implicit_plot(Elin==E0+1e-3,(x,.4,.91),(v,-.3,.3),color='gray') \\\n",
    "              for E0 in srange(Emin,-0.1,0.005)])\n",
    "plt+plt_lin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "For larger ones, there is a growing discrepancy:\n",
    "\n",
    " - for the nonlinear system, above certain energy, there are open trajectories\n",
    " - motion in an linearized system is always an ellipse. The period does not depend on the amplitude."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Period of oscillation around potential minimum\n",
    "\n",
    "We take a particle with  energy by $dE$ larger than minimum of the potential:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dE = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E0 = U(xmin)+dE\n",
    "E0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to analyze eqaution of motion, we need to solve: $$U(x)=E_0$$ and obtain the extreme values of position to the left and right from potentiall minimum. In our case it con be done analytically. \n",
    "\n",
    "However, note that it requires to neglect imaginary part, which is approximated small value of the order of [machine epsilon](https://en.wikipedia.org/wiki/Machine_epsilon)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, x1,x2 = sorted( [s_.rhs().n().real() for s_ in solve(U(x)==E0,x)])\n",
    "x1,x2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, having $x_{1,2}$ we can integrate numerically equation for period (\\ref{eq:1d_TE}) and obtain:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period = 2*sqrt(m/2.)*\\\n",
    "  integral_numerical(1/sqrt(E0-U(x)) , x1,x2, algorithm='qags')[0]\n",
    "period"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let is plot this situation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U(x) = x^3-x^2\n",
    "xmax,xmin = sorted([s_.rhs() for s_ in solve(U.diff(x)==0,x)])\n",
    "plot(U(x),(x,-0.4,1.1),figsize=4,gridlines=[None,[U(xmin),E0]])+\\\n",
    " point([xmin,U(xmin)],color='red',size=40)+\\\n",
    " point([xmax,U(xmax)],color='red',size=40)+\\\n",
    " point([x1,U(x1)],color='green',size=40)+\\\n",
    " point([x2,U(x2)],color='green',size=40)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write a function which computes the period based on previous steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def T(E0):\n",
    "    m = 1\n",
    "    _, x1,x2 = sorted( [s_.rhs().n().real() \\\n",
    "                        for s_ in solve(U(x)==E0,x)])\n",
    "    integral, error = \\\n",
    "      integral_numerical(1/sqrt(E0-U(x)), x1,x2, algorithm='qags')\n",
    "    period = 2*sqrt(m/2.) * integral\n",
    "    return period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period_num = T(U(xmin)+dE)\n",
    "period_num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We known exactly the period of oscillation for small neiborhood of potential minimum. It is simply period of the harmonic oscillator with frequecy $\\omega=\\sqrt{U''(x_{min})}$ (for $m=1$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = sqrt(U(x).diff(x,2).subs(x==xmin.n()))\n",
    "period_harm = 2*pi.n()/omega\n",
    "period_harm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How does period depent on energy?\n",
    "\n",
    "We can easily compute plot $T(E_0)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TonE = [(E_,T(E_)) for E_ in srange(U(xmin)+1e-6,-1e-5,0.001)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line(TonE, figsize=(6,2),gridlines=[None,[period_harm]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How does the tracjectory look like? \n",
    "\n",
    "Let us integrate numerically equation of motion:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_lst = srange(0, period_num, 0.01, include_endpoint=True) \n",
    "sol = desolve_odeint([v,-U.diff(x)],[x2,.0],t_lst,[x,v])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have in phase space $(x,\\dot x)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line(sol[::10,:],marker='o',figsize=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Position oscillates with time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line(zip(t_lst,sol[:,0]),figsize=(6,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Time to reach the hill\n",
    "\n",
    "The top of the potential hill  is at the beginning of the coordinate system $(x,E)$. Then we examine the limit $E\\to0$ boundary\n",
    "\n",
    "Near zero, we can approximate the potential by  the reverse parabola. Then the time to reach the hill from a certain point (for example $x=1$) reads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('E')\n",
    "assume(E>0)\n",
    "integrate(-1/sqrt(E+x^2),x,1,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "This result is divergent for $E\\to0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "limit( arcsinh(1/x),x=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It means that time needed to climb a hill with *just* enough kinetic energy is **infinite**. It is valid only for potential hills which have zero derivative at the top. On the other hand for potential barriers which do not have this property, for example, $U(x)=-|x|$, the particle can reach the top with just enought energy in  finite time.\n",
    "\n",
    "Let analyze it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U1(x) = -abs(x)\n",
    "E0 = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we can plot velocity and potential:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot([U1(x), sqrt(2*(E0 - U1(x)))],(x,-1,1),figsize=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the time of travel from $x=-1$ to $x=0$ is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:1d_TE}\n",
    "t=\\sqrt{\\frac{m}{2}} \\; \\int_{-x1}^{x1}{\\frac{dx}{(\\sqrt(E-U(x)}}\n",
    "\\end{equation}\n",
    "\n",
    "which in this case is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sqrt(m/2.)*integrate(1/sqrt((E0- U1(x))),x,-1,0).n()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of potentials which behave like $|x|^\\alpha$, for $\\alpha>1$ we can calculate time  of travel if we particle total energy is by $dE$ larger than potential barrier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def t_hill(E0):\n",
    "    m = 1\n",
    "\n",
    "    x2, = [s_.rhs().n().real() for s_ in solve(U(x)==E0,x) if s_.rhs().n().imag().abs()<1e-6]\n",
    "    \n",
    "    integral, error = \\\n",
    "      integral_numerical(1/sqrt(E0-U(x)), 0,x2, algorithm='qags')\n",
    "    m = 1\n",
    "    period = 2*sqrt(m/2.) * integral\n",
    "    return period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_hill(9.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "t_E = [(E_,t_hill(E_)) for E_ in np.logspace(-6,1,120)]\n",
    "line(t_E,axes_labels=['$E$','$t$'], figsize=(6,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_E = [(log(E_),t_hill(E_)) for E_ in np.logspace(-6,16,120)]\n",
    "line(t_E,axes_labels=['$\\log_{10} E$','$t$'], figsize=(6,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "#### Experiment with Sage!\n",
    "Examine in a similar way the system corresponding to the movement in the $U(x) = -\\cos(x)$ potential - this is a physical pendulum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\newpage"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.6",
   "language": "",
   "name": "sagemath"
  },
  "language": "python",
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "nbTranslate": {
   "displayLangs": [
    "en",
    "pl"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "pl",
   "targetLang": "en",
   "useGoogleTranslate": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}