Homework № 7

Kevin Mack (Worked with Alex DeWitt), Clarkson University 12/01/2017

Problem Statement:

In this assignment the task is create a model of the n-body problem, which is a famous problem

in physics. In this case, the system is designated to be celestial bodies who‘s force of action is

the gravitational force exerted by each body on each other. This is a problem that is notoriously

diﬃcult to solve analytically (it was actually shown by Poincarè that there is no analytical solution

given by algebraic expressions and integrals), and it can also exhibit chaotic behavior. Due to

these facts it is often solved with a numerical solution, also making it an appropriate problem to

apply equation-free methods in order to extract the dynamics from the data. A database of the

movement of celestial bodies would be an ideal place to gather data for equation-free modeling,

and is available from NASA and other astronomical societies. In this work no speciﬁc system is

modeled, however a general case of the n-body problem will be simulated using the ODE solver

in Matlab, which will then be compared with equation-free models.

Solution to the N-Body Problem:

The n-body problem can be generalized with an equation of motion that can be applied to each

body of the system. When all bodies are considered, the system will provide the motion of each

body with the consideration of the gravitational forces from every other body. The equations of

motion are given by,

j=1

j6=i



− q



− q

∂U

∂q

(1)

where q

is the position vector of the i

body, m

is the mass of the i

body, and G is the

universal gravitational constant. Using (1) it can be seen that in three spatial dimensions, the

number of diﬀerential equations used to solve the system is on the order of 6n (where n is the

number of bodies), which can become increasing computationally expensive for sets of large

particles. The ODE‘s, written in code as n_body.m (listing 1), is solved using Matlab‘s built-in

ODE solver ode113.m.

Listing 1: Matlab code – n_body.m

1 f u n c t i o n dx = n_body(~ , x )

2 %GM = 1 ; % Normalized G* mass

3 GM = 1 e2 ; % Mass*G that " works w e l l " f o r our i n i t i a l d i s t a n c e s and ←-

v e l o c i t i e s

Assignment № 7 Page 1

4 %x = re sha p e ( x , [ 2 , 6 ] ) ;

5 N = s i z e ( x , 1 ) / 6 ;

6 x = res h ape ( x , [ 6 , N] ) ‘ ; % re sha p e from v e c t o r to matrix

7 %x = [ x ( 1 : 6 ) ‘ ; x ( 7 : end ) ‘ ] ;

8 dx ( : , 1 ) = x ( : , 4 ) ; % s e t v e l o c i t y i n x

9 dx ( : , 2 ) = x ( : , 5 ) ; % s e t v e l o c i t y i n y

10 dx ( : , 3 ) = x ( : , 6 ) ; % s e t v e l o c i t y i n z

11 f o r i = 1: s i z e ( x , 1)

12 dx ( i , 4 : 6 ) = 0 ;

13 f o r j = 1 : s i z e ( x , 1)

14 i f ( j~= i )

15 df = x ( j , 1 : 3 ) - x ( i , 1 : 3 ) ; % ge t d i f f e r e n c e between each p a i r ←-

o f p a r t i c l e s

16 dx ( i , 4 : 6 ) = dx ( i , 4 : 6 ) + GM*( d f ) /( d f * df . ‘ ) . ^ ( 3 / 2 ) ; % compute←-

g r a v i t a t i o n a l e f f e c t s

17 end

18 end

19 end

20 dx = re s hap e ( dx ‘ , [N*6 , 1 ] ) ; % ODE must output a column v e c t o r

21 end

In this code, ode113 is used instead of ode45 due to its improvement in accuracy. When

using ode45 to solve these equations, the round-oﬀ error in the implementation of the Runge-

Kutta method was causing the system to lose energy artiﬁcially over time. The ode113 method

is variable-step, variable order (VSVO) Adams-Bashforth-Moulton solver of orders 1 to 13, and

has been noted for its ability to be used in celestial mechanics problems (which is similar, if

not equivalent to this problem). Using this solver and the code in n_body.m, we generate the

micro-scale model that is necessary to perform equation-free modeling.

Listing 2: Matlab code – Micro-scale simulation

