Tamaño de paso adaptativo

En matemáticas y análisis numérico , se utiliza un tamaño de paso adaptativo en algunos métodos para la solución numérica de ecuaciones diferenciales ordinarias (incluido el caso especial de la integración numérica ) con el fin de controlar los errores del método y garantizar propiedades de estabilidad como la estabilidad A. El uso de un tamaño de paso adaptativo es de particular importancia cuando existe una gran variación en el tamaño de la derivada. Por ejemplo, al modelar el movimiento de un satélite alrededor de la Tierra como una órbita Kepler estándar , un método de paso de tiempo fijo como el método de Euler puede ser suficiente. Sin embargo, la situación se complica si se desea modelar el movimiento de una nave espacial teniendo en cuenta tanto la Tierra como la Luna, como en el problema de los tres cuerpos . En este caso, surgen escenarios donde se pueden usar pasos de tiempo grandes cuando la nave espacial está lejos de la Tierra y la Luna, pero si la nave espacial se acerca a una colisión con uno de los cuerpos planetarios, entonces se necesitan pasos de tiempo pequeños. El método de Romberg y Runge-Kutta-Fehlberg son ejemplos de métodos de integración numérica que utilizan un tamaño de paso adaptativo.

Para simplificar, el siguiente ejemplo utiliza el método de integración más sencillo, el método de Euler ; en la práctica, se prefieren los métodos de orden superior, como los métodos de Runge-Kutta, debido a sus propiedades superiores de convergencia y estabilidad.

Consideremos el problema de valor inicial.

