Rigidez y múltiples escalas de tiempo

Cuando la trayectoria parece lenta, pero el integrador exige pasos diminutos

En el capítulo 25 vimos que un integrador adaptativo puede reducir el paso cuando una solución cambia con mayor rapidez. La relación entre la velocidad de la dinámica y el tamaño del paso es intuitiva. Si una señal oscila cada vez más rápido, hace falta observarla con mayor resolución temporal para reproducir su fase y su amplitud.

La rigidez plantea un problema distinto. Un sistema puede contener un modo muy rápido que desaparece casi por completo después de un transitorio breve, mientras la trayectoria de interés sigue evolucionando lentamente. Aun así, un método explícito puede verse obligado a usar pasos comparables con la escala rápida durante todo el cálculo. En ese caso el paso queda limitado por la estabilidad numérica del método, aunque la resolución de la señal visible permitiría pasos mucho mayores.

La situación aparece en la cinética química y en la farmacocinética, entre muchos otros sistemas que combinan procesos con tiempos característicos muy distintos. El artículo clásico de Curtiss y Hirschfelder (1952) ya estudiaba ecuaciones de este tipo y contribuyó a fijar el uso moderno del término inglés stiff.

No existe, sin embargo, una definición escalar única de rigidez que sirva para todos los problemas. En sistemas lineales simples puede relacionarse con los valores propios y las escalas de relajación. En problemas no lineales o dependientes del tiempo, la rigidez varía a lo largo de la trayectoria y depende del método numérico elegido. Hairer y Wanner (1996) la tratan como una propiedad que se reconoce por sus consecuencias para la integración. Ninguna razón numérica aislada basta para definirla.

El capítulo estudia el mecanismo en un ejemplo completamente soluble y compara un método explícito con métodos implícitos. Con esa base precisaremos qué significa que un integrador sea adecuado para un problema rígido.

Un sistema con dos escalas temporales

Consideremos el sistema

\dot{x}=-x, \qquad \varepsilon \dot{y}=x-y, \qquad \varepsilon=10^{-3}, \qquad x(0)=1, \qquad y(0)=0.

La variable $x$ evoluciona en una escala temporal del orden de la unidad, mientras que $y$ tiende a acercarse a $x$ con una constante temporal $\varepsilon=0{,}001$ . En forma matricial,

\frac{d}{dt} \begin{pmatrix} x\\y \end{pmatrix} = A \begin{pmatrix} x\\y \end{pmatrix}, \qquad A= \begin{pmatrix} -1 & 0\\ 1000 & -1000 \end{pmatrix}.

Los valores propios son $\lambda_1=-1$ y $\lambda_2=-1000$ , con tiempos característicos de relajación $\tau_1=1$ y $\tau_2=10^{-3}$ . En este ejemplo lineal podemos definir la razón de escalas

\kappa=\frac{|\lambda_{\max}|}{|\lambda_{\min}|}=1000.

La cifra ayuda a reconocer una separación amplia de tiempos, aunque no basta como definición universal de rigidez. En particular, una gran razón entre frecuencias de un sistema oscilatorio puede tener consecuencias numéricas distintas de las de una gran razón entre tasas de decaimiento.

La solución exacta de la primera ecuación es $x(t)=e^{-t}$ . Al sustituirla en la segunda se obtiene

y(t)=\frac{e^{-t}-e^{-t/\varepsilon}}{1-\varepsilon}.

El término $e^{-t/\varepsilon}=e^{-1000t}$ desaparece con extrema rapidez. Cae por debajo de $10^{-8}$ antes de $t=0{,}019$ . Desde entonces la solución queda dominada por el modo lento (figura 26.1).

Figura 26.1. Solución exacta del sistema con dos escalas temporales. En el panel (a), después de la capa inicial la variable y sigue de cerca al modo lento. El panel (b) amplía el intervalo entre 0 y 0,01, donde y sube desde cero en unas pocas milésimas. El modo rápido seguirá condicionando a ciertos integradores explícitos aunque ya sea casi invisible en la trayectoria. Ampliar figura ↗

Separación de escalas y rigidez