1 %% n - body problem s i m u l a t i o n

3 N = 5 0 ; % number o f bo d ie s

4 pdim = 3 ; % number o f s p a t i a l d i mensio n s

5 vdim = 3 ; % number o f v e l o c i t y dimens i o ns

6 vvar = 1 e0 ; % va r i a n c e f o r i n i t i a l v e l o c i t y of p a r t i c l e s

7 max = 1 e2 ; % max d i s t a n c e f o r i n i t a l p o s i t i o n s o f p a r t i c l e s

9 p o s i t i o n s = -max+(max+max) . * rand (N, pdim ) ;

10 v e l o c i t i e s = s q r t ( vvar ) * randn (N, vdim ) ;

11 i c = [ p o s i t i o n s , v e l o c i t i e s ] ;

12 i c = resh a pe ( i c . ‘ , [ 1 , N* 6 ] ) ;

14 % s e t time parameters

15 t_ st ar t = 0;

Assignment № 7 Page 2

16 t_end = 1 e4 ;

17 % t_step = 5 e1 ;

18 t = [ t_ s t ar t t_end ] ;

20 t i c

21 [ T,Y] = ode113 (@n_body , t , [ i c ] ) ; % Sol ve ODE‘ s

22 toc

24 % g e t output i n c o r r e c t form t o p l o t time - s e r i e s

25 Y_plot = ze r o s ( pdim+vdim , N, s i z e (T, 1 ) ) ;

26 f o r i = 1: s i z e (T, 1 )

27 Y_plot ( : , : , i ) = r e sha p e (Y( i , : ) . ‘ , [ pdim+vdim , N] ) ;

28 end

30 %% Plot output

31 f i g u r e

32 hold a l l

33 g r i d on

34 g r i d minor

35 box on

37 f o r i = 1: s i z e (T, 1 )

38 i f ( i ~= 1 && i~= s i z e (T, 1 ) )

39 p l o t 3 ( Y_plot ( 1 , : , i ) , Y_plot ( 2 , : , i ) , Y_plot ( 3 , : , i ) , ‘ b . ‘ , ‘←-

MarkerSize ‘ , 3) % p l o t p a r t i c l e path with blu e marker

40 e l s e i f ( i == 1 )

41 p l o t 3 ( Y_plot ( 1 , : , i ) , Y_plot ( 2 , : , i ) , Y_plot ( 3 , : , i ) , ‘ g . ‘ , ‘ ←-

MarkerSize ‘ , 10) % p l o t s t a r t i n g po i nt with green marker

42 e l s e

43 p l o t 3 ( Y_plot ( 1 , : , i ) , Y_plot ( 2 , : , i ) , Y_plot ( 3 , : , i ) , ‘ r . ‘ , ‘ ←-

MarkerSize ‘ , 10)% p l o t end p oi nt with red marker

44 end

45 end

46 x l a b e l ( ‘ x ‘ )

47 y l a b e l ( ‘ y ‘ )

48 z l a b e l ( ‘ z ‘ )

49 t i t l e ( [ ‘ n - body s i m u l a t i o n f o r n = ‘ num2str (N) ] )

Since the equations cannot be solved analytically, as stated above, a numerical solution in

Matlab has been obtained. The solver will compute the position and velocity vectors for each

particle in the system. Several diﬀerent values for the number of particles will be considered,

and several graphs will show the output from each solution. In the solver, the initial conditions

for position of each particle in the system are generated from a uniform random distribution,

while the initial velocities are generated from a Gaussian distribution. Further, each mass is of

equal magnitude, and all have been normalized so that m

G = 1 (where m

is the mass of the i

Assignment № 7 Page 3

body, G is the universal gravitational constant). This is to limit the complexity of the simulation

and keep the computation times as small as possible.

(a) n-body simulation for n = 3 (b) n-body simulation for n = 3

Figure 1: Solutions to the n-body problem for n = 3, and n = 5 bodies. Each mass in the system

is equal, with only the initial conditions being diﬀerent. Two trials of each case are shown, and

it is clear that the initial conditions play a large role in the manifested dynamics of the system.

