# Ecuaciones diferenciales

## Preparación

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



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

```


```julia
using DifferentialEquations

```


```julia
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:



```julia
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](https://en.wikipedia.org/wiki/Ordinary_differential_equation##Reduction_of_order).

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 




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

```

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



```julia
u₀ = 1/2; tspan_u = (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



```julia
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



```julia
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í](https://diffeq.sciml.ai/stable/tutorials/ode_example/##Choosing-a-Solver-Algorithm), 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í](https://en.wikipedia.org/wiki/Approximation_error).

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



```julia
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

```

### 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:

$$\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}$$

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

El movimiento generado por esta ecuaciones se ve similar a:




```julia
html"
<p align='center'>

</p>
"

```


```julia
## 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:



```julia
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

```


```julia
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

```


```julia
double_pendulum_problem = ODEProblem(double_pendulum, initial, tspan);

```


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

```


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

```

### 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{align}
\frac{d^2x}{dt^2}&=-\frac{\partial V}{\partial x}\\
\frac{d^2y}{dt^2}&=-\frac{\partial V}{\partial y}
\end{align}
$$

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.




```julia
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:



```julia
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

```

Resolvemos el problema:



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

```


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

```

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



```julia
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

```

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



```julia
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.



```julia
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:



```julia
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.



```julia
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

```

## Ecuaciones diferenciales Estocásticas



## Bibliografía
* [Documentación de SciML](https://sciml.ai/) 
* [Tutoriales de SciML](https://github.com/SciML/SciMLTutorials.jl)
* 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
* 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.
* 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.



449f81c57dd5db087202392ce8873a9af449990f
ipynb
Jupyter Notebook
_build/html/_sources/Ecuaciones_diferenciales.ipynb
galexbh/Intro-Julia-2021
ec97166de13c5ea793d754750ee6191af26760cb
MIT
2021-01-11T23:22:13.000Z
2021-03-21T19:43:09.000Z
_build/html/_sources/Ecuaciones_diferenciales.ipynb
galexbh/Intro-Julia-2021
ec97166de13c5ea793d754750ee6191af26760cb
MIT
2021-01-13T18:16:13.000Z
2021-01-15T22:18:25.000Z
_build/html/_sources/Ecuaciones_diferenciales.ipynb
galexbh/Intro-Julia-2021
ec97166de13c5ea793d754750ee6191af26760cb
MIT
2021-01-12T17:54:45.000Z
2021-04-13T06:21:14.000Z
Qwen/Qwen-72B
1. YES
2. YES