Estabilidad numérica

01 · Punto de partida

Manso el sistema, loco el método

El decaimiento dM/dt = −λ·M es lo más manso que hay: cae suave a cero. Pero Euler con dt grande oscila (salta arriba y abajo) y con dt mayor diverge (cada salto más grande). El sistema no cambió: el método se volvió inestable.

La estabilidad numérica no es precisión: es no-explosión. Un dt preciso pero inestable no existe.

• ¿λ·dt está bajo el límite del método?
• ¿La solución oscila donde la exacta es monótona?
• ¿El error se achica o explota al avanzar?
• ¿Achicar dt o cambiar de método lo arregla?

02 · Definición

¿Qué es la estabilidad numérica?

Que el error no se amplifique paso a paso:

El factor (1 − λ·dt) decide: positivo converge, negativo alterna, mayor que 1 en módulo explota.

03 · Anatomía

Factor, límite y remedio

La inestabilidad en tres líneas:

Factor de amplificación

Lo que multiplica cada paso: (1 − λ·dt) en Euler. Su módulo debe ser < 1.

Límite del método

Euler explícito: λ·dt < 2. RK4 aguanta más; implícitos casi todo.

Síntoma

Oscilación donde la exacta es monótona: el primer aviso. Después viene la explosión.

Remedio

Achicar dt bajo el límite o cambiar de método (RK, implícito). Afinar el error no sirve si diverge.

04 · Inestabilidad en el curso

Dónde puede explotar

Cada rigidez esconde un límite:

1. 1
   Decaimiento rápido.

   λ grande con dt ingenuo: el caso de manual. Achicar dt.

2. 2
   Resorte rígido.

   k grande oscila con Euler explícito: semi-implícito o RK.

3. 3
   SIR con β alto.

   Contagio veloz pide dt chico en días: si I oscila, es numérico.

4. 4
   Reacciones rápidas.

   Química con escalas separadas (rígida): implícitos.

5. 5
   Difusión fina.

   Malla espacial + temporal: el límite ata dt a dx².

05 · Representación en Python

Factor y régimen en dos funciones

En Python, el factor predice y la simulación confirma: si λ·dt > 2, ni hace falta correr.

def factor(lam=1.0, dt=0.5):
    return 1 - lam * dt


def regimen(lam=1.0, dt=0.5):
    f = factor(lam, dt)
    if abs(f) >= 1:
        return f"inestable ({f:+.2f})"
    if f < 0:
        return f"oscila decayendo ({f:+.2f})"
    return f"estable suave ({f:+.2f})"


for dt in [0.5, 1.5, 2.5]:
    print(dt, regimen(1.0, dt))

+0,5 estable · −0,5 oscila · −1,5 inestable: el factor firma antes de simular.

Euler que explota

La historia confirma al factor:

def euler_decae(m0=100.0, lam=1.0, t_final=6.0, dt=0.5):
    m, t, hist = m0, 0.0, [(0.0, m0)]
    while t < t_final - 1e-9:
        m, t = m * (1 - lam * min(dt, t_final - t)), t + min(dt, t_final - t)
        hist.append((round(t, 1), round(m, 1)))
    return hist


print(euler_decae(dt=2.5)[-3:])

06 · Exploración

Laboratorio: del decaimiento a la explosión

Decaimiento λ con tu dt hasta T = 10 desde M₀ = 100. Dorada: exacta que cae suave. Teal con nodos: Euler. El medidor λ·dt dice el régimen: verde estable, ámbar oscila, rojo diverge. Cruzá el 2 y observá la explosión alternada.

La exacta nunca oscila: si Euler alterna, el culpable es dt, no el sistema. En divergencia el eje se autoescala y la exacta se aplasta abajo.

Preguntas para explorar

1. Con λ = 1 y dt = 0,5, ¿Euler oscila? ¿Qué factor hay?
2. Con λ = 2 y dt = 1,5, ¿qué régimen hay? ¿Y con dt = 0,5?
3. Si λ·dt = 1,9 estable pero feo, ¿lo aceptás? ¿Qué harías?

07 · Comprensión

Confusiones frecuentes

08 · Práctica guiada

Ejercicios con Python

Ejercicio 1: predecir con el factor

Con λ = 2, clasificá dt = 0,5, 1 y 1,5 sin simular. Después verificá el primero.

for dt in [0.5, 1.0, 1.5]:
    f = 1 - 2 * dt
    print(dt, round(f, 2))

Ejercicio 2: dt máximo estable

¿Cuál es el mayor dt estable con λ = 0,5? ¿Y con λ = 4 (reacción rápida)?

for lam in [0.5, 4.0]:
    print(lam, "→ dt <", round(2 / lam, 3))

Ejercicio 3: assert que salva

Escribí un Euler que se niegue a correr si λ·dt ≥ 2. Probalo con λ = 1 y dt = 2,5.

def euler_seguro(m0, lam, t_final, dt):
    assert lam * dt < 2, f"inestable: λ·dt = {lam * dt}"
    m, t = m0, 0.0
    while t < t_final - 1e-9:
        h = min(dt, t_final - t)
        m, t = m * (1 - lam * h), t + h
    return round(m, 2)


try:
    print(euler_seguro(100, 1.0, 6, 2.5))
except AssertionError as e:
    print("frenado:", e)

09 · Síntesis

Ideas para recordar

En el próximo tema compararemos soluciones: analítica frente a numérica.