From Figure 1 it appears that the nbody simulator is working correctly for the n = 3 and

n = 5 cases. It is now important to test situations where the number of bodies is signiﬁcantly

higher. Tests for n = 50 and n = 100 are given in Figures 2 and 3. These tests caused a large

rise in computational complexity, this is due to the fact that the number of equations that are

solved by the ODE solver in Matlab is on the order of 6n (n being the number of bodies). As the

Assignment № 7 Page 4

(a) Angle 1, n-body simulation for n = 50 (b) Angle 2, n-body simulation for n = 50

(e) Angle 3, n-body simulation for n = 50 (long time) (f) Angle 4, n-body simulation for n = 50 (long time)

Figure 2: Solutions to the n-body problem for n = 50 bodies. Each mass in the system is equal,

with only the initial conditions being diﬀerent. Again, the starting positions of each mass is

labeled with a green dot, with the ﬁnal position shown as a red dot. Sub-ﬁgures (a) and (b) are

from the same simulation, with diﬀerent angles. Sub-ﬁgures (c),(d),(e), and (f) go together and

show a longer time evolution.

Assignment № 7 Page 5

(a) Angle 1, n-body simulation for n = 100

(b) Angle 2 (c) Angle 3

Figure 3: The images are similar to the previous n-body simulations, with only the number of

particles being diﬀerent. The system, as in the other ﬁgures, is heavily dependent on the initial

conditions of particle position and velocity (on the micro-scale).

Assignment № 7 Page 6

Equation Free Modeling:

Due to the fact that this model becomes computationally expensive for sets of large particles

or over long periods of time, it is appropriate to consider an equation free modeling approach.

The equation-free modeling algorithm consists of two main aspects, restricting and lifting. In

restriction, a selection of macroscopic variables is made in order to generate a model of the

system. The restriction operator converts the microscopic (ﬁne) model into the macroscopic

(course) model. The restriction operator generally completes the task of projecting the prin-

cipal components of the microscopic variables, onto the chosen macroscopic variables (using

methods such as SVD/PCA/POD). This is described in (2),

U = Mu, (2)

where U is the system of variables for the macroscopic behavior, M is the restriction op-

erator, and u is the system of variables for the microscopic behavior. It is important to note

that U ∈ R

and u ∈ R

, where generally M << m. Assuming that the microscopic system

can be adequately modeled by the macroscopic variables U, then the the remaining statistical

characteristics at the microscopic scale should be properly approximated by functionals of the

those selected for the restriction operation. If we say that the microscopic dynamics can be

transformed to a new set of variables, this gives a system of diﬀerential equations,

= g(v

, v

), (3)

where v = (v

, v

) is partitioned into M and m − M components. Then the variable v

corresponds to with the macroscopic variable U so that,

= g(v

, r(v

)) = F (U). (4)

Thus, the functions g, r, and F can be constructed systematically (not in closed form) to give

the process its equation-free approach.

In this work the cumulative distribution is used as the restriction operator. This portion of

the code computes the cumulative distributions of radial distance and radial velocity. The Lifting

operation is performed to bring the simulation from a lower-dimensional macroscale simulation

to a higher dimensional microscale model. The lifting operator generates initial conditions from

the macroscale variables for the individual particles in the simulation. The overall process can be

summarized in order by performing a Lift, evolving the microscopic variables in ﬁne-scale time

(Evolve), and then Restricting to the to the macroscopic variables for coarse-time evolution. The

code for the three steps mentioned previously (Lift, Evolve, Restrict) can be found in listing 3

below. In this case, the lifting operator will be performing Dynamic Mode Decomposition (DMD)

analysis (listing 4, dmd.m) on the chosen variables. This will also allow us to evolve the system

in time.

Listing 3: Matlab code – Lift -> Evolve -> Restrict

Assignment № 7 Page 7

1 %Alex DeWitt & Kevin Mack

2 %

3 %% n - body problem s i m u l a t i o n

4 c l e a r

6 N = 5 0 ; % number o f bo d ie s

7 pdim = 3 ; % number o f s p a t i a l d i m e nsions

8 vdim = 3 ; % number o f v e l o c i t y di mensio n s

