Ecuaciones diferenciales¶

Preparación¶

Proseguimos a activar el ambiente de trabajo y cargar los paquetes a utilizar.

import Pkg; Pkg.activate("."); Pkg.instantiate()

 Activating environment at `~/work/Intro-Julia-2021/Intro-Julia-2021/Project.toml`

using DifferentialEquations

using Plots

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:

begin
	size = 50
	Plots.default(size = (2600,2000),titlefontsize = size, tickfontsize = size, legendfontsize = size, guidefontsize = size, legendtitlefontsize = size, lw = 5)
end

Modelos de ecuaciones diferenciales ordinarias (ODEs)¶

Preliminares¶

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.

Una vez realizado, la forma que obtenemos, para ecuaciones autónomas (no dependientes de manera explícita del tiempo), es:

\[ x'(t) = F(x(t)) \]

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:

\[F: \mathbb{R}^n \rightarrow \mathbb{R}^n\]

Esto define un sistema de \(n\) ecuaciones a resolver. A continuación veremos algunos ejemplos.

Ejemplo introductorio: La estructura general¶

Se ejemplifica el procedimiento general de resolución de ecuaciones diferenciales utilizando el paquete DifferentialEquations utilizando el ejemplo:

\[\frac{du}{dx} = f(u, p, t) = \alpha u\]

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.)).

Los pasos se presentan a continuación

1. Definición de la función¶

f(u,p,t) = 1.01*u

f (generic function with 1 method)

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)\)).

2. Definición de las condiciones iniciales y de frontera¶

u₀ = 1/2; tspan_u = (0.0,1.0)

(0.0, 1.0)

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\)).

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)\}\).

3. Planteamiento del problema¶

prob_u = ODEProblem(f,u₀,tspan_u);

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.

Ejemplos de estos serán mostrados luego.

4. Resolución y exploración¶

sol_u = solve(prob_u, Tsit5(), reltol=1e-8, abstol=1e-8);

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.

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í, y de éstos se habla en mayor profundidad luego.

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í.

Posterior a resolver la ecuación diferencial, podemos explorar su comportamiento de manera gráfica:

begin
	plot(sol_u,linewidth=5, title="Solución aproximada al problema",
     xaxis="Time (t)", yaxis="u(t) (in μm)",label="Aproximación")
	plot!(sol_u.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="Solución exacta")
end

_images/Ecuaciones_diferenciales_17_0.svg

Péndulo caótico¶

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:

\[\begin{split}\frac{d}{dt} \begin{pmatrix} \alpha \\ l_\alpha \\ \beta \\ l_\beta \end{pmatrix}= \begin{pmatrix} 2\frac{l_\alpha - (1+\cos\beta)l_\beta}{3-\cos 2\beta} \\ -2\sin\alpha - \sin(\alpha + \beta) \\ 2\frac{-(1+\cos\beta)l_\alpha + (3+2\cos\beta)l_\beta}{3-\cos2\beta}\\ -\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} \end{pmatrix}\end{split}\]

Resolveremos este sistema utilizando los métodos del paquetes DifferentialEquations.jl.

El movimiento generado por esta ecuaciones se ve similar a:

html"
<p align='center'>
<img width='216' height='190' src='https://upload.wikimedia.org/wikipedia/commons/6/65/Trajektorie_eines_Doppelpendels.gif'>
</p>
"

## El prefijo const es opcional y no significa VALOR constante, si no tipo constante.
begin
	const m₁, m₂, L₁, L₂, g = 1, 2, 1, 2, 9.81   
	initial = [0, π/3, 0, 3pi/5]
	tspan = (0., 50.)
end;

Se define una función auxiliar para transformar de coordenadas polares a cartesianas:

function polar2cart(sol; dt=0.02, l1=L₁, l2=L₂, vars=(2,4))
    u = sol.t[1]:dt:sol.t[end]
    p1 = l1*map(x->x[vars[1]], sol.(u))
    p2 = l2*map(y->y[vars[2]], sol.(u))
    x1 = l1*sin.(p1)
    y1 = l1*-cos.(p1)
    (u, (x1 + l2*sin.(p2),
     y1 - l2*cos.(p2)))
end

polar2cart (generic function with 1 method)

function double_pendulum(xdot,x,p,t)
    xdot[1] = x[2]
    xdot[2] = - ((g*(2*m₁+m₂)*sin(x[1]) + m₂*(g*sin(x[1]-2*x[3]) + 
				2*(L₂*x[4]^2+L₁*x[2]^2*cos(x[1]-x[3]))*sin(x[1]-x[3])))/
		        (2*L₁*(m₁+m₂-m₂*cos(x[1]-x[3])^2)))
    xdot[3] = x[4]
    xdot[4] = (((m₁+m₂)*(L₁*x[2]^2+g*cos(x[1])) + 
			   L₂*m₂*x[4]^2*cos(x[1]-x[3]))*sin(x[1]-x[3]))/
			   (L₂*(m₁+m₂-m₂*cos(x[1]-x[3])^2))
end

double_pendulum (generic function with 1 method)

double_pendulum_problem = ODEProblem(double_pendulum, initial, tspan);

sol = solve(double_pendulum_problem, Vern7(), abs_tol=1e-10, dt=0.05)

retcode: Success
Interpolation: specialized 7th order lazy interpolation
t: 302-element Array{Float64,1}:
  0.0
  0.05
  0.1223079052234777
  0.21763564439678446
  0.32592132827555614
  0.45428546205171927
  0.609957416358271
  0.7734779293497827
  0.9540060574595153
  1.1799033713407105
  1.403442437872348
  1.5857800225971301
  1.7547534996346459
  ⋮
 48.10013873938749
 48.26086940221816
 48.51276681735055
 48.68829170708072
 48.90897270610238
 49.074198457665524
 49.26762992059672
 49.41331132773069
 49.58656471302393
 49.73635419325158
 49.929069338760755
 50.0
u: 302-element Array{Array{Float64,1},1}:
 [0.0, 1.0471975511965976, 0.0, 1.8849555921538759]
 [0.05276815671595484, 1.071438957072351, 0.09384176137401032, 1.8607344940613866]
 [0.13361722361756873, 1.1748571429557286, 0.2248435889699277, 1.749602555594999]
 [0.2537554611755441, 1.3400852896072526, 0.38108445951036796, 1.518757313193717]
 [0.40410119950098694, 1.4024410831505718, 0.5301603161398095, 1.2396279286879397]
 [0.5728649474072226, 1.169954418422348, 0.6718502052458136, 0.9854204230486661]
 [0.7088857697515267, 0.5369386864730226, 0.8084759913872458, 0.7786857878244005]
 [0.7353785982784866, -0.19338512009518283, 0.9169953717631636, 0.5291567670987339]
 [0.6472090513598847, -0.715269349725902, 0.9734453010510711, 0.0543580496299323]
 [0.46935732113951517, -0.7327102879675962, 0.887866370720516, -0.86279081696396]
 [0.37225084287475957, -0.08343101301168705, 0.5743822024889788, -1.9236009700729089]
 [0.3419732736206664, -0.6055412037246396, 0.19590319502068296, -1.9975724150479082]
 [0.11351472321247819, -2.0721701874347174, -0.08424980592934868, -1.281496345624497]
 ⋮
 [-0.42641322931863473, -0.3495803765321254, -1.0156807674608763, -0.6419671254245816]
 [-0.4973969181679087, -0.5006444011316504, -1.0574689947805727, 0.11204124902755815]
 [-0.630914582747322, -0.5147537467276151, -0.8911329999018179, 1.1679719640654815]
 [-0.67975290420246, 0.13973286964103626, -0.6514876665177031, 1.4454143327848226]
 [-0.4549919175434134, 1.9015781454592204, -0.36771741048108453, 1.0727040751077042]
 [-0.08586635119481108, 2.28552090644888, -0.19809628081445166, 1.112713727902364]
 [0.2498746728460613, 1.030480739331191, 0.08577212263454681, 1.8451390168903863]
 [0.3320815110176066, 0.25499621194857525, 0.3793211991189911, 2.071661061898589]
 [0.39238335257281043, 0.5691064464087552, 0.6955158132179681, 1.4770421393438187]
 [0.5050439638383725, 0.8813133553208242, 0.8649801258453232, 0.7945405264152658]
 [0.6711792406408479, 0.7375180832629122, 0.9459177363546021, 0.09015063400030744]
 [0.7166955252930275, 0.5339053276716273, 0.9454344718959518, -0.09727163755146666]

begin
	ts, ps = polar2cart(sol, l1=L₁, l2=L₂, dt=0.01)
	plot(ps...)
end

_images/Ecuaciones_diferenciales_26_0.svg

3. Sistema de Hénon-Heiles¶

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:

\[\begin{split} \begin{align} \frac{d^2x}{dt^2}&=-\frac{\partial V}{\partial x}\\ \frac{d^2y}{dt^2}&=-\frac{\partial V}{\partial y} \end{align} \end{split}\]

donde

\[V(x,y)={\frac {1}{2}}(x^{2}+y^{2})+\lambda \left(x^{2}y-{\frac {y^{3}}{3}}\right).\]

es conocido como el potencial de Hénon–Heiles. De éste puede derivarse su Hamiltoniano:

\[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).\]

Esta cantidad representa un invariante del sistema dinámico: Una cantidad conservada. En este caso, es la energía total del sistema.

begin
	## Parámetros
	initial₂ = [0.,0.1,0.5,0]
	tspan₂ = (0,100.)
	## V: Potencial, T: Energía cinética total, E: Energía total
	V(x,y) = 1//2 * (x^2 + y^2 + 2x^2*y - 2//3 * y^3)
	E(x,y,dx,dy) = V(x,y) + 1//2 * (dx^2 + dy^2);
end;

Definimos el modelo en una función:

function Hénon_Heiles(du,u,p,t)
    x  = u[1]
    y  = u[2]
    dx = u[3]
    dy = u[4]
    du[1] = dx
    du[2] = dy
    du[3] = -x - 2x*y
    du[4] = y^2 - y -x^2
end

Hénon_Heiles (generic function with 1 method)

Resolvemos el problema:

begin
	prob₂ = ODEProblem(Hénon_Heiles, initial₂, tspan₂)
	sol₂ = solve(prob₂, Vern9(), abs_tol=1e-16, rel_tol=1e-16);
end

retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 92-element Array{Float64,1}:
   0.0
   0.002767153900836259
   0.019390834923504494
   0.12119935187168689
   0.530301790748649
   1.1815820951240696
   1.9076818589199944
   2.760621588805973
   3.605356397694905
   4.619986523154658
   5.645461607705661
   6.628534205499573
   7.693582537682585
   ⋮
  87.998350712816
  89.12549586531028
  90.34475976436744
  91.46207401166643
  92.80356604989571
  93.85953574845666
  95.06904060013457
  96.26535461031364
  97.42922082465732
  98.66129338374192
  99.77739539466218
 100.0
u: 92-element Array{Array{Float64,1},1}:
 [0.0, 0.1, 0.5, 0.0]
 [0.0013835748315707893, 0.09999965542762242, 0.4999977028602053, -0.00024904536251688065]
 [0.00969468838029096, 0.09998307727738429, 0.4998872044881803, -0.0017456951735989768]
 [0.06042185852124362, 0.09933514655181921, 0.49560211524963727, -0.011034331488800994]
 [0.25058314954625777, 0.08601880052006405, 0.41888733977799825, -0.05742234854837401]
 [0.4444070725050881, 0.011729989621118292, 0.15861982132709898, -0.17862147581375334]
 [0.4473066382282281, -0.16320285100673101, -0.13126053699732168, -0.278919260652333]
 [0.2584572479002514, -0.35063663280889906, -0.2790546850571461, -0.09910177915852485]
 [0.004597212445906916, -0.2765685142768536, -0.312548979741221, 0.26654022206668787]
 [-0.2725985862807251, 0.0888755547420482, -0.17762619800611465, 0.36446229045208195]
 [-0.2359649523333103, 0.33539268583553616, 0.2639673601222507, 0.0963845261733765]
 [0.1265413408705088, 0.30963375665407095, 0.371111092458688, -0.1399949538666176]
 [0.3331864229707778, 0.02777322267904361, -0.008512405641007949, -0.3759318157658049]
 ⋮
 [0.4186977554232182, 0.02225898651098427, 0.21193927276905847, -0.1756847863354898]
 [0.40581590801727596, -0.26179610944305776, -0.18430403170494647, -0.25779590723676926]
 [0.09926832114771249, -0.32419636125803547, -0.28889656169215455, 0.21128012800300267]
 [-0.2114517878067161, 0.0730906156176862, -0.2152584827559103, 0.3958776870516688]
 [-0.1651820081856587, 0.39689185721893394, 0.3024690065429725, 0.05584713594737721]
 [0.1983089199690154, 0.31747189605713433, 0.2712578251378906, -0.20508331912354352]
 [0.2589198172932789, -0.09448897836357471, -0.1423573171660115, -0.4186496601369496]
 [0.012425073399517497, -0.4034830628308482, -0.22915816033577166, 0.01615845936515269]
 [-0.23244506354682215, -0.08009238032616647, -0.15414702520990145, 0.42836918054426903]
 [-0.18927473635920564, 0.34544617499081803, 0.2560388635086127, 0.20350117082002062]
 [0.1741544372800989, 0.41603304066553226, 0.26940699405279905, -0.07855686343548968]
 [0.2254464439669529, 0.3916330728294025, 0.18834406694149203, -0.14125776488921374]

plot(sol₂, vars=(1,2), title = "Órbita del sistema de Hénon-Heiles", xaxis = "x", yaxis = "y", leg=false)

_images/Ecuaciones_diferenciales_33_0.svg

Parece estar correctamente resuelto pero… examinando la evolución de la energía total/Hamiltoniano podemos encontrar lo siguiente:

begin
	energy = map(x->E(x...), sol₂.u)
	plot(sol₂.t, energy .- energy[1], title = "Cambio de la energía en el tiempo", xaxis = "Número de iteraciones", yaxis = "Cambio en energía")
end

_images/Ecuaciones_diferenciales_35_0.svg

¡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.

Para evitar eso, podemos utilizar algo conocido como un integrador simplético que considere la estructura de sistema Hamiltoniano que tiene nuestras ecuaciones.

Éste lo implementamos a continuación:

function HH_acceleration!(dv,v,u,p,t)
    x,y  = u
    dx,dy = dv
    dv[1] = -x - 2x*y
    dv[2] = y^2 - y -x^2
end;

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.

begin
	initial_positions = [0.0,0.1]
	initial_velocities = [0.5,0.0]
end;

Resolvemos el problema ahora como una ecuación de segundo orden pero con estructura simpléctica detectada:

begin
	prob₃ = SecondOrderODEProblem(HH_acceleration!, 
								 initial_velocities,
								 initial_positions,
								 tspan₂)
	sol₃ = solve(prob₃, KahanLi8(), dt=1/10)
end;

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.

begin
	energy₂ = map(x->E(x[3], x[4], x[1], x[2]), sol₃.u)
	plot(sol₃.t, energy₂ .- energy₂[1], title = "Cambio de la energía en el tiempo", xaxis = "Número de iteraciones", yaxis = "Cambio en energía")
end