Paso adaptativo

En este segundo artículo vamos a ver una de las formas de refinar nuestro motor de inferencia numérico. Anteriormente hablamos de métodos numéricos aplicables al cálculo de las sucesivas posiciones / velocidades ( "Phase Space" ) de nuestras partículas. Vimos que necesitábamos evaluar ecuaciones diferenciales de primer orden y que lo haríamos a una determinada velocidad de animación en el tiempo. Simpre siguiendo un cierto paso de nuestro sistema. Con el paso nos referimos al tiempo que dura una secuencia de animación, una iteración...pero no al tiempo real que tarda el sistema en calcularla sinó al tiempo de "escena" que se supone que pasa entre un cuadro (frame) y otro.

Recordemos un poquito el diagrama básico de nuestro sistema de simulación física:

Recordemos para que sirve cada cosa:
 

Después de cada ciclo empezamos de nuevo. Recordad que siempre tenemos un tiempo inicial de simulación y uno final. Como un bucle, iteraremos el sistema hasta que el tiempo inicial más un número N de pasos igualen o rebasen al tiempo final.

Ajustemos pués...

Y ¿qué es eso de ajustar el paso?, ¿qué quiere decir hacerlo adaptativo?...y ¿para qué después de todo?. Seguid leyendo y os responderé a todas esas dudas.

¿Es lo mismo simular una pelota de golf botando entre hierbas caóticas que animar una pelotita cayendo sobre suelo llano?...obviamente no pero ¿por qué?. En el primer caso la simulación será increiblemente más compleja. Tendremos que evaluar cientos de miles de colisiones y rozamientos de la pelota con la hierba, variar su velocidad y trayectoria a cada momento y hacerlo muy poquito a poco para no dar una sensación visual errónea. En cambio en el segundo caso podemos casi predecir lo que ocurrirá a cada momento...caerá debido a la gravedad acelerándose, rebotará invirtiendo su velocidad...y volverá a subir y así a cada momento. Los cambios en la velocidad serán obvios, tan sólo invertir el signo al chocar con el suelo, y el movimiento es todo el rato vertical, de arriba a abajo. Nada de trayectorias extrañas que se modifican cada milisegundo.

¿Vais viendo a donde quiero llegar?...En el primer caso tenemos que simular a velocidad muy lenta para no equivocarnos. Las derivadas varian casi contínuamente y si tomamos intervalos de tiempo muy largos para evaluarlas cometeremos un error grave. En el segundo caso las derivadas varian mínimamente pues los cambios en la trayectoria son insignificantes, sólo subimos y bajamos en vertical!!!, así que podemos acelerar la animación tomando pasos más grandes seguros de no equivocarnos.

Por lo tanto lo que haremos será adaptar el paso que utiliza nuestro sistema a cada momento según estemos en un "momento" de la simulación más o menos complejo. Si es fácil, corremos más y aumentamos el paso. Si es difícil disminuimos el paso yendo más lentos pero también manteniendo la cota máxima de error permitida.

De esta forma siempre iremos a la máxima velocidad posible cometiendo un error inferior al máximo. Si mantenemos el paso estático e immutable como antes iremos siempre a igual velocidad, demasiado deprisa a veces y demasiado lentos otras...nunca contentos!!

Recordad que nos interesa simular eventos físicamente correctos pero tampoco queremos eternizar la simulación con horas y horas de cálculo...y no hablemos del tiempo real!!!

¿Cómo lo implemento en mi sistema?

Hay varias formas de hacerlo, mejores y peores como siempre. Yo os comentaré una muy sencilla que ya dá resultados satisfactorios y visibles incluso a simple vista. De hecho he probado con varias simulaciones y usando un simple péndulo oscilante la velocidad de la animación se incrementaba notablemente ya que el paso se aumentaba muchas veces manteniendo el error ( hacia el render a tiempo real...es decir calcular y dibujar, calcular y dibujar, y así secuencialmente ).

La idea es sencilla y además es totalmente independiente del solver numérico empleado. Funciona igual usando el método de Euler, que el punto medio o Runge Kutta-4. Obviamente a mejor solver mejor sistema global pero a efectos de implementación del paso adaptativo, las líneas de código son prácticamente las mismas. A cada iteración del sistema haremos dos cálculos distintos para todas las partículas. De su diferencia averiguaremos cuál es el error que estamos cometiendo y según la magnitud de este, aumentaremos, disminuiremos o dejaremos igual el paso de simulación.

El algoritmo podría expresarse así:

1. Hacemos un "Back-Up" de la situación inicial del sistema (posición y velocidad de todas las partículas).
2. Evaluamos el nuevo Phase-Space de todas las partículas con el solver usando nuestro paso de valor H.
3. Guardamos los resultados en un array temporal.
4. Ahora, y usando el Back-Up inicial, volvemos a evaluar el sistema pero de la siguiente forma:
• Una primera evaluación usando un paso de H / 2.
• Los resultados obtenidos se vuelven a iterar durante otro paso de valor H / 2.
5. Guardamos los resultados en otro array temporal.
6. Restamos partícula a partícula las posiciones / velocidades y buscamos si alguna de las diferencias supera el error máximo permitido en el sistema.
7. Si es así, disminuiremos el paso y volveremos a iterar desde el Back-Up inicial utilizándolo.
8. Sinó miraremos si podemos aumentarlo y cuanto !!. Si lo aumentamos podemos iterar de nuevo desde el Back-Up o usar el paso recién aumentado para la próxima vez ( así nos ahorramos tiempo ahora y aprovechamos el nuevo paso luego ).

Para la primera iteración (me refiero a la inicial de todas, la primera) se parte de un paso H predefinido en el programa como por ejemplo de valor 1 / 25 ( 0.04 ). Luego el sistema lo va ajustando iteración tras iteración.

En cuanto al error...pues hombre valores típicos son errores máximos de 0.001, 0.0001, 0.00001 y en adelante. A menos error permitido más lenta irá la simulación pues más serán los ajustes. Hay que vigilar mucho aquí. No os paseis con las restricciones porque sinó el sistema os irá a velocidad de caracol y según lo que esteis simulando no tiene porqué ser necesaria tanta TANTÍSIMA perfección!!! ( como véis en temas de Informática siempre se llega a un compromiso...;)

Bién. Hasta aquí sabemos calcular el error en el sistema a cada iteración pero no sabemos ajustar el paso en consecuencia. La formulita mágica es la siguiente:

Bueno veamos que es cada cosa. El triangulito es el error aproximado que estamos cometiendo. Como podeis observar, es la diferencia entre los dos estados del sistema, el correspondiente a un paso de valor H y el referido a dos pasos consecutivos de valor H / 2.

En cuanto a la fórmula de abajo aquí es dónde llega la novedad. El paso nuevo a calcular ( es decir el paso óptimo, mayor o menor al anterior ) es lo que deseamos conocer y es h0. El valor h1 se corresponde con el paso que hemos utilizado para calcular el error (el que yo he estado nombrando como H ). En cuanto a los triangulitos:

Si el denominador es mayor que el numerador, la ecuación nos estará diciendo cuanto hay que disminuir el paso para recortar el error hasta el valor deseado. En caso contrario, la ecuación nos dirá en que medida podemos incrementar el paso sin peligro alguno para la siguiente iteración. Vaya!, que h0 es el paso óptimo a utilizar para conseguir un error controlado y una velocidad siempre acorde con él.

Si a cada iteración implementamos este paso adaptativo estamos cargando a nuestro sistema con más cálculos, obvio, pero la recompensa aumenta con creces el rendimiento final del sistema que se acelera cuando puede y se equivoca siempre lo mínimo!!!

Un abrazo para todos/as. Nos vemos detectando colisiones en el siguiente artículo !
 
 

ÚLTIMA REVISIÓN EN JULIO DE 1999