${\displaystyle y'(t)=f(t,y(t)),\qquad y(a)=y_{a}}$

donde y y f pueden denotar vectores (en cuyo caso esta ecuación representa un sistema de EDO acopladas en varias variables).

Se nos da la función f ( t , y ) y las condiciones iniciales ( a , y _a ), y estamos interesados ​​en encontrar la solución en t  = b . Sea y ( b ) la solución exacta en b , y sea y _b la solución que calculamos. Escribimos ${\displaystyle y_{b}+\varepsilon =y(b)}$, dónde${\displaystyle \varepsilon }$es el error en la solución numérica.

Para una secuencia ( t _n ) de valores de t , con t _n = a + nh , el método de Euler proporciona aproximaciones a los valores correspondientes de y ( t _n ) como

${\displaystyle y_{n+1}^{(0)}=y_{n}+hf(t_{n},y_{n})}$

El error de truncamiento local de esta aproximación se define por

${\displaystyle \tau _{n+1}^{(0)}=y(t_{n+1})-y_{n+1}^{(0)}}$

y por el teorema de Taylor , se puede demostrar que (siempre que f sea suficientemente suave) el error de truncamiento local es proporcional al cuadrado del tamaño del paso:

${\displaystyle \tau _{n+1}^{(0)}=ch^{2}}$

donde c es alguna constante de proporcionalidad.

Hemos marcado esta solución y su error con un${\displaystyle (0)}$.

El valor de c no nos es conocido. Apliquemos ahora de nuevo el método de Euler con un tamaño de paso diferente para generar una segunda aproximación a y ( t _{n +1} ). Obtenemos una segunda solución, que etiquetamos con un${\displaystyle (1)}$Consideremos que el nuevo tamaño del paso es la mitad del tamaño original y apliquemos dos pasos del método de Euler. Esta segunda solución es presumiblemente más precisa. Dado que debemos aplicar el método de Euler dos veces, el error local es (en el peor de los casos) el doble del error original.

${\displaystyle y_{n+{\frac {1}{2}}}=y_{n}+{\frac {h}{2}}f(t_{n},y_{n})}$

${\displaystyle y_{n+1}^{(1)}=y_{n+{\frac {1}{2}}}+{\frac {h}{2}}f(t_{n+{\frac {1}{2}}},y_{n+{\frac {1}{2}}})}$

${\displaystyle \tau _{n+1}^{(1)}=c\left({\frac {h}{2}}\right)^{2}+c\left({\frac {h}{2}}\right)^{2}=2c\left({\frac {h}{2}}\right)^{2}={\frac {1}{2}}ch^{2}={\frac {1}{2}}\tau _{n+1}^{(0)}}$

${\displaystyle y_{n+1}^{(1)}+\tau _{n+1}^{(1)}=y(t+h)}$

Aquí, asumimos un factor de error.${\displaystyle c}$es constante en el intervalo${\displaystyle [t,t+h]}$En realidad, su tasa de cambio es proporcional a${\displaystyle y^{(3)}(t)}$Al restar las soluciones se obtiene la estimación del error:

${\displaystyle y_{n+1}^{(1)}-y_{n+1}^{(0)}=\tau _{n+1}^{(1)}}$

Esta estimación de error local tiene una precisión de tercer orden.

La estimación del error local se puede utilizar para decidir el tamaño del paso.${\displaystyle h}$debe modificarse para lograr la precisión deseada. Por ejemplo, si una tolerancia local de${\displaystyle {\text{tol}}}$Si está permitido, podríamos dejar que h evolucione así:

${\displaystyle h\rightarrow 0.9\times h\times \min \left(\max \left(\left({\frac {\text{tol}}{2\left|\tau _{n+1}^{(1)}\right|}}\right)^{1/2},0.3\right),2\right)}$

El${\displaystyle 0.9}$es un factor de seguridad para garantizar el éxito en el siguiente intento. El mínimo y el máximo sirven para evitar cambios extremos con respecto al tamaño del paso anterior. En principio, esto debería dar un error de aproximadamente${\displaystyle 0.9\times {\text{tol}}}$en el siguiente intento. Si${\displaystyle |\tau _{n+1}^{(1)}|<{\text{tol}}}$Consideramos que el paso ha sido exitoso y la estimación del error se utiliza para mejorar la solución:

${\displaystyle y_{n+1}^{(2)}=y_{n+1}^{(1)}+\tau _{n+1}^{(1)}}$

Esta solución tiene una precisión de tercer orden en el ámbito local (y de segundo orden en el ámbito global), pero como no dispone de una estimación de error, no contribuye a reducir el número de pasos. Esta técnica se denomina extrapolación de Richardson .

Comenzando con un tamaño de paso inicial de${\displaystyle h=ba}$, esta teoría facilita nuestra integración controlable de la EDO desde el punto${\displaystyle a}$a${\displaystyle b}$, utilizando un número óptimo de pasos dada una tolerancia de error local. Una desventaja es que el tamaño del paso puede volverse prohibitivamente pequeño, especialmente cuando se utiliza el método de Euler de orden bajo .

Se pueden desarrollar métodos similares para métodos de orden superior, como el método de Runge-Kutta de cuarto orden. Además, se puede lograr una tolerancia de error global escalando el error local al ámbito global.

Entre los métodos de tamaño de paso adaptativo que utilizan una estimación de error denominada "incrustada" se incluyen los métodos de Bogacki-Shampine , Runge-Kutta-Fehlberg , Cash-Karp y Dormand-Prince . Estos métodos se consideran más eficientes desde el punto de vista computacional, pero presentan menor precisión en sus estimaciones de error.

Para ilustrar las ideas del método integrado, considere el siguiente esquema que actualiza${\displaystyle y_{n}}$:

${\ Displaystyle y_ {n+1} = y_ {n} + h_ {n} \ psi (t_ {n}, y_ {n}, h_ {n})}$

${\displaystyle t_{n+1}=t_{n}+h_{n}}$

El siguiente paso${\displaystyle h_{n}}$Se predice a partir de la información previa.${\displaystyle h_{n}=g(t_{n},y_{n},h_{n-1})}$.

Para el método RK embebido, el cálculo de${\displaystyle \psi }$incluye un método RK de orden inferior${\displaystyle {\tilde {\psi }}}$El error entonces se puede escribir simplemente como

${\displaystyle {\textrm {err}}_{n}(h)={\tilde {y}}_{n+1}-y_{n+1}=h({\tilde {\psi }}(t_{n},y_{n},h_{n})-\psi (t_{n},y_{n},h_{n}))}$

${\displaystyle {\textrm {err}}_{n}}$es el error no normalizado. Para normalizarlo, lo comparamos con una tolerancia definida por el usuario, que consta de la tolerancia absoluta y la tolerancia relativa:

${\displaystyle {\textrm {tol}}_{n}={\textrm {Atol}}+{\textrm {Rtol}}\cdot \max(|y_{n}|,|y_{n-1}|)}$

${\displaystyle E_{n}={\textrm {norm}}({\textrm {err}}_{n}/{\textrm {tol}}_{n})}$

Luego comparamos el error normalizado.${\displaystyle E_{n}}$contra 1 para obtener el pronosticado${\displaystyle h_{n}}$:

${\displaystyle h_{n}=h_{n-1}(1/E_{n})^{1/(q+1)}}$

El parámetro q es el orden correspondiente al método RK.${\displaystyle {\tilde {\psi }}}$, que tiene un orden inferior. La fórmula de predicción anterior es plausible en el sentido de que aumenta el paso si el error local estimado es menor que la tolerancia y lo reduce en caso contrario.

La descripción anterior es una simplificación comparada con los procedimientos utilizados en el control del tamaño de paso para los solucionadores explícitos de RK. Un tratamiento más detallado se puede encontrar en el libro de texto de Hairer. ^{[ 1 ]} El solucionador de EDO en muchos lenguajes de programación utiliza este procedimiento como estrategia predeterminada para el control adaptativo del tamaño de paso, que agrega otros parámetros de ingeniería para hacer que el sistema sea más estable.

• Cuadratura adaptativa
• Diferenciación numérica adaptativa

• William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery, Recetas numéricas en C , Segunda edición, CAMBRIDGE UNIVERSITY PRESS, 1992. ISBN 0-521-43108-5
• Kendall E. Atkinson, Análisis numérico , Segunda edición, John Wiley & Sons, 1989. ISBN 0-471-62489-6