Cómo hacer simulaciones de partículas en MatLab. III: Órbita gravitacional (1 Part. – 2D)

< II. Tiro parabólico (1 part. – 2D)

 

Hasta ahora en los dos tutoriales anteriores hemos visto en qué consiste teóricamente la simulación de un sistema dinámico y un ejemplo de simulación de una sola partícula en un espacio de dos dimensiones (2D). En esa última simulación la fuerza sobre la partícula (la gravedad en ese caso) era constante. En este y los siguientes tutoriales continuaré poniendo ejemplos de otras combinaciones de número de partículas, dimensiones espaciales, y tipos de fuerzas. Como esto puede sonar un poco lioso, he hecho una tabla con todas las combinaciones posibles. El número que aparece en las casillas corresponde al número del tutorial (el tutorial I no aparece porque era solo teoría).

 

F = const. F (x) F (v) F (x,v)

1D

2D

II (1)

III (1)

3D

*(1): Una sola partícula. (N): Un sistema de N partículas independientes. (C): Colisiones.

 

En esta ocasión, como se ve en la tabla,  veremos cómo hacer una simulación donde la fuerza sobre la partícula dependa de la posición de la misma. En general, recordemos, la fuerza puede depender de las dos variables de estado del sistema: posición y velocidad. De nuevo en este ejemplo se usará un sistema con una sola partícula, y también en el plano (2D).

 

Ley de gravitación universal

Como ejemplo de fuerza que depende de la posición, vamos a usar la gravedad. En la simulación del segundo tutorial considerábamos la gravedad como una constante, y se puede considerar razonable esa simplificación si nuestra simulación es de un sistema que se encuentra cerca de la superficie de la Tierra y es relativamente de pequeño tamaño. Pero si quisiéramos simular el comportamiento de un ascensor que llegue hasta la estratosfera, o el de un satélite que orbita la Tierra, considerar un valor constante del campo gravitatorio terrestre sería mucho simplificar.

En realidad, según dicta la Ley de la gravitación universal, dos cuerpos con masa se ejercen una fuerza entre sí (del mismo valor y sentido opuesto) que depende del cuadrado la distancia (d12) entre sus centros de masas y de la masa (m1 y m2) de cada cuerpo [1].

 

Gravedad

 

La fuerza que ejerce la partícula 2 sobre la 1, como se ve en la imagen, es un vector (la flecha roja) que tiene la dirección de la recta que une los dos puntos (la misma que el vector distancia, flecha gris), con el sentido hacia la otra partícula (la gravedad no repele, sino atrae), y su módulo calculado según la Ley de Gravitación Universal que aparece en la imagen. Es fácil comprobar que si calculamos la fuerza que hace la partícula 1 sobre la 2, será el mismo vector pero de sentido opuesto (porque ahora el vector distancia será el opuesto) y de exactamente el mismo módulo. La G es una constante de la física que vale 6.67·10-11 N·m2/kg2.

 

NOTA IMPORTANTE: En las simulaciones voy a considerar el movimiento de una sola partícula, y por lo tanto también solo la fuerza sobre esa partícula. En un caso más realista habría que considerar el movimiento de las dos partículas. Es una simplificación que hago para ser coherente con el uso de una sola partícula en estos tutoriales de momento. De todas formas en las simulaciones la masa de la “partícula estática” que he considerado es muchas veces mayor que la partícula “dinámica”, por lo que el efecto producido sobre ella es menospreciable.

 

Órbita de una partícula alrededor de otra estática

clear
figure(1)

%variables
x=[];
v=[];
F=[];

%parametros
m=1;
M=1e20;
G=6.67e-11;
h=0.01;

%condiciones iniciales
x=[-1e3;0];
v=[0;-2e3];

%variable externa
F=-G*(m*M/norm(x)^2)*(x/norm(x));

for step=1:500

hold off;
scatter(x(1),x(2), ‘black’,‘filled’);
hold on
scatter(0,0,‘black’,‘filled’);

title([‘step : ‘ num2str(step)]);
axis([-1.5e3 1.5e3 -1.5e3 1.5e3]);
pause(0.01);

%guardar valor anterior
xa=x;
va=v;

%paso integración
F=-G*(m*M/norm(xa)^2)*(xa/norm(xa));
vpm = va + (h/2)*(F/m);
xpm = xa + (h/2)*va;

F=-G*(m*M/norm(xpm)^2)*(xpm/norm(xpm));
v = va + h*(F/m);
x = xa + h*vpm;

end

 

El código tiene una estructura muy similar al del tutorial anterior, así que me limitaré a explicar solo las novedades. Consulta la entrada anterior si alguna función no aparece explicada aquí.

Parámetros: La masa m es la que vamos a animar y la masa M es la que va a quedar estática en el centro de la imagen. La masa M es mucho mayor que la m para compensar la simplificación de simular solo una partícula. G es la constante gravitatoria.

Condiciones iniciales: La posición y velocidad iniciales determinarán como será la trayectoria de la partícula. Yo he intentado conseguir una trayectoria lo más circular posible, eligiendo una velocidad tangente a la recta que une las partículas. Pero pueden probarse otros valores que producen órbitas más elípticas, o hiperbólicas (la gravedad de M hace que salga fuera de órbita).

–  Variable externa: Una pequeña diferencia con respecto a la simulación anterior es que ahora la fuerza sobre la partícula depende su posición. Como la posición de la partícula estática es el origen de coordenadas, el vector distancia es simplemente el vector posición x. Hay que tener en cuenta que como el vector distancia de la fórmula iba de la partícula que animamos hacia la estática, el vector x tiene signo opuesto (porque va del origen a la partícula animada). Por eso en el cálculo de la fuerza hay un signo menos al principio.

La función norm() calcula el módulo de un vector. Veréis que la fórmula para calcular la fuerza no es exactamente igual que la que mostré en la imagen. Eso es porque nos interesa calcular el vector fuerza (no solo el módulo). Hay que multiplicar la fórmula del módulo por un vector que sea unitario (que mida 1) y que tenga la dirección que queremos (la dirección del vector distancia, en nuestro caso x). Ese vector por el que hay que multiplicar el módulo de la fuerza es (x/norm(x)). Dividiendo un vector por su modulo siempre obtienes uno unitario de la misma dirección y sentido.

hold on/ hold off: Estas funciones permiten mantener un dibujo en una figura cada vez que se genere uno nuevo (hold on), o dejar de mantener el dibujo existente provocando que el nuevo lo borre (hold off). Los he utilizado para poder dibujar las dos partículas sin que el dibujo de la segunda me borre el de la primera.

scatter(): Otra diferencia con respecto al programa anterior es que en lugar de utilizar plot() para dibujar he utilizado esta nueva función. El motivo es que con esta función pueden editarse mejor la creación de puntos (distintos colores, tamaños, relleno de color…). Las opciones elegidas para esta simulación han sido pintar los dos puntos de negro y hacer que sean puntos con relleno. Para más opciones puedes consultar la ayuda online del programa.

– Paso de integración: Como en el tutorial anterior calculamos el valor siguiente de la posición y la velocidad a partir del método del punto medio. La diferencia es que ahora para calcular la velocidad del punto medio (vpm) la fuerza debe actualizarse primero al valor del punto anterior (no es constante como antes). Una vez calculada la posición del punto medio (xpm), se actualiza la fuerza, la posición y la velocidad.

El resultado final es algo como lo que se ve en la siguiente imagen animada. En la captura del GIF se pierde algo de fluidez, pero el programa en un ordenador normal se ve más rápido. Si tu simulación se ve lenta puedes aumentar el paso de integración h, pero no mucho porque perderás precisión y puede que estabilidad.

 

dguinstation.worspress.com

Haz clic si no ves la animación para abrirla otra ventana

 

Versión coloreada con trayectoria

La versión anterior queda un poco simple, teniendo en cuenta que podemos elegir los colores que queramos y añadir otras cosas. Lo que puede hacerse para darle un poco más de color a la simulación lo siguiente:

Pintar el fondo: Podemos elegir un color de fondo. Por ejemplo yo he elegido el negro para que parezca el espacio. Para ello usamos la siguiente línea, después de dibujar los puntos:

set(gca,’Color’,‘black’);

Editar los puntos: Para cambiar el color solo basta con cambiar las opciones de la función scatter. En lugar de negro (black), yo he utilizado azul para la tierra (blue) y amarillo para el sol (yellow). El tamaño hay que escribirlo después de las coordenadas del punto, por ejemplo:

scatter(0,0,100,‘yellow’,‘filled’);

– Dibujar la trayectoria: Para ello debemos guardar todos los puntos calculados hasta el momento en cada iteración. Los guardaremos en una matriz que declaramos al principio como el resto de variables (yo la he llamado es de “estela”). Esa matriz tendrá 2 filas y N columnas, porque cada columna será las dos coordenadas de un punto calculado hasta el momento. Podemos dibujar los puntos de la trayectoria después del punto móvil. Hay que comprobar que ya haya al menos un punto anterior calculado (que es no esté vacío). Usaremos el plot ya que es lo adecuado para hacer líneas entre puntos.

