Un Modelo básico: Integrando ecuaciones de movimiento con el método de Euler

En los anteriores posts, expliqué las razones para usar MVC y la estructura general de un programa simple usando este patrón de diseño. Ahora voy a explicar la versión mas básica de un modelo: un simulador de partículas. Bajo ninguna circunstancia este modelo es un modelo aplicable directamente para un juego, pero es una primer iteración para poder implementar algunas cosas básicas como Steering Behaviors (que voy a explicar en otro post).

Este modelo considera partículas puntuales que se mueven libremente o bajo la acción de una fuerza, en un espacio de 2 dimensiones. Con este modelo simple se pueden implementar muchos tipos de juegos 2D: Juegos de plataforma, Arcanoid, Galaga, juegos de carrera o los juegos que sugirió Juan Pablo: Choplifter, Desert strike y Cannon Fodder. Que el Modelo sea 2D no implica usar gráficos 2D (el View), incluso se pueden usar gráficos 3D. Ejemplos clásicos de esto es el Wolfenstein 3D, donde en realidad el mapa es 2D, no existe la altura en el modelo, y sin embargo la representación es 3D. Un juego como Diablo, también tiene gráficos 3D siendo el modelo 2D. Con esto quiero aclarar que tener un modelo 2D no implica para nada utilizar gráficos 2D.

Derivadas

Antes que nada voy a explicar muy básicamente lo que es una derivada. Se puede encontrar una definición mas avanzada en la Wikipedia.

La derivada de la posición respecto al tiempo es el cociente entre la variación de la posición y la variación de tiempo. Si un auto avanza 10 km en 1 hora, podemos aproximar la derivada como 10km/1h=10km/h.

Pero la derivada es más que eso, la derivada implica que la variación en el tiempo sea muy pequeña. En nuestro caso, tendríamos que ver cuanto avanzó el auto en, por ejemplo, un segundo. Si el auto avanzó 3 metros en un segundo, entonces la derivada de la posición respecto al tiempo es más precisamente 3 m/1s=3m/s=10.8 m/s.

Se puede refinar el experimento, midiendo el desplazamiento para intervalos de tiempo aun mas pequeños. La derivada es el valor que obtenemos cuando el intervalo de tiempo es infinitamente pequeño:

En esta ecuación el triángulo es una forma de escribir "La variación de". O sea que la ecuación puede ser leída como "La variación de la posición, dividido la variación del tiempo". Cuando la variación del tiempo es infinitamente pequeña, el cociente se llama derivada. El caso de la derivada de la posición respecto al tiempo,  es denominada "la velocidad". Pero uno puede derivar la velocidad, o sea, dividir un cambio de velocidad por el intervalo de tiempo en que ocurre este cambio. Esta última derivada se llama "aceleración".

Cuando derivamos dos veces, se suele decir que obtenemos "la derivada segunda". Entonces, la derivada segunda de la posición es la aceleración.

Una forma de escribir la derivada es poniendo un punto sobre la variable que se deriva. Y la derivada segunda se escribe con dos puntos sobre la variable que se deriva. Pero el significado es el mismo: el cociente entre la variación de una magnitud y la variación de tiempo durante el cual ocurre el cambio.

Ecuaciones de movimiento

Supongamos que tenemos un objeto que se mueve con una fuerza FX, FY aplicada. La ecuación de movimiento mas general que se puede escribir es:

Donde m es la masa del objeto. X e Y son las coordenadas de la partícula,  y los dos puntos sobre la variable X es la derivada segunda de la posición respecto al tiempo (o sea, derivar la posición dos veces respecto al tiempo). Como dijimos anteriormente, esta derivada segunda se conoce como aceleración. "Fuerza igual masa por aceleración" es la llamada segunda ley de Newton.

Esta fórmula debe ser resuelta para obtener ecuaciones de movimiento que podamos usar. En el caso de una fuerza constante, las ecuaciones pueden ser resueltas en forma exacta:

 

Para quien sepa derivadas, es fácil ver que si derivamos la expresión de X dos veces, obtenemos la segunda ley de Newton con una fuerza constante (a0 es la fuerza dividida la masa).

Estas ecuaciones consideran solo la coordenada X. En las ecuaciones, el tiempo aparece en forma explicita como t, y los subindices 0 en la posición, velocidad y aceleración indican que se trata de la posición, la velocidad y la aceleración iniciales. En esta ecuación, la aceleración es constante, por tratarse de una fuerza constante. La segunda ecuación es la velocidad en función del tiempo, que en algunas situaciones es útil.

Las ecuaciones para las coordenada en Y son idénticas, pero reemplazando X por Y y utilizando las velocidades en la dirección Y en lugar de usar las velocidades en X y la aceleración en Y en lugar de la aceleración en X.

Estas ecuaciones tienen la limitación de que solo son válidas mientras la fuerza se mantiene constante. De cambiar la fuerza, cambia la aceleración, y las ecuaciones dejan de ser válidas. Si en general la fuerza varía en forma complicada con el tiempo, se hace imposible resolver las ecuaciones de movimiento en forma exacta.

El caso de fuerza que varía con el tiempo debe ser resuelto para cada caso particular. Las ecuaciones que resultan son en general mas complicadas. En un juego, necesitamos utilizar algún método general para resolver las ecuaciones independientemente de como varíe la fuerza.

Para hacer esto, lo mas simple es usar el llamado Método de Euler. El método de Euler es un método numérico para resolver ecuaciones diferenciales de primer orden, o sea que solo tienen derivadas primeras. La primer ecuación que mostré es una ecuación de segundo orden, pero eso puede ser resuelto fácilmente introduciendo la velocidad. La velocidad, como dijimos, es la primer derivada de la posición:

Utilizando esta definición, es posible convertir nuestra ecuación de movimiento en varias ecuaciones de primer orden, ya que la aceleracion es la primer derivada de la velocidad:

Con estas cuatro ecuaciones, tenemos solo derivadas primeras, con lo cual podemos aplicar el método de Euler.

El método de Euler consiste simplemente en volver a la definición de la derivada, y "olvidarse" que la derivada requiere que la variación del tiempo sea infinitamente pequeña. En lugar de eso, se utiliza una variacion de tiempo finita.

Usando esa aproximación, las ecuaciones anteriores resultan ser:

La ventaja de estas ecuaciones es que podemos resolverlas algebraicamente para obtener la posición de la partícula en un tiempo t, un instante de  tiempo anterior t- Dt. Simplifiquemos un poco la notación, para que sea mas evidente como resolverlas:

• Llamemos a la posición a tiempo t utilizando un subíndice: Xt es la posición a tiempo t.
• Llamemos h a la variación del tiempo (también llamado paso de tiempo).
  
• La variación de la posición no es mas que la diferencia entre dos posiciones. O sea:
  
  
  
• Lo mismo puede decirse de la velocidad. La velocidad es la diferencia de la velocidad en dos tiempos distintos:
  
  

Teniendo eso en cuenta, resultan las siguientes ecuaciones (considerando solo las coordenadas en X):

Lo bueno de estas ecuaciones, es que solo requieren una multiplicación y una suma. Además, son absolutamente generales, pueden ser usadas para cualquier tipo de fuerza(aunque en algunas situaciones puede fallar).

La desventaja, es que es el método mas simple, y el mas impreciso. En general, la precisión de este método depende del tamaño de h (el intervalo de tiempo). Si usamos un h mas pequeño, el resultado va a ser mas preciso, pero van a ser necesarios mas cálculos para encontrar el valor de la posición.

El algoritmo en Python

Veamos como se traducen las ecuaciones anteriores para el caso de un modelo con una sola partícula, en Python.

class ParticulaSimple():
   def __init__(self, x0, y0, vx0, vy0, masa):
      self.x=x0
      self.y=y0
      self.vx=vx0
      self.vy=vy0 
      self.t=0
      self.masa=masa
      
   def get_force(self):
      '''
      Devuelve la fuerza a tiempo self.t.
      '''
      ...
      return force_x, force_y

   def update(self, h):
      t, x_prev, y_prev= self.t, self.x, self.y,
      vx_prev, vy_prev=self.vx, self.vy
      masa =  self.masa
      #Obtener la fuerza a tiempo t
      fuerza_x, fuerza_y=self.get_force()

      #Actualizar el tiempo, posición 
      #y velocidad usando el método de Euler.
      #Actualiza a tiempo t+h
      self.t = t+h
      self.x +=vx_prev*h
      self.y +=vy_prev*h
      self.vx += (fuerza_x*1.0/masa) *h
      self.vy += (fuerza_y*1.0/masa) *h

El código es bastante simple. En caso de ser más de una partícula, se debe iterar sobre todas las partículas. Usando librerías como numpy todo esto puede hacerse en una sola instrucción (multiplicando matrices por vectores). Pero eso no importa ahora.

Una parte clave del algoritmo es la función get_force, que devuelve las fuerzas en la dirección x e y. Esta función puede depender del tiempo, y en situaciones mas complejas puede depender de la velocidad (por ejemplo, la fricción es proporcional a la velocidad) o también puede depender de la posición (por ejemplo, la fuerza que hace un resorte depende linealmente del estiramiento del resorte).

