Nbody Class

This module provides the basic functions to build objects for the simulation of a N-body problem.

Build a Nbody instance

To build a Nbody instance, we simply call the constructor.

Simple constructor

import nbodypy
system = nbodypy.Nbody()

Here the documentation of the constructor.

Nbody.__init__(N=2, mass=None, init=None, dim=4, rotating=None, jsonFile=None, zdzInit=None)

Constructor

Parameters:

• N – int The number of bodies. default 2.
• mass – numpy.array of masses of the bodies. default None that implies each body has a mass equal to one.
• init – numpy.array of initial conditions (ex: [x1, y1, dx1, dy1,x2,y2,dx2,dy2...] in dimension 4 (2D)). default None, the bodies are on a circle.
• dim – number of dimension (space and velocity). default 4 (planar motion)
• rotating – rotation matrix (dim/2 times dim/2) for the frame. Default None (no rotation)
• jsonFile – name of the json file to construct the object. default none.

Let us show a more complete example.

Complete example

import nbodypy
import numpy

zinit = numpy.array([0.165649854848924,  3.0549363634996e-151 ,    0.209141794901608,
                  8.86909514289322e-29 ,     2.03225904939571, -9.31938237331704e-29,
                  -0.0828249274244622 ,    0.160070428973807,    -0.104570897450804,
                  -2.04906075007578 ,    -1.01612952469786 ,     1.15731474026243,
                  -0.0828249274244618 ,   -0.160070428973807 ,   -0.104570897450804,
                  2.04906075007578,     -1.01612952469785,     -1.15731474026243])

Arotating = numpy.array([[0.0, 5.13873206001209, 0.0],
                      [-5.13873206001209, 0.0, 0.0],
                      [0.0,0.0, 0.0]])

System = nbodypy.Nbody(N=3, dim=6, init = zinit, rotating=Arotating)

With a Json File

For the same example, one can build an instance of Nbody class with a Json file.

example.json

{
   "N" : "3",
   "dim" : "6",
   "init" : "[0.165649854848924,  3.0549363634996e-151 ,    0.209141794901608,8.86909514289322e-29 ,     2.03225904939571, -9.31938237331704e-29,-0.0828249274244622 ,    0.160070428973807,    -0.104570897450804,-2.04906075007578 ,    -1.01612952469786 ,     1.15731474026243,-0.0828249274244618 ,   -0.160070428973807 ,   -0.104570897450804,2.04906075007578,     -1.01612952469785,     -1.15731474026243]",
   "rotating" : "[[0.0, 5.13873206001209, 0.0],[-5.13873206001209, 0.0, 0.0],[0.0,0.0, 0.0]]"
}

Complete example

import nbodypy
import numpy

System = nbodypy.Nbody(jsonFile="example.json")

Access to the Parameters

Get the Number of Bodies

Nbody.get_N()

Get the number of bodies :return: (int) Return the size of the system

import nbodypy
system = nbodypy.Nbody()
Nbodies = system.get_N()

Get the Dimension

Nbody.get_dim()

Get the dimensions (position+velocity of one body) :return: (int) Return the size of the system dynamics for each body

The dimension parametrer corresponds to the dimension of the phase space for one body.

import nbodypy
system = nbodypy.Nbody()
Nbodies = system.get_dim()

Get the Current State

Nbody.get_z()

Get the position of one body

Returns:Return the current state (numpy.array of dimension N*dim)

To get the current state of the system, a numpy.array of dimension the number of bodies times the dimension.

import nbodypy
system = nbodypy.Nbody()
Nbodies = system.get_z()

One can simply use system.z.

Get the Position and the Velocity

Nbody.get_r(I)

Get the position of one body

Parameters:I – the number of the considered body (starting at 0) Returns:Return the position of the Ith Body

Nbody.get_v(I)

Get the velocity of one body

Parameters:I – the number of the considered body (starting at 0) Returns:Return the velocity of the Ith Body

We provide two methods to get the position and the velocity of one chosen body. Of course, this depends on the dimension of the system.

import nbodypy
import numpy

zinit = numpy.array([0, 0.316046051487446, 0.864077674448677, 0,  0.100941149773621, -0.158023025743723, -0.432038837224337, 2.02906597904133, -0.100941149773621 , -0.158023025743723,  -0.43203883722434, -2.02906597904133])

system = nbodypy.Nbody(N=3,init=zinit)
R1 = system.get_r(0) # get the position of the first body
V3 = system.get_v(2) # get the velocity of the third body

Integrate a solution

The dynamics

This part of the Nbody class depends on the function:

nbodypy.dynamics.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

Parameters:

• N – (int) number of bodies
• dim – (int) dimension of the phase space for one body
• z – (numpy array) state at which we want to compute the vector field. Dimension : dim*N
• t – (float) time at which we want to compute the vector field. Here the vector field is time independant
• rotating – (numpy array) rotation matrix for the frame. dimension dim/2 times dim/2

Returns:

(numpy array) dz/dt the value of the vector field at state z (dimension N*dim)

Normally, you do not have to use this function, but if you want to, you have to import it as follows:

import nbodypy.dynamics

The integrate method

The method to integrate a solution is:

Nbody.integrate(t, fileName=None)

Method to integrate a solution for a given N-body problem

Parameters:

• t – (numpy array) t of all the times for which we want to get a point
• fileName – (string) optional parameter to save the solution in a external file. The first column of the file contains the time steps. The others contain the value of the phase state of each body.

Returns:

the integrate solution at each times in a numpy array of size of t

For instance, the following exemple:

Example of integration

import nbodypy
import numpy

zinit = numpy.array([0, 0.316046051487446, 0.864077674448677, 0,  0.100941149773621, -0.158023025743723, -0.432038837224337, 2.02906597904133, -0.100941149773621 , -0.158023025743723,  -0.43203883722434, -2.02906597904133])

system = nbodypy.Nbody(N=3,init=zinit)
t = numpy.linspace(0, 1.0/3, 100)
solution = system.integrate(t)

Save the solution in a file

One can save the integration of a solution in a file. To do that, one just has to use the optional parameter fileName.

Example of integration

import nbodypy
import numpy

zinit = numpy.array([0, 0.316046051487446, 0.864077674448677, 0,  0.100941149773621, -0.158023025743723, -0.432038837224337, 2.02906597904133, -0.100941149773621 , -0.158023025743723,  -0.43203883722434, -2.02906597904133])

system = nbodypy.Nbody(N=3,init=zinit)
t = numpy.linspace(0, 1.0/3, 100)
solution = system.integrate(t,fileName="exportSol.dat")

This can be use without getting the solution in a variable:

Example of integration

system.integrate(t,fileName="exportSol.dat")

The file contains in its first column all the value of the time at each step, and the others contain the phase space states of each body. For example, in the planar case for \(N\) bodies:

\[t_i, x_{i,1}^1, x_{i,1}^2, v_{i,1}^1, v_{i,1}^2,\dots, x_{i,N}^1, x_{i,N}^2, v_{i,N}^1, v_{i,N}^2\]

Move the Initial Condition

You can also “move” the initial position system.z. With the move membre you can integrate the system during a time t, and modify the current state system.z.

Nbody.move(t)

Move the state integrating self.z during the time t :param t: time (float) for integration :return: move the self.z to the point after integration during t

Example of moving

import nbodypy
import numpy

zinit = numpy.array([0, 0.316046051487446, 0.864077674448677, 0,  0.100941149773621, -0.158023025743723, -0.432038837224337, 2.02906597904133, -0.100941149773621 , -0.158023025743723,  -0.43203883722434, -2.02906597904133])

system = nbodypy.Nbody(N=3,init=zinit)
system.move(0.3)

The Linearized Dynamics

This part of the Nbody class depends on the function:

nbodypy.dynamics.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

Parameters:

• N – (int) number of bodies
• dim – (int) dimension of the phase space for one body
• zdz – (numpy array) state at which we want to compute the vector field. Dimenstion 2*dim*N
• t – (float) time at which we want to compute the vector field. Here the vector field is time independant
• rotating – (numpy array) rotation matrix for the frame. dimension dim/2 times dim/2

Returns:

(numpy array) dz/dt the value of the vector field at state z (dimension N*dim)

Normally, you do not have to use this function, but if you want to, you have to import it as follows:

import nbodypy.dynamics

The Dintegrate Method

The method to integrate a solution and the linearized system around this solution is:

Nbody.Dintegrate(t, fileName=None)

Method to integrate a solution for a given N-body problem and the linearized system around the solution

Parameters:

• t – (numpy array) t of all the times for which we want to get a point
• fileName – (string) optional parameter to save the solution in a external file. The first column of the file contains the time steps. The others contain the value of the phase state of each body.

Returns:

the integrate solution at each times in a numpy array of size of t

For instance, the following exemple:

Example of Dintegration

import nbodypy
import numpy

zinit = numpy.array([0, 0.316046051487446, 0.864077674448677, 0,  0.100941149773621, -0.158023025743723, -0.432038837224337, 2.02906597904133, -0.100941149773621 , -0.158023025743723,  -0.43203883722434, -2.02906597904133])

system = nbodypy.Nbody(N=3,init=zinit)
t = numpy.linspace(0, 1.0/3, 100)

e1 = numpy.zeros(3*4, dtype=float)
e1[0] = 1.0
system.zdz = numpy.hstack((system.z, e1))
solution = system.Dintegrate(t)

Compute the Monodromy Matrix

The monodromy matrix gives us information about the influence of perturbations around the periodic orbit. Let us consider a periodic solution \(x^{*}(\cdot)\) of period \(T\), the flow from the initial point \(x_{0}^{*}=x^{*}(0)\) is denoted by \(\varphi(t,x_{0}^{*})\).

The matrix

\[M=\dfrac{\partial\varphi(T,x_{0}^{*})}{\partial x_{0}}\]

determines whether initial perturbations \(x_{0}^{*}\) from the periodic orbit \(x^{*}(\cdot)\) decay or grow. This matrix is called the monodromy matrix.

Integrating the linearized system

One way to compute it is to solve a certain dynamical system.

Let us consider the dynamical system \(\dot x=f(x)\) and a periodic solution \(x^{*}\). The monodromy matrix \(M\) is the matrix \(\Phi(T)\) where \(\Phi\) is a solution of the linearized system

\[\begin{split}\left\{ \begin{array}{l} \dfrac{\partial\Phi}{\partial t}=\dfrac{\partial f(x^*)}{\partial x}\Phi,\\ \Phi(0)=Id. \end{array} \right.\end{split}\]

The function to compute the monodromy matrix following this way is:

Nbody.monodromy(t)

Compute the monodromy matrix starting à self.z position pour a periodic solution of period t. The computation is done integrating the linearized equation around the periodic solution.

Parameters:t – (float) time for integration. The period of the periodic solution Returns:(numpy.array) the monodromy matrix. Dim (N*dim, N*dim)

This can be used as follow:

Example of computing the monodromy matrix

import nbodypy
import numpy

zinit = numpy.array([0, 0.316046051487446, 0.864077674448677, 0,  0.100941149773621, -0.158023025743723, -0.432038837224337, 2.02906597904133, -0.100941149773621 , -0.158023025743723,  -0.43203883722434, -2.02906597904133])

system = nbodypy.Nbody(N=3,init=zinit)
M = system.monodromy(1.0)

By Finite Differences

There is another function to compute the monodromy matrix using finite differences. It is faster but less precise.

Nbody.monodromyDF(t, dx=1e-09)

Compute the monodromy matrix starting à self.z position pour a periodic solution of period t. The computation is done by finite differences.

Parameters:t – (float) time for integration. The period of the periodic solution Returns:(numpy.array) the monodromy matrix. Dim (N*dim, N*dim)

Example of computing the monodromy matrix

import nbodypy
import numpy

zinit = numpy.array([0, 0.316046051487446, 0.864077674448677, 0,  0.100941149773621, -0.158023025743723, -0.432038837224337, 2.02906597904133, -0.100941149773621 , -0.158023025743723,  -0.43203883722434, -2.02906597904133])

system = nbodypy.Nbody(N=3,init=zinit)
M = system.monodromyDF(1.0)