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 la estabilidad A. 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 pasos de tiempo fijos como el método de Euler puede ser suficiente. Sin embargo, las cosas son más difíciles si uno 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 uno puede tomar grandes pasos de tiempo cuando la nave espacial está lejos de la Tierra y la Luna, pero si la nave espacial se acerca a colisionar con uno de los cuerpos planetarios, entonces se necesitan pequeños pasos de tiempo. El método de Romberg y el de 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 métodos de orden superior como los métodos de Runge-Kutta debido a sus propiedades superiores de convergencia y estabilidad.

Consideremos 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 nos interesa hallar 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 , donde es 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 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)}$

No conocemos el valor de c . 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 a . Supongamos que el nuevo tamaño de paso es la mitad del tamaño de paso original y apliquemos dos pasos del método de Euler. Esta segunda solución es presumiblemente más precisa. Dado que 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í, suponemos que el factor de error es constante en el intervalo . En realidad, su tasa de cambio es proporcional a . Al restar las soluciones, obtenemos 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 de 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 de la siguiente manera: ${\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)}$

El 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 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)}}$

Esta solución es, en realidad, precisa en tercer orden en el ámbito local (en segundo orden en el ámbito global), pero como no hay una estimación de error para ella, no ayuda a reducir la cantidad de pasos. Esta técnica se denomina extrapolación de Richardson .

Partiendo de un tamaño de paso inicial de , esta teoría facilita nuestra integración controlable de la EDO desde el punto hasta , utilizando una cantidad óptima de pasos dada una tolerancia de error local. Una desventaja 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 Runge-Kutta de cuarto orden. Además, se puede lograr una tolerancia de error global escalando el error local a un alcance global.

Estimaciones de errores incorporados

Los métodos de tamaño de paso adaptativo que utilizan una estimación de error denominada "integrada" incluyen los métodos Bogacki-Shampine , Runge-Kutta-Fehlberg , Cash-Karp y Dormand-Prince . Se considera que estos métodos son más eficientes desde el punto de vista computacional, pero tienen una 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 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 de incluye un método RK de orden inferior . El error 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 contra 1 para obtener el valor predicho : ${\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 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. ${\displaystyle {\tilde {\psi }}}$

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

Véase 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.

Lectura adicional

• 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