9 vvar = 5 e0 ; % v a r i a n c e f o r i n i t i a l v e l o c i t y o f p a r t i c l e s

10 max = 1 e3 ; % max d i s t a n c e f o r i n i t a l p o s i t i o n s o f p a r t i c l e s

12 p o s i t i o n s = -max+(max+max) . * rand (N, pdim ) ;

13 v e l o c i t i e s = s q r t ( vvar ) * randn (N, vdim ) ;

14 i c = [ p o s i t i o n s , v e l o c i t i e s ] . ‘ ;

15 output_600_t = z e r o s (10+1 ,6 ,N) ;

16 output_600_t ( 1 , : , : ) = i c ;

17 num_iters = 1 0 0 ;

18 s t a r t = 1 ;

19 q u a r t e r = f l o o r ( num_iters /4) ;

20 h a l f = f l o o r ( num_iters /2) ;

21 t hre e _qu art e r = f l o o r (3 * num_iters /4) ;

23 f o r i t e r a t o r = 1 : num_iters

24 i c = output_600_t ( i t e r a t o r , : , : ) ;

25 i c = r esh a pe ( i c , [ 1 , N* 6 ] ) ;

26 % s e t time parameters

27 t_ st ar t = 0;

28 t_end = 5 e2 ;

29 t = [ t_ s t ar t : 1 : t_end ] ;

31 t i c

32 [ T,Y] = ode113 (@n_body , t , [ i c ] ) ; % Sol ve ODE‘ s

33 toc

34 Tc = s i z e (T, 1 ) ;

35 % g e t output i n c o r r e c t form t o p l o t time - s e r i e s

36 Y_plot = ze r o s ( pdim+vdim , N, Tc) ;

37 f o r i = 1: Tc

38 Y_plot ( : , : , i ) = r e sha p e (Y( i , : ) . ‘ , [ pdim+vdim , N] ) ;

39 end

41 Y_r = z e r o s (Tc ,N) ;

42 hboxs_r = l i n s p a c e ( 0 ,4 0 0 0 , 1 0 ) ;

43 f o r i = 1: Tc

44 Y_r( i , : ) = Y_plot ( 1 , : , i ) .^2+Y_plot ( 2 , : , i ) .^2+Y_plot ( 3 , : , i ) . ^ 2 ;

45 Y_r( i , : ) = s q r t (Y_r( i , : ) ) ;

46 r_h ( i , : ) = h i s t (Y_r( i , : ) , hboxs_r ) ;

47 end

Assignment № 7 Page 8

50 Y_v = ze r o s ( Tc ,N) ;

51 hboxs_v = 0 : 1 : 1 0 ;

52 f o r i = 1: Tc

53 Y_v( i , : ) = Y_plot ( 4 , : , i ) .^2+Y_plot ( 5 , : , i ) .^2+Y_plot ( 6 , : , i ) . ^ 2 ;

54 Y_v( i , : ) = s q r t (Y_v( i , : ) ) ;

55 v_h( i , : ) = h i s t (Y_v( i , : ) , hboxs_v ) ;

56 end

58 % DMD

59 X = [ r_h ,v_h ] . ‘ ;

60 dt = ( t_end - t _ s t a r t ) /(Tc - 1 ) ;

61 t = [T;T( 2 : 1 0 0 )+T( end ) ] ;

62 [ ~ ,X_dmd, modes ] = dmd( t ,X, dt ) ;

63 X_dmd = r e a l (X_dmd) ;

65 %%%%%%%% L i f t i n g %%%%%%%%%%

66 xx = l i n s p a c e (0 ,5 0 0 0 , 1 0 0 0 ) ;

67 dT = 5 0 00/10 0 0;

68 [ pos_pdf ] = c sa p s ( hboxs_r ,X_dmd( 1 : 1 0 , t ( end ) ) , . 9 9 , xx ) ; % smooth with s p l i n e

69 pos_pdf ( pos_pdf < 0) = 0 ;

70 pos_cdf = cumtrapz ( xx , pos_pdf ) ;

71 pos_cdf = pos_cdf . / pos_cdf ( end ) ;

