02-ODEs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "# System of autonomous differential equations\n",
    "\n",
    "\n",
    "In the previous part  we presented two  examples of systems described by  ordinary differential equations. The first class of systems is described by classical mechanics and its Newton equations. Another class of systems is the phenomenological equation describing growth processes, chemical kinetics processes, population dynamics in biological systems, etc. These two classes of systems are particular cases of a set  of autonomous ordinary differential equations in  the form \n",
    "\n",
    "$$\\frac{dx_1}{dt} = F_1(x_1, x_2, x_3, ..., x_n)$$\n",
    "\n",
    "$$\\frac{dx_2}{dt} =   F_2(x_1, x_2, x_3, ..., x_n)$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "$$\\frac{dx_n}{dt} = F_n(x_1,   x_2, x_3, ..., x_n) $$\n",
    "\n",
    "This set of equations can be written in shorter vector notation  as \n",
    "\n",
    "$$\\frac{d\\vec x}{dt} = \\vec F(\\vec x), \\quad \\quad \\vec x(0)  = \\vec x_0  \\quad \\quad \\quad \\quad (1) $$\n",
    "\n",
    "where vectors $\\vec x$ and $\\vec F$ read \n",
    "\n",
    "$$\\vec x = [x_1, x_2, x_3, ...., x_n], \\quad \\quad \\vec F = [F_1, F_2, F_3, ..., F_n]$$\n",
    "\n",
    "and a set of initial conditions is: \n",
    "\n",
    "$$\\vec x(0) = [x_1(0), x_2(0), x_3(0), ... , x_n(0)] = \\vec x_0 = [x_1^{(0)}, x_2^{(0)}, x_3^{(0)}, ... ,  x_n^{(0)}]  $$\n",
    "\n",
    "The subscript $n$ tell us how many differential equations are \"hidden\" in the above vector notation. In other words, we consider the set of $n$ differential equations characterized by $n$ scalar functions $F_i(x_1, x_2, x_3, ..., x_n), (i=1,2,3, ..., n)$. Note that we are considering $F_i$ functions that do not explicitly depend on time. In this case, we say that it is the autonomous system of $n$ ordinary differential equations. In addition, the vector $\\vec F$  can be regarded as a vector field associated with the system of differential equations or a vector field generated by the system of differential equations. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "pl"
   },
   "source": [
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Existence and uniqueness of solutions\n",
    "\n",
    "In order to describe phenomena using ordinary differential equations with corresponding initial conditions at e.g. $t=0$, we need to know solutions for times $t>0$ (time evolutions). \n",
    " \n",
    "**EXAMPLE 1:**\n",
    "\n",
    "The differential equation \n",
    "\n",
    "$$\\frac{dx}{dt}=-2x, \\quad \\quad \\quad x(0) = x_0$$\n",
    "\n",
    "is a linear differential equation with all possible values $x\\in (-\\infty, \\infty)$. It is one of the simplest and frequently exploited equation. One can check that its solution is given by the function \n",
    "\n",
    "$$x(t) = x_0   e^{-2t}$$\n",
    "\n",
    "and fulfills the initial condition $x (0) = x_0$. This exponential function is defined for all times \n",
    "$t \\in (- \\infty, \\infty)$.  It is the only solution and unique solution.  Below, we depict solutions for selected initial conditions.   Let us note that all solutions tend to the same final value $x = 0$  when $t \\to \\infty$. It is called a stationary state (stationary solution, steady-state, ...). \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('t')\n",
    "g(t,a) = a*exp(-2*t)\n",
    "p1=plot(g(t,1),(t,0,2),figsize=(6, 3), legend_label=\"x(0)=1\", color='blue' )\n",
    "p2=plot(g(t,0),(t,0,2),figsize=(6, 3), legend_label=\"x(0)=0\", color='red' )\n",
    "p3=plot(g(t,-1),(t,0,2),figsize=(6, 3), legend_label=\"x(0)=-1\", color='green' )\n",
    "show(p1+p2+p3, axes_labels=[r'$t$',r'$x(t)$'], frame=True, axes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "**EXAMPLE 2:**\n",
    "\n",
    "Equation\n",
    "$$\\frac{dx}{dt}= 3 x^2, \\quad \\quad \\quad x(0) = x_0$$\n",
    "is a nonlinear differential equation. The right side of this equation is determined for all values $x \\in (-\\infty, \\infty)$. It can be  solved by the variable separation method. We obtain the function \n",
    "\n",
    "$$x(t) = \\frac{x_0}{1-3 x_0 t} $$\n",
    "\n",
    "which is the solution and obeys the initial condition $x(0) = x_0$. This function is not defined for all finite values of time  $t \\in  (-\\infty,   \\infty)$. There are restrictions on its  value. But this is the only solution and unique solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('t')\n",
    "#detect_poles - wykrywanie i rysowanie biegunów\n",
    "g=plot(-4.0/(1 +12*t),t,0,5,detect_poles='show',legend_label=r'$x(0)= -4$', color='blue')\n",
    "g+=plot(lambda t: 0.0,t,0,5,legend_label=r'$x(0)=0$',color='red')\n",
    "g+=plot(1.0/(1-3*t),t,0,0.33,detect_poles='show', legend_label=r'$x(0)=1$',color='green')\n",
    "g.show(axes_labels=[r'$t$',r'$x$'],tick_formatter='latex',xmin=0,xmax=0.5,ymin=-4.1,ymax=8, figsize=(7,4), frame=True, axes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "The solutions with  the negative initial condition $x(0) < 0$ (blue curve) are determined for all $t > 0$ . The same applies to the solution for the initial condition $x(0)=0$ (red curve). However, the solution with a positive initial condition $x(0) > 0$ diverges at finite time $t$ (green curve). \n",
    " \n",
    "**EXAMPLE  3:**\n",
    "\n",
    "The differential equation \n",
    "\n",
    "$$\\frac{dx}{dt}=  2 \\sqrt x, \\quad \\quad \\quad x(0) = x_0 \\ge 0$$\n",
    "\n",
    "is a nonlinear differential equation.  The right hand side of this equations is determined for the non-negative values of $x$. One can solve this equation by the method of variable separation.  The  result is:  \n",
    "\n",
    "$$x(t) = (t +  \\sqrt x_0)^2 $$\n",
    "\n",
    "This function is defined for al  $t> 0$.   If the initial condition $x(0) > 0$ then it is the only solution. \n",
    "If, however, the initial condition is  $x (0) = 0$ then there are 2 solutions! One of the possible solution is the function \n",
    "\n",
    "$$x(t) = x_1(t)  = t^2, \\quad \\quad \\quad x(0) = 0$$\n",
    "\n",
    "The second solution  reads: \n",
    "\n",
    "$$x(t) =x_2(t)= 0, \\quad \\quad \\quad x(0)=0$$ \n",
    "\n",
    "In this case, there is not a unique solution. If this differential equation describes a real process, it would be \n",
    "impossible to predict the future evolution of the system! \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('t')\n",
    "p1=plot(t*t,(t,0,1),figsize=(6, 3), legend_label=\"x(0)=1\", color='blue' )\n",
    "p2=plot(0,(t,0,1),figsize=(6, 3), legend_label=\"x(0)=0\", color='red' )\n",
    "show(p1+p2, frame=True, axes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Why in the last example  there are two solutions, not only one solution? It is  behaviour of the function $\\sqrt{x}$ which is in the right hand side of this equation. For the uniqness the function should satisfied the so-called Lipschitz condition. \n",
    "\n",
    "\n",
    "We say that the function $f(x)$ satisfies the Lipschitz condition on the open set $U$ if there is such a constant $L > 0$ that\n",
    "\n",
    "$$|f(x) -f(y)| \\le L|x- y|$$\n",
    "\n",
    "for all $x, y \\in U$.\n",
    "The Lipschitz condition can be written as\n",
    "\n",
    "$$|f(x+h) -f(x)| \\le L h \\quad \\quad \\mbox{or as } \\quad \\quad \\frac{f(x+h) - f(x)}{h}| \\le L$$\n",
    "\n",
    "It follows that if $f(x)$ has a finite value of derivative (the derivative at a given point is bounded), then it satisfies  the Lipschitz condition. There are of course non-fifferentialbe functions that satisfy the Lipschitz condition.\n",
    "\n",
    "**Picard's theorem:**\n",
    "If the function $f(x)$ is continuous in $U$ and satisfies the Lipschtz condition in $U$ then the differential equation\n",
    "\n",
    "$$\\frac{dx}{dt} = f(x), \\quad \\quad \\quad x(0) = x_0$$\n",
    "\n",
    "has exactly one solution in $U$.\n",
    "There are several versions of this theorem, but for us this simplest version is enough.\n",
    "\n",
    "In the case of a system of differential equations, the Lipschitz condition has the form\n",
    "\n",
    "$$|F_i(x_1, x_2, x_3, ..., x_n) - F_i(y_1, y_2, y_3, ..., y_n)| \\le L  \\sum_{k=1}^n|x_k-y_k|$$\n",
    "\n",
    "This inequality has to be fulfilled for all functions $F_i$ and Picard's theorem is formulated in a similar way.  The Lipschitz condition is satisfied if partial derivatives are baunded,\n",
    "\n",
    "$$\\lvert\\frac{\\partial F_i}{\\partial x_k}\\rvert \\le K$$\n",
    "\n",
    "for positive $K$.\n",
    "\n",
    "Now we can answer why in the third example the solution is non-uniques: the function $f(x) = 2\\sqrt x$  does not satisfy the Lipschitz condition at the point $x=0$  because the derivative\n",
    "\n",
    "$$\\frac{df(x)}{dx} = \\frac{1}{\\sqrt x}$$\n",
    "\n",
    "at this point $x=0$ is divergent. At points $x>0$ the derivative is bounded  and Picard's theorem is satisfied. Therefore, solutions are unique for all initial conditions  $x(0) > 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "pl"
   },
   "source": [
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "\n",
    "## PHASE SPACE\n",
    "\n",
    "\n",
    "Once again, we rewrite $(1)$ differential equations in notation:\n",
    "\n",
    "$$\\frac{dx_1}{dt} = F_1(x_1, x_2, x_3, ..., x_n)$$ \n",
    "\n",
    "$$\\frac{dx_2}{dt} =  F_2(x_1, x_2, x_3, ..., x_n)$$ \n",
    "\n",
    "$$  .....................................$$ \n",
    "\n",
    "$$\\frac{dx_n}{dt} = F_n(x_1,  x_2, x_3, ..., x_n) $$\n",
    "\n",
    "The above system of differential equations defines a certain dynamic system (the mathematical definition of a dynamic system can be very abstract, but for our needs, what we have written is enough). The set of all possible $\\{x_1, x_2, x_3, ..., x_n\\}$ values ​​creates a set which we call the phase space of the $(2)$ system. The dimension of this space is $n$, which is the same as there are differential equations.\n",
    "Depending on the context, we will use different notations for the same equations.\n",
    "Examples:\n",
    "\n",
    "1. One differential equation. Usually we will use the notation\n",
    "\n",
    "$$\\frac{dx}{dt} = \\dot x = f(x)$$\n",
    "\n",
    "The phase space is 1-dimensional.\n",
    "\n",
    "2. Two differential equations. Usually we will use the notation \n",
    "\n",
    "$$\\frac{dx}{dt} = \\dot x = f(x, y)$$\n",
    "\n",
    "$$\\frac{dy}{dt} = \\dot y= g(x, y)$$\n",
    "\n",
    "The phase space is 2-dimensional.\n",
    "\n",
    "3. Three differential equations. Usually we will use the notation \n",
    "\n",
    "$$\\frac{dx}{dt} = \\dot x = f(x, y, z)$$ \n",
    "\n",
    "$$\\frac{dy}{dt} = \\dot y= g(x, y, z)$$ \n",
    "\n",
    "$$\\frac{dz}{dt} = \\dot z= h(x, y, z)$$\n",
    "\n",
    "The phase space is 3-dimensional.\n",
    "\n",
    "4. The Newton's equation for a particle moving only along the $OX$ axis and  which is driven by the force $F(x)$  dependent only on the particle coordinate: \n",
    "\n",
    "$$m \\ddot x= F(x)$$\n",
    "\n",
    "where $m$ is the mass of the particle. This is a 2-nd order differential equation and it is equivalent to a system of 2 first-order differential equations:\n",
    "$$\\dot x = v$$ \n",
    "\n",
    "$$\\dot v = \\frac{1}{m} F(x)$$\n",
    "\n",
    "The phase space is 2-dimensional: it is a set of all possible positions and velocities of the particle, $\\{x, v\\}$. Despite its simplicity, this model is extremely important. It provides a starting point for understanding many crucial  aspects of dynamical systems.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Geometric properties of phase space\n",
    "\n",
    "###  Phase curves\n",
    "\n",
    "To avoid abstract definitions at this stage, as an example we consider a 2-dimensional dynamical system \n",
    "defined by a set of differential equations: \n",
    "\n",
    "$$ \\dot x = f(x, y), \\quad \\quad x(0) = x_0 $$\n",
    "\n",
    "$$ \\dot y= g(x, y),\\quad \\quad y(0) = y_0 $$\n",
    "\n",
    "The phase space of this system is 2-dimensional. It can be a plane or part of it. But it can be a more complicated 2-dimensional set. For example, it can be a sphere (similar to the surface of a ball), it can be a torus (similar to a bicycle inner tube). These can be even more complicated 2-dimensional objects. But for our purposes it is enough to consider the plane.  We introduce a Cartesian coordinate system with the OX and OY axes on the plane. The initial condition $\\{x_0=x(0), y_0=y(0)\\}$ is the point on this plane.  We solve the above system of differential equations numerically using the definition of serivative and applying  the simplest scheme:\n",
    "\n",
    "$$\\frac{x(t+h) - x(t)}{h} = f(x(t), y(t))$$ \n",
    "\n",
    "$$\\frac{y(t+h) - y(t)}{h} = g(x(t), y(t))$$\n",
    "\n",
    "We will re-write these equation in the form\n",
    "\n",
    "$$x(t+h) = x(t) + f(x(t), y(t)) h$$ \n",
    "\n",
    "$$y(t+h) = y(t) + g(x(t), y(t)) h$$\n",
    "\n",
    "I. Numerical calculations must be started with the initial condition at the moment $t=0$, i.e. in the first step we calculate\n",
    "\n",
    "$$x_1 =x(h) = x(0) + f(x(0), y(0)) h$$ \n",
    "\n",
    "$$y_1 = y(h) = y(0) + g(x(0), y(0)) h$$\n",
    "\n",
    "On the plane we get a point with coordinates $\\{x_1, y_1\\}$. Let's mark it on the plane. We now have 2 points:\n",
    "\n",
    "$$\\{x_0, y_0\\}, \\quad \\quad \\{x_1, y_1\\}$$\n",
    "\n",
    "II. In the next step we choose $t=h$ time:\n",
    "\n",
    "$$x_2 =x(h+h) = x(2h) =  x(h) + f(x(h), y(h)) h$$ \n",
    "\n",
    "$$y_2 = y(h+h) = y(2h) =  y(h) + g(x(h), y(h)) h$$\n",
    "\n",
    "We will use the above notation: $x_1 =  x(h),  y_1 = y(h)$ and write these equations in the form\n",
    "\n",
    "$$x_2 =  x_1 + f(x_1, y_1) h$$\n",
    "\n",
    "$$y_2 =  y_1 + g(x_1, y_1) h$$\n",
    "\n",
    "On the plane we get a point with coordinates $\\{x_2, y_2\\}$. Let's mark it on the plane. We now have 3 points:\n",
    "\n",
    "$$\\{x_0, y_0\\}, \\quad \\quad \\{x_1, y_1\\},  \\quad \\quad \\{x_2, y_2\\}$$\n",
    "\n",
    "III. We can see immediately that in step 3 we get the equations\n",
    "\n",
    "$$x_3 =  x_2 + f(x_2, y_2) h$$ \n",
    "\n",
    "$$y_3 =  y_2 + g(x_2, y_2) h$$\n",
    "\n",
    "and we get the point with coordinates $\\{x_3, y_3\\}$.\n",
    "\n",
    "\n",
    "IV. We notice that in the n-th step we get equations\n",
    "\n",
    "$$x_n =  x_{n-1} + f(x_{n-1}, y_{n-1}) h$$ \n",
    "\n",
    "$$y_n =  y_{n-1} + g(x_{n-1}, y_{n-1}) h$$\n",
    "\n",
    "V. More often you write what you get in the next step, i.e. n + 1:\n",
    "\n",
    "$$x_{n+1} =  x_n + f(x_n, y_n) h$$ \n",
    "\n",
    "$$y_{n+1} =  y_n + g(x_n, y_n) h$$\n",
    "\n",
    "We obtain recurrence equations that allow us to determine the evolution of the system, i.e. the numerical solution of the system. On the $XY$ plane we get a collections of points with coordinates\n",
    "\n",
    "$$\\{x_n, y_n\\}$$\n",
    "\n",
    "If the time step  $h$ is small enough, we connect this series of points with a solid line and we get a curve on the plane. This curve is called the phase curve of the dynamical system. It is one of the most important geometric object related to the system of differential equations. Having drawn such a phase curve, we can conclude about the evolution of the system and the characteristics of the system behavior. Below we present two examples: phase curves for a harmonic oscillator and a damped harmonic oscillator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "#### A. HARMONIC OSCILLATOR\n",
    "\n",
    "\n",
    "For the harmonic oscillator, the force acting on the particle has the form \n",
    "\n",
    "$$ F = -kx$$\n",
    "\n",
    "where $k$ is characterized by \"springiness\" of the spring. The Newton equation  then reads:\n",
    "\n",
    "$$m\\ddot x = -kx, \\quad \\quad \\mbox{lub w postaci} \\quad \\quad \\ddot x = -(k/m) x = -\\omega^2 x$$\n",
    "\n",
    "where $\\omega^2 = k/m$. This second order al differentiequation is equivalent to 2 first order equations:\n",
    "\n",
    "$$\\dot x = y, \\quad \\quad x(0) = x_0$$ \n",
    "\n",
    "$$ \\dot y = -\\omega^2 x, \\quad \\quad y(0) = y_0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "#### A. DAMPED  HARMONIC OSCILLATOR\n",
    "\n",
    "The previous example of the harmonic oscillator is rather idealized bacause it is  assumed that tha particle moves in  a vacuum. In real situation, it moves in  in air, water or other liquid. Then the particle is affected by an additional force, namely the friction (damping)  force $F_d$. The friction force is proportional to the velocity of the particle and opposite to the direction of motion, i.e., \n",
    "\n",
    "$$F_d = -\\gamma_0 v = -\\gamma_0 \\dot x$$\n",
    "\n",
    "where $\\gamma_0$ is called the friction (or damping) coefficient.\n",
    "In such a case, the Newton equation takes the form\n",
    "\n",
    "$$m\\ddot  x = -kx - \\gamma_0 \\dot x, \\quad \\quad \\mbox{or in the form} \\quad \\quad  \\ddot x = -\\frac{k}{m} x - \\frac{\\gamma_0}{m} x = -\\omega^2 x - \\gamma \\dot x$$\n",
    "\n",
    "where $\\omega^2 = k/m$ and $\\gamma = \\gamma_0/m$. The above equation is equivalent to 2 first-order differential equations:\n",
    "\n",
    "$$\\dot x = y, \\quad \\quad x(0) = x_0$$ \n",
    "\n",
    "$$ \\dot y = -\\gamma y -\\omega^2 x, \\quad \\quad y(0) = y_0$$\n",
    "\n",
    "Of course, when $\\gamma = 0$, then we get the harmonic oscillator equation without friction (undamped motion).\n",
    "Below the phase curves for these 2 examples are depicted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "#### Harmonic oscillator without friction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('x y')\n",
    "\n",
    "def schemat_eulera2D(vec, ics, Tlist):\n",
    "    i = 0\n",
    "    dx = [ics[0]]\n",
    "    dy = [ics[1]]\n",
    "    h = Tlist[i+1] - Tlist[i]\n",
    "\n",
    "    iks(x,y) = vec[0R]*h\n",
    "    igrek(x,y) = vec[1R]*h\n",
    "    \n",
    "    for time in Tlist:\n",
    "        dx.append(dx[i] + iks(dx[i],dy[i]))\n",
    "        dy.append(dy[i] + igrek(dx[i],dy[i]))\n",
    "        i += 1\n",
    "        \n",
    "    return zip(dx,dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(layout={'top':[['_omega','x0','y0']],'bottom':[['T','h']]})\n",
    "def _(title=['a','b'], h=selector(['0.005','0.01','0.05','0.1','0.5','1'], default='0.1', buttons=True),x0=input_box(2,label=r'$x_0$', width=10), y0=input_box(4,label=r'$y_0$', width=10), T=input_box(0, width=10), _omega=input_box(1,label=r'$\\omega$', width=10)):\n",
    "    global oscylator_nietlumiony, background\n",
    "    \n",
    "    if T == 0:\n",
    "        T = 2*pi/_omega\n",
    "        \n",
    "    listT = srange(0,T,float(h), include_endpoint=True)\n",
    "    background = desolve_odeint(vector([y,-_omega^2*x]), [x0, y0], srange(0,T+0.1,0.1,include_endpoint=True), [x,y])\n",
    "    oscylator_nietlumiony = schemat_eulera2D([y,-_omega^2*x], [x0, y0], listT)\n",
    "\n",
    "    print r'Parametry modelu'\n",
    "    html(r'$\\omega=%s, x_0=%s, y_0=%s$'%(_omega,x0,y0))\n",
    "    print r'Parametry symulacji'\n",
    "    html(r'$h=%s, T=%s$'%(h,T))\n",
    "    print '\\nDla T=0 lista generowana jest automatycznie dla jednego okresu własnego oscylatora'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact\n",
    "def _(krok=slider(1, len(oscylator_nietlumiony), 1, default=1, label=r'krok')):\n",
    "    \n",
    "    buf = zip(*oscylator_nietlumiony)\n",
    "    minx, maxx, miny, maxy = min(buf[0]), max(buf[0]), min(buf[1]), max(buf[1])\n",
    "    kroki = oscylator_nietlumiony[:krok]\n",
    "    kroki_plot = list_plot(kroki, figsize=(4,4), axes_labels=[r'$x$',r'$y$'], size=30, xmin=minx, xmax=maxx, ymin=miny, ymax=maxy)\n",
    "    \n",
    "    txt_plot = text(r'$[x_0,y_0]$',kroki[0],vertical_alignment='bottom',horizontal_alignment='left')\n",
    "    for i in range(1,len(kroki)):\n",
    "        txt_plot += text(r'$[x_{%d},y_{%d}]$'%(i,i),kroki[i],vertical_alignment='bottom',horizontal_alignment='left')\n",
    "    \n",
    "    full_plot = list_plot(oscylator_nietlumiony, plotjoined=True, figsize=(4,4), axes_labels=[r'$x$',r'$y$'])\n",
    "    full_plot += list_plot(background.tolist(), plotjoined=True, color='grey', alpha=0.5)\n",
    "    print(\"krzywe fazowe dla oscylatora bez tarcia\")\n",
    "    (full_plot+kroki_plot).show()\n",
    "    (kroki_plot+txt_plot).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "In the case of an ideal (not damped) oscillator, the phase curves are closed and it means that the motion is periodic.  The period of such periodic motion is the time needed for the particle starting from a point e.g. $\\{x_0, y_0\\}$ and moving along a phase curve to reach the same point $\\{x_0, y_0\\}$ again.\n",
    "In the case of a damped oscillator, the phase curve is a spiral that winds up to $\\{0, 0\\}$. Spiral motion means that both $x$ and $v=y$ decrease over time and for long times the location of $x$ and the speed of $v$ go to zero, i.e. the particle slows down and finally stops. This is damped motion: the vibration amplitude decreases over time. One can observe more  complex phase curves but  the general rule is: when $x$ grows, it means the particle's position increases. When $y$ decreases, it means that the particle velocity decreases. When $x$ decreases, the position coordinate of the particle decreases. As $y$ increases, the particle velocity increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(layout={'top':[['_omega','_gamma','x0','y0']],'bottom':[['T','h']]})\n",
    "def _(title=['a','b'], h=selector(['0.05','0.01','0.1','0.5','1'], default='0.1', buttons=True),x0=input_box(2,label=r'$x_0$', width=10), y0=input_box(4,label=r'$y_0$', width=10), T=input_box(0, width=10), _omega=input_box(1,label=r'$\\omega$', width=10), _gamma=input_box(0.5,label=r'$\\gamma$', width=10)):\n",
    "    global oscylator_tlumiony, glob_gamma, glob_omega, background2\n",
    "    glob_gamma = _gamma\n",
    "    glob_omega = _omega\n",
    "    \n",
    "    if T == 0: \n",
    "        T = 2*pi/_omega\n",
    "    \n",
    "    listT = srange(0,T,float(h),include_endpoint=True)\n",
    "    background2 = desolve_odeint(vector([y,-_omega^2*x-_gamma*y]), [x0, y0], srange(0,2*pi/_omega+0.1,0.1,include_endpoint=True), [x,y])\n",
    "    oscylator_tlumiony = schemat_eulera2D([y,-_omega^2*x-_gamma*y], [x0, y0], listT)\n",
    "\n",
    "    print r'Parametry modelu'\n",
    "    html(r'$\\gamma=%s, \\omega=%s, x_0=%s, y_0=%s$'%(_gamma,_omega,x0,y0))\n",
    "    print r'Parametry symulacji'\n",
    "    html(r'$h=%s, T=%s$'%(h,T))\n",
    "    \n",
    "    print '\\nDla T=0 lista generowana jest automatycznie dla jednego okresu własnego oscylatora'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact\n",
    "def _(krok=slider(1, len(oscylator_tlumiony), 1, default=1, label=r'krok')):\n",
    "    \n",
    "    buf = zip(*oscylator_tlumiony)\n",
    "    minx, maxx, miny, maxy = min(buf[0]), max(buf[0]), min(buf[1]), max(buf[1])\n",
    "    kroki = oscylator_tlumiony[:krok]\n",
    "    kroki_plot = list_plot(kroki, figsize=(4,4), axes_labels=[r'$x$',r'$y$'], size=30, xmin=minx, xmax=maxx, ymin=miny, ymax=maxy)\n",
    "    \n",
    "    txt_plot = text(r'$[x_0,y_0]$',kroki[0],vertical_alignment='bottom',horizontal_alignment='left')\n",
    "    for i in range(1,len(kroki)):\n",
    "        txt_plot += text(r'$[x_{%d},y_{%d}]$'%(i,i),kroki[i],vertical_alignment='bottom',horizontal_alignment='left')\n",
    "    \n",
    "    full_plot = list_plot(oscylator_tlumiony, plotjoined=True, figsize=(4,4), axes_labels=[r'$x$',r'$y$'])\n",
    "    full_plot += list_plot(background2.tolist(), plotjoined=True, color='grey', alpha=0.5)\n",
    "    \n",
    "    print(\"krzywe fazowe dla oscylatora tłumionego\")\n",
    "    (full_plot+kroki_plot).show()\n",
    "    (kroki_plot+txt_plot).show()\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "### VECTOR FIELD\n",
    "\n",
    "Right sides of the system of differential equations \n",
    "\n",
    "$$ \\dot x = f(x, y), \\quad \\quad x(0) = x_0 $$ \n",
    "\n",
    "$$ \\dot y= g(x, y),\\quad \\quad y(0) = y_0 $$ \n",
    "\n",
    "can be treated as components of a certain vector field: \n",
    "\n",
    "$$\\vec F = [F_x, F_y] = [f(x, y), g(x, y)]$$ \n",
    "\n",
    "At each point  $\\{x, y\\}$ of the coordinate plane we draw a vector with components $[f(x, y), g(x, y)]$. \n",
    " In this way we get a vector field. Well, but what information  can be obtained  from this vector field. Let's do the following numerical exercise: We start from the initial condition $\\{x_0, y_0\\}$ and draw a vector with $[f(x_0, y_0), g(x_0, y_0)]$ components at this point:  \n",
    "\n",
    "$$\\mbox{at the point }  \\quad \\{x_0, y_0\\}   \\quad \\mbox{we draw a vector of components} \\quad [f(x_0, y_0), g(x_0, y_0)]$$ \n",
    "\n",
    "As before, we numerically solve the system of differential equations and calculate $\\{x_1, y_1\\}$: \n",
    "\n",
    "$$\\mbox{at the point }  \\quad \\{x_1, y_1\\}   \\quad \\mbox{we draw a vector of components} \\quad [f(x_1, y_1), g(x_1, y_1)]$$ \n",
    "\n",
    "Then we calculate $\\{x_2, y_2\\}$: \n",
    "\n",
    "$$\\mbox{at the point }  \\quad \\{x_2, y_2\\}   \\quad \\mbox{we draw a vector of components} \\quad [f(x_2, y_2), g(x_2, y_2)]$$ \n",
    "\n",
    "In the nth step of the iteration, we calculate $\\{x_n, y_n\\}$:\n",
    "\n",
    "$$\\mbox{at the point }  \\quad \\{x_n, y_n\\}   \\quad \\mbox{we draw a vector of components} \\quad [f(x_n, y_n), g(x_n, y_n)]$$\n",
    "\n",
    "Because all of the above points $\\{x_i, y_i\\}$ are located on the phase curve, the vector $[f(x_i, y_i), g(x_i, y_i)]$ is  attached to the phase curve. We note that these vectors are tangent to the phase curve. If $\\{x_i, y_i\\}$ has interpretations of the position of the particle, then the vector $[f(x_i, y_i), g(x_i, y_i)]$ \n",
    "can be interpreted as components of the volocity vector because $\\dot x = f(x, y)$ and $\\dot y = g(x,  y)$. We know that $\\dot x = v_x$ is the x-component, while $\\dot y = v_y$ is the y-component of the velocity. In other words, the resulting vector field is the velocity field of the fictitious particle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(layout={'top':[['_omega','_gamma','x0','y0']],'bottom':[['T','h']]})\n",
    "def _(title=['a','b'], h=selector(['0.05','0.01','0.1','0.5','1'], default='0.1', buttons=True),x0=input_box(2,label=r'$x_0$', width=10), y0=input_box(4,label=r'$y_0$', width=10), T=input_box(20, width=10), _omega=input_box(1,label=r'$\\omega$', width=10), _gamma=input_box(0.5,label=r'$\\gamma$', width=10)):\n",
    "    global oscylator_tlumiony, glob_gamma, glob_omega\n",
    "    glob_gamma = _gamma\n",
    "    glob_omega = _omega\n",
    "    \n",
    "    listT = srange(0,T,float(h))\n",
    "    oscylator_tlumiony = desolve_odeint(vector([y,-_omega^2*x-_gamma*y]), [x0, y0], listT, [x,y])   \n",
    "\n",
    "    print r'Parametry modelu'\n",
    "    html(r'$\\gamma=%s, \\omega=%s, x_0=%s, y_0=%s$'%(_gamma,_omega,x0,y0))\n",
    "    print r'Parametry symulacji'\n",
    "    html(r'$h=%s, T=%s$'%(h,T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vf = lambda u,a,b: (u[0]+u[1],u[1]-a*u[0]-b*u[1])\n",
    "\n",
    "@interact\n",
    "def _(krok=slider(1, len(oscylator_tlumiony), 1, default=1, label=r'krok')):\n",
    "    kroki = oscylator_tlumiony[:krok]\n",
    "    kroki_plot = list_plot(kroki.tolist(), figsize=(4,4), axes_labels=[r'$x$',r'$y$'], size=30, xmin=-4.5, xmax=4.5, ymin=-4.5, ymax=4.5)\n",
    "    \n",
    "    pole_wektorowe = arrow(kroki[0],vf(kroki[0],glob_omega^2,glob_gamma),color='red',xmax=vf(kroki[0],glob_omega^2,glob_gamma)[0])\n",
    "    for krok in kroki[1:]:\n",
    "        pole_wektorowe += arrow(krok,vf(krok,glob_omega^2,glob_gamma),color='red', width=.4)\n",
    "        \n",
    "    txt_plot = text(r'$[x_0,y_0]$',kroki[0],vertical_alignment='bottom',horizontal_alignment='left')\n",
    "    for i in range(1,len(kroki)):\n",
    "        txt_plot += text(r'$[x_{%d},y_{%d}]$'%(i,i),kroki[i],vertical_alignment='bottom',horizontal_alignment='left')\n",
    "        \n",
    "    shadowplot = list_plot(oscylator_tlumiony.tolist(), plotjoined=True, figsize=(4,4), axes_labels=[r'$x$',r'$y$'], alpha=0.2) \n",
    "    full_plot = list_plot(oscylator_tlumiony.tolist(), plotjoined=True, figsize=(4,4), axes_labels=[r'$x$',r'$y$']) + plot_vector_field([y,-glob_omega^2*x-glob_gamma*y], (x, -4.5, 4.5), (y, -4.5, 4.5), plot_points=20, color='lime')\n",
    "        \n",
    "    print(\"krzywe fazowe dla oscylatora tłumionego\")\n",
    "    (full_plot+kroki_plot).show()\n",
    "    (shadowplot+kroki_plot+txt_plot+pole_wektorowe).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Below you will find a cell in which we encourage everyone to experiment with different models.  Have fun..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#########\n",
    "# Model #\n",
    "#########\n",
    "\n",
    "# zmienne\n",
    "var('x y')\n",
    "\n",
    "# parametry\n",
    "# UWAGA: jeżeli Twój model będzie zależny od innych parametrów \n",
    "#        tu właśnie musisz je wszystkie wyspecyfikować\n",
    "alpha = 1\n",
    "\n",
    "# warunki początkowe\n",
    "x0 = 1\n",
    "y0 = 1\n",
    "\n",
    "# model\n",
    "dx = y\n",
    "dy = -alpha*x - y\n",
    "\n",
    "# czas (T) i krok (h) symulacji\n",
    "T = 12\n",
    "h = 0.1\n",
    "\n",
    "\n",
    "###################################################\n",
    "# Symulacje + wizualizacja                        #\n",
    "###################################################\n",
    "listT = srange(0,T,float(h), include_endpoint=True)\n",
    "numeryka = desolve_odeint(vector([dx, dy]), [x0, y0], listT, [x,y])\n",
    "\n",
    "przestrzen_fazowa = list_plot(numeryka.tolist(), plotjoined=True, figsize=(4,4), axes_labels=[r'$x$',r'$y$'])\n",
    "pole_wektorowe = plot_vector_field([dx,dy], (x, numeryka[:,0].min(),\\\n",
    "                                             numeryka[:,0].max()), \\\n",
    "                                   (y, numeryka[:,1].min(), numeryka[:,1].max()), plot_points=10, color='lime')\n",
    "\n",
    "show(przestrzen_fazowa+pole_wektorowe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\newpage"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.7",
   "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
}