La presencia de tiempos muy diferentes es una señal frecuente de rigidez, aunque ambos conceptos no deben identificarse sin condiciones adicionales. Imaginemos dos osciladores conservativos, uno con frecuencia 1 y otro con frecuencia 1000. Si hay que reproducir ambas oscilaciones, un paso del orden de $10^{-3}$ puede ser necesario simplemente para resolver la señal rápida. En ese caso el costo elevado corresponde a contenido dinámico que se quiere observar.

En nuestro sistema ocurre algo distinto. El modo con tasa 1000 es fuertemente disipativo. Después de un transitorio muy breve su amplitud es despreciable y la evolución relevante tiene una escala cercana a 1. Aun así, Euler explícito seguirá restringido por el modo rápido, porque su región de estabilidad no incluye pasos grandes para $\lambda=-1000$ .

Este contraste da una caracterización operacional útil. Un problema es rígido para un método dado cuando la estabilidad obliga a usar pasos mucho menores que los que exigiría la precisión de las componentes lentas de interés (Butcher, 2016; Hairer y Wanner, 1996). La formulación deja sin resolver algunas sutilezas de la teoría, pero identifica el conflicto numérico central.

La rigidez tampoco es una propiedad de las ecuaciones consideradas fuera de contexto, porque puede depender del intervalo de integración. Si solo estudiáramos el transitorio inicial, harían falta pasos pequeños por precisión aunque se usara un método implícito. Si integramos durante mucho tiempo después de que el transitorio desapareció, la restricción de estabilidad del método explícito se vuelve mucho más costosa en relación con la información que se busca.

Euler explícito y el modo que ya no vemos

Para la ecuación lineal $u'=\lambda u$ , Euler explícito produce $u_{n+1}=(1+h\lambda)u_n$ . La solución exacta decae cuando $\operatorname{Re}(\lambda)<0$ . Para que el modo numérico también disminuya en magnitud se necesita $|1+h\lambda|<1$ , lo que para $\lambda$ real y negativo equivale a

0<h<\frac{2}{|\lambda|}.

En nuestro sistema, el modo lento admite $h<2$ , mientras que el modo rápido exige $h<0{,}002$ .

Supongamos que elegimos $h=0{,}01$ . Para el modo lento, $1+h\lambda_1=0{,}99$ , un factor del todo razonable. Para el modo rápido, en cambio, $1+h\lambda_2=1-10=-9$ . Cada paso multiplica su amplitud por $-9$ . La solución del modelo elimina ese modo casi al instante, mientras que el algoritmo lo convierte en una oscilación explosiva. Desde $(1,0)$ , los primeros valores de $y$ son $10$ , $-80{,}1$ y $730{,}7$ . Después de diez pasos, $y$ vale unos $-3{,}5\times10^9$ , cuando la solución exacta en $t=0{,}1$ es cercana a $0{,}9$ (figura 26.2).

Figura 26.2. Euler implícito y Euler explícito con el mismo paso h = 0,01. En el panel (a), Euler implícito sigue la solución exacta de y después del primer paso. En el panel (b), Euler explícito es inestable porque el modo rápido queda fuera de su región de estabilidad, así que la variable alterna de signo y crece en magnitud. Las dos trayectorias son discretizaciones del mismo modelo continuo. Ampliar figura ↗

Un error de este tipo puede confundirse con un comportamiento del sistema si no se examina la estabilidad del integrador. La alternancia creciente es una trayectoria del sistema discreto que produce Euler con un paso inadecuado. El modelo diferencial no la predice.

El costo de un paso pequeño

Euler explícito puede usarse con un paso menor, por ejemplo $h=0{,}001$ . El factor rápido pasa a ser $1-1000h=0$ , de modo que en este caso particular el modo rápido desaparece numéricamente en un paso. Para integrar hasta $T=1$ hacen falta mil pasos. La aproximación obtenida da $x(1)\approx0{,}3676954$ , frente al valor exacto $e^{-1}\approx0{,}3678794$ , con un error absoluto de unos $1{,}84\times10^{-4}$ .

El procedimiento funciona, pero su costo es alto. El modo lento cambia en una escala de orden uno y, aun así, hay que avanzar en incrementos de una milésima para impedir que un modo ya extinguido genere una inestabilidad numérica. La dificultad es de eficiencia, porque la solución es perfectamente calculable. Con ese paso, además, el primer valor de $y$ es 1, mientras que la solución exacta en $t=0{,}001$ vale $0{,}632$ . El paso que asegura la estabilidad no basta para describir con precisión la capa inicial.

