Runge-Kutta, tolerancias y paso adaptativo

Por qué una sola pendiente puede ser insuficiente

En el capítulo 24 aproximamos una ecuación diferencial mediante el método de Euler explícito. La regla $y_{n+1}=y_n+h f(t_n,y_n)$ usa una sola evaluación del campo, la pendiente al comienzo del intervalo. Esa simplicidad es su principal ventaja pedagógica y también su limitación, porque durante un paso finito la trayectoria verdadera puede curvarse y su pendiente cambiar de manera apreciable.

Una respuesta natural consiste en evaluar el campo en varios puntos dentro del mismo paso. Los métodos de Runge-Kutta desarrollan esa idea de forma sistemática. Combinan varias pendientes calculadas en estados intermedios que el propio método construye, sin reconstruir la solución exacta entre $t_n$ y $t_{n+1}$ . Los coeficientes de esas combinaciones se eligen para que el desarrollo de la actualización coincida con el de Taylor de la solución hasta un orden determinado (Butcher, 2016; Hairer et al., 1993).

Así se obtiene alta precisión sin calcular derivadas de orden superior de $f$ . Un método de Taylor de cuarto orden necesitaría derivadas sucesivas de la solución, mientras que el Runge-Kutta clásico de cuarto orden usa cuatro evaluaciones de $f$ y alcanza un error global $O(h^4)$ bajo condiciones regulares.

Un programa que entrega una trayectoria sin más explicación suele ocultar decisiones que es necesario separar. Hay que saber qué integrador se usó y cómo se estimó el error numérico, además de cómo se decidió el tamaño del paso. Todas esas decisiones pertenecen al experimento computacional y ninguna es una propiedad empírica del sistema modelado.

La estructura general de un método de Runge-Kutta

Consideremos el problema

\frac{dy}{dt}=f(t,y), \qquad y(t_0)=y_0.

Un método explícito de Runge-Kutta con $s$ etapas calcula pendientes intermedias

k_i=f\left(t_n+c_i h,\ y_n+h\sum_{j<i}a_{ij}k_j\right)

y actualiza mediante

y_{n+1}=y_n+h\sum_{i=1}^{s}b_i k_i.

Los números $a_{ij}$ , $b_i$ y $c_i$ definen el método y suelen registrarse en una tabla de Butcher. Esta notación permite entender distintos integradores como elecciones dentro de una misma arquitectura formal (Butcher, 2016).

El método de Euler es el Runge-Kutta explícito más sencillo, con una etapa, $c_1=0$ y $b_1=1$ . La familia incluye, por tanto, a Euler como su caso de primer orden.

Un ejemplo intermedio es el método del punto medio explícito,

k_1=f(t_n,y_n), \qquad k_2=f\!\left(t_n+\frac{h}{2},\ y_n+\frac{h}{2}k_1\right), \qquad y_{n+1}=y_n+h k_2.

La pendiente $k_2$ se evalúa en un estado estimado para la mitad del intervalo. Bajo hipótesis de suavidad el error global es de segundo orden, de modo que dividir el paso por dos debería reducir el error, asintóticamente, por un factor cercano a cuatro.

El Runge-Kutta clásico de cuarto orden, abreviado RK4, emplea cuatro etapas,

k_1=f(t_n,y_n), \qquad k_2=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1\right),

k_3=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2\right), \qquad k_4=f(t_n+h,y_n+h k_3),

y actualiza con

\boxed{ y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4). }

El cuadro 25.1 muestra su tabla de Butcher.

Cuadro 25.1. Tabla de Butcher del método RK4. La última fila contiene los pesos $b_i$ .

$c_i$	$a_{i1}$	$a_{i2}$	$a_{i3}$	$a_{i4}$
$0$
$1/2$	$1/2$
$1/2$	$0$	$1/2$
$1$	$0$	$0$	$1$
	$1/6$	$1/3$	$1/3$	$1/6$

Estos pesos no son un simple promedio empírico de cuatro pendientes. Las etapas están acopladas, porque $k_2$ , $k_3$ y $k_4$ se evalúan sobre estados intermedios producidos por las pendientes anteriores. Las condiciones de orden garantizan que, para campos suficientemente suaves, el desarrollo resultante coincide con el de Taylor hasta el orden buscado (Butcher, 2016; Iserles, 2008).

