Sistemas, estabilidad y rigidez

01 · Más de una variable dinámica

Un estado que evoluciona como vector

Muchos modelos contienen variables acopladas: posición y velocidad, varias especies químicas, poblaciones que interactúan o temperaturas en distintos puntos. Se escriben como:

El vector 𝐲=[y₁,…,yₘ] representa el estado completo. La función 𝐟 devuelve simultáneamente la velocidad de todas sus componentes.

02 · Una representación uniforme

Convertir ecuaciones de orden superior

Toda ecuación de orden m puede transformarse en un sistema de m ecuaciones de primer orden. Para el oscilador amortiguado:

definimos y₁=x e y₂=x′:

Esta conversión permite aplicar Euler, Runge-Kutta y otros integradores sin diseñar un algoritmo distinto para cada orden.

function osciladorAmortiguado(omega, zeta, fuerza) {
  return (t, [posicion, velocidad]) => [
    velocidad,
    fuerza(t) - 2 * zeta * omega * velocidad - omega ** 2 * posicion
  ];
}

const modelo = osciladorAmortiguado(2, 0.15, t => Math.sin(t));
console.log("derivada en t=0:", modelo(0, [1, 0]));
console.log("derivada en t=1:", modelo(1, [0.4, -0.2]));

03 · Modos del sistema

El modelo lineal 𝐲′=𝐀𝐲

En un sistema lineal autónomo, una matriz constante determina la evolución:

Los autovalores λᵢ de 𝐀 describen modos característicos:

• Re(λᵢ)<0: el modo decae;
• Re(λᵢ)>0: el modo crece;
• Im(λᵢ)≠0: aparece oscilación;
• módulos muy diferentes: coexisten escalas temporales separadas.

Si la matriz puede diagonalizarse, la solución es una combinación de términos eλᵢt. Esta conexión lleva el análisis escalar de estabilidad a sistemas.

const desplazar = (estado, escala, derivada) =>
  estado.map((valor, i) => valor + escala * derivada[i]);

function pasoRK4Vector(f, t, estado, h) {
  const k1 = f(t, estado);
  const k2 = f(t + h / 2, desplazar(estado, h / 2, k1));
  const k3 = f(t + h / 2, desplazar(estado, h / 2, k2));
  const k4 = f(t + h, desplazar(estado, h, k3));
  return estado.map((valor, i) =>
    valor + h * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6);
}

const dosEscalas = (t, [u, v]) => [-u, -20 * v];
const siguiente = pasoRK4Vector(dosEscalas, 0, [1, 1], 0.05);
console.log("estado después de un paso:", siguiente);

Los arreglos auxiliares evitan mezclar componentes de etapas diferentes. Cada llamada a f recibe un estado vectorial coherente.

04 · Linealización local

El Jacobiano en problemas no lineales

Para 𝐲′=𝐟(t,𝐲), el Jacobiano respecto del estado es:

Cerca de un estado, una perturbación δ𝐲 evoluciona aproximadamente según:

Los autovalores del Jacobiano revelan velocidades locales de crecimiento, decaimiento u oscilación. En un problema no lineal cambian con el tiempo y con el estado, por lo que la restricción de paso también puede variar durante la integración.

05 · El método frente al decaimiento

Estabilidad absoluta

Aplicamos un método de un paso a la ecuación de prueba y′=λy. El resultado adopta la forma:

R(z) es la función de estabilidad. Si la solución exacta decae, queremos que el modo numérico tampoco crezca:

La región de estabilidad absoluta es el conjunto de valores complejos z que satisfacen esa condición. Un paso estable exige que cada producto hλᵢ relevante quede dentro de la región.

06 · Comparar integradores

Funciones y regiones de estabilidad

Mayor orden y mayor región estable son propiedades distintas. RK2 mejora la precisión de Euler, pero no amplía su intervalo estable sobre el eje real negativo.

const factores = {
  euler: z => 1 + z,
  rk2: z => 1 + z + z ** 2 / 2,
  rk4: z => 1 + z + z ** 2 / 2 + z ** 3 / 6 + z ** 4 / 24,
  eulerImplicito: z => 1 / (1 - z)
};

for (const z of [-1, -2, -3, -10]) {
  const linea = Object.entries(factores).map(([nombre, R]) =>
    `${nombre}: |R|=${Math.abs(R(z)).toFixed(3)}`);
  console.log(`z=${z} - ${linea.join("; ")}`);
}

07 · Escalas muy separadas

Qué significa que un problema sea rígido

Un problema es rígido cuando la estabilidad obliga a un método explícito a usar pasos mucho menores que los necesarios para representar la parte de la solución que nos interesa.

La rigidez aparece con frecuencia en:

• cinética química con reacciones rápidas y lentas;
• circuitos con constantes de tiempo muy diferentes;
• difusión espacial discretizada en mallas finas;
• modelos biológicos con transitorios rápidos y evolución prolongada.

No significa simplemente que la solución cambie rápido. Un transitorio puede desaparecer pronto, pero su autovalor muy negativo continúa restringiendo el paso de un método explícito.

08 · Un ejemplo mínimo

Dos modos de decaimiento

Consideremos:

Las soluciones son u=e−t y v=e−100t. La componente rápida es prácticamente cero después de un tiempo breve; luego solo interesa seguir el lento decaimiento de u.

Sin embargo, Euler explícito necesita:

para que el modo rápido no crezca artificialmente. Integrar hasta t=10 requiere más de 500 pasos por estabilidad, aunque la componente visible u cambia en una escala del orden de uno.

09 · Cuando el cálculo contradice al modelo

Oscilaciones y crecimiento artificial

Para el modo rápido v′=−100v, Euler produce:

• con h=0.005, el factor es 0.5: decae sin cambiar de signo;
• con h=0.015, el factor es −0.5: decae oscilando;
• con h=0.025, el factor es −1.5: oscila y crece.

La solución exacta es positiva y decreciente en los tres casos. La inestabilidad es un artefacto del discretizador, no una propiedad del sistema físico.

10 · Evaluar en el nuevo estado

Euler implícito

Euler hacia atrás utiliza la pendiente al final del paso:

El nuevo estado aparece en ambos lados. Para y′=λy puede despejarse:

Si Re(λ)<0, el factor permanece dentro del disco unidad para cualquier h>0. Euler implícito es A-estable y amortigua fuertemente modos muy rápidos.

11 · Sistemas lineales por paso

Euler implícito para 𝐲′=𝐀𝐲

La fórmula vectorial es:

Por tanto, cada paso resuelve un sistema lineal. Para el ejemplo diagonal de dos escalas el despeje es directo:

function pasoEulerImplicitoDiagonal(estado, lambdas, h) {
  if (estado.length !== lambdas.length) throw new Error("Dimensión incompatible");
  return estado.map((valor, i) => valor / (1 - h * lambdas[i]));
}

let estado = [1, 1];
const lambdas = [-1, -100];
const h = 0.1;
for (let n = 0; n < 10; n++) {
  estado = pasoEulerImplicitoDiagonal(estado, lambdas, h);
}

console.log("aproximación en t=1:", estado);
console.log("solución exacta:", [Math.exp(-1), Math.exp(-100)]);

El paso h=0.1 sería inestable para Euler explícito sobre el modo rápido. Euler implícito lo amortigua, aunque conserva solo precisión de primer orden para el modo lento.

12 · El costo oculto

Problemas implícitos no lineales

Si 𝐟 es no lineal, debemos encontrar 𝐲ₙ₊₁ que anule:

Newton aplica iteraciones:

Cada iteración puede requerir formar o aproximar el Jacobiano y resolver un sistema lineal. En problemas grandes se reutilizan factorizaciones, matrices dispersas y precondicionadores.

El paso debe rechazarse si Newton no converge. Aceptar silenciosamente una iteración incompleta mezcla error algebraico con error de integración.

13 · Diagnóstico y elección práctica

Cómo reconocer y tratar la rigidez

• Un método explícito reduce el paso repetidamente pese a que la solución visible es suave.
• Al aumentar levemente h aparecen oscilaciones o explosiones no físicas.
• El Jacobiano tiene autovalores estables con magnitudes muy separadas.
• RK4 necesita demasiadas evaluaciones para cubrir un intervalo largo.

Ante estas señales:

1. confirmá que el modelo y las unidades sean correctos;
2. estimá escalas o autovalores relevantes;
3. compará refinamientos y monitoreá invariantes físicos;
4. considerá métodos implícitos A-estables o especializados;
5. usá control adaptativo de error y registrá pasos rechazados.

Estabilidad permite calcular; precisión determina si el resultado calculado es útil. Ambas condiciones deben cumplirse.

14 · Exploración interactiva

Laboratorio: dos escalas y cuatro métodos

El sistema combina u′=−u con v′=−μv. Las curvas continuas son numéricas y las discontinuas son exactas. Ajustá μ, el método y la cantidad de pasos.

Preguntas para explorar

1. Con RK4 y μ=50, reducí los pasos hasta cruzar el límite aproximado z=−2.785.
2. Compará Euler y RK2: ¿por qué tienen el mismo umbral real aunque RK2 sea más preciso?
3. Elegí Euler implícito y aumentá μ. ¿Qué ocurre con el factor del modo rápido?
4. Una vez estable, aumentá los pasos y observá la precisión de la componente lenta.

15 · Comprensión

Confusiones frecuentes

16 · Práctica guiada

Ejercicios

Ejercicio 1: convertir a sistema

Transformá x″+4x′+5x=cos(t) en un sistema de primer orden.

Ejercicio 2: paso estable de Euler

Para un sistema con autovalores −2 y −80, determiná el máximo paso estable de Euler explícito.

Ejercicio 3: factor implícito

Calculá el factor de Euler implícito para λ=−100 y h=0.1.

Ejercicio 4: reconocer rigidez

Un modelo evoluciona suavemente durante 20 segundos, pero RK4 solo funciona con h<10⁻⁴ por un transitorio que desaparece en 0.01 segundos. ¿Qué sugiere?

Ejercicio 5: matriz implícita

Para 𝐀=diag(−1,−50) y h=0.2, escribí la matriz que Euler implícito debe resolver.

17 · Síntesis

Ideas para recordar

En el próximo tema reuniremos estos conceptos en un proyecto integrador: una biblioteca numérica reutilizable, verificable y extensible.