if (size(es,2)>0)
plot(es(1,:),es(2,:));
end

La parte es(1,:) significa que queremos la fila 1 de es y  todas las columnas (:), es decir, que obtendremos toda la primera fila (todas las coordenadas horizontales). Lo mismo para las coordenadas verticales.

Por último es necesario actualizar la matriz es cada vez que se calcula un punto nuevo, guardando así el anterior. Así que una vez que hemos calculado un nuevo punto hay que hacer:

es = [xa es];

Eso guarda el punto anterior xa al principio de la matriz. Pero realmente el orden no importa porque luego el plot lo dibujara todo de golpe.

El resultado es algo como la siguiente simulación:

 

Tutorial_02c

Haz clic si no ves la animación para abrirla otra ventana

 

Para el que quiera un reto

Se pueden hacer muchas más cosas que animo a probar. Por ejemplo simplemente variando las condiciones iniciales un poco esta última simulación con la trayectoria puede dar lugar a dibujos muy bonitos, como este:

 

2c

Órbita con condiciones iniciales diferentes

 

También pueden hacerse otros cambios con la estela. Por ejemplo yo he dibujado una estela de un número finito de puntos, con un degradado de color del azul al blanco:

 

Tutorial_02b

Haz clic si no ves la animación para abrirla otra ventana

 

En la siguientes entradas probaré otras combinaciones posibles de la tabla que puse al principio. Aún no lo tengo decidido pero seguramente será un ejemplo de péndulo con fricción, por lo que se tratará de una partícula, en 2D, y con la fuerza que depende de la posición y la velocidad.

 

 

[1] Que la gravedad es una fuerza que actúa “a pares” es algo que se olvida a menudo. Por ejemplo, no es del todo correcto decir que la Tierra gire alrededor del Sol. Lo que Newton demostró es que el sistema Tierra-Sol tiene un centro de masas (casi en el centro de Sol, porque tiene mucha más masa que la Tierra), y que tanto el Sol como la Tierra orbitan ese centro de masas. Pero como nuestras simulaciones son de partículas (que son puntos), no debemos preocuparnos por encontrar centros de masas del sistema. El centro de masas de una partícula está en el mismo lugar que ocupa la partícula.

IV. Péndulo (1 Part. – 2D) >


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

22 comentarios

  1. Gracias por compartir esta explicación tan clara! Sin embargo tengo la duda con la trayectoria. Dices que creas una matriz “es”. Esta parte del proceso no me queda aún muy clara. Podrías por favor detallar este paso? Mil gracias.

    1. Al principio (donde pone %variables por ejemplo) puedes inicializar la variable haciendo:

      es=[ ];

      Para dibujar hay que hacer lo siguiente (lo puedes poner después de los otros dibujos):

      if (size(es,2)>0)
      plot(es(1,:),es(2,:));
      end

      La primera vez no dibujará nada porque aun no hay puntos anteriores. Pero lo importante es como se actualiza ‘es’. Lo puedes pones al final del for:

      es = [xa es];

      Como ‘xa’ es un vector columna con dos números (coordenada ‘x’ e ‘y’), lo que hace la línea anterior es simplemente ampliar la matriz es con una nueva posición. Por ejemplo supongamos que el primer punto es [1;2], entonces:

      es =

      1
      2

      si el segundo punto es [3;4], entonces:

      es =

      3 1
      4 2

      Y así se va ampliando ‘es’ con posiciones ‘anteriores’. Espero haberte ayudado.

  2. Gracias por responder! Sin duda para este caso ya está claro. Pero y si quisieramos incluir más cuerpos, se necesitarían dos matrices x y y que guarden la respectiva cordenada de cada cuerpo? Si es así como le decimos que guarde cada coordenada en dichas matrices por cada timestep? Gracias de antemano si puedes aclararme.

    1. Yo lo que hago es guardar las dos coordenadas juntas en un vector, en lugar de guardarlas por separado. Por eso los vectores de posicion y velocidad tiene dos componentes. Si tienes más partículas las posiciones y velocidades se pueden guardar en una matriz 2xN donde N es el número de partículas. Tendría la misma forma que la matriz ‘es’ que dije antes pero ahora cada columna sería la posicion o velocidad de una particula.

      Para la estela con más de una partícula puedes definir diferentes matrices, por ejemplo ‘es1’ y ‘es2’ si hay dos partículas. También podrían guardarse en una matriz 3D, pero puede que eso sea más complicado.

  3. Gracias por tu aclaración! :)

  4. Muy bueno, me funciona. Mil gracias! C:

¿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: