Cómo hacer simulaciones de partículas en Matlab. I: Introducción

 

Desde que publiqué mis ejemplos de simulaciones en Matlab he ido recibiendo muchos comentarios de personas que necesitaban ayuda un poco más detallada que la que mostraba en esas entradas. Aquellos post no eran tutoriales, aunque explicara más o menos en que se basaban las simulaciones, así que no sirven para saber detalladamente cómo hay que hacerlas, pero sí para saber algunas cosa que pueden hacerse en Matlab. En esta serie de entradas lo que me propongo es hacer una guía muy básica sobre las simulaciones, en principio, solo de partículas. Más adelante quizás lo amplíe a cuerdas o superficies.

Los ejemplos que utilice en estos tutoriales deben interpretarse como una plantilla muy básica. Los programas se pueden complicar y ampliar tanto como desee el programador, y ni mucho menos explicaré las infinitas posibilidades que hay de simular partículas. Lo que mostré en la entrada de simulación de partículas eran algunos ejemplos que se me ocurrieron a mí, como los fuegos artificiales o la colisión con esferas desde dentro, pero las posibilidades son ilimitadas. De todas formas, si tienes alguna idea nueva y no sabes como enfocar la solución, deja un comentario y quizás pueda ayudarte.

Antes de empezar también quiero comentar que doy por supuesto que se sabe usar Matlab a un nivel básico (cálculos con vectores y matrices, uso de variables…). Utilizar la ayuda del programa puede ser muy útil cuando se empieza (escribiendo help en la ventana de comandos) o algún manual como este [PDF].

Antes de empezar a escribir el código tenemos que saber qué estamos simulando y en qué consiste una simulación. Este primer post será todo teoría y no hay nada de código, pero es importante para saber qué estamos haciendo.

 

Simular un sistema dinámico

Lo que queremos es hacer una representación gráfica de un sistema cuyas propiedades van cambiando con el tiempo. Un sistema puede ser mecánico, eléctrico, biológico, termodinámico, económico… Cada sistema tiene unas variables que definen el estado en el que se encuentra en cada momento ese sistema (las variables de estado). En un sistema mecánico pueden ser la velocidad y la posición, en uno eléctrico la corriente y el voltaje, en uno biológico la concentración de una substancia en una disolución, en uno termodinámico la presión y la temperatura, etc. Esas son las propiedades que se representan gráficamente, de una forma u otra. Cuando se simula temperatura, puede usarse  una escala de colores. Cuando se simulan partículas, pueden dibujarse éstas para indicar su posición, y hacer cambiar esa posición según su velocidad.

En nuestro caso, lo que nos interesa estudiar es el comportamiento mecánico de las partículas. No nos interesa, por ejemplo, su comportamiento eléctrico o termodinámico, solo su movimiento (cinemática) y la causa de ese movimiento (dinámica). La cinemática y la dinámica son las componentes principales de la mecánica, y están relacionadas según las leyes de Newton.

 

Las leyes de la mecánica

Según la Segunda Ley de Newton, la variación del momento cinético de una partícula en un instante es igual a la suma de fuerzas sobre la partícula en ese instante. El momento cinético (p) es el producto de la masa (m) de la partícula por su velocidad (v). Como la velocidad es un vector, el momento cinético también lo es. También consideramos que la masa de la partícula es constante. Así, la Segunda Ley de Newton queda (las negritas son vectores):

 

image

 

La derivada (d()/dt) es la herramienta matemática para describir la variación instantánea que comentaba en la definición de la Segunda Ley de Newton. No es importante saber derivar para hacer la simulación, pero sí entender que es una variación de una función respecto al tiempo, en un periodo de tiempo muy pequeño. La fuerza puede depender de la posición de la partícula (como en un muelle) o de su velocidad (si hay rozamiento con el aire, por ejemplo). También podemos suponer que la fuerza sobre la partícula es una constante (como un campo gravitatorio que sea uniforme), en ese caso no depende ni de la posición ni de la velocidad. Todo esto dependerá de la complejidad que le queramos dar a nuestra simulación.  

 

Discretización temporal

