{ "cells": [ { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "# Appendix: a gentle introduction to differential equations\n", "\n", "\n", "\n", "## What is the differential equation?\n", "\n", "\n", "The differential equation is the relationship between the function sought and its function\n", "derivative.\n", "\n", "In the general case, we are talking about the equation of the $n$th rank if we have\n", "relations:\n", "\n", "$$F(y^{(n)}(x),y^{(n-1)}(x),\\dots,y(x),x) = 0,$$\n", "\n", "where:\n", "\n", "- $y^{(n)}(x) = \\frac{d^n f(x)}{dx^n}.$\n", "- $y(x)$ is the function you are looking for, or a dependent variable,\n", "- $x$ is called an independent variable\n", "\n", "We may also have a situation that we have $m$ $n$ equations on $m$\n", "$y_i.$ function A special case is the $m$ system of the first equations\n", "$m$ degree of function. It turns out that it is possible from the $n$ equation of this degree,\n", "create an equivalent $n$ system of first order equations.\n", "\n", "## Example: Newton equation for one particle in one dimension\n", "\n", "\n", "The particle's motion is described by:\n", "\n", "$$m a = F$$\n", "\n", "Acceleration is the second derivative of the position over time, and the force is in\n", "generality some function of $x$ position and time. So we have:\n", "\n", "$$m \\ddot x = F(\\dot x, x, t)$$\n", "\n", "Let's introduce the new $v = \\frac{dx(t)}{dt}$ function now. Substituting in\n", "the previous equation can be written:\n", "\n", "$$\\begin{cases} \\dot x = v \\\\ \\dot v = - \\frac{F(\\dot x,x,t)}{m} x. \\end{cases}$$\n", "\n", "We see that we have received two systems from one equation of the second degree\n", "first degree equations.\n", "\n", "First order equations are often presented in a form in which after\n", "the right side of the equal sign stands for the derivative and the left for the expression\n", "depending on the function:\n", "\n", "$$\\underbrace{\\frac{dx}{dt}}_{\\text{derivative }} = \\underbrace{f(x,t)}_{\\text{Right Hand Side}}$$\n", "\n", "## Geometric interpretation of differential equations.\n", "\n", "\n", "Consider the system of two equations:\n", "\n", "$$\\begin{cases} \\dot x = f(x,y) \\\\ \\dot y = g(x,y) \\end{cases}.$$\n", "\n", "This is the so-called two-dimensional autonomous system of differential equations\n", "ordinary. Autonomy means the independence of right parties from time\n", "(ie, independent variable). An example of such a system can be traffic\n", "particles in one dimension with forces independent of time.\n", "\n", "Equations from the above system can be approximated by substituting derivatives\n", "differential quotient:\n", "\n", "$$\\begin{cases} \\frac{x(t+h)-x(t)}{h} = f(x,y) \\\\ \\frac{y(t+h)-y(t)}{h} = g(x,y) \\end{cases},$$\n", "\n", "multiplying each equation by $h$ and transferring a member from value\n", "dependent variables at the moment $t$ to the right page we get:\n", "\n", "$$\\begin{cases} x(t+h) = x(t) + h \\cdot f(x,y) \\\\ y(t+h) = y(t) +h \\cdot g(x,y). \\end{cases}$$\n", "\n", "According to the definition of a derivative, in the $h\\to\\infty$ border of the $f(x,y)$ expression\n", "can be taken at the time between $t$ and $t+h$. Let's assume for ease,\n", "that we will take a moment $t$.\n", "\n", "
\n", "\n", "Comment: this choice leads to the so-called overt algorithm if\n", "we would take a moment eg $t+h$ it would be an entangled algorithm and execution\n", "the step would be related to the solution of algebraic equations.\n", "\n", "
\n", "\n", "This system has an interesting interpretation:\n", "\n", "- firstly, notice that the pair of functions determines the vector field on\n", "    plane\n", "- secondly, these equations give us a recipe like the value of the function in\n", "    of the moment $t$ get the value from the \"next\" moment $t+h$ which can be\n", "    useful for recreating the $(x(t),y(t))$ curve.\n", "\n", "## Vector field\n", "\n", "\n", "A vector field is a function that gives each point of space\n", "assigns a certain vector size. If the space will be e.g.\n", "$\\mathbb{R}^2$ is a function that will consist of two functions\n", "scalar:\n", "\n", "$$\\vec F(x,y) = \\begin{cases}f(x,y) \\\\ g(x,y) \\end{cases}$$\n", "\n", "This vector field can be visualized by drawing arrows for a certain one\n", "number of points on the plane. An example of widespread use of such\n", "vector field is wind speed field.\n", "\n", "In Sage, we can draw a vector field with `plot _vector_ field`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "var('x y')\n", "f(x,y) = -0.4*x - 0.9*y\n", "g(x,y) = 0.9*x - 0.4*y\n", "plt_v = plot_vector_field((f(x,y),g(x,y)),(x,-1.2,1.2),(y,-1.2,1.2),figsize=6,aspect_ratio=1)\n", "plt_v" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "## Graphical solution of the system of two differential equations\n", "\n", "\n", "Using the derived approximated formulas to allow calculating\n", "solution of the system of differential equations at the moment $t+h$ knowing them in\n", "At the moment of Mathath2 we can try to sketch a solution based on the chart\n", "vector field. It is enough to move in small steps according to\n", "the local direction of the arrows.\n", "\n", "Let's try to do it with the help of the algorithm:\n", "\n", "1. we take the starting point in $t$\n", "2. we calculate the point at the moment $t+h$\n", "3. we draw the end point on the graph\n", "4. we take the final point as the starting point\n", "5. we go back to 1." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0,y0 = (1,0)\n", "h = 0.2" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "By executing this cell many times we get the next steps of the algrorytm:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x1,y1 = x0+h*f(x0,y0),y0+h*g(x0,y0) \n", "plt_v = plt_v + point((x0,y0)) + arrow2d( (x0,y0), (x0+h*f(x0,y0),y0+h*g(x0,y0) ),width=1,arrowsize=2,arrowshorten=-10,aspect_ratio=1 )\n", "x0,y0 = x1,y1 \n", "plt_v" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h = 0.2\n", "x0,y0 = (1,0)\n", "plts = [plot_vector_field((f(x,y),g(x,y)),\\\n", " (x,-1.2,1.2),(y,-1.2,1.2),\\\n", " figsize=(3,3),aspect_ratio=1)]\n", "for i in range(5):\n", " x1,y1 = x0+h*f(x0,y0),y0+h*g(x0,y0) \n", " plt_v = plt_v + point((x0,y0)) + arrow( (x0,y0), (x0+h*f(x0,y0),y0+h*g(x0,y0)) ,width=1,arrowsize=2,arrowshorten=-10 )\n", " x0,y0 = x1,y1 \n", " plts.append(plt_v)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plts[-1].show()\n", "# animate(plts).show()" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "We have the following conclusions:\n", "\n", "1. The solution of the system of 2 first order equations is the curve w\n", "    $\\mathbb{R}^2.$ space\n", "2. The curve depends on the selection of the starting point.\n", "3. Two solutions coming from different starting points may be\n", "    go down to one point, but **can not intersect!**\n", "4. Because we have an unlimited selection of starting points and\n", "    there is (3) the solution of the system of two equations is two-parameter\n", "    family of flat curves.\n", "\n", "The differential equation (or system of equations) with the initial condition is called\n", "in the mathematics of Cauchy. Point (3) is known as\n", "Piccard's theorem on the existence and uniqueness of solutions to the problem\n", "Cauchy and it's worth noting that it imposes some restrictions on\n", "variability of the right sides of the system of equations.\n", "\n", "## Analytical solutions of differential equations\n", "\n", "Differential equations can be analyzed using the graphical method a\n", "numeric values ​​can be obtained with any accuracy using\n", "approximate method. These methods do not limit the form in any way\n", "right sides of the layout.\n", "\n", "Is it possible to obtain an analytical formula for the family of functions being\n", "the solution of the differential equation?\n", "\n", "This is difficult in the general case, however there are several forms of equations\n", "differentials in which we can always find an analytical solution.\n", "One of such cases is one separable equation of the first\n", "degree. Separability means that the right side is the product of the function\n", "$x$ and $t$:\n", "\n", "$$\\frac{dx}{dt} = f(x,t) = a(x)\\cdot b(t).$$\n", "In this case, we can write the equation, treating the derivative as\n", "product of differentials:\n", "\n", "$$\\frac{dx}{dt} = a(x)\\cdot b(t)$$\n", "and integrate the above expression on both sides. Because the left side is not\n", "it explicitly includes time integration after $x$ we are doing as if $x$ was\n", "independent variable.\n", "\n", "### Example:\n", "\n", "$$\\frac{dx}{dt} = - k x$$\n", "\n", "$$\\frac{dx}{x} =-k dt$$\n", "\n", "$$\\log(x(t)) =-k t + C$$\n", "assuming that $x>0$.\n", "When we solve $x$ we have:\n", "\n", "$$x(t) = e^{-kt +C}$$\n", "Let's see how the integration constant depends on the initial condition. Let\n", "$x(0)=x_0$, we have:\n", "\n", "$$x(t=0) = e^{-k0+C} =e^{C}.$$\n", "So we can save the solution with the initial condition $x(0)=x_0$ as:\n", "\n", "$$x(t) =x_0 e^{-kt}.$$\n", "Let's check if this solution agrees with the obtained approximate method:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "L = []\n", "k = 1.0\n", "dt = 0.01\n", "x0 = 1.2\n", "X = x0\n", "czas = 0\n", "xt = [X]\n", "ts = [0]\n", "for i in range(500):\n", " X = X + dt*(-k*X)\n", " czas = czas + dt\n", " if not i%10:\n", " xt.append(X)\n", " ts.append(czas)\n", " \n", "var('t')\n", "p1 = plot( x0*exp(-k*t) ,(t,0,5),color='red',figsize=(5,2) )\n", "p2 = point(zip(ts,xt))\n", "p1 + p2" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "## Solving ODEs using `desolve_odeint`\n", "\n", "\n", "There are several algorithms built in to the Sage system that significantly\n", "more accurately and more efficiently solve differential equations. Without\n", "going into the details of their implementation, it is worth learning them\n", "use.\n", "\n", "A good choice in general case is the function `desolve_ system:`\n", "\n", "`desolve_odeint (right sides of differential equations, initial conditions, times, searched)`\n", "\n", "\n", "For our example, we have the use of this procedure looks like the following\n", "way:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = -k*x\n", "ic = 1.2\n", "t = srange(0,5.01,0.5)\n", "sol = desolve_odeint(f,ic,t,x)\n", "p = points(zip(t,sol[:,0]),size=40,color='green')\n", "(p1 + p).show()\n", "print k, t" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "The solution is passed in the form of a matrix (in fact, type\n", "np.array from the numpy package) in which for every $n$ equations each row contains\n", "$n$ variable values in subsequent time periods.\n", "In our case, we have one equation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sol.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "type(sol)" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "### Example: harmonic oscillator\n", "\n", "A system of two differential equations corresponding to the motion of a particle in\n", "Potential (1d)\n", "\n", "$$U(x) = \\frac{1}{2} k x^2$$\n", "\n", "Newton's equation:\n", "\n", "$$m \\ddot x = m a = F = -U'(x) = -k x$$\n", "\n", "what can you save:\n", "\n", "$$\\begin{cases} \\dot x = v \\\\ \\dot v = - k x \\end{cases}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "var('t')\n", "var('x, v')\n", "k = 1.2\n", "times = srange(0.0, 11.0, 0.025, include_endpoint=True) \n", "sol = desolve_odeint([v, -k*x], [1,0], times, [x,v])" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "The solution is a numpy array (see [Introduction to\n", "numpy] (https://sage2.icse.us.edu.pl/home/pub/114/)), which can be\n", "conveniently and efficiently searched by the \"slicing\" technique, e.g." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "lang": "en" }, "outputs": [], "source": [ "sol [::200, :]" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "The parametric dependence of $(x(t),v(t))$ can be presented on the plane\n", "(X, v):\n", "\n", "The dependencies on time, speed and position are given by functions\n", "periodic:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "var('x v')\n", "k = 1.2\n", "sol = desolve_odeint([v, -k*x], [3.1,0], times, [x,v])\n", "px = line(zip(times,sol[:,0]),figsize=(5,2))\n", "pv = line(zip(times,sol[:,1]),figsize=(5,2),color='red')\n", "px+pv" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "Because this system is known as a harmonic oscillator and we know that\n", "solution for the initial condition $x(0)=1$, $v(0)=0$ is in the form:\n", "\n", "$$x(t) = \\cos(\\sqrt{k}t), v(t) = -\\sin(\\sqrt{k}t).$$\n", "\n", "therefore, we can compare the result of the approximate method and the solution\n", "Analytical.\n", "\n", "An analytical solution can also be obtained using the Sage function\n", "desolve, which solves the differential equations symbolically:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sage: var('t k')\n", "sage: assume(k>0)\n", "sage: x = function('x')(t)\n", "sage: de = diff(x,t,2) == -k*x\n", "sage: desolve(de, x,ivar=t)" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "Even if we know the form of the differential equation solution, then we can\n", "always use desolve, this is the correct application of the condition\n", "initial. Take, for example, the harmonic oscillator in which at the moment\n", "the initial $x(0)=x_0$ and $v(0)=v_0$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "var('t k')\n", "assume(k>0)\n", "x = function('x')(t)\n", "de = diff(x,t,2) == -k*x\n", "var('v0,x0')\n", "show( desolve(de, x,ics=[0,x0,v0],ivar=t))" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "Let's compare the numerical and analytic solution for the condition\n", "initial $x_0,v_0=0,1$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sage: var('t x v')\n", "sage: k=1.22\n", "sage: sol = desolve_odeint([v, -k*(x)], [1.,0], times, [x,v])\n", "sage: px = line(zip(times,sol[:,0]),figsize=(5,2))\n", "sage: px+plot(cos(sqrt(k)*t),(t,0,10),color='green')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sage: var('t')\n", "sage: pv = line(zip(times,sol[:,1]),figsize=(5,2),color='red')\n", "sage: pv+plot(-sqrt(k)*sin(sqrt(k)*t),(t,0,10),color='green')" ] }, { "cell_type": "markdown", "metadata": { "lang": "en" }, "source": [ "### Example 2: mathematical pendulum:\n", "\n", "Newton's equation:\n", "\n", "$$m \\ddot x = m a = F = -U'(x) = -k \\sin(x)$$\n", "\n", "can be written in a form of a system od two ODEs:\n", "\n", "$$\\begin{cases} \\dot x = v \\\\ \\dot v = - k \\sin(x) \\end{cases}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\newpage" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 8.7", "language": "", "name": "sagemath" }, "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 }