Ecuaciones_diferenciales.ipynb

{
    "cells": [
    {
 	    "cell_type": "markdown", 
 	    "source": [ 
 		"# Ecuaciones diferenciales" 
 	    ], 
 	    "metadata": {} 
	},
        {
            "cell_type": "markdown",
            "source": [
                "## Preparación\n",
                "\n",
                "Proseguimos a activar el ambiente de trabajo y cargar los paquetes a utilizar.\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "import Pkg; Pkg.activate(\".\"); Pkg.instantiate()\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 67
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "using DifferentialEquations\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 68
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "using Plots\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 69
        },
        {
            "cell_type": "markdown",
            "source": [
                "La siguiente configuración por defecto se coloca para la versión específica de Julia, Plots y backend GR utilizada en construir este documento. Puede ser óptimo cambiarlo a otra para uso personal:\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\tsize = 50\n",
                "\tPlots.default(size = (2600,2000),titlefontsize = size, tickfontsize = size, legendfontsize = size, guidefontsize = size, legendtitlefontsize = size, lw = 5)\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 70
        },
        {
            "cell_type": "markdown",
            "source": [
                "## Modelos de ecuaciones diferenciales ordinarias (ODEs) \n",
                "### Preliminares\n",
                "Gracias a la gran eficiencia de métodos iterativos y de álgebra lineal, la forma preferencial de pensar en ecuaciones diferenciales es en su forma de *sistema de ecuaciones de primer órden*. Esto siempre es posible de obtener mediante [reducción de órden](https://en.wikipedia.org/wiki/Ordinary_differential_equation##Reduction_of_order).\n",
                "\n",
                "Una vez realizado, la forma que obtenemos, para ecuaciones autónomas (no dependientes de manera explícita del tiempo), es:\n",
                "\n",
                "$$\n",
                "x'(t) = F(x(t))\n",
                "$$\n",
                "\n",
                "donde $x: \\mathbb{R} \\rightarrow \\mathbb{R}^n$ es una trayectoria en el espacio para la cual deseamos resolver y obedece la dinámica del campo vectorial:\n",
                "\n",
                "$$F: \\mathbb{R}^n \\rightarrow \\mathbb{R}^n$$\n",
                "\n",
                "Esto define un sistema de $n$ ecuaciones a resolver. A continuación veremos algunos ejemplos.\n",
                "\n",
                "### Ejemplo introductorio: La estructura general\n",
                "\n",
                "Se ejemplifica el procedimiento general de resolución de ecuaciones diferenciales utilizando el paquete `DifferentialEquations` utilizando el ejemplo:\n",
                "\n",
                "$$\\frac{du}{dx} = f(u, p, t) = \\alpha u$$\n",
                "\n",
                "donde $u$ es el observable dinámico a encontrar, y la forma general $f(u, p, t)$ describe la evolución temporal de $u$ y puede depender de dicho observable: $u$, el tiempo paramétrico que describe su dinámica: $t$, y parámetros potencialmente exteriores que caracterizan el contexto específico del fenómeno (ejemplo: constantes físicas (masa, aceleración debido a gravedad, permitividad, etc.), indicadores socio-económicos (tasas constantes, natalidad, migración, etc.)).\n",
                "\n",
                "Los pasos se presentan a continuación \n",
                "\n",
                "#### 1. Definición de la función \n",
                "\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "f(u,p,t) = 1.01*u\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 71
        },
        {
            "cell_type": "markdown",
            "source": [
                "Como descrito anteriormente, $f$ depende de $u, p, t$, pero solamente es obligatoria la dependencia de $u$ (salvo el caso trivial donde $f$ es constante y en dado caso no necesitamos resolución numérica), mientras que, en este caso, no tenemos parámetros movibles $p$ ni tampoco el tiempo $t$ explícito, aunque por supuesto $u$ se asume ser función del tiempo ($u = u(t)$).\n"
            ],
            "metadata": {}
        },
        {
            "cell_type": "markdown",
            "source": [
                "#### 2. Definición de las condiciones iniciales y de frontera\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "u₀ = 1/2; tspan_u = (0.0,1.0)\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 72
        },
        {
            "cell_type": "markdown",
            "source": [
                "Aquí $u_0$ es el valor inicial tomado por la variable dinámica $u$. El tiempo inicial se define en el intervalo `tspan_u`, que es una tupla de un valor inicial y un valor final ($t_0$, $t_f$). \n",
                "\n",
                "Es este intervalo el cual será particionado en una cadena discreta de puntos $\\{t_0, t_1, t_2, \\ldots, t_{n-1}, t_f\\}$ y obtener los valores de $u$ en dichos tiempos: $\\{u(t_0), u(t_1), u(t_2), \\ldots, u(t_{n-1}), u(t_f)\\}$.\n",
                "\n",
                "#### 3. Planteamiento del problema\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "prob_u = ODEProblem(f,u₀,tspan_u);\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 73
        },
        {
            "cell_type": "markdown",
            "source": [
                "Este es un ejemplo de un problema planteado como una ecuación diferencial ordinaria de primer orden. En su lugar, pudimos haber definido una de orden superior, una estocástica, integrodiferencial, etc.\n",
                "\n",
                "Ejemplos de estos serán mostrados luego.\n",
                "\n",
                "#### 4. Resolución y exploración\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "sol_u = solve(prob_u, Tsit5(), reltol=1e-8, abstol=1e-8);\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 74
        },
        {
            "cell_type": "markdown",
            "source": [
                "La función `solve` es una de las funciones más complejas (en cuanto a completitud y flexibilidad) del paquete, siendo la misma redefinida mediante multiple-dispatch para poder reaccionar sin aparente dificultad a cualquier tipo de problema planteado y permitiendo argumentos que personalicen el proceso de solución para cada una de ellas.\n",
                "\n",
                "Como ejemplo, aquí tenemos que el problema `prob_u` se resuelve mediante un **solver** llamado `Tsit5`. La lista de solvers incluidos en el paquete puede encontrarse en su documentación [aquí](https://diffeq.sciml.ai/stable/tutorials/ode_example/##Choosing-a-Solver-Algorithm), y de éstos se habla en mayor profundidad luego. \n",
                "\n",
                "Además, se ilustra que podemos elegir un error relativo y error absoluto que queremos aceptar para que el solver lo garantice. Leer más sobre los tipos de errores de aproximación [aquí](https://en.wikipedia.org/wiki/Approximation_error).\n",
                "\n",
                "Posterior a resolver la ecuación diferencial, podemos explorar su comportamiento de manera gráfica:\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\tplot(sol_u,linewidth=5, title=\"Solución aproximada al problema\",\n",
                "     xaxis=\"Time (t)\", yaxis=\"u(t) (in μm)\",label=\"Aproximación\")\n",
                "\tplot!(sol_u.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label=\"Solución exacta\")\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 75
        },
        {
            "cell_type": "markdown",
            "source": [
                "### Péndulo caótico\n",
                "El péndulo doble es un sistema muy famosamente estudiado por exhibir un comportamiento **caótico** (cuya definición matemática es rigurosa). Sus ecuaciones del movimiento son las siguientes:\n",
                "\n",
                "$$\\frac{d}{dt}\n",
                "\\begin{pmatrix}\n",
                "\\alpha \\\\ l_\\alpha \\\\ \\beta \\\\ l_\\beta\n",
                "\\end{pmatrix}=\n",
                "\\begin{pmatrix}\n",
                "2\\frac{l_\\alpha - (1+\\cos\\beta)l_\\beta}{3-\\cos 2\\beta} \\\\\n",
                "-2\\sin\\alpha - \\sin(\\alpha + \\beta) \\\\\n",
                "2\\frac{-(1+\\cos\\beta)l_\\alpha + (3+2\\cos\\beta)l_\\beta}{3-\\cos2\\beta}\\\\\n",
                "-\\sin(\\alpha+\\beta) - 2\\sin(\\beta)\\frac{(l_\\alpha-l_\\beta)l_\\beta}{3-\\cos2\\beta} + 2\\sin(2\\beta)\\frac{l_\\alpha^2-2(1+\\cos\\beta)l_\\alpha l_\\beta + (3+2\\cos\\beta)l_\\beta^2}{(3-\\cos2\\beta)^2}\n",
                "\\end{pmatrix}$$\n",
                "\n",
                "Resolveremos este sistema utilizando los métodos del paquetes `DifferentialEquations.jl`. \n",
                "\n",
                "El movimiento generado por esta ecuaciones se ve similar a:\n",
                "\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "html\"\n",
                "<p align='center'>\n",
                "<img width='216' height='190' src='https://upload.wikimedia.org/wikipedia/commons/6/65/Trajektorie_eines_Doppelpendels.gif'>\n",
                "</p>\n",
                "\"\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 76
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "## El prefijo const es opcional y no significa VALOR constante, si no tipo constante.\n",
                "begin\n",
                "\tconst m₁, m₂, L₁, L₂, g = 1, 2, 1, 2, 9.81   \n",
                "\tinitial = [0, π/3, 0, 3pi/5]\n",
                "\ttspan = (0., 50.)\n",
                "end;\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 77
        },
        {
            "cell_type": "markdown",
            "source": [
                "Se define una función auxiliar para transformar de coordenadas polares a cartesianas:\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "function polar2cart(sol; dt=0.02, l1=L₁, l2=L₂, vars=(2,4))\n",
                "    u = sol.t[1]:dt:sol.t[end]\n",
                "",
                "    p1 = l1*map(x->x[vars[1]], sol.(u))\n",
                "    p2 = l2*map(y->y[vars[2]], sol.(u))\n",
                "",
                "    x1 = l1*sin.(p1)\n",
                "    y1 = l1*-cos.(p1)\n",
                "    (u, (x1 + l2*sin.(p2),\n",
                "     y1 - l2*cos.(p2)))\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 78
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "function double_pendulum(xdot,x,p,t)\n",
                "    xdot[1] = x[2]\n",
                "    xdot[2] = - ((g*(2*m₁+m₂)*sin(x[1]) + m₂*(g*sin(x[1]-2*x[3]) + \n",
                "\t\t\t\t2*(L₂*x[4]^2+L₁*x[2]^2*cos(x[1]-x[3]))*sin(x[1]-x[3])))/\n",
                "\t\t        (2*L₁*(m₁+m₂-m₂*cos(x[1]-x[3])^2)))\n",
                "    xdot[3] = x[4]\n",
                "    xdot[4] = (((m₁+m₂)*(L₁*x[2]^2+g*cos(x[1])) + \n",
                "\t\t\t   L₂*m₂*x[4]^2*cos(x[1]-x[3]))*sin(x[1]-x[3]))/\n",
                "\t\t\t   (L₂*(m₁+m₂-m₂*cos(x[1]-x[3])^2))\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 79
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "double_pendulum_problem = ODEProblem(double_pendulum, initial, tspan);\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 80
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "sol = solve(double_pendulum_problem, Vern7(), abs_tol=1e-10, dt=0.05)\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 81
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\tts, ps = polar2cart(sol, l1=L₁, l2=L₂, dt=0.01)\n",
                "\tplot(ps...)\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 82
        },
        {
            "cell_type": "markdown",
            "source": [
                "### 3. Sistema de Hénon-Heiles \n",
                "Es un sistema dinámico que modela el movimiento de una estrella orbitando alrededor de su centro galáctico mientras yace restringido en un plano. Éste es un ejemplo de un sistema Hamiltoniano, y tiene la forma:\n",
                "\n",
                "$$\n",
                "\\begin{align}\n",
                "\\frac{d^2x}{dt^2}&=-\\frac{\\partial V}{\\partial x}\\\\\n",
                "\\frac{d^2y}{dt^2}&=-\\frac{\\partial V}{\\partial y}\n",
                "\\end{align}\n",
                "$$\n",
                "\n",
                "donde\n",
                "\n",
                "$$V(x,y)={\\frac {1}{2}}(x^{2}+y^{2})+\\lambda \\left(x^{2}y-{\\frac {y^{3}}{3}}\\right).$$\n",
                "\n",
                "es conocido como el **potencial de Hénon–Heiles**. De éste puede derivarse su Hamiltoniano:\n",
                "\n",
                "$$H={\\frac {1}{2}}(p_{x}^{2}+p_{y}^{2})+{\\frac {1}{2}}(x^{2}+y^{2})+\\lambda \\left(x^{2}y-{\\frac {y^{3}}{3}}\\right).$$\n",
                "\n",
                "Esta cantidad representa un invariante del sistema dinámico: Una cantidad conservada. En este caso, es la energía total del sistema.\n",
                "\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\t## Parámetros\n",
                "\tinitial₂ = [0.,0.1,0.5,0]\n",
                "\ttspan₂ = (0,100.)\n",
                "\t## V: Potencial, T: Energía cinética total, E: Energía total\n",
                "\tV(x,y) = 1//2 * (x^2 + y^2 + 2x^2*y - 2//3 * y^3)\n",
                "\tE(x,y,dx,dy) = V(x,y) + 1//2 * (dx^2 + dy^2);\n",
                "end;\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 83
        },
        {
            "cell_type": "markdown",
            "source": [
                "Definimos el modelo en una función:\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "function Hénon_Heiles(du,u,p,t)\n",
                "    x  = u[1]\n",
                "    y  = u[2]\n",
                "    dx = u[3]\n",
                "    dy = u[4]\n",
                "    du[1] = dx\n",
                "    du[2] = dy\n",
                "    du[3] = -x - 2x*y\n",
                "    du[4] = y^2 - y -x^2\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 84
        },
        {
            "cell_type": "markdown",
            "source": [
                "Resolvemos el problema:\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\tprob₂ = ODEProblem(Hénon_Heiles, initial₂, tspan₂)\n",
                "\tsol₂ = solve(prob₂, Vern9(), abs_tol=1e-16, rel_tol=1e-16);\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 85
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "plot(sol₂, vars=(1,2), title = \"Órbita del sistema de Hénon-Heiles\", xaxis = \"x\", yaxis = \"y\", leg=false)\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 86
        },
        {
            "cell_type": "markdown",
            "source": [
                "Parece estar correctamente resuelto pero... examinando la evolución de la energía total/Hamiltoniano podemos encontrar lo siguiente:\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\tenergy = map(x->E(x...), sol₂.u)\n",
                "\tplot(sol₂.t, energy .- energy[1], title = \"Cambio de la energía en el tiempo\", xaxis = \"Número de iteraciones\", yaxis = \"Cambio en energía\")\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 87
        },
        {
            "cell_type": "markdown",
            "source": [
                "¡La energía total cambia! Eso quiere decir que la ley de conservación que esperamos que se cumpla físicamente no parece cumplirse en nuestra simulación. Esto delata un error en la resolución de la ecuación, específicamente por el método utilizado.\n",
                "\n",
                "Para evitar eso, podemos utilizar algo conocido como un **integrador simplético** que considere la estructura de sistema Hamiltoniano que tiene nuestras ecuaciones.\n",
                "\n",
                "Éste lo implementamos a continuación:\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "function HH_acceleration!(dv,v,u,p,t)\n",
                "    x,y  = u\n",
                "    dx,dy = dv\n",
                "    dv[1] = -x - 2x*y\n",
                "    dv[2] = y^2 - y -x^2\n",
                "end;\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 88
        },
        {
            "cell_type": "markdown",
            "source": [
                "Debemos ahora definir la condición inicial por separado, pues nuestro sistema, al ser Hamiltoniano, tiene una segmentación natural en esos pares de variables.\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\tinitial_positions = [0.0,0.1]\n",
                "\tinitial_velocities = [0.5,0.0]\n",
                "end;\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 89
        },
        {
            "cell_type": "markdown",
            "source": [
                "Resolvemos el problema ahora como una ecuación de segundo orden pero con estructura simpléctica detectada:\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\tprob₃ = SecondOrderODEProblem(HH_acceleration!, \n",
                "\t\t\t\t\t\t\t\t initial_velocities,\n",
                "\t\t\t\t\t\t\t\t initial_positions,\n",
                "\t\t\t\t\t\t\t\t tspan₂)\n",
                "\tsol₃ = solve(prob₃, KahanLi8(), dt=1/10)\n",
                "end;\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 90
        },
        {
            "cell_type": "markdown",
            "source": [
                "Ahora podemos observar cómo la energía se mantiene muy cercana a cero, oscilando solamente por errores de precisión numérica pero sin existir una tendencia de crecimiento anómala.\n"
            ],
            "metadata": {}
        },
        {
            "outputs": [],
            "cell_type": "code",
            "source": [
                "begin\n",
                "\tenergy₂ = map(x->E(x[3], x[4], x[1], x[2]), sol₃.u)\n",
                "\tplot(sol₃.t, energy₂ .- energy₂[1], title = \"Cambio de la energía en el tiempo\", xaxis = \"Número de iteraciones\", yaxis = \"Cambio en energía\")\n",
                "end\n",
                "",
                ""
            ],
            "metadata": {},
            "execution_count": 91
        },
        {
            "cell_type": "markdown",
            "source": [
                "## Ecuaciones diferenciales Estocásticas\n"
            ],
            "metadata": {}
        },
        {
            "cell_type": "markdown",
            "source": [
                "\n",
                "## Bibliografía\n",
                "* [Documentación de SciML](https://sciml.ai/) \n",
                "* [Tutoriales de SciML](https://github.com/SciML/SciMLTutorials.jl)\n",
                "* Hénon, Michel (1983), \\\"Numerical exploration of Hamiltonian Systems\\\", in Iooss, G. (ed.), Chaotic Behaviour of Deterministic Systems, Elsevier Science Ltd, pp. 53–170, ISBN 044486542X\n",
                "* Aguirre, Jacobo; Vallejo, Juan C.; Sanjuán, Miguel A. F. (2001-11-27). \\\"Wada basins and chaotic invariant sets in the Hénon-Heiles system\\\". Physical Review E. American Physical Society (APS). 64 (6): 066208. doi:10.1103/physreve.64.066208. ISSN 1063-651X.\n",
                "* Steven Johnson. 18.335J Introduction to Numerical Methods . Spring 2019. Massachusetts Institute of Technology: MIT OpenCourseWare, https://ocw.mit.edu. License: Creative Commons BY-NC-SA.\n",
                "\n"
            ],
            "metadata": {}
        }
    ],
    "nbformat_minor": 2,
    "metadata": {
        "language_info": {
            "file_extension": ".jl",
            "mimetype": "application/julia",
            "name": "julia",
            "version": "1.5.1"
        },
        "kernelspec": {
            "name": "julia-1.5",
            "display_name": "Julia 1.5.1",
            "language": "julia"
        }
    },
    "nbformat": 4
}