73 pos_r = z e r o s (N, 1 ) ;

74 uni = rand (N, 1 ) ;

76 f o r i = 1:N

77 pos_r ( i ) = fi n d ( pos_cdf>uni ( i ) , 1 ) . * dT;

78 end

80 th et a = rand (N, 1 ) *2* pi ;

81 phi = rand (N, 1 ) * p i ;

83 o j = z e r o s ( 6 ,N) ;

84 o j ( 1 , : ) = pos_r . * s i n ( phi ) . * cos ( th e t a ) ;

85 o j ( 2 , : ) = pos_r . * s i n ( phi ) . * s i n ( t he ta ) ;

86 o j ( 3 , : ) = pos_r . * c os ( phi ) ;

88 cdf_dummy = @cumtrapz ;

89 %r a d i a l v e l o c i t y

90 xx = l i n s p a c e (0 , 1 0 , 1 0 00 ) ;

91 dT = 1 0/1 0 00;

92 [ v_pdf ] = c s ap s ( hboxs_v ,X_dmd( 1 1 : 2 1 , t ( end ) ) , . 9 9 , xx ) ;%smooth with s p l i n e

93 v_pdf ( v_pdf < 0) = 0 ;

94 v_cdf = cdf_dummy ( xx , v_pdf ) ;

Assignment № 7 Page 9

95 v_cdf = v_cdf . / v_cdf ( end ) ;

97 v_r = z e r o s (N, 1 ) ;

98 uni = rand (N, 1 ) ;

99 f o r i = 1:N

100 v_r ( i ) = fi n d ( v_cdf>uni ( i ) , 1 ) . * dT;

101 end

102

103 % p l o t h is t ogr a ms and graphs

i t e r a t o r == t h ree _ qu a rte r | | i t e r a t o r == num_iters )

105 bins_pos = f l o o r ( s q r t ( s i z e ( pos_r , 1) ) ) ;

106 bins_v = f l o o r ( s q r t ( s i z e (v_r , 1) ) ) ;

107 f i g u r e

108 hold on

109 histogram ( pos_r , bins_pos )

110 t i t l e ( { ‘ Histogram o f Radial D ista n ces , i t e r a t i o n ‘ i t e r a t o r })

111

112 f i g u r e

113 hold on

114 histogram ( pos_r , bins_v )

115 t i t l e ( { ‘ Histogram o f Radial V e l o c i t i e s , i t e r a t i o n ‘ i t e r a t o r })

116

117 f i g u r e

118 hold on

119 gr i d on

120 gr i d minor

121 box on

122 pl o t ( v_pdf )

123 t i t l e ( { ‘CDF o f Ra d i a l V el o c i ty , i t e r a t i o n ‘ i t e r a t o r })

124 end

125

126 %r a d i a l v e l o c i t y i s i n the d i r e c t i o n as r

127 o j ( 4 , : ) = v_r . * s i n ( phi ) . * co s ( t h e t a ) ;

128 o j ( 5 , : ) = v_r . * s i n ( phi ) . * s i n ( t h e t a ) ;

129 o j ( 6 , : ) = v_r . * c os ( phi ) ;

130

131 %t a n g e n t i a l v e l o c i t y i s p r o p o r t i o n a l to r

132 th et a = rand (N, 1 ) *2* p i ;

133 phi = p hi + pi / 2 ;

134

135 %compute c i r c u l a r o r b i t s

136 v_t = s q r t (100*N. / pos_r ) ;%

137 o j ( 4 , : ) = o j ( 4 , : ) ‘ + v_t . * s i n ( phi ) . * c os ( t h e t a ) ;

138 o j ( 5 , : ) = o j ( 5 , : ) ‘ + v_t . * s i n ( phi ) . * s i n ( t h e t a ) ;

139 o j ( 6 , : ) = o j ( 6 , : ) ‘ + v_t . * c o s ( phi ) ;

140

Assignment № 7 Page 10

141 %%%%%%%%o j = new Y v e c t o r%%%%%%%%%%%

142

143 output_600_t ( i t e r a t o r + 1 , : , : ) = o j ;

144 end