Si $\varepsilon$ bajara a $10^{-6}$ , la restricción sería tres órdenes de magnitud más severa. En sistemas grandes pueden coexistir muchas escalas de relajación. Reducir el paso puede entonces volver impracticable un cálculo conceptualmente sencillo.

Euler implícito y evaluación del campo en el estado futuro

Euler implícito modifica la regla y evalúa el campo al final del paso, $z_{n+1}=z_n+h f(t_{n+1},z_{n+1}).$

Para el sistema lineal $\dot z=Az$ , esa regla da $z_{n+1}=z_n+hAz_{n+1}$ , es decir, $(I-hA)z_{n+1}=z_n$ . Cada paso exige resolver un sistema lineal, $z_{n+1}=(I-hA)^{-1}z_n.$

En la ecuación de prueba $u'=\lambda u$ , el factor de amplificación es

R(z)=\frac{1}{1-z}, \qquad z=h\lambda,

y cumple $|R(z)|<1$ siempre que $\operatorname{Re}(z)<0$ . Euler implícito es, por tanto, A-estable, porque su región de estabilidad contiene todo el semiplano izquierdo. Dahlquist (1963) formalizó ese concepto en la teoría moderna de estabilidad de los métodos numéricos.

Para el modo rápido de nuestro ejemplo con $h=0{,}01$ , $z=-10$ y $R(-10)=1/11\approx0{,}0909$ . La solución exacta reduciría ese modo en cada intervalo por el factor $e^{-10}\approx4{,}54\times10^{-5}$ . Euler implícito reproduce con poca precisión la magnitud de esa relajación dentro del primer paso, aunque la amortigua con fuerza. Si el interés está en la dinámica lenta posterior, esa propiedad puede importar más que reproducir el transitorio rápido con miles de pasos.

A-estabilidad y L-estabilidad

La A-estabilidad no resuelve todos los problemas de rigidez. La regla trapezoidal lo muestra con claridad,

z_{n+1}=z_n+\frac{h}{2}\left[f(t_n,z_n)+f(t_{n+1},z_{n+1})\right].

Aplicada a $u'=\lambda u$ , su función de estabilidad es

R(z)=\frac{1+z/2}{1-z/2}.

También es A-estable, pero para $z=-10$ da $R(-10)=-2/3$ . El método no diverge, aunque transforma un modo que debería desaparecer casi por completo en una oscilación numérica que decae lentamente.

Para analizar este comportamiento se usa una condición más fuerte. Un método A-estable es L-estable si además $\lim_{z\to-\infty}R(z)=0$ (Iserles, 2008). Euler implícito cumple esa condición. La regla trapezoidal no la cumple, porque

\lim_{z\to-\infty} \frac{1+z/2}{1-z/2}=-1,

de modo que puede conservar modos extremadamente rápidos como oscilaciones alternantes.

Con $h\lambda=-10$ , reducir el modo rápido por debajo de una milésima de su amplitud inicial exige 18 pasos de la regla trapezoidal y 3 de Euler implícito (figura 26.3). La diferencia crece con la rigidez. Si $\varepsilon=10^{-6}$ y $h=0{,}01$ , el factor trapezoidal vale $-0{,}9996$ y el modo rápido conserva el 96 % de su amplitud después de cien pasos.

Figura 26.3. Amortiguación del modo rápido con hλ = −10. La solución exacta lo elimina casi por completo en un paso y Euler implícito lo amortigua en pocos pasos. La regla trapezoidal permanece estable, pero conserva durante mucho más tiempo una oscilación alternante. Ampliar figura ↗

La comparación deja una conclusión metodológica. Un método puede evitar el crecimiento numérico y, al mismo tiempo, describir mal la desaparición de las escalas rápidas, de modo que la estabilidad no asegura la fidelidad cualitativa.

Costo algebraico de los métodos implícitos

Los métodos explícitos calculan el estado siguiente directamente a partir de información disponible, mientras que los implícitos incluyen el estado futuro dentro de la ecuación que debe resolverse. Esa diferencia tiene un costo. En un sistema no lineal, la regla $z_{n+1}=z_n+h f(t_{n+1},z_{n+1})$ no puede despejarse, en general, mediante una inversión matricial fija. En cada paso hay que resolver

