🕰️ Módulo 6: Ecuaciones Diferenciales Ordinarias (EDOs)
#

⏳ EDOs: Modelando la Evolución Temporal de un Sistema
#

Una Ecuación Diferencial Ordinaria (EDO) relaciona una función desconocida con sus derivadas y es la herramienta principal para modelar la evolución temporal de sistemas físicos (ej. cinemática, circuitos, transferencia de calor). Resuelven problemas de valor inicial, donde se busca $y(t)$ conociendo el estado del sistema en $t_0$, $y(t_0)$.

$$\frac{dy}{dt} = f(t, y)$$

Resolver una EDO numéricamente significa determinar el valor de la función $y(t)$ en una serie de puntos discretos, $y_0, y_1, y_2, \dots$, mediante un proceso iterativo que utiliza la tasa de cambio actual para proyectar el estado futuro:

$$\text{Nuevo Estado} = \text{Estado Actual} + (\text{Tasa de Cambio} \times \text{Paso de Tiempo } h)$$

#

🛠️ Métodos Clásicos de Solución (Paso Simple)
#

1. Método de Euler
#

El método de Euler es el algoritmo más simple, basado en la aproximación lineal de la Serie de Taylor de primer orden. Es el punto de partida para entender la integración de EDOs.

Euler Explícito (Forward Euler)
#

Este método utiliza la pendiente calculada en el punto inicial ($t_i, y_i$) para proyectar el nuevo punto ($y_{i+1}$). Es el más sencillo de implementar.

$$\mathbf{y_{i+1} = y_i + f(t_i, y_i) \cdot h}$$

Ventajas: Cálculo directo.
Desventajas: Condicionalmente estable; requiere pasos $h$ muy pequeños para EDOs con tasas de decaimiento rápidas (problemas stiff o rígidos).

Euler Implícito (Backward Euler)
#

Este método utiliza la pendiente calculada en el punto final ($t_{i+1}, y_{i+1}$) para determinar el nuevo valor.

$$\mathbf{y_{i+1} = y_i + f(t_{i+1}, y_{i+1}) \cdot h}$$

Diferencia Crucial: La variable $y_{i+1}$ aparece en ambos lados de la ecuación, dentro de la función $f$. Esto significa que $y_{i+1}$ no se puede despejar directamente, sino que se debe resolver una ecuación algebraica (lineal o no lineal) en cada paso.
Ventajas: Es incondicionalmente estable. Es excelente para problemas stiff, permitiendo el uso de pasos $h$ grandes.
Desventajas: Requiere resolver un sistema de ecuaciones en cada iteración, lo que lo hace computacionalmente más costoso por paso que el explícito.

2. Métodos de Runge-Kutta de Cuarto Orden (RK4)
#

El método de RK4 es el más utilizado para EDOs no rígidas debido a su alta precisión (error global $O(h^4)$). Utiliza un promedio ponderado de cuatro estimaciones de la pendiente dentro del intervalo $[t_i, t_{i+1}]$:

$$k_1 = f(t_i, y_i)$$$$k_2 = f(t_i + \frac{1}{2}h, y_i + \frac{1}{2}k_1 h)$$$$k_3 = f(t_i + \frac{1}{2}h, y_i + \frac{1}{2}k_2 h)$$

$$k_4 = f(t_i + h, y_i + k_3 h)$$$$\mathbf{y_{i+1} = y_i + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) h}$$

⚙️ Sistemas de EDOs de Orden Superior
#

Cualquier EDO de orden $m$ o sistema acoplado debe transformarse primero en un sistema de $m$ EDOs acopladas de primer orden.

Ejemplo: EDO de segundo orden (Vibración)

$$\frac{d^2x}{dt^2} + 0.5\frac{dx}{dt} + 16x = 0$$

Transformación:

Definir $y_1 = x$ (posición) y $y_2 = \frac{dx}{dt}$ (velocidad).
Sistema de primer orden: $\begin{cases} \frac{dy_1}{dt} = y_2 \\ \frac{dy_2}{dt} = -0.5y_2 - 16y_1 \end{cases}$

Los métodos numéricos (Euler, RK4) se aplican a cada ecuación del sistema simultáneamente.

💻 Ejemplo de Programación Manual (Euler Explícito)
#

El siguiente ejemplo implementa el Método de Euler Explícito sin utilizar librerías de terceros (solo math para funciones básicas) para resolver $\frac{dy}{dt} = 0.5y$, con $y(0)=1$ en $t \in [0, 4]$.

import math

# 1. Definición de la EDO (f(t, y))
def f_edo(t, y):
    """Define la tasa de cambio dy/dt = 0.5 * y"""
    return 0.5 * y

# 2. Parámetros Iniciales y de Simulación
t_inicial = 0.0
t_final = 4.0
y_inicial = 1.0
h = 0.5  # Tamaño del paso

# Inicialización de la solución
t_valores = [t_inicial]
y_valores = [y_inicial]

# 3. Implementación del Método de Euler Explícito
n_pasos = int((t_final - t_inicial) / h)
y_actual = y_inicial
t_actual = t_inicial

print(f"--- 🧪 Método de Euler Explícito (h={h}) ---")
print(f"Paso 0: t={t_actual:.2f}, y={y_actual:.4f}")

for i in range(1, n_pasos + 1):
    # Cálculo de la pendiente (f(ti, yi))
    pendiente = f_edo(t_actual, y_actual)
    
    # Nuevo valor de y (Euler Explícito)
    y_siguiente = y_actual + pendiente * h
    
    # Actualización para el siguiente paso
    t_actual = t_inicial + i * h
    y_actual = y_siguiente
    
    t_valores.append(t_actual)
    y_valores.append(y_actual)
    
    print(f"Paso {i}: t={t_actual:.2f}, y={y_actual:.4f}")

# 4. Comprobación (Valor exacto en t=4: y(4) = 1 * e^(0.5*4) = e^2 ≈ 7.3891)

📚 Herramientas de Python: `scipy.integrate.solve_ivp`
#

Para resolver problemas reales, se utiliza scipy.integrate.solve_ivp. Esta función emplea métodos de paso adaptativo (como RK45 o BDF), ajustando automáticamente $h$ para maximizar la eficiencia y mantener el error bajo tolerancia, superando la limitación de la estabilidad de los métodos manuales de paso fijo.

import numpy as np
from scipy.integrate import solve_ivp

# Definición del Sistema de EDOs (Oscilador Armónico: d2x/dt2 = -9.8x)
def sistema_edos(t, y):
    # y = [posicion, velocidad]
    dy1_dt = y[1] 
    dy2_dt = -9.8 * y[0] 
    return [dy1_dt, dy2_dt]

# Condiciones Iniciales [x(0), x'(0)]
y0 = [1.0, 0.0] 
t_span = [0, 10] 

# Resolución con solve_ivp (RK45 por defecto, paso adaptativo)
solucion = solve_ivp(
    sistema_edos, 
    t_span, 
    y0, 
    method='RK45'
)

# Resultados
print("\n--- 🚀 Simulación con solve_ivp ---")
print(f"Posición final (t=10s): {solucion.y[0][-1]:.4f}")

📝 Ejercicios Propuestos
#

A. EDOs de Primer Orden (Método de Euler)
#

Crecimiento Poblacional Simple (Euler Manual): Resuelve $\frac{dy}{dt} = 0.1y$, con $y(0) = 100$, en $t \in [0, 5]$.
- Tarea: Implementa manualmente el Método de Euler Explícito con un paso $h = 1.0$.
Decaimiento Radioactivo (Euler Manual): El decaimiento radioactivo sigue $\frac{dA}{dt} = -0.05A$, con $A(0) = 50$.
- Tarea: Usa el Método de Euler Explícito con $h=0.5$ para estimar la cantidad de sustancia $A$ en $t=2.0$.

B. EDOs de Primer Orden (Método RK4 y `solve_ivp`)
#

Carga/Descarga de un Condensador (Circuito RC): La carga $Q(t)$ en un circuito RC sigue $\frac{dQ}{dt} = \frac{V_s}{R} - \frac{Q}{RC}$. (Usar $V_s=10 \text{ V}$, $R=2 \Omega$, $C=0.1 \text{ F}$).
- Tarea: Utiliza scipy.integrate.solve_ivp para encontrar la carga $Q(t)$ en el intervalo $[0, 5]$ segundos, con $Q(0)=0$. ¿Cuál es la carga en $t=5$?
Caída Libre con Resistencia del Aire: La velocidad $v$ de un objeto que cae es: $\frac{dv}{dt} = g - \frac{c}{m}v$. (Usar $g=9.8 \text{ m/s}^2$, $m=5 \text{ kg}$, $c=0.2 \text{ kg/s}$).
- Tarea: Resuelve la EDO con solve_ivp para encontrar la velocidad en $[0, 20]$ segundos, asumiendo $v(0)=0$. ¿Cuál es la velocidad terminal?

C. Sistemas de EDOs de Orden Superior
#

Oscilador Amortiguado: La posición $x$ de un sistema de masa-resorte con amortiguamiento es: $\frac{d^2x}{dt^2} + 0.5\frac{dx}{dt} + 16x = 0$.
- Tarea: Transforma esta EDO a un sistema de dos EDOs de primer orden.
- Tarea: Resuelve el sistema usando scipy.integrate.solve_ivp en el intervalo $[0, 15]$ con condiciones iniciales $x(0)=1$ y $x'(0)=0$.
Sistema Presa-Depredador (Lotka-Volterra): Simula las poblaciones de conejos ($x$) y zorros ($y$).
$$\frac{dx}{dt} = x(0.5 - 0.01y)$$

$$\frac{dy}{dt} = y(-0.2 + 0.005x)$$
- Tarea: Utiliza solve_ivp para simular las poblaciones en $[0, 100]$ con $x(0)=50$ y $y(0)=20$.
Diferencia de Estabilidad (Análisis Conceptual): Para la EDO $\frac{dy}{dt} = -20y$ (un problema rígido), si intentas usar el Método de Euler Explícito con $h=0.1$, la solución numérica explota.
- Pregunta de Análisis: ¿Por qué un ingeniero preferiría un método Implícito (o un método adaptativo como BDF en solve_ivp) en lugar de Euler Explícito para este tipo de problema?
Péndulo Simple No Lineal: La ecuación del péndulo es $\frac{d^2\theta}{dt^2} + \frac{g}{L}\sin(\theta) = 0$. (Usar $g=9.8 \text{ m/s}^2$, $L=1 \text{ m}$).
- Tarea: Transforma a un sistema de primer orden.
- Tarea: Resuelve usando solve_ivp en $[0, 10]$ con condiciones iniciales $\theta(0)=1.5 \text{ radianes}$ (casi 90°) y $\theta'(0)=0$. Grafica el ángulo $\theta(t)$.