Indice | Previo: ProblemaNCuerpos.SolucionNumerica | Siguiente: Problema2Cuerpos
6.4.5. Constantes de movimiento y teorema del virial¶
Una interesante primera "aplicación" de la solución numérica al problema de los N cuerpos vista en estas secciones, es la de verificar "experimentalmente" los resultados analíticos descritos en la Sección ¿Solución analítica? y en la Sección Energía y virial.
Para ello estudiaremos un sistema de 5 partículas con masas, posiciones y velocidades generadas al azar. La solución numérica a las e.d.m.r. del sistema, obtenida con los métodos vistos en esta sección, nos permitirá obtener las listas de sus posiciones y velocidades para distintos valores del tiempo. Con estas listas podremos calcular y graficar los valores de las constantes de movimiento, momento lineal, momento angular, energía y de las cantidades críticas para el teorema del virial.
Comencemos pues por generar las condiciones iniciales del sistema usando, entre otras cosas, la rutina sistema_a_Y del Alg. (sistema_a_Y) y las rutinas de generación de números aleatorios que usamos en la Sección Centro de masa:
#Importa todos las rutinas y constantes de otras libretas
from pymcel.export import *
#Número de partículas
N=5
#Generación de las condiciones para cada partícula
from numpy.random import uniform,seed
seed(7)
#Condiciones iniciales
sistema=[]
for i in range(N):
particula=dict(
m=uniform(0.0,1.0),
r=uniform(-1.0,1.0,size=3),
v=uniform(-1.0,1.0,size=3)
)
sistema+=[particula]
N,mus,Y0s=sistema_a_Y(sistema)
#Tiempos
from numpy import linspace
Nt=100
ts=linspace(0.0,10.0,Nt)
print(f"N = {N}")
print(f"mus = {mus}")
print(f"Y0s = {Y0s}")
Ahora podemos resolver las ecuaciones de movimiento y extraer las posiciones y velocidades de las partículas:
#Solución
from scipy.integrate import odeint
solucion=odeint(edm_ncuerpos,Y0s,ts,args=(N,mus))
#Extracción de las posiciones y velocidades
rs,vs=solucion_a_estado(solucion,N,Nt)
Una gráfica del movimiento de las partículas en tres dimensiones se puede obtener usando la rutina plot_ncuerpos_3d que definimos en el Alg. (plot_ncuerpos_3d):
fig=plot_ncuerpos_3d(rs,vs,marker='.');
Ahora podemos calcular las constantes de movimiento. En este caso, sin embargo, la dificultad algorítmica estriba en que las posiciones y velocidades de las partículas están guardadas en las matrices rs y vs que no son triviales de manipular.
Así, por ejemplo el momento lineal inicial de la partícula 0 esta dado por:
p_0_0=mus[0]*vs[0,0,:]
Pero si queremos el momento lineal de esa partícula en cualquier tiempo será:
p_0=mus[0]*vs[0,:,:]
La cuadratura de momento lineal total $C_{P_\mathrm{CM}}$ la podemos obtener si sumamos uno a uno los momentos lineales en cada tiempo de todas las partículas del sistema:
from numpy import zeros
C_PCM=zeros((Nt,3))
for i in range(N):
C_PCM=C_PCM+mus[i]*vs[i,:,:]
Y como vemos el valor del momento lineal es el mismo, que es lo que esperabamos de acuerdo con la teoría.
Por otro lado la cuadratura de momento angular será:
from numpy import zeros,cross
C_L=zeros((Nt,3))
for i in range(N):
C_L=C_L+mus[i]*cross(rs[i,:,:],vs[i,:,:])
Que de nuevo resulta constante como esperabamos.
Finalmente la cuadratura de energía se puede calcular usando la fórmula para la energía potencial la dada por la Ec. (ncuerpos_potencial):
from numpy import zeros
from numpy.linalg import norm
C_E=zeros(Nt)
K=zeros(Nt)
U=zeros(Nt)
for i in range(N):
K=K+0.5*mus[i]*norm(vs[i,:,:],axis=1)**2
for j in range(N):
if i==j:continue
rij=norm(rs[i,:,:]-rs[j,:,:],axis=1)
U+=-0.5*mus[i]*mus[j]/rij
C_E=K+U
Como vemos el valor de la energía es negativo, lo que podría implicar que el sistema es ligado (tal y como sugieren las trayectorias de las partículas.) Sin embargo, como mencionamos en la Sección Virial la condición $E<0$ es necesaria más no suficiente. Para saber si el sistema es ligado debemos evaluar los promedios a largo plazo de las energía cinética y potencial y compararlas con la energía total:
E=C_E[0]
Kmean=K.mean()
Umean=U.mean()
print(f"-E = {-E}")
print(f"<K> = {Kmean}")
print(f"-<U>/2 = {-Umean/2}")
Como vemos la identidad $\langle K\rangle=-E=-\langle U\rangle/2$ se cumple aproximadamente para la ventana de tiempo en la que estudiamos el sistema y podríamos sospechar que es estable.
6.4.6. Una algoritmo general¶
Usando lo visto en esta y en las secciones anteriores, podemos construir ahora un algoritmo general que nos servirá en lo sucesivo para partiendo de la descripción de un sistema de N cuerpos obtener las posiciones y velocidades de las partículas que lo constituyen, resolviendo numéricamente las e.d.m.r. para un determinado conjunto de tiempo provistos.
El algoritmo presentado a continuación define esa rutina. Para ello usa las rutinas definidas en los Algs. (edm_ncuerpos, sistema_a_Y, solucion_a_estado) y los procedimientos para el cálculo de las constantes de movimiento presentados en la Sección Constantes de movimiento y teorema del virial.
def ncuerpos_solucion(sistema,ts):
#Condiciones iniciales
from pymcel.export import sistema_a_Y
N,mus,Y0s=sistema_a_Y(sistema)
#Masa total
M=sum(mus)
#Número de tiempos
Nt=len(ts)
#Solución
from scipy.integrate import odeint
solucion=odeint(edm_ncuerpos,Y0s,ts,args=(N,mus))
#Extracción de las posiciones y velocidades
from pymcel.export import solucion_a_estado
rs,vs=solucion_a_estado(solucion,N,Nt)
#Calcula las constantes de movimiento
from numpy import zeros
PCM=zeros(3)
for i in range(N):
PCM=PCM+mus[i]*vs[i,0,:]
#Posición del CM como función del tiempo
RCM=zeros((Nt,3))
for i in range(N):
RCM=RCM+mus[i]*rs[i,:,:]
RCM/=M
#Momento angular
from numpy import zeros,cross
L=zeros(3)
for i in range(N):
L=L+mus[i]*cross(rs[i,0,:],vs[i,0,:])
#Energía total
from numpy.linalg import norm
K=zeros(Nt)
U=zeros(Nt)
for i in range(N):
K=K+0.5*mus[i]*norm(vs[i,:,:],axis=1)**2
for j in range(N):
if i==j:continue
rij=norm(rs[i,:,:]-rs[j,:,:],axis=1)
U+=-0.5*mus[i]*mus[j]/rij
E=K[0]+U[0]
#Constantes
constantes=dict(M=M,
RCM=RCM,PCM=PCM,
L=L,K=K,U=U,E=E)
#Posiciones y velocidades relativas al centro de masa
from numpy import subtract
rps=rs-RCM
vps=subtract(vs,PCM/M)
#Devuelve las posiciones y velocidades
return rs,vs,rps,vps,constantes
La rutina se invocaría así:
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,ts)
Una representación gráfica de las posiciones de las partículas en el sistema de referencia del centro de masa:
fig=plot_ncuerpos_3d(rps,vps,marker='.')
Indice | Previo: ProblemaNCuerpos.SolucionNumerica | Siguiente: Problema2Cuerpos