Simulación física. Métodos numéricos

Bién está claro que deseamos animar objetos deformables. Para poder hacerlo necesitamos un motor de cálculo numérico muy robusto. Este será el motor encargado de calcular el nuevo valor para todas las posiciones y velocidades de cada una de nuestras partículas.

¿Qué has dicho Oscar?, ¿partículas?....¿y eso qué es?. Pues bién. Crearemos modelos deformables a partir de partículas conectadas entre sí con muelles. De esta manera conseguiremos interesantísimos comportamientos como la simulación del agua de mar, del fuego, del polvo, bandadas de pájaros, estampidas, todo es posible!!!!

Pero antes de empezar a disfrutar como cosacos tenemos que "zamparnos" algo de matemática.

Como os iba diciendo, generar fuerzas que serán aplicadas a partículas requiere evaluar ecuaciones. De hecho son ecuaciones diferenciales de primer orden en las que hay derivadas implicadas. ¿Quién es el guapo que evaluará esas derivadas?...pues el motor de inferencia dinámica, o sea, el motor interno de cálculo numérico.

Evaluación de las ecuaciones diferenciales

Una ecuación diferencial describe la relación entre una función cualquiera y el valor de sus derivadas. Para verificar la ecuación deberán cumplirse algunas condiciones. En la mayoría de los casos, y en el nuestro por supuesto, hablaremos de lo que se conoce como "problemas de valor inicial" dónde el comportamiento del sistema puede resumirse como una función derivada del tipo:

Aquí f es una función conocida que depende de x y de t; x es el estado del sistema y x’ su derivada. Esto es una ecuación diferencial ordinaria o EDO.

Normalmente sabremos algo así...

...es decir el valor del estado para un cierto tiempo, y a partir de aquí el sistema tendrá que ir variando su comportamiento según las fuerzas externas que se le apliquen.

Intentaremos obtener la derivada de una función en distintos intervalos discretos ya que no podemos hacerlo de forma contínua. Aplicaremos por tanto una solución numérica para resolver este tipo de problemas.

Tenemos una función derivada f que nos evalua los cambios en el estado del sistema x, que son lo que denominamos Incx, y lo hace a lo largo de un intervalo de tiempo llamado Inct. Entonces lo que haremos será incrementar el estado, o sea x, con el valor calculado Incx, obteniendo así el nuevo estado.

Los métodos numéricos evaluan una o más derivadas en un mismo intervalo de tiempo. Esto los hace distintos en cuanto a eficiencia y robustez.

Os comentaré un par de ellos para ir entrando en materia y poderlos aplicar a posteriori a nuestra animación de objetos deformables.

Método de Euler

Este método es el más simple de todos ellos. Si partimos de una condición inicial como:

entenderemos que un paso de tiempo después tendremos

dónde h es lo que se llama paso (delta o incremento).

En el método de Euler se evalua un nuevo estado del sistema partiendo de los valores conocidos del estado anterior...

De esta forma sabemos cual es el valor de nuestro sistema un tiempo h después del tiempo anterior to, gracias a conocer el estado anterior xo y la derivada x’(to).

Me parece que alguno se estará poniendo vizco de leer y no entender nada...ya sé que es muy oscuro al principio pero hay que digerirlo a base de cocotazos con la pared y una vez lo tienes en el estómago te lo pasas pipa simulando!!!

Tranquilos que luego hablaremos de ejemplillos y lo iréis viendo bastante más claro.

Este método es altamente inestable y poco eficaz. Puede "colgarnos" literalmente el ordenador en muchas ocasiones (no por falta de recursos sinó por problemas de tendencias a infinito y demás...) y suele evaluar valores de derivada incorrectos con lo cual lo que vemos en pantalla dista muchísimo de ser físicamente correcto o incluso visualmente aceptable.

De hecho esto es de esperar si hacéis un poco de memória e intentáis recordar la famosa fórmula de la série de Taylor para una función cualquiera. Si recordais, con Taylor podiamos aproximar una función cualquiera por una série infinita de polinomios de orden creciente. Algo así:

dónde el valor de nuestra función un pasito más lejos del punto t0, es decir en t0+h, se aproxima tal y cómo os muestro.

La fórmula de Euler proviene del truncamiento de esta fórmula a partir del segundo término con lo cual el error cometido en la aproximación de la función es notable y por tanto la solución es incorrecta.

Habra casos donde el error sea negligible, ínfimo, pero en muchos otros no será así creedme.

Existen métodos mucho más robustos, que evaluan 4 o 5 derivadas por iteración, que aún haciendo el coste computacional de nuestra aplicación más alto, nos dan unos resultados físicamente casi correctos y visualmente geniales.

Hablemos pues del método que encarecidamente os recomiendo usar para evaluar derivadas. El método Runge-Kutta de orden 4.

Método de Runge-Kutta, de orden 4 (RK4)

Dos matemáticos alemanes, Runge y Kutta como podéis imaginar, desarrollaron diversos algoritmos para resolver ecuaciones diferenciales eficientemente y manteniendo siempre un error mínimo.

