No2 - ODEs y NODEs - Datos, Ecuaciones Diferenciales e Inteligencia Artificial

Fecha: 13/04/2026

Estas notas cubren la definición de Ecuaciones Diferenciales Ordinarias (ODEs), su aplicación en modelado dinámico, la conversión de otros tipos de ecuaciones a ODEs, y cómo se utilizan para la estimación de parámetros a partir de datos observacionales, culminando con la introducción a las Neural ODEs (NODEs).

Ecuaciones diferenciales ordinarias (ODEs)¶

Las ODEs describen la evolución de un sistema en función de una única variable independiente, típicamente el tiempo $t$ .

Definición y Parámetros¶

Una ODE queda definida de la siguiente manera:

Variable de estado: $x(t) \in \mathbb{R}^n$ representa el estado del sistema en el instante $t$ .
Dinámica: Está determinada por una función $f$ , conocida como campo vectorial, y un conjunto de parámetros del modelo $\theta \in \mathbb{R}^p$ .
Forma general:
$\frac{dx}{dt} = f(x, t, \theta)$
(1)
donde $f: \mathbb{R}^n \times \mathbb{R} \times \mathbb{R}^p \rightarrow \mathbb{R}^n$ .
Condición inicial: $x(t_0) = x_0$ . Define el estado de partida del sistema (condición inicial) en el tiempo inicial $t_0$ .

Salvo casos muy particulares, las soluciones a ODEs no son analiticias. En la mayoria de los casos, vamos a recurrir a métodos numéricos para resolver dichas ecuaciones (ver Clase 4).

Reducción de Orden: El Oscilador Armónico

Una ecuación diferencial de segundo orden, como la del oscilador armónico forzado dada por

\frac{d^2x}{dt^2} + \omega^2x = f(x,t),

(2)

no se presenta inicialmente en la forma estándar de una ODE de primer orden ( $\dot{x} = f(x,t)$ ). Sin embargo, es posible transformar cualquier ecuación de orden superior en un sistema de ODEs de primer orden mediante la definición de variables de estado adicionales.

Conversión a sistema de primer orden. Definimos la velocidad como una nueva variable de estado $v = \frac{dx}{dt}$ . Esto nos permite descomponer la ecuación original en un sistema de dos ecuaciones de primer orden:

$\frac{dx}{dt} = v$
$\frac{dv}{dt} = -\omega^2x + f(x,t)$ con sus respectivas condiciones iniciales.

Representación matricial. Para el vector de estado $\mathbf{u}(t) = \begin{pmatrix} x \\ v \end{pmatrix}$ , el sistema se expresa de forma compacta:

\frac{d}{dt}\begin{pmatrix} x \\ v \end{pmatrix} = \begin{pmatrix} v \\ -\omega^2x+f \end{pmatrix}=\begin{pmatrix} 0 & 1 \\ -\omega^2 & 0\end{pmatrix}\begin{pmatrix}x \\ v \end{pmatrix}+\begin{pmatrix}0 \\ f \end{pmatrix}

(3)

La idea central es que siempre podemos llevar una ecuación de orden más alto a un sistema de ODEs de primer orden aumentando la dimensionalidad del vector de estado.

Ecuaciones Diferenciales Parciales (PDEs) y el Método de Líneas

El enfoque de reducir problemas a sistemas matriciales también puede aplicarse a Ecuaciones Diferenciales Parciales (PDEs) mediante la discretización espacial, una técnica conocida como el Método de Líneas (Ascher (2008)).

Ejemplo: Ecuación de Difusión. Consideremos un campo $u(x,t)$ , con $x \in \mathbb{R}$ , que evoluciona según la ecuación de difusión:

\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}

(4)

Para resolver este sistema, también debemos especificar condiciones iniciales ( $u(x,t)=u_0(x)$ ) y de borde (por ejemplo, $u(0,t)=0$ y $u(1,t)=1$ ).

Discretización espacial. Para transformar esta PDE en un sistema de ODEs, aproximamos la segunda derivada espacial utilizando diferencias finitas centradas con un paso espacial $h$ :

\frac{\partial^2 u}{\partial x^2} \approx \frac{u(x+h) - 2u(x) + u(x-h)}{h^2}

(5)

Al aplicar esto, la derivada espacial desaparece y nos queda una ecuación que depende únicamente de una derivada temporal, convirtiéndose efectivamente en una ODE.

Representación matricial. Definiendo un vector de estado $\mathbf{u}$ con los valores discretizados, obtenemos el siguiente sistema matricial:

\frac{d}{dt} \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{pmatrix} = \frac{D}{h^2} \begin{pmatrix} -2 & 1 & 0 & \dots \\ 1 & -2 & 1 & \dots \\ 0 & 1 & -2 & \dots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{pmatrix}

(6)

Donde la matriz de coeficientes es tridiagonal, con -2 en la diagonal principal y 1 en la sub y supra diagonal.

Ejemplos¶

El Modelo Dinámico: Lotka-Volterra (Depredador-Presa)

Este modelo describe la dinámica temporal del estado de un sistema compuesto por dos poblaciones: conejos ( $x$ ) y lobos ( $y$ ).

Intuición de la dinámica:

Si los conejos están solos, se reproducen exponencialmente: $\frac{dx}{dt}=\alpha x$ .
Si los lobos están solos, mueren exponencialmente: $\frac{dy}{dt}=-\beta y$ .
Para modelar la interacción (los lobos se comen a los conejos), se agrega un término no lineal $xy$ .

Las ecuaciones del modelo estan entonces dadas por

\frac{dx}{dt} = \alpha x - \gamma xy

(7)

\frac{dy}{dt} = -\beta y + \eta xy

(8)

Componentes del sistema:

Vector de estado: $\begin{pmatrix} x \\ y \end{pmatrix}(t; \theta)$ . Depende del tiempo $t$ y está parametrizado por $\theta$ . Para cada conjunto de parámetros, las curvas temporales (trayectorias) serán diferentes. Suelen ser oscilatorias.
Condición inicial: $x(t_0)=x_0$ e $y(t_0)=y_0$ .
Vector de parámetros: $\theta = \begin{pmatrix} \alpha \\ \gamma \\ \beta \\ \eta \end{pmatrix} \in \mathbb{R}^4$ . Como hay 4 parámetros y términos cruzados, el sistema es no lineal.
- $\alpha > 0$ : tasa de crecimiento de conejos.
- $\gamma > 0$ : tasa de depredación.
- $\beta > 0$ : tasa de mortalidad de lobos.
- $\eta > 0$ : tasa de crecimiento de lobos por interacción.

Inferencia estadística¶

En la realidad, si uno tuviera conocimiento absoluto de la dinámica, conocería la trayectoria perfecta. Sin embargo, nunca se observan estas trayectorias puras. En su lugar, se observan datos, usualmente en tiempos discretos, que se asemejan a la trayectoria subyacente, pero posiblemente contaminados con ruido aleatorio.

Ecuación del modelo observacional. Vamos a considerar un modelo para los datos de la forma

x_i^{\text{obs}} = x(t_i; \theta) + \varepsilon_i

(9)

Donde:

$x_i^{\text{obs}}$ : Observación en el tiempo $t_i$ .
$x(t_i; \theta)$ : Estado del modelo dinámico.
$\varepsilon_i$ : Ruido observacional.

Características del ruido observacional $\varepsilon_i$ . Ejemplos comununes incluyen:

Caso estándar: Se asume que los ruidos son independientes e idénticamente distribuidos (i.i.d.) de forma Gaussiana $\varepsilon_i \sim N(0, \sigma^2)$ , con valor medio nulo y varianza constante para cada sitio muestreado. Se supone que este ruido no depende del valor de $x$ , aunque no siempre es cierto.
Ruido correlacionado: Común en series de tiempo, donde la correlación entre dos errores es distinta de cero: $\mathbb{E}[\varepsilon_i, \varepsilon_j] \neq 0$ para $i \neq j$ .
- Nota estadística: Si dos distribuciones son Gaussianas y tienen correlación 0, son independientes. Si no son Gaussianas, tener correlación 0 no implica independencia.

Ajuste de Trayectorias

El problema central es: dadas las observaciones $x_i^{\text{obs}}$ e $y_i^{\text{obs}}$ , ¿cómo estimamos los parámetros $\theta$ que mejor describen la dinámica subyacente? El objetivo es convertir esto en un problema de optimización, buscando minimizar una función de costo $\mathcal{L}(\theta; \text{DATOS})$ que compare las observaciones reales con las trayectorias generadas por el modelo $x(t; \theta)$ . Esto se va a realizar mediante el método de cuadrádos mínimos no lineales, el cual en el contexto de este curso llamaremos ajuste de trayectorias (trajectory matching) Ramsay & Hooker (2017):

Ajuste de Trayectorias: A diferencia del fiteo lineal de una recta, acá se compara con una función que depende de $\theta$ de manera no trivial.

\min_{\theta} \sum_{i=1}^{N} (x_i^{\text{obs}} - x(t_i; \theta))^2

(10)

Más adelante veremos que minimizar el cuadrado de los errores asume inherentemente un ruido de naturaleza Gaussiana.

Ejemplo (continuado)¶

Inferencia en el sistema Lotka-Volterra: Si salimos al campo y medimos las poblaciones de presas ( $x_i^{\text{obs}}$ ) y depredadores ( $y_i^{\text{obs}}$ ) a lo largo del tiempo, nuestra función de costo (o pérdida) a minimizar será:

\min_{\theta} \mathcal{L}(\theta) = \sum_{i=1}^{N} \left[ (x_i^{\text{obs}} - x(t_i; \theta))^2 + (y_i^{\text{obs}} - y(t_i; \theta))^2 \right]

(11)

Un optimizador resolverá numéricamente las ODEs de Lotka-Volterra en cada iteración, ajustando sistemáticamente el vector $\theta = (\alpha, \gamma, \beta, \eta)^\top \in \mathbb{R}^4$ hasta que la trayectoria predicha pase lo más cerca posible de nuestros datos ruidosos.

Ecuaciones diferenciales ordinarias neuronales (NODEs)¶

Las Neural Ordinary Differential Equations (NODEs), introducidas formalmente por Chen et al. (2018), representan un cambio de paradigma al fusionar el aprendizaje profundo con los sistemas dinámicos continuos.

1. Generalización de Modelos Clásicos¶

Una forma intuitiva de entender las NODEs es viéndolas como una generalización de modelos dinámicos preexistentes. Por ejemplo, en el modelo de Lotka-Volterra, podemos reemplazar o aumentar los términos de interacción utilizando redes neuronales:

\frac{dx}{dt} = \alpha x - \text{NN}_1(x, y)

(14)

\frac{dy}{dt} = -\beta y + \text{NN}_2(x, y)

(15)

En este caso, los parámetros internos de $\text{NN}_1$ y $\text{NN}_2$ (pesos y sesgos) pasan a formar parte del vector global de parámetros a estimar, $\theta$ , potencialmente junto con $\alpha$ y $\beta$ .

2. Definición Formal¶

Si generalizamos por completo (sin suponer ningún término de crecimiento o muerte específico predefinido), obtenemos una NODE. En su forma más abstracta, una NODE parametriza la derivada del estado continuo de un sistema directamente a través de una red neuronal:

\frac{du}{dt} = \text{NN}(u; \theta)

(16)

Donde $u \in \mathbb{R}^n$ representa el estado del sistema en un tiempo dado, y $\theta$ engloba todos los parámetros entrenables de la red neuronal.

3. Propiedades y Consideraciones Clave¶

Aproximación universal: Al estar basadas en redes neuronales, las NODEs heredan la capacidad de ser aproximadores universales Goodfellow et al. (2016). Tienen la flexibilidad necesaria para aprender y representar una gama casi ilimitada de dinámicas continuas a partir de datos empíricos.
Intratabilidad numérica: Advertencia: Al llevar el modelo a este nivel de complejidad no lineal, un problema frecuente en la optimización es caer en regiones del espacio de parámetros donde el sistema se vuelve matemáticamente inestable o demasiado costoso de resolver para los solvers de las ecuaciones diferenciales.

Implementación Computacional en Julia¶

A continuación, implementamos el modelo Lotka-Volterra en Julia.

Nota: El código completo de este ejemplo se puede referenciar acá: 01_LV_forward. Ahí, además de resolver el sistema, simulamos que tenemos datos empíricos con ruido y vemos cómo cambia el paisaje de la función de pérdida (Loss Landscape). Esos detalles quedan disponibles en el enlace para quienes deseen profundizar.

Nosotros vamos a contar brevemente los componentes principales de la resolución numérica.

En Julia, usamos ! en el nombre de la función para indicar que modifica sus argumentos in-place (es decir, no tenemos un return)

using DifferentialEquations
using Plots
using Statistics
using Random

# Definición del sistema de Ecuaciones Diferenciales Ordinarias
function lotka_volterra!(du, u, p, t)
    x, y = u            # x = presas, y = depredadores (esto es lo mismo que hacer u[1], u[2])
    α, β, δ, γ = p      # Desempaquetamos el vector p
    du[1] = α * x - β * x * y
    du[2] = δ * x * y - γ * y
end

# Definimos los parámetros de la ecuación
α = 1.0     # Nacimiento de presas
β = 0.1     # Tasa de depredación
δ = 0.075   # Reproducción del depredador
γ = 1.5     # Muerte del depredador
p_true = [α, β, δ, γ]

# Condiciones iniciales y tiempo
u0 = [10.0, 5.0]
tspan = (0.0, 30.0)

# ODEProblem es una función de Julia para resolver una ecuación diferencial ordinaria, 
# por eso le mandamos la función, la condición inicial y un tiempo para resolver.
prob = ODEProblem(lotka_volterra!, u0, tspan, p_true) # Acá definimos el problema matemático

# Acá resolvemos el problema, elegimos con qué solver (en nuestro caso Tsit5) y cada 
# cuánto se va a guardar, nivel de tolerancia, nivel de error. Se define la parte numérica.
sol = solve(prob, Tsit5(), saveat=0.1) 

# En sol está la solución al problema
print(sol)

References¶

Ascher, U. M. (2008). Numerical methods for evolutionary differential equations. SIAM.
Ramsay, J., & Hooker, G. (2017). Dynamic data analysis. Springer New York, New York, NY. Doi, 10, 978–1.
Chen, R. T., Rubanova, Y., Bettencourt, J., & Duvenaud, D. K. (2018). Neural ordinary differential equations. Advances in Neural Information Processing Systems, 31.
Goodfellow, I., Bengio, Y., Courville, A., & Bengio, Y. (2016). Deep learning (Vol. 1). MIT press Cambridge.