012-1d_potential_well.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Particle in one-dimensional potential well"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Period of oscillations  in potential well "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Dynamics of a particle of mass $m$ moving in one dimension $OX$ is described the Newton equation \n",
    "\n",
    "$$m\\ddot x =m\\dot v = F(x) = -U'(x),$$ \n",
    " \n",
    "where $F(x)$ is a force acting on tha particle and $U(x)$ is potential energy of the particle. We assume that there are no dissipative forces. Therefore, as we know from the previous section, the total energy $E$ of the particle is conserved, i.e.,  \n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:cons_E1d}\n",
    "\\frac{mv^2}{2} + U(x) = E = const.\n",
    "\\end{equation}\n",
    "\n",
    "If we treat the velocity $v=\\frac{dx}{dt}$ as an  independent variable, it is implicit equation of the orbit (the phase curve) in **phase space**  $(x,v)$ of the  particle with  energy $E$. This equation can be rewritten in the form \n",
    "\n",
    "$$\\frac{m}{2}\\left(\\frac{dx}{dt}\\right)^2+U(x)=E$$\n",
    "\n",
    "It is a first-order differential equation and can be solved by the method of variable separation. The result reads \n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:t_cons_E1d}\n",
    " t=\\sqrt{\\frac{m}{2}} \\; \\int_{a}^{b}{\\frac{dx}{\\sqrt{E-U(x)}}}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Here, $a$ is the initial position $x(0)=a$ and $b$ is the final position $x(t)=b$ of the particle. The time $t$ is time for moving the particle from the point $a$ to the point $b$ under the condition that $E \\ge U(x)$ for all $x \\in (a, b)$. \n",
    "\n",
    "This is a useful formula which allows to calculate period of oscillations of the particle in a potential well.  \n",
    "\n",
    "In this section we will consider a class of potentials in the form  \n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:Uxn}\n",
    "U(x)=A |x|^n,\n",
    "\\end{equation}\n",
    "\n",
    "where $n$ is a positive real number.\n",
    "\n",
    "These potentials are similar: they are bounded from below,   have only one minimum at $x=0$ and tends to infinity when $x\\to \\pm \\infty$. In such potential the particle motion is bounded and the particle oscillates between two symmetrical positions $x_0$ and $-x_0$ which in general depend on the total energy $E$ and are determined by the equation  \n",
    "\n",
    "$$U(\\pm x_0) = E$$\n",
    "\n",
    "Because of symmetry, the period $T$ of oscillations can be determined by the equation \n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:T}\n",
    " T=4 \\; \\sqrt{\\frac{m}{2}} \\; \\int_{0}^{x_0}{\\frac{dx}{\\sqrt{E-U(x)}}}\n",
    "\\end{equation}\n",
    "\n",
    "This class of potentials is simple but nevertheless analysis of $T$ in dependence of the index $n$ and the total energy $E$ is very interesting and instructive. \n",
    "\n",
    "We will use computer algebra and numerical methods to investigate properties of motion is such potential wells.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load('cas_utils.sage')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = var('t')\n",
    "m = var('m')\n",
    "A = var('A')\n",
    "assume(A > 0)\n",
    "assume(m > 0)\n",
    "y = function('y')(t)\n",
    "de = m*diff(y,t,2) + 2*A*y == 0\n",
    "showmath( desolve(de, y,ivar=t) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is an analytical solution of the Newton equation in the case when $n=2$ (a harmonic potential). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Particle in potential $x^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "For $n=2$ the system is  a harmonic oscillator:\n",
    "\n",
    "$$U(x)=A x^2.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reset()\n",
    "var('m A x E')\n",
    "forget()\n",
    "assume(A > 0)\n",
    "assume(E > 0)\n",
    "assume(E,'real')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "To obtain the integration limit $x_0$  in the formula for the period of oscillations, we must solve the equation:\n",
    "\n",
    "$$U(x)=E$$\n",
    "\n",
    "So for the $Ax^2$ potential, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U(A,x) = A*x^2\n",
    "xextr = solve (U(A=A,x=x)==E,x)\n",
    "showmath(xextr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "These formulas describe the values of the oscillator's extremal positions for a given energy. Let's put them into the formula for $T$:\n",
    ":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period = 2*sqrt(m/2)*integrate( 1/sqrt(E-U(A,x)),(x,x.subs(xextr[0]),x.subs(xextr[1])))\n",
    "period = period.canonicalize_radical()\n",
    "showmath(period)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "We see that the period $T$ does not depend on energy of the oscillator. It means that it does not depend on the initial conditions because the total energy of the particle depends on them. In turn, it means that it does not depend on the distance  between the points $-x_0$ and $x_0$. It seems to be unusual behavior: time to travel  from $-1$ to $1$ and back is the same as time to travel from $-10000$ to $10000$ and back. In the second case the distance is much,  much longer but time is the same. This is an exceptional property valid only for the harmonic potential! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Particle in $|x|^n$ potential"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "If $n\\neq2$, the general formula for the period can be written as:\n",
    "\n",
    "$$T=4 \\sqrt{\\frac{m}{2}} \\; \\int_0^{x_0}\\frac{dx}{\\sqrt{E- A x^n}}$$\n",
    "\n",
    "or in the equivalent form: \n",
    "\n",
    "$$T=4 \\sqrt{\\frac{m}{2}}\\frac{1}{\\sqrt{E}}\\int_0^{x_0}\\frac{dx}{\\sqrt{1-Ax^n/E}}$$\n",
    "\n",
    "This  integral can be transformed  to a dimensionless form  by substitution \n",
    "\n",
    "$$\\frac{A}{E}x^n=y^n.$$\n",
    "\n",
    "It  is in fact a  linear relationship between $x$ and $y$:\n",
    "\n",
    "$$\\left(\\frac{A}{E}\\right)^{\\frac{1}{n}}x=y.$$\n",
    "\n",
    "Therefore, we can change the integration variable to $y$.  To do this,  we use SAGE to transform the expression under integral  in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('dx dy A E x y')\n",
    "var('n',domain='integer')\n",
    "assume(n>=0)\n",
    "assume(A>0)\n",
    "assume(E>0)\n",
    "ex1 = dx/sqrt(1-A/E*x^n)\n",
    "showmath(ex1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "and we substitute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex2 = ex1.subs({x:(E/A)^(1/n)*y,dx:dy*(E/A)^(1/n)})\n",
    "showmath( ex2.canonicalize_radical().full_simplify() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Let's take out the expression that depends on the parameters $A$ and $E$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr2 = (ex2/dy*sqrt(-y^n + 1)).full_simplify()\n",
    "showmath( expr2.canonicalize_radical() )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prefactor = expr2*sqrt(m/2)*4*1/sqrt(E)\n",
    "showmath( prefactor.canonicalize_radical() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Finally, we obtain:\n",
    "\n",
    "$$T=4 \\sqrt{\\frac{m}{2}}\\frac{1}{A^{1/n}} E^{\\frac{1}{n}-\\frac{1}{2}}\\int_0^{y_0}\\frac{dy}{\\sqrt{1- y^n}}$$\n",
    "\n",
    "\n",
    "For $n=2$, dependence on $E$ disappears, as we already have seen in the previous case.\n",
    "\n",
    "We still need to calculate the upper limit  $y_0$ of integration.  In the integral, the upper limit is the position in which the total energy is the potential energy:\n",
    "\n",
    "$$U(x)=E$$\n",
    "\n",
    "In this case $$Ax^n=E.$$\n",
    "\n",
    "By changing the variables we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solve( (A*x^n == E).subs({x:(E/A)^(1/n)*y}), y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "That is, the integration limit is $y_0=1.$\n",
    "\n",
    "Therefore the period of oscillations is given by the relation:\n",
    "\n",
    "$$T=4 \\sqrt{\\frac{m}{2}}\\frac{1}{A^{1/n}} E^{\\frac{1}{n}-\\frac{1}{2}}\\int_0^{1}\\frac{dy}{\\sqrt{1- y^n}}$$\n",
    "\n",
    "We note that only for the case $n=2$, the period $T$ does not depend on E (i.e. on initial conditions, i.e. on the distance). In other cases it depends on the total energy $E$, i.e. it depends on initial conditions, i.e. it depends on the distance  between the points $-x_0$ and $x_0$. \n",
    "\n",
    "The above equation shows how much time the particle needs to travel the distance for one oscillation in dependence on $E$ and in consequence on the distance:   If energy $E$ is higher then the distance $4x_0$ is longer. \n",
    "\n",
    "The scaled integral can be expressed by the beta function of Euler http://en.wikipedia.org/wiki/Beta_function.\n",
    "We can calculate it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('a')\n",
    "assume(a,'integer')\n",
    "assume(a>0)\n",
    "print( assumptions() )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate(1/sqrt(1-x^(a)),(x,0,1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "We get a formula containing the beta function. It can be evaluated numerically for any values ​​of the $a$ parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(beta(1/2,1/a)/a).subs({a:2}).n()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Let's examine this formula numerically. You can use the $\\beta$ function, or numerically estimate the integral. This second approach allows you to explore any potential, not just $U(x)=ax^n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def beta2(a,b):\n",
    "    return gamma(a)*gamma(b)/gamma(a+b)\n",
    "\n",
    "a_list = srange(0.02,5,0.1)\n",
    "a_list2 = [1/4,1/3,1/2,1,2,3,4,5]\n",
    "\n",
    "integr_list = [ integral_numerical( 1/sqrt(1-x^a_) ,0,1, algorithm='qng',rule=2)[0] \\\n",
    "               for a_ in a_list ]\n",
    "integr_list_analytical = [ beta2(1/2, 1/a_)/a_ for a_ in a_list2 ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we obtain some analytically simple formulas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath( integr_list_analytical )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not  we can compare those analytical numbers with numerical results, for example on the plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt_num = list_plot(zip( a_list,integr_list), plotjoined=True )\n",
    "plt_anal = list_plot(zip( a_list2,integr_list_analytical),color='red')\n",
    "(plt_num + plt_anal).show(ymin=0,figsize=(6,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Having an analytical solution, you can examine the asymptotics for large $n$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('x')\n",
    "asympt = limit( beta2(1/2, 1/x)/x,x=oo )\n",
    "asympt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt_asympt = plot(asympt,(x,0,5),linestyle='dashed',color='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Let's add a few points for which the integral takes exact values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = zip(a_list2[:5],integr_list_analytical[:5])\n",
    "showmath(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_point_labels(l):\n",
    "    p=[]\n",
    "    for x,y in l:\n",
    "        p.append(  text( \"$(\"+latex(x)+\", \"+latex(y)+\")$\" ,(x+0.1,y+0.2) , fontsize=14,horizontal_alignment='left',color='gray') )\n",
    "        p.append( point ( (x,y),size=75,color='red' ) )\n",
    "    return sum(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "some_points = plot_point_labels(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt_all = plt_num+plt_anal+plt_asympt+some_points\n",
    "plt_all.show(figsize=(6,3),ymin=0,ymax=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Numerical convergence\n",
    "\n",
    "The integral  \n",
    "$$\\int_0^1  \\frac{dx}{\\sqrt{1-x^n}}$$ seems to be divergent for $n:$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath( numerical_integral( 1/sqrt(1-x^(0.25)) , 0, 1)  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "However, choosing the right algorithm gives the correct result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_ = 1/4. # exponent in integral\n",
    "integral_numerical( 1/sqrt(1-abs(x)^a_) , 0, 1, algorithm='qags')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "lets check it out with an exact formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(beta(1/2,1/a)/a).subs({a:a_}).n()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, we  see that carefull numerical integration gives finite result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## The dependence of  period on  energy for different $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('E x n')\n",
    "def draw_E(n,figsize=(6,2.5)): \n",
    "    p = []\n",
    "    p2 = []\n",
    "    p.append( plot(abs(x)^n,(x,-2,2),\\\n",
    "                   ymin=0,ymax=4,legend_label=r\"$U(x)=|x|^{%s}$\" % n )   )\n",
    "    p.append( plot( (x)^2,(x,-2,2),\\\n",
    "                   color='gray',legend_label=r\"$U(x)=x^{2}$\",\\\n",
    "                   axes_labels=[r\"$x$\",r\"$U(x)$\"]  ))\n",
    "    \n",
    "    p2.append(  plot(  4/sqrt(2)*(beta(1/2, 1/n)/n)* E^(1/n-1/2),\\\n",
    "                     (E,0.00,3),ymin=0,ymax=7,axes_labels=[r\"$E$\",r\"$T$\"]  ) )\n",
    "    p2.append(  plot(  4/sqrt(2)*(beta(1/2, 1/2)/2),\\\n",
    "                     (E,0.00,3) ,color='gray' ) )\n",
    "    \n",
    "    show( sum(p), figsize=figsize )\n",
    "    show( sum(p2), figsize=figsize )\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbtest"
    ]
   },
   "outputs": [],
   "source": [
    "@interact\n",
    "def _(n=slider([1/4,1/2,2/3,3/4,1,3/2,2,3])):\n",
    "    draw_E(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "if 'PDF' in os.environ.keys():\n",
    "    draw_E(1/2)\n",
    "    draw_E(1)\n",
    "    draw_E(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the dependence of period $T$ on energy (i.e. amplitude) $T(E)$ for different values of $n$. In figure belowe we see that if $n>2$ then oscillations are faster as energy grows. On the other hand if $n<1$, oscillations are getting slower with growing energy.\n",
    "\n",
    "Another interesting observation is that for potentials with $n>>1$ and $n<<1$ oscillations will become arbitrarily slow and fast, respectively when $E\\to0$. For $n>1$ the potential well is *shallow* at the bottom and *steep* far away from the minimum and for $n<1$ the opposite is true."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_s  = [1/3,2/3, 1, 2, 3] \n",
    "plot( [4/sqrt(2)*(beta(1/2, 1/n_)/n_)* E^(1/n_-1/2) \\\n",
    "                 for n_ in  n_s],\\\n",
    "    (E,0.00,2),axes_labels=[r\"$E$\",r\"$T$\"],\\\n",
    "     legend_label=['$n='+str(latex(n_))+'$' for n_ in n_s],\\\n",
    "     ymax=7, figsize=(6,3) )\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Numerical integration of equations of motion \n",
    "\n",
    "Here we will investigate numerically period  $T$  and compare with the analytical formula.\n",
    "First, let's have a closer look how the potential and force behave for different index $n$:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_Uf(n_):\n",
    "    U(x) = abs(x)^(n_)\n",
    "    plt = plot( [U(x),-diff( U(x),x)],(x,-1,1),\\\n",
    "               detect_poles='show',ymin=-3,ymax=3,\n",
    "              legend_label=[r\"$U(x)$\",'$f(x)$'])\n",
    "    plt.axes_labels([r\"$x$\",r\"$U(x)=|x|^%s$\"%latex(n_)])\n",
    "    show(plt,figsize=(6,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbtest"
    ]
   },
   "outputs": [],
   "source": [
    "@interact\n",
    "def _(n_=slider(0.1,2,0.1)):\n",
    "    plot_Uf(n_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'PDF' in os.environ.keys():\n",
    "    plot_Uf(1)\n",
    "    plot_Uf(2)\n",
    "    plot_Uf(1/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that for $n \\ge 1$ the force and potential are continuous. If $n=1$ then force has finite jump (discontinuity). Both those cases should not be a problem for numerical integration. \n",
    "\n",
    "However for $n<1$ we have infinite force at $x=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problems with numerical integration\n",
    "\n",
    "There is a possibility that  if the ODE integrator comes too close, it will blow out!\n",
    "\n",
    "We can fix this problem by softening the  potential singularity by adding small number: \n",
    "$$ |x| \\to |x| + \\epsilon. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('x',domain='real')\n",
    "var('v t')\n",
    "eps = 1e-6\n",
    "U(x) = (abs(x)+eps)^(1/2)\n",
    "showmath( U.diff(x).expand().simplify() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "to make sure that Sage will not leave $x/|x|$ unsimplified we can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w0 = SR.wild(0)\n",
    "w1 = SR.wild(1)\n",
    "f = -U.diff(x).subs({w0*w1/abs(w1):w0*sign(w1)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath( f(x) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode_pot = [v,f(x)]\n",
    "\n",
    "t_lst = srange(0,10,0.01)\n",
    "sol = desolve_odeint(ode_pot,[1,.0],t_lst,[x,v])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = line(zip(t_lst, sol[:,0])) + line(zip(t_lst, sol[:,1]), color='red')\n",
    "p.axes_labels(['$t$','$x(t),v(t)$'])\n",
    "p + plot(1,(x,0,10),linestyle='dashed',color='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can evaluate the period $T$ from the trajectory obtained via numerical solutions. For this purpose one might need an interpolation of numerical table returned by `desolve_odeint`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def find_period(x,t):\n",
    "    zero_list=[]\n",
    "    x0 = x[0]\n",
    "    for i in range(1,len(x)):\n",
    "        if x[i]*x[i-1] < 0:\n",
    "            zero_list.append( - (t[i-1]*x[i] - t[i]*x[i-1])/(x[i-1] - x[i]) )\n",
    "    lnp = np.array(zero_list)\n",
    "    return 2*(  (lnp-np.roll(lnp,1))[1:] ).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('x1 x2 t1 t2 a b ')\n",
    "showmath(  (-b/a).subs(  solve([a*t1+b==x1,a*t2+b==x2],[a,b], solution_dict=True)[0] ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "We find numerically a period of trajectory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = find_period( sol[:,0],t_lst)\n",
    "T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbtest"
    ]
   },
   "outputs": [],
   "source": [
    "assert abs(T-7.54250)<1e-4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Exact results  for comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for n=2 2*pi/sqrt(2)==(2*pi/sqrt(2)).n()\n",
    "table( [[\"n\",\"T\"]]+[ [n_,((4/sqrt(2)*(beta(1/2, 1/n_)/n_)* E^(1/n_-1/2)).subs({E:1})).n()] \n",
    "                    for n_ in [1/4,1/3,1/2,2/3,1,2,3,4,5] ] )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Using the formula for the period to reproduce the trajectory of movement\n",
    "\n",
    "We take $m=1$ and $A=1$ (then $x=E$), then we can reproduce the trajectory reversing the formula for $T(E)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('x')\n",
    "U(A,x) = A*x^2\n",
    "A = 1/2\n",
    "E = 1\n",
    "m = 1.\n",
    "x1=0.1\n",
    "showmath( solve(E-U(A,x), x) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_lst  = [ (sqrt(m/2.)*integrate( 1/sqrt(E-U(A,x)),(x,0,x1)).n(),x1)  \\\n",
    "          for x1 in srange(0,sqrt(2.)+1e-10,1e-2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point(t_lst ,color='red')+\\\n",
    " plot(sqrt(2)*sin(x),(x,0,pi),figsize=(6,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, if we known the dependence of $T(E)$ then we can calculate exactly the potential! \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\newpage"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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
}