Ahora sabemos que, para cada instante de tiempo, la posición cambia según lo indique la velocidad (la segunda ecuación de las anteriores), y la velocidad cambia según indican la fuerza sobre la partícula y la masa (la Segunda Ley de Newton). Pero si queremos que un ordenador calcule TODAS las posiciones y velocidades por las que pasa nuestra partícula entre dos instantes, tenemos un problema. Imaginaos que en un experimento real, observamos lo que pasa en un intervalo de tiempo, pongamos, entre el segundo 0 y el segundo 5. ¿Cuántos instantes de tiempo hay en medio? La respuesta es infinitos. Por tanto, ¿cuántas velocidades y posiciones intermedias hay entre los dos extremos de ese intervalo? También infinitas. Lo que hay que tener claro con las simulaciones es que el ordenador siempre calcula un número finito de estados del sistema. Es decir, hay que considerar el tiempo como una variable discreta en lugar de continua.

En una simulación, por lo tanto, no calculamos las infinitas posiciones de la partícula, aunque en ocasiones podamos hacerlo solucionando las ecuaciones diferenciales anteriores. Eso nos daría una función de la posición respecto al tiempo, que luego podríamos representar. Pero eso no es una simulación dinámica. Nuestra simulación consiste en el cálculo de las posiciones y velocidades para un número finito de instantes de tiempo de forma recurrente, es decir, calculando los valores del instante siguiente a partir del instante anterior. Tiene sentido hablar a partir de ahora de instantes “siguientes” y “anteriores” porque el tiempo es una variable discreta, no continua.

Para hacer esos cálculos de las variables de estado siguientes a partir de las anteriores (velocidades y posiciones) hacemos uso de los métodos numéricos (en contraposición al método analítico que nos daría la solución explicita de la posición respecto al tiempo). Como hemos visto, ahora el tiempo no es una variable continua, sino discreta, por lo que tampoco veremos un cambio del resto de variables (velocidad y posición) de forma continua. En lugar de eso, solo tendremos un número finito de valores en nuestro intervalo de tiempo en el que dure la simulación. Es como si hiciéramos fotografías cada cierto tiempo para ver donde se encuentran las partículas y qué velocidad tienen, y a partir de esas posiciones y velocidades calculásemos la de la siguiente fotografía. La siguiente imagen ayudará a entenderlo en el caso de considerar solo dos dimensiones espaciales. 

 

simulacion_paso

 

Como no pasamos por los infinitos instantes de tiempo, sino que vamos dando “saltos”, debemos saber cómo de grandes son esos saltos. Al tamaño de esos intervalos de tiempo le llamaremos h. Si h es de 1 segundo, por ejemplo, significa que calcularemos las posiciones y velocidades de 1 segundo en 1 segundo. Cuanto más pequeño sea ese paso, más instantes estaremos considerando y más precisión conseguiremos. Pero eso también significará más tiempo de cálculo y otros posibles problemas de estabilidad de los métodos numéricos que utilicemos.

Durante ese intervalo de tiempo h se podría decir que estamos ciegos. Solo vemos lo que pasa en los instantes inicial y final, pero no en medio. Imaginemos que en el instante inicial (t1) conocemos la posición ([x1,y1]) y la velocidad (v1). De momento no quiero explicar cómo calcular las velocidades, simplemente estoy usando el ejemplo del cálculo de la posición para comentar algo muy importante. Así que pensemos que de alguna forma (ya lo veremos en adelante) también sé la velocidad del instante t2, es decir, conocemos v2. Sabemos por tanto la velocidad inicial y final, pero no las infinitas velocidades que tiene la partícula durante todo el intervalo h. Aquí viene la gracia de los métodos numéricos: hay que suponer que alguna variable que cambia de forma continua con el tiempo, es constante al menos en un intervalo pequeño de tiempo.

Obviamente eso no será verdad, pero es la simplificación que hacemos para poder calcular el cambio de posición y velocidad. Hay que suponer que esos cambios se producen a un ritmo (la derivada) constante durante el intervalo h. Es decir, durante ese intervalo, la velocidad (a la que le he puesto el nombre de v1,2) se supone constante. Y lo mismo se supone para la aceleración, que es la derivada de la velocidad.

 

Ecuaciones de estado discretas

Podemos reescribir las dos ecuaciones escritas anteriormente (la Segunda Ley de Newton y la definición de velocidad), pero ahora considerando la discretización temporal. Recapitulemos un poco. Ahora avanzar en el tiempo consiste en dar saltos de duración h (el incremento de tiempo t2-t1). Para calcular la posición siguiente supondremos que la velocidad en el intervalo h es constante. Y para calcular la nueva velocidad supondremos una aceleración también constante en ese intervalo. Si escribimos las derivadas (que son relaciones entre variaciones infinitesimales de dos magnitudes) en forma de incrementos (valor final – valor inicial) en el intervalo h, obtenemos las siguientes ecuaciones:

 

image

 

Conociendo la posición y velocidad iniciales, siempre podremos calcular las “siguientes” posición y velocidad. Para ello necesitaremos dos cosas. En el cálculo de la velocidad, necesitamos saber qué aceleración constante considerar en el intervalo de tiempo h. En general, esa aceleración dependerá (aparte de la masa) de la posición y la velocidad. La segunda cosa que necesitamos, esta vez para calcular la posición siguiente, es qué velocidad constante vamos a considerar para el intervalo h.

La única diferencia entre un método numérico y otro para calcular la posición final radica precisamente en el criterio que escojamos para decidir esas constantes. Podríamos escoger por ejemplo el valor de velocidad en el instante inicial del intervalo, v1 (método de Euler explícito), o podríamos escoger el valor en el instante final v2 (método de Euler implícito). O quizás podríamos escoger la media aritmética de las dos velocidades (Método del valor medio).

El método que utilizaré en adelante es el que se conoce como método del Punto Medio, y consiste en escoger para el intervalo la velocidad que corresponde al punto medio que se encuentra entre la posición inicial y la final. Esa posición y velocidad intermedias se calculan con el método de Euler Explícito suponiendo que pasa un tiempo h/2 (la mitad del intervalo). La posición del punto medio xpm es:

 

image

 

Y la velocidad vpm que corresponde a ese punto medio es:

 

image

 

La posición y la velocidad calculadas para el punto medio con el método de Euler Explícito serán las que usaremos para el cálculo de la aceleración y la velocidad en el intervalo h. Así podemos calcular el punto final del intervalo (x2, v2).

 

image

 

En resumen,

– Simular un sistema dinámico consiste en representar gráficamente las variables de estado de un sistema que cambia con el tiempo. Nuestro sistema será una partícula (o un conjunto independiente de partículas). Si nos interesa el comportamiento mecánico del sistema, necesitamos calcular las posiciones y velocidades de las partículas.

– Según la Segunda Ley de Newton, lo que determina la variación de la velocidad en el tiempo de la partícula (la aceleración) es la fuerza que actúa sobre la partícula dividida entre su masa. La velocidad por su parte es por definición la variación de la posición respecto al tiempo.

– Aunque en el mundo real las variables de estado cambian de forma continua, un ordenador solo puede representar un número finito de estados. Es necesario discretizar el tiempo, es decir, considerar que solo calculamos un número finito de estados cada cierto tiempo h. El calculo se hace de forma recurrente: calculando el estado "siguiente” a partir del “anterior”.

– Para cada estado conocido, debe hacerse una predicción del estado siguiente del sistema (la siguiente posición y velocidad). Para ello se usan los métodos numéricos. El método de integración que usaré en adelante es el del Punto Medio, que consiste en predecir un estado intermedio mediante el método explícito, y luego utilizar esas variables de estado intermedias para hacer la predicción del estado final del intervalo.

 

Hasta aquí la primera parte de esta serie de entradas, de momento muy teórica. A partir de la siguiente entrega veremos como implementar una simulación sencilla con el programa Matlab.

 

 

II. Tiro parabólico >


Esta entrada pertenece a la colección Animación con MatLab

4 comentarios

  1. Información Bitacoras.com

    Valora en Bitacoras.com:   Desde que publiqué mis ejemplos de simulaciones en Matlab he ido recibiendo muchos comentarios de personas que necesitaban ayuda un poco más detallada que la que mostraba en esas entradas. Aquellos post no eran tutoriales, ..…

  2. […] Colección: Tutorial de Simulación de partículas en MatLab (Parte I, Parte II, Parte III) Debido casi a la petición popular, este año comencé una serie de entradas […]

  3. hola, estamos intentando hacer un proyecto con partículas que interaccionan de manera diferente según la temperatura a la que se encuentran,
    ¿podría ayudarnos?
    gracias

¿Qué te parece?

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

A %d blogueros les gusta esto: