Tamaño de paso adaptable

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 integración numérica ) con el fin de controlar los errores del método y garantizar propiedades de estabilidad como A- estabilidad . El uso de un tamaño de paso adaptativo es de particular importancia cuando hay 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 estándar de Kepler , un método de paso de tiempo fijo como el método de Euler puede ser suficiente. Sin embargo, las cosas son más difíciles 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 . Allí surgen escenarios en los que se pueden tomar grandes pasos de tiempo cuando la nave espacial está lejos de la Tierra y la Luna, pero si la nave espacial está a punto de colisionar con uno de los cuerpos planetarios, entonces se necesitan pequeños pasos de tiempo. 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.

Ejemplo

Para simplificar, el siguiente ejemplo utiliza el método de integración más simple, 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.

Considere el problema del 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 , ¿dónde está el error en la solución numérica? ${\displaystyle y_{b}+\varepsilon =y(b)}$ ${\displaystyle \varepsilon }$

Para una secuencia ( t _n ) de valores de t , con t _n = a + nh , el método de Euler da 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 mediante 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 una constante de proporcionalidad.

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

No conocemos el valor de c . Apliquemos ahora el método de Euler nuevamente 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 . Tome el nuevo tamaño del paso como la mitad del tamaño del paso original y aplique dos pasos del método de Euler. Esta segunda solución es presumiblemente más precisa. Como tenemos que aplicar el método de Euler dos veces, el error local es (en el peor de los casos) el doble del error original. ${\displaystyle (1)}$

${\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 que el factor de error es constante durante el intervalo . En realidad, su tasa de cambio es proporcional a . Restar soluciones da la estimación del error: ${\displaystyle c}$ ${\displaystyle [t,t+h]}$ ${\displaystyle y^{(3)}(t)}$

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

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

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

${\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)}$

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

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

En realidad, esta solución tiene una precisión de tercer orden en el ámbito local (segundo orden en el ámbito global), pero como no hay una estimación de error para ella, esto no ayuda a reducir el número de pasos. Esta técnica se llama extrapolación de Richardson .

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

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 a errores globales escalando el error local al alcance global.

Estimaciones de error integradas

Los métodos adaptativos de tamaño de paso que utilizan la denominada estimación de error "incrustada" incluyen los métodos de Bogacki-Shampine , Runge-Kutta-Fehlberg , Cash-Karp y Dormand-Prince . Se considera que estos métodos son más eficientes computacionalmente, pero tienen menor precisión en sus estimaciones de error.

Para ilustrar las ideas del método integrado, considere el siguiente esquema que se 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 se predice a partir de la información anterior . ${\displaystyle h_{n}}$ ${\displaystyle h_{n}=g(t_{n},y_{n},h_{n-1})}$

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

${\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 con 1 para obtener lo previsto : ${\displaystyle E_{n}}$ ${\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 , que tiene 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. ${\displaystyle {\tilde {\psi }}}$

La descripción proporcionada anteriormente es un procedimiento simplificado utilizado en el control de tamaño de pasos para solucionadores RK explícitos. Se puede encontrar un tratamiento más detallado en el libro de texto de Hairer. ^[1] El solucionador ODE en muchos lenguajes de programación utiliza este procedimiento como estrategia predeterminada para el control adaptativo del tamaño de los pasos, que agrega otros parámetros de ingeniería para hacer que el sistema sea más estable.

Ver también

• Cuadratura adaptativa
• Diferenciación numérica adaptativa

Referencias

1. E. Hairer, SP Norsett G. Wanner, "Resolución de ecuaciones diferenciales ordinarias I: problemas no rígidos", sec. II.

Otras lecturas

• 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