Tres Métodos, Tres Regiones
Para la ecuación de prueba dy/dx = λy, tres métodos explícitos de EDO tienen las siguientes regiones de estabilidad en el plano complejo hλ:
Método de Euler (primer orden): la región de estabilidad es el disco |1 + hλ| ≤ 1, un círculo de radio 1 centrado en (-1, 0). Los hλ reales negativos deben estar en [-2, 0].
Runge-Kutta 2 (método del punto medio) (segundo orden): la región de estabilidad es |1 + hλ + (hλ)²/2| ≤ 1. Más grande que el disco de Euler, pero aún acotada.
Runge-Kutta 4 (cuarto orden): la región de estabilidad satisface |1 + hλ + (hλ)²/2 + (hλ)³/6 + (hλ)⁴/24| ≤ 1. Los hλ reales negativos se extienden a aproximadamente -2.785. La región es sustancialmente más grande que la de Euler.
Backward Euler (implícito): la región de estabilidad es el plano complejo entero excepto el disco |1 - hλ|⁻¹ > 1, equivalentemente |1/(1-hλ)| ≤ 1. Para λ en el semiplano izquierdo (Re(λ) < 0), esto es incondicionalmente estable — sin restricción en h de la estabilidad.
La Función de Amplificación
Para cualquier método de Runge-Kutta, el factor de amplificación por paso R(hλ) es una aproximación polinomial a e^(hλ):
- Euler: R(z) = 1 + z (truncado en grado 1)
- RK2: R(z) = 1 + z + z²/2 (truncado en grado 2)
- RK4: R(z) = 1 + z + z²/2 + z³/6 + z⁴/24 (truncado en grado 4)
La región de estabilidad es {z ∈ ℂ : |R(z)| ≤ 1}. La amplificación de la solución verdadera: |e^z| = e^(Re(z)). Para Re(z) < 0 (EDO estable), la solución verdadera decae. El método numérico es estable si |R(z)| ≤ 1 — coincidiendo con el comportamiento decayente.
Valores Propios Puramente Imaginarios: Sistemas Oscilatorios
Muchos sistemas físicos tienen valores propios puramente imaginarios: λ = iω (oscilaciones sin amortiguamiento). El sistema masa-resorte, mecánica orbital, dinámicas de péndulo.
Para λ = iω: hλ = ihω está en el eje imaginario.
Estabilidad de Euler en el eje imaginario: |1 + ihω|² = 1 + (hω)² > 1 para cualquier h > 0. Euler es inestable para cualquier tamaño de paso en valores propios puramente imaginarios. La 'oscilación' calculada crece sin límite.
Estabilidad de RK4 en el eje imaginario: la región de estabilidad se extiende a aproximadamente |hω| ≤ 2.83 en el eje imaginario. Para h suficientemente pequeño, RK4 maneja oscilaciones sin amortiguamiento. Euler no puede.
Esta geometría es por qué Euler falla en sistemas conservativos (masa-resorte, órbitas, ecuaciones de onda) incluso con h pequeño, mientras que RK4 los maneja bien.
La Geometría de Problemas Rígidos
Un sistema de EDO rígido tiene valores propios con magnitudes muy diferentes. La relación de rigidez: κ = max|Re(λᵢ)| / min|Re(λᵢ)| >> 1.
Por qué la rigidez es costosa para solucionadores explícitos:
La estabilidad requiere h·max|λᵢ| ≤ C (donde C depende del método). El valor propio más negativo establece el límite.
La precisión para la dinámica lenta requiere h·min|λᵢ| ≥ ε (resolver adecuadamente el modo más lento).
Si κ es grande, estos dos requisitos fuerzan un h diminuto: suficientemente pequeño para la estabilidad del modo rápido, lo suficientemente grande para muestrear el modo lento. El número de pasos escala con κ.
Imagen geométrica en el espectro de valores propios: los valores propios de la Jacobiana ∂f/∂y forman un conjunto de puntos en el plano complejo. La región de estabilidad de un solucionador explícito debe contener todos los puntos h·λᵢ. Si los valores propios van de -1 a -1000, la región de estabilidad debe cubrir un rango de 1000 a lo largo del eje real — requiriendo un h muy pequeño.
Solucionadores implícitos: la región de estabilidad de backward Euler cubre todo el semiplano izquierdo. Todos los valores propios con Re(λ) < 0 están automáticamente dentro de la región de estabilidad independientemente de h. La restricción en h proviene solo de la precisión, no de la estabilidad.
Relación de Rigidez & Costo
Considera una red de reacciones químicas con reacciones rápidas (escala de tiempo 10⁻⁶ s) y reacciones lentas (escala de tiempo 1 s).
Relación de rigidez: κ = 10⁶ / 1 = 10⁶.
Con RK4 (límite de estabilidad h·|λ_max| ≤ 2.785): h_max = 2.785 / 10⁶ ≈ 2.8 × 10⁻⁶ s.
Para integrar sobre 10 s de tiempo de reacción: pasos = 10 / (2.8 × 10⁻⁶) ≈ 3.6 × 10⁶.
Con backward Euler (incondicionalmente estable): h puede elegirse para la precisión de las reacciones lentas. h = 10⁻² s (100 muestras sobre escala de 1 s). Pasos = 10 / 10⁻² = 1000.
Relación de costo: explícito 3.6 millones de pasos vs implícito 1000 pasos — un factor de 3600. Cada paso implícito requiere resolver un sistema lineal (costo por paso es más alto), pero el costo total es mucho más bajo para problemas muy rígidos.
Por Qué los Tubos n-Dimensionales No Son Lo Que Piensas
En 2D, un 'tubo' de radio ε alrededor de una curva C es el conjunto de puntos dentro de distancia ε de C. La sección transversal es un círculo de radio ε. El volumen del tubo crece proporcionalmente a su longitud.
En n dimensiones, la geometría del tubo cambia fundamentalmente, debido al fenómeno del Capítulo 9:
La paradoja de la esquina n-dimensional: en el espacio n-dimensional, casi todo el volumen de un hipercubo n-dimensional está en las esquinas — no en la región central. A medida que n aumenta, la fracción de volumen dentro de distancia ε del centro tiende a cero para cualquier ε fijo.
Aplicado a tubos de solución de EDO:
En 2D: si la solución verdadera pasa por el centro de un tubo, la mayoría de los puntos cercanos están cerca de la curva. Las pequeñas perturbaciones te mantienen cerca de la solución verdadera.
En altas dimensiones: la mayoría de los puntos dentro de la caja delimitadora del tubo están en realidad lejos de la curva de solución verdadera. El 'volumen' del tubo está dominado por las esquinas — regiones que están lejos del centro en múltiples dimensiones simultáneamente.
Consecuencia para la simulación: con 28 EDOs acopladas (problema de intercepción de la Marina de Hamming), una perturbación de tamaño ε en cada dimensión puede producir un desplazamiento total de ε√28 ≈ 5.3ε de la solución verdadera. El tubo debe entenderse en términos de la norma L2 en todas las dimensiones, no solo el desplazamiento máximo en ninguna dimensión.
Estabilidad en altas dimensiones: un sistema donde cada componente decae independientemente (cada valor propio tiene parte real negativa) aún puede mostrar desplazamientos combinados grandes porque los errores de los componentes se suman en la norma L2. El tubo 28-dimensional no es solo 28 tubos 1-dimensionales independientes — la geometría los acopla.
De la Geometría al Diseño
Los conocimientos geométricos de los Capítulos 18-20 se unen como un conjunto de principios de diseño para la simulación numérica:
Selección de tamaño de paso: h debe colocar h·λ dentro de la región de estabilidad para cada valor propio. Para sistemas rígidos, los métodos implícitos eliminan la restricción de estabilidad, dejando solo requisitos de precisión.
Acumulación de errores en altas dimensiones: el error global es un vector en el espacio n-dimensional. Su norma crece como √n veces el error por componente. Las simulaciones de alta dimensionalidad necesitan requisitos de precisión por paso más ajustados.
Retroalimentación como estabilizador: si la simulación incorpora retroalimentación (la salida calculada influye en las entradas posteriores, como en un sistema de guía), la retroalimentación convergente amortigua los errores. La simulación puede tolerar entradas imprecisas para cantidades dentro del bucle de retroalimentación.
La inestabilidad como señal: para problemas con campos de dirección divergentes, la inestabilidad puede explotarse: la dirección de divergencia lleva información sobre el error de la condición inicial, permitiendo ajuste correctivo.