S2_TD4_python_notebook_en_student.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "F1WI4b8ShUv1"
   },
   "source": [
    "# <span style=\"color:blue\"> <center>Teacher / TD4 : 4TPU279U $-$ Bachelor 1st year $-$ spring 2023</center></span>\n",
    "# <center>Introduction to python programming</center>\n",
    "# <hr style=\"border:1px solid black\"><center>  Finite differences and integration of movement equations  </center><hr style=\"border:1px solid black\">\n",
    "</br>\n",
    "\n",
    "<div style=\"text-align: right\"> Credits: L. Truflandier and Simon Villain-Guillot </div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following topics will be cover in this notebook:\n",
    "- Derivatives and finite differences\n",
    "- Finite differences and Euler method\n",
    "- Application to dynamical equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <hr style=\"border:1px solid black\">  Dérivées et différences finies  <hr style=\"border:1px solid black\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forward and backward finite differences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dérivative $f'(x)$ for a function $f : \\mathbb{R}\\rightarrow\\mathbb{R}$ at $x$ is formally given by:\n",
    "\n",
    "$$f'(x) = \\lim_{h\\rightarrow 0}\\frac{f(x+h)-f(x)}{h}$$\n",
    "or\n",
    "$$f'(x) = \\lim_{h\\rightarrow 0}\\frac{f(x)-f(x-h)}{h}$$\n",
    "\n",
    "From a geometrical point of view, the derivative corresponds to the slope of the tangent at $f$ in position $x$.\n",
    "\n",
    "> Example: from the analytical derivative, represent the tangents to $f(x)=3x^2-2x^3$ at $x_1=0$ and $x_2=1$ with $x\\in[-1.0,2.0]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. We can start by defining $f(x)$ and $f'(x)$ as Python functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 3*x**2 - 2*x**3\n",
    "\n",
    "def fp(x):\n",
    "    return 6*x - 6*x**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Discretize along the $(Ox)$ axis with 128 grid points. Then calculate $y=f(x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linspace\n",
    "\n",
    "x   = linspace(-1.0,2.0,128)\n",
    "y   = f(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. We define the abscissa points $x_1=0$ and $x_2=1$. Then we calculate the derivatives and equations of the tangents at these points, ie.\n",
    "\n",
    "$$y'_1=f'(x_1)(x-x_1) + f(x_1)$$ $$y'_2=f'(x_2)(x-x_2) + f(x_2)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xp1 = 0.0\n",
    "xp2 = 1.0\n",
    "\n",
    "yp1 = fp( xp1 ) * (x - xp1) + f( xp1 )\n",
    "yp2 = fp( xp2 ) * (x - xp2) + f( xp2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. On représente les courbes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "\n",
    "plt.plot(x,y,   linestyle='-',color='black',label='$f(x)=x^2$')\n",
    "\n",
    "plt.plot(x,yp1,linestyle='--', color='red',  label='tengente at $x_1=0$')\n",
    "plt.plot(x,yp2,linestyle='--', color='green',label='tengente at $x_2=1$')\n",
    "\n",
    "plt.plot(xp1,f(xp1),marker='o',color='red')\n",
    "plt.plot(xp2,f(xp2),marker='o',color='green')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlim(-1.0,+2.0)\n",
    "plt.ylim(-0.5,+1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the analytical derivative of $f$ is unknown or the values of $f(x_i)$ are only known for a discrete set of points $\\{x_i\\}^N_{i=1}$ (e.g. the results of a series of experiments), then we need to evaluate \n",
    "\n",
    "$$f'(x) = \\lim_{h\\rightarrow 0}\\frac{f(x+h)-f(x)}{h}$$ \n",
    "\n",
    "by numerical approximation. By fixing a value for $h$, knowing $f(x+h)$ and $f(x)$ then $f'(x)$ at the neighborhood of $x+h$ is evaluated by finite *progressive* difference, as follows: \n",
    "\n",
    "$$f'(x) \\simeq \\frac{f(x+h)-f(x)}{h}\\quad\\quad\\textcolor{green}{(1)}$$ \n",
    "\n",
    "or by *regressive* finite difference, as follows: \n",
    "\n",
    "$$f'(x) \\simeq \\frac{f(x)-f(x-h)}{h}\\quad\\quad\\textcolor{green}{(2)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 4.0 <hr style=\"border:1px solid grey\">\n",
    "\n",
    "For the function $f(x)=3x^2-2x^3$, with a step $h=0.5$, represent the derivatives by progressive finite differences at the point $x_1=0$. Use the previous example of graphical representation for this exercise. Compare with analytical results. Repeat for $h=0.1$ and $h=0.01$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Taylor development and error on *progressive* and *regressive* finite differences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\require{color}$\n",
    "\n",
    "By recalling that a taylor development (series) applied to a function $f$ which is $n$ times differentiable in $a$ is given by:\n",
    "\n",
    "$$f(x) = \\sum_{k=0}^n\\frac{(x-a)^k}{k!}f^{(k)}(a)$$\n",
    "\n",
    "Using a 3rd-order Taylor expansion along $h$ for $x$ around $x+h$ gives us[<sup>1</sup>](#fn1):\n",
    "\n",
    "$$f(x+h) = f(x) + hf'(x) + \\frac{h^2}{2!}f''(x) + \\frac{h^3}{3!}f'''(x) + h^3\\epsilon(h)\\quad\\quad\\textcolor{green}{(3)}$$\n",
    "\n",
    "From this relationship, we can deduce the error from the *progressive* finite difference approximation: \n",
    "\n",
    "$$f'(x) = \\frac{f(x+h)-f(x)}{h} - \\left( \\textcolor{blue}{\\frac{h}{2!}f''(x)} + \\frac{h^2}{3!}f'''(x) + h^2\\epsilon(h) \\right)$$\n",
    "\n",
    "In the above expression, the dominant term is $\\frac{h}{2!}f''(x) $ in which case the error is linear with respect to $h$ and we obtain:\n",
    "\n",
    "$$f'(x) = \\frac{f(x+h)-f(x)}{h} + O(h)$$\n",
    "\n",
    "In the same way, Talor's development can be used to approximate *regressive* finite differences:\n",
    "\n",
    "$$f(x-h) = f(x) - hf'(x) + \\frac{h^2}{2!}f''(x) - \\frac{h^3}{3!}f'''(x) - h^3\\epsilon(h)\\quad\\quad\\textcolor{green}{(4)}$$\n",
    "\n",
    "It follows:\n",
    "\n",
    "$$f'(x) = \\frac{f(x)-f(x-h)}{h} - \\left(\\textcolor{blue}{- \\frac{h}{2!}f''(x)} + \\frac{h^2}{3!}f'''(x) + h^2\\epsilon(h)\\right)$$\n",
    "\n",
    "The error on $f(x-h)$ is also linear with respect to $h$.\n",
    "\n",
    "<span id=\"fn1\"> $^1$ $h^3\\epsilon(h)$ represents the truncation error of the Taylor expansion with  $\\lim_{h\\rightarrow 0}\\epsilon(h)= 0$. We use the abbreviation $O(h^p)$ for $h^p\\epsilon(h)$ .</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 4.1 <hr style=\"border:1px solid grey\">\n",
    "\n",
    "For $f(x)=3x^2-2x^3$ calculate the analytical value of the derivative at $x_1 = 0$ and $x_2=1$. Display the result. Calculate the value of the numerical derivative using the *progressive* finite difference formula with $h=0.01$. Finally, calculate the order correction $h$ and refine the numerical result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 4.2 <hr style=\"border:1px solid grey\">\n",
    "\n",
    "Repeat ***Exercise 4.1*** using the *regressive* finite difference method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Central finite difference approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By subtracting the *progressive* Taylor expansion of the ***equation*** $\\textcolor{green}{(3)}$ from the *regressive* Taylor expansion of the ***equation*** $\\textcolor{green}{(4)}$, we obtain a new definition for the finite difference, called the *centered* finite difference approximation:\n",
    "\n",
    "$$f'(x) = \\frac{f(x+h)-f(x-h)}{2h} - \\left(\\textcolor{blue}{\\frac{h^2}{3!}f'''(x)} + h^2\\epsilon(h)\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see that the error is second-order with respect to $h$. For a given value of $h$, this new approximation is expected to be more accurate than the *progressive* and *regressive* versions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 4.3 <hr style=\"border:1px solid grey\">\n",
    "\n",
    "For $f(x)=3x^2-2x^3$ calculate the value of the numerical derivative at $x_1 = 0$ and $x_2=1$ using the finite difference formula *centered* with $h=0.01$. Then calculate the order correction $h^2$ and refine the numerical result.[<sup>1</sup>](#fn1)\n",
    "\n",
    "<span id=\"fn1\"> $^1$ *You should find that the numerical result is identical to the analytical one. Prove this to yourself using a pen and paper. Indeed, this is a general result: any Taylor development of order $n$ applied to a polynomial of same order are identical in every respect.*</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <hr style=\"border:1px solid black\">  Finite differences and Euler's method <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the finite-difference methods for calculating derivatives described above, we can also develop a method for solving differential equations. Here, we will only consider first-order equations of the type:  \n",
    "\n",
    "$$\\frac{dy}{dt}=f(t)\\quad\\quad\\textcolor{green}{(5)}$$ \n",
    "\n",
    ">From the ***equation*** $\\textcolor{green}{(1)}$, neglecting orders higher than $h$ we can approximate $y'(t)$ by:$$\\frac{dy}{dt}=f(t)\\simeq\\frac{y(t+\\delta t) - y(t)}{\\delta t}$$\n",
    "where $\\delta t\\equiv h$ corresponds to the time step. We can then proceed to a step-by-step integration of the following ***equation*** $\\textcolor{green}{(5)}$:\n",
    "Knowing the initial conditions $y_0=y(t_0)$ and the value of $f(t_i)$ for each value $t_i\\in[t_0,t_\\infty]$ we can predict the value of $y$ at step $t_{i+1}=t_i+h$. The ***equation*** $\\textcolor{green}{(5)}$ can be rewritten as follows:\n",
    "$$y(t_{i+1}) = y(t_i) + hf(t_i)\\quad\\quad\\textcolor{green}{(6)}$$ \n",
    "\n",
    "This integration method is called the explicit Euler method.\n",
    "\n",
    "> Note that there's another way of finding the above recurrence relation, by explaining the integration in successive time steps. In other words, from the ***equation*** $\\textcolor{green}{(5)}$, we obtain: $$dy=f(t)dt\\quad\\quad\\textcolor{green}{(7)}$$ The numerical integration of the ***equation*** $\\textcolor{green}{(7)}$ over the interval $(t_{i+1}-t_i)$ is written: $$\\int^{t_{i+1}}_{t_{i}}dy = \\int^{t_{i+1}}_{t_{i}}f(t)dt$$\n",
    "where the right-hand side of the equation can be approximated by the area of the rectangle\n",
    "with sides $h=t_{i+1} - t_{i}$ and $f(t_i)$.\n",
    "This gives: $$y(t_{i+1}) - y(t_{i}) \\simeq hf(t_i)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shown below are rectangles of area $hf(t_i)$ for $f=3t^2-2t^3$ and $t\\in[-1,2]$ with a step size of $h=0.2$. Note the error associated with the rectangles method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linspace, arange\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "tmin = -1\n",
    "tmax = +2\n",
    "h    = +0.2\n",
    "\n",
    "t = arange(tmin, tmax+h, h)\n",
    "f = 3*t**2 - 2*t**3\n",
    "plt.plot(t,f,marker='o',linestyle='-',color='blue')\n",
    "plt.fill_between(t, f, alpha=0.2, interpolate=True)\n",
    "\n",
    "for i in range(len(t)-1):\n",
    "    t_rect = [t[i], t[i], t[i+1], t[i+1], t[i]]\n",
    "    f_rect = [0   , f[i], f[i]  , 0     , 0   ]\n",
    "    plt.plot(t_rect,f_rect,marker='',linestyle='-',color='red')\n",
    "    \n",
    "plt.xlabel('$t_i$')\n",
    "plt.ylabel('$f(t_i)$')\n",
    "plt.xticks(arange(tmin, tmax+h, h*2))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 4.4 <hr style=\"border:1px solid grey\">\n",
    "\n",
    "Using the previous program, calculate the area of the function $f(x)=3x^2-2x^3$ for $h=0.2$ and $t\\in[-1,2]$. Compare with the theoretical value of $3/2$, displaying the error. Repeat for $h=0.1,0.01,0.001,0.0001,0.00001$ (you can use a loop)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example of application of Euler's method\n",
    "\n",
    "Take, for example, the following equation:\n",
    "$$\\frac{dy}{dt}=-y$$\n",
    "with initial condition $y_0=100$ and analytical solution $y(t)=y_0\\exp(-t)$.\n",
    "> Let's discretize $t$ for $t\\in[0,5]$ with a step $h=0.1$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import arange, linspace, exp, zeros\n",
    "\n",
    "h = 0.1\n",
    "t = arange(0,5+h,h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Let's define the right-hand side of $\\textcolor{green}{(5)}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return -x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Let's calculate step by step $y(t_{i+1})=y(t_i) + hf(t_i)$ with $y(t_0)=100$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = 100\n",
    "yi =  y0\n",
    "for i in range(1,len(t)):\n",
    "    yi = yi + h*f(yi)\n",
    "    #print('index %s : ti = %.3f, yi = %.3f'%(i,t[i],yi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Backup each value of $y(t_i)$ in an array `y_euler`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = 100\n",
    "y_euler    = zeros(len(t))\n",
    "y_euler[0] = y0\n",
    "for i in range(len(t)-1):\n",
    "    y_euler[i+1] = y_euler[i] + h*f(y_euler[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Let's represent $y(t)$ obtained by Euler's method $(y_\\textrm{Euler})$ and by the analytical solution $(y_\\textrm{exact})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "\n",
    "y_exact= y0*exp(-t)\n",
    "\n",
    "plt.plot(t,y_euler, linestyle='-', label='Euler')\n",
    "plt.plot(t,y_exact, linestyle='--',label='Analytic')\n",
    "plt.xlabel('$t$ (in sec.)')\n",
    "plt.ylabel('$y(t)$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Show the deviation of the numerical solution from the analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "\n",
    "y_err = y_exact - y_euler\n",
    "\n",
    "plt.plot(t,y_err, linestyle='-', label='Euler')\n",
    "plt.xlabel('$t$ (in sec.)')\n",
    "plt.ylabel('$\\Delta y(t)$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 4.5 <hr style=\"border:1px solid grey\">\n",
    "\n",
    "Take the numerical solution of the equation $$\\frac{dy}{dt}=-y,\\quad y_0=100$$ for time steps $h=\\{0.1,0.05,0.01\\}$ and plot the error against the analytical results. What do we observe? Is it expected?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 4.6 <hr style=\"border:1px solid grey\">\n",
    "\n",
    "From the *centered* finite difference approximation the numerical solution of the differential equation $$\\frac{dy}{dt}=-y,\\quad y_0=100$$ is written:\n",
    "$$y(t+h) = y(t-h) + 2hf(t)$$\n",
    "that is,\n",
    "$$y(t_{i+1}) = y(t_{i-1}) + 2hf(t_i)$$\n",
    "1. Write the program to solve this problem. Take $h=0.01$, $y_0 =100$ and $y_1 =99$.\n",
    "2. Plot the analytical and numerical results on the same graph.\n",
    "3. On a second graph, plot the error of finite differences *centered* and finite differences *progressive*.\n",
    "4. Repeat with the following initial conditions: $y_0 =100$ and $y_1 =99.005$. Conclude."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:red\">**solution :**</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import arange, linspace, exp, zeros\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return -x\n",
    "\n",
    "y0 = 100\n",
    "\n",
    "h  = 0.01\n",
    "t  = arange(0,5+h,h)\n",
    "y_euler_c    = zeros(len(t))\n",
    "y_euler_c[0] = y0\n",
    "y_euler_c[1] = 99\n",
    "for i in range(1,len(t)-1):\n",
    "    y_euler_c[i+1] = y_euler_c[i-1] + 2*h*f(y_euler_c[i])\n",
    "    \n",
    "y_exact = y0*exp(-t)    \n",
    "y_err_c = y_exact - y_euler_c\n",
    "\n",
    "plt.plot(t,y_euler_c, linestyle='-', label='Euler central FD')\n",
    "plt.plot(t,y_exact, linestyle='--',  label='Analytic')\n",
    "plt.xlabel('$t$ (in sec.)')\n",
    "plt.ylabel('$y(t)$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <hr style=\"border:1px solid black\">  Application to the equations of dynamics <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The aim here is to construct the trajectories of material points subjected to various forces (spring, weight, friction, gravity) numerically, i.e. without analytically solving the differential equations of dynamics:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\vec f &=m \\vec a =m \\frac {d \\vec v}{dt}\\quad\\quad\\textcolor{green}{(8)}\\\\\n",
    "\\vec v &= \\frac {d \\vec r}{dt}\n",
    "\\end{align}$$\n",
    "\n",
    "with $\\vec{v}$ the velocity vector with coordinates $(v_x,v_y,v_z)$ and $\\vec{r}$ the position vector with coordinates $(x,y,z)$. Using Euler's method, cf. ***equation*** $\\textcolor{green}{(6)}$, the calculation of $\\vec{v}$ and $\\vec{r}$ can be obtained by the following recurrence:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\vec{r}(t_{i+1}) &= \\vec{r}(t_i) + h\\ \\vec{v}(t_i)\\\\\n",
    "\\vec{v}(t_{i+1}) &= \\vec{v}(t_i) + \\frac{h}{m}\\vec{f}(t_i)\n",
    "\\end{align}$$\n",
    "\n",
    "with $h=t_{i+1}-t_{i}$ the integration step of the equations of motion. In the case of unidirectional motion, we have $\\vec{r}=x\\vec{e}_x$, $\\vec{v}=v\\vec{e}_x$ and $\\vec{f}=f\\vec{e}_x$.\n",
    "This gives us one equation for position:\n",
    "\n",
    "$$\\begin{align}\n",
    "x(t_{i+1})   &= x(t_i) + h\\;v(t_i)\\\\   \n",
    "\\end{align}$$\n",
    "\n",
    "and one equation for speed :\n",
    "\n",
    "$$\\begin{align}\n",
    "v(t_{i+1}) &= v(t_i)+\\frac{h}m f\\left(x(t_i)\\right) \\\\ \n",
    "\\end{align}$$\n",
    "\n",
    "From these equations, we can see that if we know $(x,v)$ at a given instant $t_i$, we can then determine $(x,v)$ at a later instant $t_{i+1}$ and thus build step by step the trajectory of a material point over time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1D harmonic oscillator motion and Euler method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> We wish to study the motion of an object of mass $m$ subjected to a restoring force (Hooke's law): $$f=-kx \\quad\\quad\\textcolor{green}{(9)}$$ with $m=1$ kg, $k=4$ N/m. Friction will be neglected. The study time will be $t_\\textrm{max}=4\\pi$ s, with a discretization step of $h=0.01$ s. Initial conditions are $x(t_0)=0.1$ and $v(t_0)=0$.  We recall the analytical solution of the equation of motion: $x(t) = x_0\\cos(\\omega t)$ with $\\omega=\\sqrt{k/m}$.\n",
    "\n",
    "> We begin by defining the right-hand side of the equation of dynamics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    k = 4.0\n",
    "    return -k*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The time axis is discretized:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import pi, arange, zeros, cos, sqrt\n",
    "\n",
    "tmin = 0.0\n",
    "tmax = 4*pi\n",
    "h    = 0.01\n",
    "t    = arange(tmin,tmax,h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> We initialize $m$, $k$, the arrays for $x$ and $v$ as well as the initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 1.0\n",
    "k = 4.0\n",
    "\n",
    "x = zeros(len(t))\n",
    "v = zeros(len(t))\n",
    "\n",
    "x0 = 0.1\n",
    "v0 = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> We apply Euler's recurrence relation to $x$ and $v$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[0] = x0\n",
    "v[0] = v0\n",
    "for i in range(len(t)-1):\n",
    "    x[i+1] = x[i] + h*v[i]\n",
    "    v[i+1] = v[i] + h*f(x[i])/m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Plot the time laws $x(t)$ and $v(t)$ and the analytical solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "\n",
    "x_analy = x0*cos( sqrt(k/m)*t )\n",
    "\n",
    "plt.plot(t,x,marker='',label='$x(t)$')\n",
    "plt.plot(t,v,marker='',label='$v(t)$')\n",
    "plt.plot(t,x_analy,marker='',linestyle='--',color='black',label='$x(t)$ analytic',)\n",
    "plt.xlabel('$t$ (in s)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Draw the phase portrait $(x(t),v(t))$ and locate the initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x,v,marker='')\n",
    "plt.plot(x0,v0,marker='o')\n",
    "plt.xlabel('$x(t)$')\n",
    "plt.ylabel('$v(t)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the ideal, frictionless case, the mechanical energy $(E_m)$ must be a constant of motion. This energy is given by the sum of kinetic energy:\n",
    "$$T(t)= \\frac{1}{2}mv(t)^2$$\n",
    "and the potential energy, which in the case of a 1D harmonic oscillator is given by:\n",
    "$$V(t)= \\frac{1}{2}kx(t)^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Plot the evolution of $E_m(t)$, $T(t)$ and $V(t)$ with\n",
    "$$T(t)= \\frac{1}{2}mv(t)^2$$\n",
    "$$V(t)= \\frac{1}{2}kx(t)^2$$\n",
    "$$E_m(t)=T(t)+V(t)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 1/2*m*v**2\n",
    "V = 1/2*k*x**2\n",
    "Em = T + V\n",
    "\n",
    "plt.plot(t,T,marker='',label='$T(t)$')\n",
    "plt.plot(t,V,marker='',label='$V(t)$')\n",
    "plt.plot(t,Em,marker='',label='$E_m(t)$')\n",
    "plt.xlabel('$t$ (in sec)')\n",
    "plt.ylabel('Energy (in J)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> From the ***equations*** $\\textcolor{green}{(8)}$ and $\\textcolor{green}{(9)}$ we can show that: $$E_m(t_{i+1}) = E_m(t_{i})\\left(1+h^2\\frac{k}{m}\\right)\\quad\\quad\\textcolor{green}{(10)}$$\n",
    "Note that the total energy is conserved if and only if $h=0$. In all other cases, energy increases at a rate of $(1+\\alpha h^2)$. Below, the evolution of mechanical energy as a function of time is represented using ***equation*** $\\textcolor{green}{(10)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Em_analy = zeros(len(t))\n",
    "\n",
    "Em_analy[0] = T[0] + V[0]\n",
    "for i in range(len(t)-1):\n",
    "    Em_analy[i+1] = Em_analy[i]*(1 + h**2*k/m) \n",
    "\n",
    "plt.plot(t,T,marker='',label='$T(t)$')\n",
    "plt.plot(t,V,marker='',label='$V(t)$')\n",
    "plt.plot(t,Em,marker='',label='$E_m(t)$')\n",
    "plt.plot(t,Em_analy,marker='',linestyle='--',label='$E_m(t)$ analytique')\n",
    "\n",
    "plt.xlabel('$t$ (in sec)')\n",
    "plt.ylabel('Energy (in J)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 4.7 <hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have seen that Euler's method of integrating the equations of motion does not conserve energy. In the 1D case, these equations are given by:\n",
    "\n",
    "$$\\begin{align}\n",
    "x(t_{i+1}) &= \\textcolor{blue}{x(t_i)} + h\\;v(t_i)\\\\   \n",
    "v(t_{i+1}) &= v(t_i)+\\frac{h}m \\textcolor{blue}{f\\left(x(t_{i})\\right)} \\\\ \n",
    "\\end{align}$$\n",
    "\n",
    "Here we propose to study the Hooke-Newton method described below:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\textcolor{blue}{x(t_{i+1})} &= x(t_i) + h\\;v(t_i)\\\\   \n",
    "v(t_{i+1}) &= v(t_i)+\\frac{h}m \\textcolor{blue}{f\\left(x(t_{i+1})\\right)} \\\\ \n",
    "\\end{align}$$\n",
    "\n",
    "for the latter, note that to predict velocity at time $t_{i+1}$ the force is now evaluated at position $x(t_{i+1})$.\n",
    "\n",
    "Resume the study of the motion of a 1D harmonic oscillator using these new equations. Use the same calculation parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iW7MFdnUhUwd"
   },
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Diaporama",
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}