Para ejemplificar muestro tres tipos de fuerzas:

Gravedad cerca de la tierra

Esta fuerza es la fuerza que experimentan todos los objetos cercanos a la superficie de la tierra. Es una fuerza constante, con dirección "hacia abajo".

...
   get_force(self):
      #La constante gravitatoria es un valor 
      #que depende solamente del planeta
      #donde se está parado. 
      #Para la tierra es g=9.8 m/s^2
      #Para la luna es g=1.6 m/s^2
      
      g=9.8
      return 0, -self.masa*g  #La fuerza es vertical
      #, apuntando hacia abajo(y<0 es abajo)

   ...

Fuerza gravitatoria general

Si el objeto en cuestión se encuentra alejado del planeta, la fuerza depende de la distancia al centro del planeta. La dirección de la fuerza es directamente apuntando al centro del planeta.

...
   get_force(self):
      #En este caso, la fuerza depende 
      #de la distancia al planeta y 
      #de la masa del mismo. 
      from math import sqrt
      #Un planeta en la coordenada (0,0)
      planetaX, planetaY= 0, 0 
      #Un planeta de la masa de la tierra
      k=3.98*10^14 

      distancia_X=planetaX-self.x
      distancia_Y=planetaY-self.y
      distancia_al_cuadrado=distancia_X*distancia_X\\
         +distancia_Y*distancia_Y

      F_magnitud=k*self.masa*1.0/distancia_al_cuadrado

      x_dir=distancia_X/sqrt(distancia_al_cuadrado)
      y_dir=distancia_Y/sqrt(distancia_al_cuadrado)

      return F_magnitud*x_dir, F_magnitud*Y_dir

   ...

Fricción de aire

Esta fuerza es una fuerza que depende de la velocidad. Para velocidades pequeñas la dependencia es lineal y puede programarse de la siguiente manera:

...
   get_force(self):
      #El coefficiente de fricción 
      #es un indice de cuanta fricción hay
      # si es 0 no hay ninguna fricción 
      #(la fuerza de fricción es siempre 0)

      fric_coef=0.1
      return -fric_coef*self.vx, -fric_coef*self.vy
   ...

Por razones que no voy a explicar ahora, esta fuerza no funciona muy bien con el método de Euler. Pero se puede usar si se mantiene el coeficiente de fricción muy chico o el paso de tiempo suficientemente pequeño. Lo mejor para un juego es probar distintos coeficientes hasta encontrar un valor que dé resultados plausibles.

Errores

Este método, como cualquier otro método de integración numérica, tiene errores. Para mostrar la diferencia entre el método exacto y el método aproximado voy a usar el ejemplo de la gravedad cercana a la tierra, que tiene solución exacta. Supongamos que un personaje de un juego de plataformas viene corriendo a 7 metros por segundo, y salta adquiriendo la misma velocidad.  La velocidad inicial es 7 m/s. vertical y 7 m/s horizontal La trayectoria del personaje va a estar dada en forma exacta por:

Comparemos la solución exacta con la solución numérica usando h=0.05

Comparación de resultados para una fuerza constante en dirección vertical.

Claramente hay un error. Generalmente, para juegos este tipo el resultado sea suficiente, ya que el movimiento es cualitativamente correcto y el algoritmo usado es rápido.

Si se usa un h mas chico, el error va a ser menor (pero van a ser necesarios más FPS del modelo). Existen otros algoritmos que resultan en valores mas precisos para iguales valores de h.

Aun más, en algunas aplicaciones el método de Euler no funciona adecuadamente. Para estas aplicaciones es necesario utilizar otros algoritmos como Verlet Velocidad.

Conclusiones

Quizá la matemática de estos algoritmos sea un poco complicada. Pero programarlos es casi trivial. Existen algunas aplicaciones en las cuales el paso de tiempo h necesario para tener un resultado adecuado es demasiado pequeño, con lo cual se deben utilizar otro tipo de métodos de integración.

En una aplicación normal, la fuerza y velocidad varían según las instrucciones del controlador. En cada paso de tiempo, se llama a la función update del modelo para actualizar las posiciones y al update de la vista para utilizar las posiciones para ubicar los sprites en pantalla. A pesar de que en el modelo utilizo unidades "reales" como metros, en la vista estas unidades se transforman en pixels de pantalla.

Por ahora, para explicar Steering Algorithms, es suficiente con este método de integración. En el próximo post voy a explicar como se crea un controller que utilize fuerzas para dirigir un vehiculo a su destino.