#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
The N-body class dynamical system
================
:Example:
>>> import nbodypy.dynamics
This is a subtitle
-------------------
"""
import nbodypy
import numpy
# the N-body system
[docs]def system(z,t,m,N,dim,rotating):
"""
The vector field for the Nbody problem
coordinates are organized as (r1,v1,r2,v2,..,rN,vN)
where r is the position in R^n where n=dim/2
and v the velocity in R^n
:param N: (int) number of bodies
:param dim: (int) dimension of the phase space for one body
:param z: (numpy array) state at which we want to compute the vector field. Dimension : dim*N
:param t: (float) time at which we want to compute the vector field. Here the vector field is time independant
:param rotating: (numpy array) rotation matrix for the frame. dimension dim/2 times dim/2
:return: (numpy array) dz/dt the value of the vector field at state z (dimension N*dim)
"""
dzdt = numpy.zeros(N*dim)
for i in range(N):
#print "Point i", i, z[(i*dim+2):(i*dim+dim)]
qi = z[(i*dim):(i*dim+dim/2)]
vi = z[(i*dim+dim/2):(i*dim+dim)]
dzdt[i*dim:(i*dim+dim/2)] = vi # dr/dt=v
for j in range(N):
qj = z[(j*dim):(j*dim+dim/2)]
vj = z[(j*dim+dim/2):(j*dim+dim)]
if(i!=j):
rij = numpy.linalg.norm(qj-qi)
dzdt[(i*dim+dim/2):(i*dim+dim)] = dzdt[(i*dim+dim/2):(i*dim+dim)]+ (m[j]*(qj-qi))/numpy.power(rij,3.0)
# if rotating frame
if(isinstance(rotating,numpy.ndarray)):
dzdt[(i*dim+dim/2):(i*dim+dim)] = dzdt[(i*dim+dim/2):(i*dim+dim)]-2.0*numpy.dot(rotating,vi)-numpy.dot(rotating,numpy.dot(rotating,qi))
return dzdt
# the N-body system resolvant equation
[docs]def Dsystem(zdz,t,m,N,dim,rotating):
"""
The vector field for the Nbody problem
coordinates are organized as (r1,v1,r2,v2,..,rN,vN) and (dr1,dv1,dr2,dv2,..,drN,dvN)
where r is the position in R^n where n=dim/2
v the velocity in R^n and dr and dv are the variation position and velocity
:param N: (int) number of bodies
:param dim: (int) dimension of the phase space for one body
:param zdz: (numpy array) state at which we want to compute the vector field. Dimenstion 2*dim*N
:param t: (float) time at which we want to compute the vector field. Here the vector field is time independant
:param rotating: (numpy array) rotation matrix for the frame. dimension dim/2 times dim/2
:return: (numpy array) dz/dt the value of the vector field at state z (dimension N*dim)
"""
# compute the dynamics
# position and velocity
z = zdz[0:N*dim].copy()
dzdt = system(z,t,m,N,dim,rotating)
# some constant parts
n = dim/2
nZero = numpy.zeros((n,n),dtype=float)
nId = numpy.eye(n,dtype=float)
zdzdt = numpy.zeros(2*N*dim)
zdzdt[0:N*dim] = dzdt.copy()
# the dynamical matrix
A = numpy.zeros((dim*N, dim*N))
for i in range(N):
qi = z[(i*dim):(i*dim+dim/2)]
vi = z[(i*dim+dim/2):(i*dim+dim)]
for j in range(N):
qj = z[(j*dim):(j*dim+dim/2)]
vj = z[(j*dim+dim/2):(j*dim+dim)]
elA = numpy.zeros((dim,dim))
# dfi_1/dq_j = 0
if(i==j):
elA[0:n,n:2*n] = nId.copy() # dfi_1/dvi
for k in range(N):
if(k!=i):
qk = z[(k*dim):(k*dim+dim/2)]
rik = numpy.linalg.norm(qk-qi)
# dfi_2/dqi
elA[n:2*n,0:n] = elA[n:2*n,0:n]-m[k]/numpy.power(rik,3.0)*nId+3.0*m[k]/numpy.power(rik,5.0)*(numpy.dot(numpy.array([qk-qi]).T,numpy.array([qk-qi])))
if(isinstance(rotating,numpy.ndarray)):
elA[n:2*n,n:2*n] = -2*rotating # dfi_2/dvi
elA[n:2*n,0:n] = elA[n:2*n,0:n]-numpy.dot(rotating,rotating) #dfi_2/dqi
else: #j!=i
rij = numpy.linalg.norm(qj-qi)
#dfi_2/dqj
elA[n:2*n,0:n] = m[j]/numpy.power(rij,3.0)*nId-3.0*m[j]/numpy.power(rij,5.0)*(numpy.dot(numpy.array([qj-qi]).T,numpy.array([qj-qi])))
A[i*dim:(i+1)*dim,j*dim:(j+1)*dim]=elA.copy()
## Validation with finite differences
# dx = 1e-09
# Adf = numpy.zeros((dim*N,dim*N))
# dzdt0 = system(z,t,m,N,dim,rotating)
# for i in range(dim*N):
# dzi = z.copy()
# dzi[i] = dzi[i]+dx
# dzim = z.copy()
# dzim[i] = dzim[i]-dx
# dzdti = system(dzi,t,m,N,dim,rotating)
# dzdtim = system(dzim,t,m,N,dim,rotating)
# Adf[i,:] = (dzdti-dzdtim)/(2.0*dx)
# print ""
# print "A=",A
# print "Adf=", Adf.T
# print "norm Adf-A = ", numpy.linalg.norm(Adf-A)
# print "norm Adf.T-A = ", numpy.linalg.norm(Adf.T-A)
zdzdt[N*dim:2*N*dim]=numpy.dot(A,zdz[N*dim:2*N*dim])
return zdzdt