El defecto local de RK4 es $O(h^5)$ y su error global es $O(h^4)$ , bajo las hipótesis usuales de regularidad y estabilidad que impiden que los errores locales se amplifiquen sin control.

Un oscilador como prueba de integración

Usemos el oscilador armónico

\dot{x}=v, \qquad \dot{v}=-x, \qquad x(0)=1, \qquad v(0)=0,

cuya solución exacta es $x(t)=\cos t$ , $v(t)=-\sin t$ . La cantidad

E=\frac{1}{2}(x^2+v^2)

permanece exactamente en $1/2$ , de modo que en el plano $(x,v)$ la trayectoria es una circunferencia de radio uno. El ejemplo permite evaluar algo más que el error en un punto final, porque también muestra si el integrador respeta la geometría de la órbita.

Para Euler explícito, la actualización es $z_{n+1}=M_Ez_n$ con

M_E=\begin{pmatrix}1 & h\\-h & 1\end{pmatrix}, \qquad M_E^{\mathsf T}M_E=(1+h^2)I.

Por tanto,

x_{n+1}^2+v_{n+1}^2=(1+h^2)(x_n^2+v_n^2), \qquad E_{n+1}=(1+h^2)E_n.

Euler aumenta el cuadrado del radio en cada paso. La espiral hacia afuera se deriva exactamente de la regla numérica y nada tiene que ver con el redondeo. En este problema, Euler explícito no reproduce la estabilidad neutral del oscilador con ningún $h>0$ .

RK4 es mucho más preciso con pasos pequeños, aunque tampoco conserva exactamente la energía del oscilador. Con $h=0{,}1$ y hasta $T=20$ , la energía que calcula RK4 es de unos $0{,}4999986$ , mientras que Euler llega a unos $3{,}6580$ y la solución exacta conserva $0{,}5$ .

Órbitas del oscilador obtenidas por solución exacta, Euler y Runge–Kutta. — Figura 25.1. Trayectoria exacta del oscilador y aproximaciones de Euler y de RK4 en el plano de fases, con h = 0,1 hasta T = 20. Euler genera una espiral hacia afuera. RK4 permanece visualmente sobre la órbita exacta, aunque no conserva exactamente la energía. Ampliar figura ↗

La comparación no implica que RK4 sea siempre la mejor elección para osciladores de larga duración. Existen integradores geométricos y simplécticos diseñados para preservar mejor ciertas estructuras (Hairer et al., 1993). Aquí usamos RK4 porque permite estudiar con claridad el orden y el costo por paso antes de introducir esas familias especializadas.

Función de estabilidad y comportamiento del oscilador

La diferencia entre Euler y RK4 también se ve a través de la función de estabilidad. Para la ecuación lineal de prueba $y'=\lambda y$ , un método de un paso produce $y_{n+1}=R(h\lambda)y_n$ , donde $R(z)$ resume cómo el integrador amplifica un modo lineal. Para Euler explícito, $R_E(z)=1+z$ . Para RK4,

R_{RK4}(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},

que coincide con el polinomio de Taylor de $e^z$ hasta cuarto orden (LeVeque, 2007).

En un problema disipativo, con $\operatorname{Re}(\lambda)<0$ , interesa que $|R(h\lambda)|$ no produzca crecimiento numérico espurio. En un oscilador lineal los valores propios son puramente imaginarios, iguales a $\pm i$ . La solución exacta multiplica cada modo por $e^{\pm ih}$ , de módulo uno. Euler produce

|R_E(ih)|=|1+ih|=\sqrt{1+h^2}>1

para todo $h>0$ , lo que reproduce algebraicamente la espiral anterior. Para RK4 se obtiene

|R_{RK4}(ih)|^2=1-\frac{h^6}{72}+\frac{h^8}{576},

que es menor que uno cuando $0<h<2\sqrt2\approx2{,}83$ . En ese intervalo RK4 amortigua levemente el oscilador. Con $h=0{,}1$ , cada paso multiplica la energía por un factor que difiere de uno en unos $1{,}39\times10^{-8}$ . Doscientos pasos dan el valor $0{,}4999986$ observado.

Esta formulación evita atribuir la diferencia a una superioridad abstracta de un método, porque ambos son polinomios distintos aplicados al mismo operador. La pregunta relevante es qué región del plano complejo y qué propiedades del problema debe reproducir fielmente el método. En el capítulo 26 esta misma perspectiva explicará por qué un método explícito puede exigir pasos extremadamente pequeños cuando coexisten escalas de relajación muy distintas.

La función de estabilidad analiza modos lineales. En un sistema no lineal sigue siendo útil para estudiar linealizaciones y restricciones numéricas locales, aunque no reemplaza el análisis completo de la dinámica.

El orden debe observarse mediante refinamiento

Una trayectoria visualmente razonable no demuestra que el método converja al ritmo esperado. Repitamos el oscilador hasta $T=4$ con distintos pasos y calculemos

e(T)=\|y_{\mathrm{num}}(T)-y_{\mathrm{exacta}}(T)\|_2.

El cuadro 25.2 muestra los resultados.

Cuadro 25.2. Error en $T=4$ para el oscilador armónico con Euler y con RK4.

$h$	Error de Euler	Error de RK4
0,400	1,135919	$8{,}5187\times10^{-4}$
0,200	0,484407	$5{,}3312\times10^{-5}$
0,100	0,220676	$3{,}3330\times10^{-6}$
0,050	0,105091	$2{,}0833\times10^{-7}$
0,025	0,051262	$1{,}3021\times10^{-8}$

En el régimen más fino, los errores de Euler se acercan al comportamiento proporcional a $h$ y los de RK4 disminuyen aproximadamente como $h^4$ . Un ajuste log-log sobre los cinco puntos da pendientes de $1{,}11$ para Euler y $4{,}00$ para RK4. La desviación de Euler recuerda que el orden es una propiedad asintótica, porque los pasos relativamente grandes pueden quedar fuera del régimen donde la potencia dominante describe bien el error.

Figura 25.2. Error global en T = 4 frente al tamaño del paso para Euler y RK4. Las líneas de referencia son proporcionales a h y a h⁴. Ambos métodos resuelven el mismo problema y se evalúan con la misma norma. Ampliar figura ↗

Un paso de RK4 exige cuatro evaluaciones de $f$ , frente a una de Euler, así que su mayor precisión no es gratuita. Aun así, un orden mayor puede reducir mucho el número de pasos necesarios para alcanzar un error dado. RK4 logra un error de $2\times10^{-7}$ con $h=0{,}05$ , es decir, con 320 evaluaciones. Si el error de Euler siguiera siendo proporcional a $h$ , alcanzar ese error exigiría un paso cercano a $10^{-7}$ y unos cuarenta millones de evaluaciones, sin contar la acumulación del redondeo. La eficiencia se evalúa, por tanto, relacionando el costo total con la precisión obtenida.

Por qué un paso fijo puede ser ineficiente

Supongamos ahora que la velocidad angular cambia con el tiempo,

\omega(t)=1+0{,}1t, \qquad \dot{x}=-\omega(t)y, \qquad \dot{y}=\omega(t)x,

con $(x(0),y(0))=(1,0)$ . Si

\phi(t)=\int_0^t \omega(s)\,ds=t+0{,}05t^2,

la solución exacta es $x(t)=\cos\phi(t)$ , $y(t)=\sin\phi(t)$ .

La trayectoria sigue siendo una circunferencia de radio uno, pero se recorre cada vez más rápido. Un paso fijo elegido para resolver bien el final del intervalo puede ser innecesariamente pequeño al comienzo, mientras que un paso adecuado para el comienzo puede resultar demasiado grande cuando la frecuencia aumenta. Un integrador adaptativo intenta resolver ese problema con pasos distintos según la dificultad local de la trayectoria, bajo un criterio de error explícito.

La adaptación no requiere que el algoritmo «comprenda» el sistema. Le basta una estimación numérica del error de cada intento de paso y una regla para aumentar o reducir $h$ .

Pares embebidos y estimación del error sin solución exacta

En un problema real rara vez se dispone de la solución exacta para medir el error. Una estrategia eficiente consiste en construir dos aproximaciones de distinto orden que compartan la mayor parte de las evaluaciones del campo. Si $y^{(p)}_{n+1}$ es una fórmula de orden alto y $y^{(q)}_{n+1}$ otra de orden menor, la diferencia

e^{\mathrm{est}}_{n+1}=y^{(p)}_{n+1}-y^{(q)}_{n+1}

estima el error local en la escala del par. En un par de órdenes consecutivos, esa diferencia está dominada por el error de la fórmula de menor orden.

Dormand y Prince (1980) desarrollaron una familia embebida 5(4) que reutiliza etapas para producir una solución de quinto orden junto con una estimación de error asociada a una fórmula de cuarto orden. La implementación RK45 de SciPy controla el error suponiendo la precisión de la fórmula de cuarto orden y avanza con el estado de quinto orden (Shampine, 1986; SciPy Developers, s. f.-a).

No hace falta memorizar los coeficientes de las siete etapas del par de Dormand y Prince. Basta entender su función, que es obtener a la vez la aproximación que se usa como solución y otra que permite decidir si el paso propuesto fue suficientemente preciso.

Tolerancia absoluta y relativa

Un error de $10^{-4}$ no significa lo mismo si una variable vale $10^3$ que si vale $10^{-6}$ . Por eso los integradores adaptativos normalizan el estimador con una combinación de tolerancia absoluta y relativa. Una escala habitual para la componente $i$ es

s_i=\mathrm{atol}_i+\mathrm{rtol}\,|y_i|,

o una variante que usa las magnitudes del estado antes y después del paso. La documentación de RK45 expresa el objetivo local con una escala de la forma atol + rtol*abs(y) y admite tolerancias absolutas distintas por componente (SciPy Developers, s. f.-a). Con esa escala se construye un error normalizado, por ejemplo mediante una norma RMS,

E= \sqrt{\frac{1}{d}\sum_{i=1}^{d} \left(\frac{e^{\mathrm{est}}_i}{s_i}\right)^2}.

El criterio habitual acepta el paso si $E\leq1$ . Si $E>1$ , el paso se rechaza y se intenta de nuevo con un $h$ menor. Si $E$ queda muy por debajo de uno, el paso siguiente puede aumentar.

La tolerancia absoluta impide exigir una precisión relativa absurda cerca de cero, mientras que la relativa adapta la escala permitida a la magnitud de la variable. Si un sistema contiene componentes con unidades o magnitudes muy distintas, una única atol escalar puede asignar prioridades poco apropiadas. En ese caso hay que justificar tolerancias por componente.

Estas tolerancias son parámetros del algoritmo y nada dicen sobre la incertidumbre de las mediciones o de los parámetros del modelo.

La elección de unidades también importa. Si una variable pasa de metros a milímetros, su valor numérico se multiplica por mil. Una tolerancia absoluta que no se transforme de manera coherente cambia entonces el criterio computacional, aunque el sistema físico sea el mismo. La tolerancia relativa es menos sensible a un cambio de unidades cuando la variable está lejos de cero, pero su combinación con atol sigue requiriendo atención.

En sistemas de varias variables, la norma que combina las componentes también es una elección. Una norma RMS reparte la contribución entre dimensiones y una norma del máximo exige que el peor componente cumpla el criterio. Ninguna es correcta en general. Si una variable decisiva, clínica o físicamente, ocupa una escala numérica pequeña, puede necesitar su propia tolerancia absoluta para que el resto del vector no la subordine.

Aquí reaparece una idea de los capítulos 11 y 16, según la cual una operación numérica puede depender de las coordenadas y de la escala elegidas. Para que la conclusión sea interpretable hay que explicar cómo esas decisiones corresponden a las magnitudes del problema.

Cómo se modifica el paso

Supongamos que el estimador se comporta localmente como $E\propto h^{q+1}$ , donde $q$ es el orden de la fórmula cuyo error controla el estimador. Para obtener aproximadamente $E_{\mathrm{nuevo}}=1$ bastaría con $h_{\mathrm{nuevo}}=hE^{-1/(q+1)}$ . En la práctica se multiplica por un factor de seguridad menor que uno y se limitan los cambios máximos de tamaño,

\boxed{ h_{\mathrm{nuevo}} =s\,h E^{-1/(q+1)}, \qquad 0<s<1. }

En el código de SciPy 1.17.1, RK45 usa $q=4$ y $s=0{,}9$ . El cambio de paso queda limitado, además, a un factor entre $0{,}2$ y $10$ .

Los controladores reales pueden incorporar información de pasos anteriores y otras salvaguardas. El exponente concreto depende del estimador y del método, de modo que la fórmula anterior explica la lógica del control sin describir cada implementación. Un paso adaptativo responde con una regla explícita a su estimación local del error, sin que eso lo haga óptimo en algún sentido particular.

El oscilador de frecuencia creciente

Integramos el sistema de frecuencia creciente, desde $t=0$ hasta $t=40$ , con el método RK45 de Dormand y Prince. Como conocemos la solución exacta, podemos evaluar tanto el comportamiento del controlador como el error global resultante. El cuadro 25.3 resume tres configuraciones calculadas con SciPy 1.17.1. El paso mediano de cada intervalo se calculó con los pasos aceptados que comienzan en él.

Cuadro 25.3. Integración con RK45 del oscilador de frecuencia creciente hasta $t=40$ .

`rtol`	`atol`	Pasos aceptados	Evaluaciones de $f$	Error global máximo	Paso mediano en $1<t<5$	Paso mediano en $35<t<39$
$10^{-3}$	$10^{-6}$	126	764	$3{,}955\times10^{-2}$	0,7310	0,2072
$10^{-6}$	$10^{-9}$	585	4262	$1{,}602\times10^{-5}$	0,1586	0,04377
$10^{-9}$	$10^{-12}$	2405	15092	$1{,}455\times10^{-8}$	0,03955	0,01107

En las tres configuraciones, el paso mediano del final es entre 3,5 y 3,6 veces menor que el del comienzo. Entre esos mismos intervalos, la frecuencia angular media pasa de $1{,}3$ a $4{,}7$ , un cociente de $3{,}6$ . El controlador no conoce esa fórmula, pero responde al error estimado de las etapas, que crece con la frecuencia.

Figura 25.3. Tamaños de paso aceptados por RK45 en el oscilador de frecuencia creciente con dos tolerancias. Al aumentar la frecuencia, el integrador reduce el paso. La tolerancia más estricta desplaza toda la secuencia hacia pasos menores. Ampliar figura ↗

El costo aumenta mucho al exigir más precisión. Entre $\mathrm{rtol}=10^{-3}$ y $\mathrm{rtol}=10^{-9}$ , las evaluaciones del campo pasan de 764 a 15092. El error global máximo disminuye varios órdenes de magnitud, pero la relación entre rtol y el error final dista de ser una identidad simple (figura 25.4).

Figura 25.4. Evaluaciones del campo frente al error global máximo observado para tres tolerancias. La figura describe este problema y este integrador, sin valor general como curva de costo de RK45. Ampliar figura ↗

Una tolerancia local y el error global

El primer caso del cuadro 25.3 merece atención. Con $\mathrm{rtol}=10^{-3}$ , el error global máximo de estado fue de unos $3{,}96\times10^{-2}$ . No hay contradicción, porque el controlador mantiene acotado un estimador local normalizado. Nada promete sobre la distancia global a la solución exacta después de cuarenta unidades de tiempo.

En un oscilador, los pequeños errores de fase se acumulan durante muchos ciclos. Distintos observables, por otra parte, pueden mostrar errores muy diferentes. En la misma simulación con $\mathrm{rtol}=10^{-3}$ , la máxima desviación del radio respecto de uno fue de unos $7{,}85\times10^{-4}$ , cincuenta veces menor que el error de estado. Una evaluación basada solo en la conservación del radio parecería excelente, aunque la fase ya se hubiera desviado de manera apreciable.

Por eso declarar una tolerancia no reemplaza la pregunta científica sobre qué propiedad necesita precisión. Un estudio puede requerir exactitud de fase, el tiempo de cruce de un umbral, la amplitud, la conservación de masa o alguna función derivada. La tolerancia del integrador controla el error según su propia norma, sin controlar cada una de esas cantidades por separado.

Tampoco garantiza la conservación exacta de invariantes ni la ausencia de sesgos sistemáticos a largo plazo. No asegura la precisión de los eventos detectados entre pasos ni la estabilidad frente a la rigidez. La fidelidad del modelo al fenómeno empírico depende de la ecuación elegida, que el integrador resuelve tal como se le entrega.

Verificación práctica de un integrador

Una estrategia reproducible exige más que elegir una tolerancia pequeña. Los problemas con solución conocida, como el oscilador armónico, permiten comparar directamente errores de estado y órdenes de convergencia. El refinamiento añade otra prueba. En un método de paso fijo, reducir $h$ debería acercar los resultados de acuerdo con su orden, siempre dentro del régimen asintótico. En un método adaptativo, reducir rtol y atol debería estabilizar los observables de interés.

Los invariantes y las identidades aportan pruebas estructurales. Si el modelo conserva masa, probabilidad u otra cantidad, su deriva puede vigilarse. El algoritmo no siempre la conserva exactamente, pero una violación grande indica un problema. Comparar integradores de familias distintas ayuda, además, a identificar resultados que dependen de una implementación concreta. La documentación debe registrar el método, la versión del software, las tolerancias, los límites de paso, el tratamiento de eventos y cualquier ajuste no predeterminado.

Estas comprobaciones no validan la ecuación como explicación empírica. Evalúan si el resultado computacional reproduce con suficiente fiabilidad la ecuación que decidimos estudiar.

Implementación mínima

Una implementación de RK4 con paso fijo puede escribirse de manera compacta.

def rk4_step(f, t, y, h):
    k1 = f(t, y)
    k2 = f(t + h/2, y + h*k1/2)
    k3 = f(t + h/2, y + h*k2/2)
    k4 = f(t + h, y + h*k3)
    return y + h*(k1 + 2*k2 + 2*k3 + k4)/6

El código reproduce literalmente las cuatro etapas matemáticas. Para un algoritmo adaptativo es preferible usar una implementación probada antes que reconstruir todos sus mecanismos de control. En SciPy, el oscilador de frecuencia creciente se integra así.

from scipy.integrate import solve_ivp


def campo(t, z):
    omega = 1.0 + 0.1 * t
    return [-omega * z[1], omega * z[0]]


sol = solve_ivp(campo, (0.0, 40.0), [1.0, 0.0], method="RK45",
                rtol=1e-6, atol=1e-9)
print(sol.t.size - 1, sol.nfev)

La llamada imprime 585 4262, los pasos aceptados y las evaluaciones del cuadro 25.3 (SciPy Developers, s. f.-b). Su brevedad no elimina las decisiones científicas. Las tolerancias, el método y el horizonte forman parte de la especificación computacional. Si el problema contiene discontinuidades, eventos o escalas temporales muy diferentes, también hay que declarar cómo se tratan. El capítulo 26 estudiará precisamente el caso en que reducir el paso de un método explícito deja de ser una solución razonable, porque el sistema es rígido.

Qué significa comparar integradores

Decir que un método es «mejor» sin especificar la tarea es insuficiente. Una comparación puede atender al error para un costo dado o al costo para una precisión dada. También puede considerar la estabilidad numérica, la conservación de invariantes, el tratamiento de eventos, la robustez frente a escalas múltiples o la facilidad para reproducir y auditar el cálculo.

RK4 supera ampliamente a Euler en el oscilador de este capítulo si el criterio es el error de estado con un paso fijo comparable. Esa conclusión no autoriza a afirmar que RK4 sea superior en cualquier ecuación. El capítulo siguiente mostrará un caso donde la estabilidad de un método explícito impone pasos muy pequeños aunque la trayectoria de interés evolucione lentamente.

La elección del integrador forma parte del modelo computacional. Las ecuaciones especifican una dinámica matemática y el método numérico especifica cómo se obtiene una aproximación finita de ella. Confundir ambas capas dificulta distinguir un comportamiento del sistema de un artefacto del cálculo.

Consolidación

RK4 usa varias pendientes dentro de cada paso, con coeficientes que cumplen condiciones de orden y producen un error global de cuarto orden bajo hipótesis regulares. Ese orden cuesta cuatro evaluaciones del campo por paso, de modo que su ventaja depende de cuánto reduce el número de pasos necesarios para la precisión buscada.

Un par embebido estima el error sin conocer la solución exacta, porque dos fórmulas de distinto orden comparten etapas y su diferencia controla el paso. Las tolerancias rtol y atol definen una escala de error numérico local. No son incertidumbres estadísticas ni cotas directas del error global.

El paso adaptativo responde a la dificultad local. En el oscilador de frecuencia creciente, el integrador redujo el paso en la misma proporción en que aumentó la frecuencia. El error debe evaluarse, sin embargo, respecto de la propiedad que importa, ya sea la fase, la amplitud, la conservación de una cantidad o el tiempo de un evento.

Un cálculo numéricamente preciso puede corresponder a un modelo empíricamente incorrecto, porque la verificación computacional y la validación científica son operaciones distintas. El capítulo 26 examinará los problemas rígidos, en los que la estabilidad limita el paso de los métodos explícitos aunque la precisión permitiría pasos mucho mayores.

Ejercicios

Ejercicio 25.1. (a) Aplique un paso del método del punto medio a $y'=-y$ con $y_0=1$ y $h=0{,}5$ . (b) Compare el resultado con $e^{-0{,}5}$ y con un paso de Euler. (c) Muestre que la función de estabilidad del punto medio es $R(z)=1+z+z^2/2$ .

Ejercicio 25.2. Deduzca la función de estabilidad de RK4 sustituyendo $f(t,y)=\lambda y$ en sus cuatro etapas.

Ejercicio 25.3. Para el oscilador armónico, (a) verifique que $|R_{RK4}(ih)|^2=1-h^6/72+h^8/576$ . (b) ¿Para qué pasos amortigua RK4 el oscilador? (c) Calcule la energía después de 200 pasos con $h=0{,}1$ y compárela con la que da Euler.

Ejercicio 25.4. Un sistema tiene dos componentes con valores $y=(1000;\ 10^{-8})$ y estimaciones de error $e=(0{,}5;\ 2\times10^{-6})$ , con $\mathrm{rtol}=10^{-3}$ y $\mathrm{atol}=10^{-6}$ para ambas. (a) Calcule el error normalizado RMS y decida si el paso se acepta. (b) ¿Qué componente decide el resultado? (c) Repita el cálculo si la segunda componente se expresa en una unidad mil veces menor, de modo que valga $10^{-5}$ con un error de $2\times10^{-3}$ . Comente el efecto de mantener la misma atol.

Ejercicio 25.5. Un controlador usa $h_{\mathrm{nuevo}}=0{,}9\,hE^{-1/5}$ . (a) Calcule el factor de cambio del paso si $E=2{,}5$ y si $E=0{,}01$ . (b) Si el estimador fuera exactamente proporcional a $h^5$ , ¿qué valor de $E$ tendría el paso repetido después de $E=2{,}5$ ? (c) ¿Para qué sirve el factor $0{,}9$ ?

Ejercicio 25.6. En el oscilador de frecuencia creciente, (a) calcule cuántas vueltas completa la trayectoria hasta $t=40$ . (b) Explique por qué el cociente entre los pasos medianos del comienzo y del final se acerca a $\omega(37)/\omega(3)$ . (c) Con $\mathrm{rtol}=10^{-3}$ , el error de estado fue de $0{,}0396$ y el del radio, de $0{,}000785$ . ¿Qué error de fase corresponde aproximadamente al error de estado?

Consultar las soluciones de los ejercicios

Soluciones de los ejercicios

Solución 25.1. (a) Se obtiene $k_1=-1$ , luego $k_2=-(1-0{,}25)=-0{,}75$ y finalmente $y_1=1+0{,}5\times(-0{,}75)=0{,}625$ . (b) El valor exacto es $e^{-0{,}5}\approx0{,}606531$ , de modo que el error del punto medio es $0{,}0185$ . Euler da $0{,}5$ , con un error de $0{,}1065$ , casi seis veces mayor. (c) Con $f=\lambda y$ se tiene $k_1=\lambda y_n$ y $k_2=\lambda(1+z/2)y_n$ , así que $y_{n+1}=y_n+hk_2=(1+z+z^2/2)\,y_n$ .