Con uno de sus métodos, el Runge-Kutta de orden 4, conseguimos mejorar en mucho lo que Euler nos ofrecía. De hecho Euler es un Runge-Kutta de orden 2. El problema es que Euler avanza el sistema un intervalo de tiempo h usando la información de la derivada sólo al principio del intervalo. Con el RK4 evaluamos cuatro derivadas distintas, una al principio del intervalo, otra al final y dos en el medio. Después lo promediamos todo y obtenemos algo que se parece mucho a lo que realmente debiera ser. En una palabra, nos estamos quedando con más términos de la série de Taylor y por tanto cometemos menos error al aproximar nuestra función.

En un RK4 hay que evaluar en cada intervalo lo siguiente:

Y gráficamente echadle ojo a este diagramilla para ver como nuestro sistema pasa de un estado Yn a otro Yn+1 según este método:

Tal y como os comentaba se evaluan4 derivadas que se promedian para obtener un valor final creible.

Pero paraaaaa!!!...Oscar paraaaa el carro que nos vas a matarrrr!!!...¿pero a que viene todo este megarollo?...vale, un pequeño mini ejemplillo para que veais que haremos con todo esto.

Ejemplillo: Queremos modelar un trozo de tela. Modelamos la tela como una malla poligonal de triangulos interconectados entre si. Bien pues ahora en cada vértice de cada triangulo ponemos una particula con su masa, posición XYZ inicial, velocidad XYZ inicial, y entonces hacemos que todos los lados de cada triangulo sean muelles con una determinada rigidez y longitud en reposo, es decir, le damos valor a las constantes de cada uno de los resortes.

Ahora aplicamos fuerzas a las partículas y, ¿qué ocurre?..pues que como bien imaginais todo empieza a moverse ya que está interconectado y si lo que hacemos es crear una fuerza de gravedad y dejamos caer el "paño" desde 10 metros, entonces tendremos que el estado del sistema variará contínuamente debido a la gravedad.

Cada partícula tiene masa y posición/velocidad inicial. Al haber gravedad hay aceleración ¿no?, y si hay aceleración varía la velocidad, y si varía ésta varía la posición, y si esto ocurre nos encontramos con un modelo 3D de un trozo de tela que se precipita desde 10 metros hasta el suelo!!!!!...tachánnnnn...ahí lo teneis!!!...

¿Y quién calculará en cada momento las variaciones en la aceleración de cada partícula debidas a las fuerzas existentes en la escena?...pues precisamente un método numérico como Euler o el RK4.

De hecho aunque el algoritmo pueda parecer intratable, no es demasiado complejo de codificar. Introducimos los valores de las variables independientes de nuestro sistema y obtenemos los nuevos valores tras "mover" todo éste durante un tiempo h.

En otro artículo hablaré de paso adaptativo para que veais como variando el valor de h en cada iteración podemos conseguir minimizar el error contínuamente. Más adelante que por ahora ya os he mareado demasiado...

A por las particulitas....

Bueno pues vamos a aplicar las matemáticas a nuestra idea de las partículas y los muelles. Como os decía podemos construir un objeto como un entramado de partículas conectadas por muelles y después aplicar fuerzas, dejando que las "leyes de la física" hagan el resto del trabajo.

Según la popular ley de Newton, una partícula se gobierna según la fórmula f = ma, que puede expresarse también como x'' = f/m. Es decir, que si conocemos la fuerza que actua sobre una partícula y la masa de ésta, ya sabemos cual es el valor de la aceleración, es decir, de la segunda derivada de la posición respecto del tiempo.

Pero esto no es una EDO de primer orden sino de segundo!!!...¿y eso como lo solucionamos?, pues la convertimos en dos de primer orden y solucionaremos cada una de ellas con uno de los métodos comentados anteriormente. Por tanto a partir de x'' = f/m, podemos deducir:

v' = f/m

x' = v

lo cual es exactamente lo mismo y nos ha simplificado el problema de una ecuación de 2º orden a dos de orden 1, resolubles ya con lo que sabemos.

Por tanto para cada partícula de nuestro objeto, de nuestro sistema de partículas, deberemos evaluar estas dos ecuaciones cada vez que iteremos de nuevo todo el sistema. Por tanto para cada partícula tenemos algo así como 6 valores muy importantes a calcular. Es lo que se llama el "phase space" de una partícula. Consiste en un vector de la forma:

   

para cada una de las n partículas que tiene el sistema. N partículas implica evaluar n vectores como este a cada iteración.

Centrémonos más en las particulillas...

Las partículas tienen masa, posición y velocidad como atributos principales. Éstas deberían reaccionar fielmente a las fuerzas que se les aplican. Por tanto una buena simulación no sólo las implica a ellas sinó también a los objetos-fuerza definidos en el sistema.

Podemos añadir algunos atributos a una partícula como un identificador, un ancla (según el valor del ancla sabemos si la partícula es libre, se mueve, o es estática, está anclada y no debe moverse ) y también un acumulador de fuerza para ir añadiendo a éste todas las fuerzas que actuan sobre ella y al final calcular la aceleración que tiene tomando este "sumatorio" como referencia.