145

146 markers = { ‘ k . ‘ , ‘ b . ‘ , ‘ g . ‘ , ‘ r . ‘ , ‘m. ‘ } ;

147 p l o t s = 0 ;

148 f i g u r e

149 hold on

150 % Plot b o di e s at d i f f e r e n t p o i n t s i n the s i m u l a t i o n

151 f o r i =1: i t e r a t o r

== num_iters )

153 p l o t 3 ( s qu ee z e ( output_600_t ( i , 1 , : ) ) , s qu ee z e ( output_600_t ( i , 2 , : ) ) , ←-

sq u ee ze ( output_600_t ( i , 3 , : ) ) , markers { p l o t s })

154 a x i s ( [ - 5 0 0 0 , 5 0 0 0 , - 5 0 0 0 , 5 0 0 0 , - 5 0 0 0 , 5 0 0 0 ] ) ;

155 end

156 end

Listing 4: Matlab code – dmd.m

1 f u n c t i o n [ t ,u_dmd, u_modes ] = dmd( t ,X, dt )

2 u = X( : , 1 ) ; % c r e a t e i n t i a l c o n d i t i o n s

3 X1 = X( : , 1 : end - 1 ) ;

4 X2 = X( : , 2 : end ) ;

6 %perform dmd

7 [U, Sigma ,V] = svd (X1 , ‘ econ ‘ ) ;

8 Sigma = Sigma + d iag (10^ -10* ones ( s i z e ( Sigma , 1 ) , 1 ) ) ;

9 Sigma_inv = di a g ( 1 . / d iag ( Sigma ) ) ;

10 S = U‘ * X2*V*Sigma_inv ;

11 [ ev ,D] = e i g ( S ) ;

12 mu = di a g (D) ;

13 omega = l o g (mu) /( dt ) ;

14 Phi = U*ev ;

16 y0 = pinv ( Phi ) *u ; %pseudo - i n v e r s e on i n i t i a l c o n d i t i o n s

17 u_modes = z e r o s ( s i z e (V, 2 ) , l e n g t h ( t ) ) ;

18 f o r i t e r = 1 : l e n g t h ( t )

19 u_modes ( : , i t e r ) = ( y0 . * exp ( omega* t ( i t e r ) ) ) ;

20 end

21 u_dmd = Phi*u_modes ; % r e t u r n the dmd modes

In conclusion, the code above provides the macro-evolution of the system, with the aid from

simulating small steps in the microscale evolution at regular time intervals. This allows for the

Assignment № 7 Page 11

modeling of the behavior of the chosen macro-scale variables (radial distance and radial veloc-

ity) over comparatively large time scales to the microscale simulation. Figures 4 and 5 represent

some of the histograms created for the restriction operator. As for the initial conditions of the

simulation, the starting positions of the particles were computed from a uniform distribution,

while the initial velocities were computed from a random normal distribution. For the lifting op-

erator, the new initial conditions for microscale simulation were again generated from random

distributions, but this time in compliance with the histograms of the macro-scale variables at that

speciﬁc time. The overall result is the evolution of the macro-scale variables across coarse time

(Figures 4 and 5). The graphs were generated through 100 iterations (meaning 100 instances

of micro-scale iterations).

Assignment № 7 Page 12

(a) CDF of radial position for iteration 1 (b) CDF of radial velocity for iteration 1

(e) CDF of radial position for iteration 100 (f) CDF of radial velocity for iteration 100

Figure 4: The ﬁgures in the left column represent the CDF of the radial position, with the right

column being the distribution for the radial velocity for the particles in the simulation. The CDF

is used in the restriction phase of the simulation.

Assignment № 7 Page 13

(a) Histogram of radial position for iteration 1 (b) Histogram of radial velocity for iteration 1

(e) Histogram of radial position for iteration 100 (f) Histogram of radial velocity for iteration 100

Figure 5: The histograms in the left column represent the radial distances of the particles at

diﬀerent points in the simulation, withe the histograms on the right representing the radial veloc-

ities. These histograms represent the evolution of the macro-scale variables over coarse-scale

time.

Assignment № 7 Page 14