Solución 25.2. Con $z=h\lambda$ , las etapas son $k_1=\lambda y_n$ , $k_2=\lambda(1+z/2)y_n$ , $k_3=\lambda(1+z/2+z^2/4)y_n$ y $k_4=\lambda(1+z+z^2/2+z^3/4)y_n$ . Al sustituirlas en la actualización se obtiene $y_{n+1}=\left[1+\frac{z}{6}\left(6+3z+z^2+\frac{z^3}{4}\right)\right]y_n$ , que es $R_{RK4}(z)y_n$ .

Solución 25.3. (a) Con $z=ih$ , la parte real de $R$ es $1-h^2/2+h^4/24$ y la imaginaria, $h-h^3/6$ . La suma de sus cuadrados da $1-h^6/72+h^8/576$ . (b) El módulo es menor que uno cuando $h^8/576<h^6/72$ , es decir, cuando $0<h<2\sqrt2\approx2{,}83$ . Con pasos mayores, RK4 amplifica el oscilador. (c) RK4 da $0{,}5\,(1-10^{-6}/72+10^{-8}/576)^{200}\approx0{,}4999986$ . Euler da $0{,}5\times1{,}01^{200}\approx3{,}658$ .

Solución 25.4. (a) Las escalas son $s_1\approx1{,}000001$ y $s_2\approx1{,}00001\times10^{-6}$ . Los cocientes $e_i/s_i$ valen $0{,}5$ y $2{,}0$ , de modo que $E=\sqrt{(0{,}25+4)/2}\approx1{,}458$ y el paso se rechaza. (b) Decide la segunda componente, que es numéricamente diminuta, pero cuyo error duplica su escala. (c) Con la misma atol, la escala de la segunda componente sigue cerca de $10^{-6}$ y su cociente pasa a unos 1980, lo que da $E\approx1400$ . El cambio de unidades vuelve el criterio mil veces más exigente para esa componente, así que la tolerancia absoluta debe convertirse junto con la variable, en este caso a $10^{-3}$ .

Solución 25.5. (a) Con $E=2{,}5$ el factor es $0{,}9\times2{,}5^{-1/5}\approx0{,}749$ , una reducción de alrededor de un cuarto. Con $E=0{,}01$ el factor es $0{,}9\times0{,}01^{-1/5}\approx2{,}26$ . (b) El nuevo estimador valdría $2{,}5\times0{,}749^5\approx0{,}59$ y el paso se aceptaría. En ese caso ideal, el resultado es siempre $0{,}9^5\approx0{,}59$ , cualquiera sea el valor inicial de $E$ . (c) Sin el factor, el paso repetido quedaría en el límite $E=1$ y cualquier desviación del comportamiento supuesto provocaría un nuevo rechazo. El factor $0{,}9$ deja un margen.

Solución 25.6. (a) Como $\phi(40)=40+80=120$ radianes, la trayectoria completa $120/(2\pi)\approx19{,}1$ vueltas. (b) Si el error de cada paso depende sobre todo del ángulo recorrido, $\omega h$ , mantenerlo cerca de la tolerancia exige que $\omega h$ sea aproximadamente constante. Entonces $h$ es inversamente proporcional a $\omega$ y el cociente de pasos se acerca a $\omega(37)/\omega(3)=4{,}7/1{,}3\approx3{,}6$ . (c) En la circunferencia unitaria, dos puntos separados por un ángulo $\delta$ distan $2\sin(\delta/2)\approx\delta$ . Como el error de radio es pequeño, el error de estado corresponde casi por completo a un desfase de unos $0{,}04$ radianes, cerca de $2{,}3$ grados.

Referencias

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

Dormand, J. R. y Prince, P. J. (1980). A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, 6(1), 19–26. https://doi.org/10.1016/0771-050X(80)90013-3

Hairer, E., Nørsett, S. P. y Wanner, G. (1993). Solving ordinary differential equations I: Nonstiff problems (2.ª ed. rev.). Springer.

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

LeVeque, R. J. (2007). Finite difference methods for ordinary and partial differential equations: Steady-state and time-dependent problems. Society for Industrial and Applied Mathematics.

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

SciPy Developers. (s. f.-b). 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

Shampine, L. F. (1986). Some practical Runge-Kutta formulas. Mathematics of Computation, 46(173), 135–150.