Nuestros objetos a simular seran sistemas de partículas. Estaran formados por una lista de partículas, una lista de muelles de conexión, un identificador, una lista de fuerzas, incluso podemos añadirles unos tiempos iniciales y finales de simulación así como el incremento de tiempo a usar ( h ) en cada iteración.

Digamos que algo como en la figura:

Por tanto para simular uno de nuestros sistemas de partículas empezamos por el "tiempo inicial" e iteramos usando el incremento h cada vez hasta llegar al "tiempo final". Por tanto todas las partículas variarán su "phase space" a cada iteración según las fuerzas que actuen sobre ellas.

Hablemos de fuerzas...

La dinámica, el movimiento implica aplicar fuerzas sobre un cuerpo. En nuestro caso tenemos objetos construidos por partículas conectadas entre sí por muelles.

Tenemos que mantener una lista de objetos-fuerza para cada sistema de partículas definido. Cada objeto-fuerza afectará a algunas partículas o incluso a todas. Cada vez que el sistema itere tendremos que aplicar cada fuerza a su grupo de partículas asociado.

Mirad la figura para verlo gráficamente:

Hablemos de gravedad...

La gravedad es muy fácil de implementar y simular. Según la 2a ley de Newton:

f = ma

y si conocemos la constante de gravitación g, podemos decir que:

fi = mi · g

dónde el valor de g para las tres coordenadas XYZ depende del lugar físico de simulación y mi es la masa de la partícula i. En principio el vector g apuntará verticalmente hacía abajo pero podemos variar esta característica para conseguir comportamientos distintos a los habituales.

Si todas las partículas estan afectadas por la gravedad tan sólo tendremos que recorrelas y añadir esta fuerza a sus acumuladores de fuerza.

Podemos tener diferentes gravedades en el mismo sistema de partículas, en la misma simulación. Tan sólo asociaremos cada una a un conjunto de partículas diferente.

Hablemos de fricción...

Una fuerza de fricción ideal tiene la forma:

fi = -kf · vi

dónde kf es el coeficiente de fricción de la fuerza en cuestión y vi es el vector de velocidad de la partícula afectada. La partícula se verá más o menos afectada según sea su velocidad.

El efecto que produce la fricción es el de restar velocidad, evitar el movimiento, hasta parar a la partícula si ésta no recibe otras influencias. El mismo viento es un buen ejemplo de este efecto.

Es recomendable aplicar una cierta fricción al sistema no sólo por los efectos que podamos conseguir sinó también para evitar problemas de estabilidad numérica.

Fuerzas de tipo muelle

¿Cómo podemos modelar una cierta afinidad entre partículas?, pués por ejemplo siguiendo la popular ley del resorte, ampliamente estudiada por casi todos nosotros en nuestros estudios secundarios en el instituto. Todo puede conectarse usando muelles para ello y podemos conseguir interesantísimos efectos si procuramos valores correctos para las constantes de cada resorte.

La fuerza que se ejerce sobre dos partículas a y b conectadas por un muelle puede escribirse así:

dónde fa y fb es la fuerza en cada partícula, l equivale al valor pa-pb (posición de a menos posición de b), l' es la derivada de l y por tanto equivale a va-vb (velocidad de a menos velocidad de b), ks es la constante del muelle y kd es la constante de pérdidas por rozamiento de éste. El valor de r se corresponde con la longitud en reposo del muelle.

Podemos observar, y de hecho la misma lógica ya nos lo indica, que la fuerza que ejerce el muelle es proporcional a la distancia entre ambas partículas. El rozamiento es proporcional a la velocidad con la que ambas se acercan/alejan. Ambas artículas reciben exactamente la misma cantidad de fuerza pero con signos contrarios.

Implementar un muelle con pérdidas es muy simple. Simplemente tendremos que crear una estructura de tipo muelle que guarde información sobre sus características internas (longitud en reposo y constantes) así como dos punteros a las dos partículas que interconecta. Para aplicarlo, buscamos las partículas, calculamos la expresión anteriormente mostrada y añadimos los valores de las fuerzas calculadas a los acumuladores de cada una.

Según el valor de las constantes conseguiremos que nuestros muelles sean más o menos rígidos lo cuál nos conduce a conseguir mucha flexibilidad en el comportamiento final del sistema.

Aquí tenéis uno de los ejemplos de modelo deformable que he programado. En la simulación de mi "PSF system" (Particle System File system) he introducido gravedad y fricción, y el modelo está creado a partir de partículas y muelles interconectados.

La secuencia debe observarse en este orden:

 1er frame 2o frame 3er frame
6o frame 5o frame 4o frame

Ahí va!!

   
     
 En este caso se ha programado también un sistema de detección de colisiones con lo cuál el pendulo rebota en el suelo, perdiendo así energía, al colisionar.

¿Lo váis viendo?, espero que sí porque de verdad que vale mucho la pena meterse en estos temas...muy duro al inicio pero super reconfortante una vez ya estás de pleno !!!

Bueno para más información conectaros a la página de un genio de la simulación dinámica de objetos deformables, David Baraff. Allí tiene su curso on-line sobre simulación dinámica que presenta cada año en SIGGRAPH con el cual yo he aprendido muchísimo.