G(w)=w-z_n-hf(t_{n+1},w)=0.

Un esquema de Newton usa la matriz $DG(w)=I-hJ_f(t_{n+1},w)$ y actualiza mediante

w^{(m+1)}=w^{(m)}- \left[DG(w^{(m)})\right]^{-1}G(w^{(m)}).

Aparecen así nuevas decisiones computacionales. Hay que decidir cómo estimar o entregar el jacobiano y qué tolerancia usar en la resolución no lineal. También hay que decidir cuándo reutilizar factorizaciones y cómo declarar que la iteración no convergió (Ascher y Petzold, 1998).

Los métodos implícitos admiten pasos mucho mayores en muchos problemas rígidos, aunque cada paso cuesta más. Lo que importa comparar es el costo total para obtener la precisión y la estabilidad requeridas (Butcher, 2016; Hairer y Wanner, 1996).

En sistemas grandes importa la estructura del jacobiano. Una matriz dispersa puede permitir resolver sistemas lineales con mucha más eficiencia que una matriz densa del mismo tamaño. Conocer la estructura matemática del modelo puede, por eso, modificar sustancialmente el costo computacional sin cambiar las ecuaciones.

Una comparación cuantitativa

Integramos el sistema hasta $T=1$ . La solución exacta vale $x(1)\approx0{,}3678794$ e $y(1)\approx0{,}3682477$ . El cuadro 26.1 resume los resultados numéricos.

Cuadro 26.1. Errores en $T=1$ de tres métodos aplicados al sistema con dos escalas temporales.

Método	$h$	Pasos	Error en $x(1)$	Error en $y(1)$
Euler explícito	0,001	1000	$1{,}84\times10^{-4}$	$1{,}84\times10^{-4}$
Euler implícito	0,01	100	$1{,}83\times10^{-3}$	$1{,}83\times10^{-3}$
Regla trapezoidal	0,01	100	$3{,}07\times10^{-6}$	$3{,}07\times10^{-6}$

El cuadro no establece un ganador universal. La regla trapezoidal es muy precisa para el estado final de este ejemplo, aunque describe mal el transitorio rápido con $h=0{,}01$ . Euler implícito produce mayor error en el modo lento porque solo es de primer orden, pero elimina con eficacia el modo rápido. Euler explícito necesita diez veces más pasos para permanecer estable.

Evaluar solo el estado final puede ocultar diferencias importantes. Si el objetivo fuera estimar el transitorio durante las primeras $0{,}005$ unidades de tiempo, un paso de $0{,}01$ sería inaceptable incluso para un método estable. La rigidez no elimina la necesidad de resolver las escalas que sean científicamente relevantes.

Sistemas lento-rápidos y reducción del modelo

Nuestro ejemplo tiene la forma

\dot{x}=f(x,y), \qquad \varepsilon\dot{y}=g(x,y), \qquad 0<\varepsilon\ll1.

Cuando $\varepsilon$ tiende a cero, la segunda ecuación sugiere la restricción $g(x,y)=0$ . En nuestro caso esa restricción es $x-y=0$ , de modo que la aproximación de orden dominante es $y\approx x$ . Así se obtiene una descripción reducida centrada en la dinámica lenta. En $t=1$ , el modelo reducido da $y=e^{-1}\approx0{,}367879$ , mientras que la solución completa vale $0{,}368248$ . La diferencia, de unos $3{,}7\times10^{-4}$ , coincide con el orden de $\varepsilon e^{-1}$ que se espera de esa aproximación.

Reducir el modelo y elegir un integrador implícito son, sin embargo, operaciones conceptualmente distintas. Un método implícito busca aproximar las ecuaciones originales completas. Una reducción lento-rápida las reemplaza por otra representación, basada en una aproximación asintótica, que puede perder la capa inicial y otras propiedades ligadas a ella.

Tampoco basta que $\varepsilon$ sea pequeño para justificar cualquier eliminación, porque se requieren condiciones de estabilidad y de regularidad de la variedad lenta. Más adelante estudiaremos con mayor profundidad los sistemas lento-rápidos y las perturbaciones singulares. Aquí basta reconocer que la rigidez puede motivar dos estrategias distintas, que son cambiar el método numérico o cambiar, con justificación, el modelo efectivo.

Rigidez en sistemas no lineales

En un sistema no lineal $\dot{z}=f(t,z)$ , los valores propios relevantes pueden cambiar a lo largo de la trayectoria, porque el jacobiano $J_f(t,z)=\partial f/\partial z$ depende del estado y del tiempo. Una misma simulación puede atravesar regiones rígidas y no rígidas. Puede comenzar con un transitorio rápido y entrar luego en una fase lenta prolongada, en la que más tarde se activen nuevas escalas. Por eso una única razón espectral calculada en un punto no describe necesariamente todo el problema.

Los solvers modernos combinan la estimación del error y el control del paso con información sobre el jacobiano, cuando está disponible. Los métodos de Radau IIA y las fórmulas BDF son herramientas habituales para ecuaciones rígidas. La documentación de SciPy recomienda métodos explícitos de Runge-Kutta para problemas no rígidos y Radau o BDF cuando un método explícito necesita un número inusualmente grande de pasos (SciPy Developers, s. f.).

Esa recomendación es heurística y no prueba la rigidez. Un solver explícito puede necesitar muchos pasos por oscilaciones genuinamente rápidas o por tolerancias muy exigentes. El diagnóstico debe relacionarse con el modelo y con la razón por la que disminuye el paso.

Comparación de un solver explícito y uno para problemas rígidos

El bloque siguiente integra el sistema con solve_ivp mediante un método explícito y dos implícitos, a los que entrega el jacobiano exacto.

import numpy as np
from scipy.integrate import solve_ivp

eps = 1e-3


def f(t, z):
    x, y = z
    return np.array([-x, (x - y) / eps])


def jac(t, z):
    return np.array([[-1.0, 0.0],
                     [1.0 / eps, -1.0 / eps]])


z0 = np.array([1.0, 0.0])
exacta = np.array([np.exp(-1.0),
                   (np.exp(-1.0) - np.exp(-1.0 / eps)) / (1.0 - eps)])

for metodo in ("RK45", "Radau", "BDF"):
    opciones = {} if metodo == "RK45" else {"jac": jac}
    sol = solve_ivp(f, (0.0, 1.0), z0, method=metodo,
                    rtol=1e-6, atol=1e-9, **opciones)
    error = np.max(np.abs(sol.y[:, -1] - exacta))
    print(metodo, sol.nfev, sol.t.size - 1, f"{error:.1e}")

En modelos grandes y rígidos, entregar un jacobiano correcto o su patrón de dispersión puede reducir mucho el costo. Un jacobiano incorrecto, en cambio, puede ser peor que dejar que el software lo aproxime. La reproducibilidad exige registrar el nombre del solver junto con las tolerancias, la versión del software, el jacobiano usado, los límites de paso y cualquier opción adicional.

La salida del software tampoco reemplaza la verificación. En este ejemplo la solución exacta permite comprobar los errores directamente. En problemas sin solución cerrada hay que recurrir al refinamiento, a la comparación de métodos y a los invariantes o balances desarrollados en los capítulos anteriores.

Un diagnóstico computacional sencillo

La rigidez también puede hacerse visible en los contadores de un solver adaptativo. El cuadro 26.2 reúne los valores que imprime el bloque anterior con SciPy 1.17.1, junto con una segunda ejecución de RK45 con tolerancias mil veces más laxas.

Cuadro 26.2. Contadores de solve_ivp para el sistema con dos escalas temporales en $0\leq t\leq1$ . Los métodos implícitos recibieron el jacobiano exacto.

Método	`rtol`	`atol`	Evaluaciones de $f$	Pasos aceptados	Error máximo en $T=1$
RK45	$10^{-6}$	$10^{-9}$	2162	325	$3{,}3\times10^{-8}$
RK45	$10^{-3}$	$10^{-6}$	2078	309	$3{,}1\times10^{-4}$
Radau	$10^{-6}$	$10^{-9}$	421	59	$5{,}2\times10^{-10}$
BDF	$10^{-6}$	$10^{-9}$	232	115	$5{,}5\times10^{-7}$

La comparación entre las dos filas de RK45 es la más reveladora. Al relajar mil veces las tolerancias, los pasos solo bajan de 325 a 309, aunque el error final crece unas diez mil veces. Lo que limita el paso de RK45 es la estabilidad frente al modo rápido, porque la precisión pedida ya permitiría pasos mucho mayores. La región de estabilidad de la fórmula de Dormand y Prince corta el eje real negativo cerca de $-3{,}31$ . Con $|\lambda_2|=1000$ , el paso no puede superar unos $3{,}3\times10^{-3}$ y se necesitan al menos unos 300 pasos para llegar a $t=1$ .

Estas cifras no establecen una clasificación universal de solvers. El costo relativo de una evaluación del campo y de una factorización lineal cambia mucho de un problema a otro. En sistemas grandes, una sola factorización puede costar más que varias evaluaciones explícitas. También influyen las tolerancias, la estructura dispersa, la implementación y el equipo usado.

RK45 no fracasa ni diverge, porque obtiene una solución precisa. Lo que ocurre es que realiza muchas más evaluaciones, ya que su estabilidad restringe el paso. Los métodos implícitos invierten más trabajo en cada paso, pero pueden atravesar la región lenta sin seguir resolviendo artificialmente la escala $10^{-3}$ . El número de pasos o de evaluaciones indica rigidez solo cuando se interpreta junto con el espectro, las tolerancias y la dinámica esperada. Una ejecución lenta, considerada aisladamente, no demuestra rigidez.

Alcance del diagnóstico de rigidez

La rigidez es una propiedad numérica y dinámica del problema de integración. No implica que el sistema sea caótico ni que sus parámetros sean inciertos. Nuestro ejemplo es lineal y exactamente soluble, sin caos ni incertidumbre paramétrica. Aun así, es rígido para los métodos explícitos convencionales cuando $\varepsilon$ es pequeño. Del mismo modo, un sistema muy no lineal puede no ser rígido en el intervalo estudiado.

Tampoco debe inferirse que un modo rápido sea irrelevante porque desaparece pronto. Puede determinar condiciones iniciales efectivas o activar umbrales. Puede ser también, precisamente, el fenómeno de interés. La decisión de ignorarlo corresponde al problema científico y el solver no puede tomarla.

Un método implícito estable puede integrar sin dificultad una ecuación mal especificada, de modo que su capacidad para avanzar con pasos largos no aporta evidencia sobre la validez empírica del modelo.

Verificación práctica de un problema potencialmente rígido

Un diagnóstico empieza por preguntar si existen escalas de tiempo muy separadas. En un sistema lineal, el espectro da una primera aproximación. En uno no lineal hay que examinar el jacobiano a lo largo de la trayectoria. También importa si los modos rápidos son disipativos, porque una frecuencia oscilatoria alta y una tasa rápida de decaimiento imponen restricciones numéricas distintas.

Después hay que averiguar si el método explícito está limitado por la estabilidad o por la precisión. Si relajar la tolerancia apenas cambia el número de pasos, como ocurrió con RK45 en el cuadro 26.2, la rigidez es una explicación plausible. La comparación con un método implícito debe mostrar entonces si el costo baja sustancialmente sin cambiar los observables relevantes, con tolerancias y criterios declarados.

Queda por comprobar si la solución conserva los balances, la positividad y otras propiedades esperadas, porque la A-estabilidad no las garantiza. Si el transitorio rápido importa científicamente, tampoco pueden usarse pasos grandes solo porque el método sea estable. Estas preguntas convierten el diagnóstico de rigidez en una evaluación conjunta del cálculo y de la pregunta científica.

Consolidación

Las múltiples escalas temporales pueden producir rigidez, aunque la rigidez aparece solo cuando las restricciones de estabilidad numérica son mucho más severas que las que exige la dinámica lenta que se quiere resolver. En el ejemplo, un modo rápido de amplitud casi nula, con $\lambda=-1000$ , obliga a Euler explícito a usar $h<0{,}002$ durante una integración cuya trayectoria lenta evoluciona en una escala de orden uno.

La estabilidad y la precisión son propiedades distintas. Euler implícito admite $h=0{,}01$ sin divergir, aunque reproduce mal el transitorio rápido y solo tiene precisión de primer orden. La A-estabilidad y la L-estabilidad tampoco son equivalentes. La regla trapezoidal es A-estable, pero conserva durante muchos pasos los modos muy rígidos que Euler implícito amortigua con fuerza.

Los métodos implícitos trasladan parte del costo desde el número de pasos hacia la resolución algebraica de cada paso, que en problemas no lineales requiere Newton y jacobianos. Cambiar de solver es, además, una operación distinta de reducir el modelo. El integrador implícito aproxima las ecuaciones originales, mientras que una reducción lento-rápida modifica la representación matemática y necesita una justificación propia.

El ejemplo central del capítulo es lineal y exactamente soluble, lo que muestra que la rigidez no implica caos ni incertidumbre empírica. El capítulo 27 examinará la aritmética de precisión finita, una fuente de error que acompaña a todos los métodos.

Ejercicios

Ejercicio 26.1. Con $\varepsilon=10^{-3}$ , aplique Euler explícito desde $(1,0)$ . (a) Calcule los tres primeros valores de $y$ con $h=0{,}01$ . (b) Repita el cálculo con $h=0{,}001$ . (c) Compare el primer valor de cada caso con la solución exacta en el tiempo correspondiente. ¿Qué muestra el caso (b) sobre la diferencia entre estabilidad y precisión?

Ejercicio 26.2. (a) Muestre que $R(z)=1/(1-z)$ cumple $|R(z)|<1$ cuando $\operatorname{Re}z<0$ . (b) Muestre que la función de estabilidad de la regla trapezoidal cumple $|R(z)|<1$ cuando $\operatorname{Re}z<0$ y $|R(is)|=1$ para todo $s$ real. (c) Calcule el límite de ambas funciones cuando $z\to-\infty$ e interprete el resultado.

Ejercicio 26.3. Suponga $\varepsilon=10^{-6}$ y $h=0{,}01$ . (a) Calcule el factor del modo rápido para la regla trapezoidal y para Euler implícito. (b) ¿Qué fracción de la amplitud del modo rápido queda después de cien pasos con cada método? (c) ¿Qué propiedad explica la diferencia?

Ejercicio 26.4. Un paso de Euler implícito con $h=0{,}5$ para $u'=-u^3$ desde $u_0=1$ exige resolver $w+0{,}5w^3=1$ . (a) Aplique tres iteraciones de Newton desde $w^{(0)}=1$ . (b) ¿Cómo disminuye el residuo $G(w)$ entre iteraciones? (c) Compare el resultado con un paso de Euler explícito y con la solución exacta $u(t)=1/\sqrt{1+2t}$ .

Ejercicio 26.5. (a) Calcule en qué instante el término $e^{-t/\varepsilon}$ cae por debajo de $10^{-8}$ con $\varepsilon=10^{-3}$ . (b) Compare en $t=1$ el valor de $y$ del modelo reducido con el de la solución completa. (c) ¿Qué parte de la solución pierde el modelo reducido?

Ejercicio 26.6. La fórmula de quinto orden de Dormand y Prince tiene la función de estabilidad $R(z)=1+z+z^2/2+z^3/6+z^4/24+z^5/120+z^6/600$ . (a) Calcule con un programa dónde corta $|R(x)|=1$ el eje real negativo. (b) Deduzca el paso máximo estable de RK45 para el sistema del capítulo y el número mínimo de pasos hasta $t=1$ . (c) ¿Qué ocurriría con ese número si $\varepsilon$ fuera $10^{-4}$ ? Compruébelo con el bloque del capítulo.

Consultar las soluciones de los ejercicios

Soluciones de los ejercicios

Solución 26.1. (a) Con $h=0{,}01$ la regla es $y_{n+1}=-9y_n+10x_n$ con $x_n=0{,}99^n$ , lo que da $10$ , $-80{,}1$ y $730{,}701$ . (b) Con $h=0{,}001$ la regla es $y_{n+1}=x_n$ , lo que da $1$ , $0{,}999$ y $0{,}998001$ . (c) La solución exacta vale $y(0{,}01)\approx0{,}991$ e $y(0{,}001)\approx0{,}632$ . El caso (b) es estable, pero su primer paso sobrepasa la capa inicial y da 1 donde la solución vale $0{,}632$ . La estabilidad no asegura la precisión del transitorio.

Solución 26.2. (a) Si $\operatorname{Re}z<0$ , entonces $|1-z|^2=(1-\operatorname{Re}z)^2+(\operatorname{Im}z)^2>1$ , así que $|R(z)|=1/|1-z|<1$ . (b) Se cumple $|1+z/2|^2-|1-z/2|^2=2\operatorname{Re}z$ , que es negativo cuando $\operatorname{Re}z<0$ . Sobre el eje imaginario la diferencia es nula y los dos módulos coinciden. (c) El límite vale 0 para Euler implícito y $-1$ para la regla trapezoidal. Euler implícito es L-estable. La regla trapezoidal conserva los modos muy rígidos como oscilaciones que alternan de signo.

Solución 26.3. (a) Con $z=-10^4$ , la regla trapezoidal da $R=-4999/5001\approx-0{,}9996$ y Euler implícito da $R=1/10001\approx10^{-4}$ . (b) Con la regla trapezoidal queda $0{,}9996^{100}\approx0{,}961$ , es decir, el 96 % de la amplitud. Con Euler implícito queda $10^{-400}$ , un valor que en doble precisión se redondea a cero. (c) La diferencia se debe a la L-estabilidad, que Euler implícito tiene y la regla trapezoidal no.

Solución 26.4. (a) Con $G(w)=w-1+0{,}5w^3$ y $G'(w)=1+1{,}5w^2$ , las iteraciones dan $0{,}8$ , $0{,}771429$ y $0{,}770917$ . (b) Los residuos valen $0{,}5$ , $0{,}056$ , $9{,}7\times10^{-4}$ y $3{,}0\times10^{-7}$ . El número de cifras correctas se duplica aproximadamente en cada iteración, como corresponde a la convergencia cuadrática de Newton. (c) Euler explícito da $u_1=0{,}5$ . La solución exacta en $t=0{,}5$ es $1/\sqrt2\approx0{,}707107$ , de modo que el error de Euler implícito es $0{,}064$ y el de Euler explícito, $0{,}207$ .

Solución 26.5. (a) El término cae por debajo de $10^{-8}$ cuando $t>\varepsilon\ln10^8\approx0{,}0184$ . (b) El modelo reducido da $y(1)=e^{-1}\approx0{,}367879$ y la solución completa, $0{,}368248$ . La diferencia es de $3{,}68\times10^{-4}$ , casi igual a $\varepsilon e^{-1}$ . (c) El modelo reducido pierde la capa inicial. Impone $y=x$ desde $t=0$ , lo que da $y(0)=1$ aunque la condición inicial es $y(0)=0$ .

Solución 26.6. (a) El corte está en $x\approx-3{,}307$ . (b) El paso máximo estable es $3{,}307/1000\approx3{,}3\times10^{-3}$ , de modo que llegar a $t=1$ exige al menos unos 303 pasos. El resultado explica por qué RK45 necesitó entre 309 y 325 pasos con tolerancias muy distintas. (c) Con $\varepsilon=10^{-4}$ , el paso máximo se divide por diez y se necesitan al menos unos 3000 pasos. El bloque del capítulo, con esa modificación, da 3043 pasos aceptados.

Referencias

Ascher, U. M. y Petzold, L. R. (1998). Computer methods for ordinary differential equations and differential-algebraic equations. Society for Industrial and Applied Mathematics.

Butcher, J. C. (2016). Numerical methods for ordinary differential equations (3.ª ed.). Wiley. https://doi.org/10.1002/9781119121534

Curtiss, C. F. y Hirschfelder, J. O. (1952). Integration of stiff equations. Proceedings of the National Academy of Sciences of the United States of America, 38(3), 235–243. https://doi.org/10.1073/pnas.38.3.235

Dahlquist, G. G. (1963). A special stability problem for linear multistep methods. BIT Numerical Mathematics, 3(1), 27–43. https://doi.org/10.1007/BF01963532

Hairer, E. y Wanner, G. (1996). Solving ordinary differential equations II: Stiff and differential-algebraic problems (2.ª ed. rev.). Springer. https://doi.org/10.1007/978-3-642-05221-7

Iserles, A. (2008). A first course in the numerical analysis of differential equations (2.ª ed.). Cambridge University Press.

SciPy Developers. (s. f.). scipy.integrate.solve_ivp. Recuperado el 19 de septiembre de 2026, de https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html