N-body simulations

Various ways to simulate the motion of N bodies subject to a pairwise potential.
While most of the time we take N=3 and d=2 with a Newtonian potential, these
functions should work for arbitrary number of particles and dimensions, as well as other potentials.

SetDirectory@NotebookDirectory[];

Newtonian gravity and initial conditions

ClearAll@NewtonianPotential;
NewtonianPotential[v_]:=-(1/(Sqrt[Total@(v^2)]));

ClearAll@F;
F[pot_][a_,a_]:=0;
F[pot_][a_,b_]:=With[{vars=Table[dummy[i],{i,Length@a}]},-Grad[pot[vars-b],vars]/.dummy[i_]:>a[[i]]]

Initial states and their properties

A state is represented by a point in phase space, i.e. a list of N coordinates and N momenta, {X’s,P’s}

Arbitrary number of particles and dimensions

ClearAll@KineticEnergy;
KineticEnergy[state_]:=Total[Last@((state^2)/2),∞];

ClearAll@RandomInitialConditions;
RandomInitialConditions[n_Integer,d_Integer]:=
With[{X0=RandomReal[{-1,1},{n,d}],P0=RandomReal[{-1,1},{n,d}]},
{#-Mean@X0&/@X0,#-Mean@P0&/@P0}
]

ApplyVariables[{state_,vars_},values_]:=state/.(First@#->Last@#&/@Transpose@{vars,values})

Ones[n_Integer]:=Table[1,{n}]

2D motion of N particles

ClearAll@AngularMomentum;
AngularMomentum[state_]:=Last@Sum[Append[First@part,0]\[Cross]Append[Last@part,0],{part,Transpose@state}]

ClearAll@SymInitialConditions2D;
SymInitialConditions2D[n_Integer][Px_,Py_]:={Table[{i,0},{i,-(n-1)/2,(n-1)/2}],Join[{{Px,Py}},Table[(((-1))^i){2Px,2Py},{i,n-2}],{{Px,Py}}]}

ClearAll@L0InitialConditions;
L0InitialConditions[n_Integer]:=
Block[{xs,ps,solvars,sol},
xs=Table[{x[i],y[i]},{i,n}];
ps=Table[{px[i],py[i]},{i,n}];
solvars={x[n],y[n],px[n],py[n],py[n-1]};
sol=First@Solve[{Mean@xs,Mean@ps,AngularMomentum@{xs,ps}}==0//LogicalExpand,solvars];
{{xs,ps}/.sol,Complement[Flatten@{xs,ps},solvars]}
]

RandomL0InitialConditions[n_Integer]:=
With[{L0init=L0InitialConditions@n},
ApplyVariables[L0init,RandomReal[{-1,1},Length@Last@L0init]]
]

L0InitialConditions[3]

Some examples of initial conditions

RandomInitialConditions[3,2]

Figure8Init=SymInitialConditions2D[3][0.3471128135672417`,0.532726851767674`]

AngularMomentum@Figure8Init

KineticEnergy@Figure8Init

Plotting and animations

These functions work for any number of particles, but only in 2D.

PlotSimulation[sim_]:=
Block[{Ti,Tf},
{Ti,Tf}=First@First[sim]["Domain"];
ParametricPlot[Evaluate@Through[sim[t]],{t,Ti,Tf},AxesLabel->{"x","y"},ImageSize->Large,PlotRange->All]
]

AnimateSimulation[sim_]:=
Module[{Ti,Tf,curve=PlotSimulation@sim,plot},
{Ti,Tf}=First@First[sim]["Domain"];
Animate[
Show[
curve,
ListPlot[{#}&/@Through[sim[t]],PlotStyle->PointSize@.025]
,AspectRatio->1/GoldenRatio]
,{t,Ti,Tf},AnimationRunning->False]
]

PlotState[sim_,t0_Real]:=
Block[{},
Show[
Join[
{ParametricPlot[Through[sim[t]],{t,0,6.325},PlotStyle->Black],ListPlot[{#}&/@Through[sim[t0]],PlotStyle->PointSize@.025]},
Table[Graphics[{Black,Arrow[{traj[t0],traj[t0]+(1/4)Normalize[traj'[t0]]}]}],{traj,sim}]
]
,ImageSize->720,AxesLabel->{"x","y"}]
]

Fixed-step integrator

ClearAll@FSIntegrator;
FSIntegrator[n_Integer,ms_,pot_,init_][dT_Real,T_Real]:=
Block[{X,P,dP},
{X[0],P[0]}=init;
dP[k_Integer]:=dP[k]=Table[Sum[F[pot][X[k][[i]],X[k][[j]]],{j,n}],{i,n}];
X[k_Integer]:=X[k]=X[k-1]+((P[k-1])/ms)dT+dP[k-1]((dT^2)/(2ms));
P[k_Integer]:=P[k]=P[k-1]+dP[k-1]dT;
Interpolation/@Transpose/@({dT Range[0,T/dT],#}&/@Transpose[X/@Range[0,T/dT]])
]

The figure 8 solution

Figure8FS=FSIntegrator[3,Ones@3,NewtonianPotential,Figure8Init][.0001,6.325/3];

AnimateSimulation@Figure8FS

Using NDSolve

We use WhenEvent to stop integrating if the particles escape further than a maximum radius

ClearAll@MergeInterpolations;
MergeInterpolations[{x_InterpolatingFunction,y_InterpolatingFunction}]:=Interpolation[Transpose@{First/@x["Grid"],Flatten/@Transpose@{x["ValuesOnGrid"],y["ValuesOnGrid"]}}];
MergeInterpolations[{x_InterpolatingFunction,y_InterpolatingFunction,z__InterpolatingFunction}]:=MergeInterpolations@{MergeInterpolations@{x,y},z};

NDSIntegrator[n_Integer,d_Integer,ms_,pot_,init_][T_Real,R_Real]:=
Block[{eqs,bc,vars,maxnorm},
eqs=Table[ms[[i]]Table[x[i,k]''[t],{k,d}]-Sum[F[pot][Table[x[i,k][t],{k,d}],Table[x[j,k][t],{k,d}]],{j,n}],{i,n}];
bc={Table[x[i,k][0],{i,n},{k,d}]-First@init,Table[ms[[i]]x[i,k]'[0],{i,n},{k,d}]-Last@init};
vars=Table[x[i,k],{i,n},{k,d}];
maxnorm=Max[Norm[#,2]&/@Through/@Through@vars@t]>R;
MergeInterpolations/@NDSolveValue[{LogicalExpand[{eqs,bc}==0],WhenEvent[Evaluate@maxnorm,"StopIntegration"]},vars,{t,0,T}]
]

The figure 8 solution

Figure8NSD=NDSIntegrator2[3,2,Ones@3,NewtonianPotential,Figure8Init][6.325,2.]

AnimateSimulation@Figure8NSD

Some random initial conditions (chaos)

Three bodies with arbitrary boundary conditions (fixed center of mass)

RandomInit=RandomInitialConditions[3,2]
NDSIntegrator[3,2,Ones@3,NewtonianPotential,RandomInit][2.,10.];
AnimateSimulation@%

Small changes in the initial conditions result in exponentially large deviations in the final state

δ={{{0,0},{0,0},{0,0}},{{(10^(-2)),0},{0,0},{0,0}}};
AnimateSimulation@NDSIntegrator[3,2,Ones@3,NewtonianPotential,RandomInit+δ][2.,10.]

We can find funny “bound-states”

RandomInit=RandomInitialConditions[3,2]
NDSIntegrator[3,2,Ones@3,NewtonianPotential,%][10.,2.];
AnimateSimulation@%

We can also use more particles

RandomInitialConditions[5,2]
NDSIntegrator[5,2,Ones@5,NewtonianPotential,%][2.,10.];
AnimateSimulation@%

Export["animation.pdf",%212]

Random initial conditions with zero angular momentum

RandomInit=RandomL0InitialConditions@3
NDSIntegrator[3,2,Ones@3,NewtonianPotential,RandomInit][10.,2.]
AnimateSimulation@%

Some other periodic orbits found in the literature

AnimateSimulation@NDSIntegrator[3,2,Ones@3,NewtonianPotential,SymInitialConditions2D[3][0.306892758965492,0.125506782829762]][6.23564136316479,10.]

AnimateSimulation@NDSIntegrator[3,2,Ones@3,NewtonianPotential,SymInitialConditions2D[3][0.464445237398184,0.396059973403921]][14.8939113169584,10.]

AnimateSimulation@NDSIntegrator[3,2,Ones@3,NewtonianPotential,SymInitialConditions2D[3][0.513938054919243,0.304736003875733]][17.328369755004,10.]

Using NBodySimulation

MMAIntegrator[n_Integer,d_Integer,ms_,init_][T_Real]:=
Block[{state},
state=Table[<|"Mass"->ms[[i]],"Position"->init[[1,i]],"Velocity"->((init[[2,i]])/(ms[[i]]))|>,{i,n}];
NBodySimulation["InverseSquare",state,T][All,"Position"]
]

The figure 8 solution

Figure8MMA=MMAIntegrator[3,2,Ones@3,Figure8Init][6.325];

AnimateSimulation@Figure8MMA

A solution with N=10 bodies

RandomL0InitialConditions@10
AnimateSimulation@MMAIntegrator[10,2,Ones@10,%][1.]

Can we find a chain?

Block[{v=.75,θ=0.7,px,py},
{px,py}=v{Cos@θ,Sin@θ};
Print@{px,py};
AnimateSimulation@MMAIntegrator[5,2,Ones@5,SymInitialConditions2D[5][px,py]][2.]
]

Conserved quantities in the figure 8

Dynamic@{a,b}

Block[{x1,x2,x3,i1,i2,i3,ϵ=(10^(-3)),di1,di2,di3,f,minv,minp,mina,minb},
{x1,x2,x3}=Figure8MMA;
i1[t_Real]:=((Norm[x1[t]-x2[t],2])^2)((Norm[x1[t]-x3[t],2])^2)((Norm[x2[t]-x3[t],2])^2);
i2[t_Real]:=((Norm[x1[t]-x2[t],2])^2)+((Norm[x1[t]-x3[t],2])^2)+((Norm[x2[t]-x3[t],2])^2);
i3[t_Real]:=((Norm[x1'[t],2])^2)+((Norm[x2'[t],2])^2)+((Norm[x3'[t],2])^2);
di1[t_Real]:=((i1[t+ϵ]-i1[t])/ϵ);
di2[t_Real]:=((i2[t+ϵ]-i2[t])/ϵ);
di3[t_Real]:=((i3[t+ϵ]-i3[t])/ϵ);
f[a_Real,b_Real]:=NIntegrate[di1[t]+a di2[t]+b di3[t],{t,0,6.325}];
FindMinimum[f[a,b],{{a,1},{b,1}}]
]

FindMinimum

Block[{x1,x2,x3},
{x1,x2,x3}=Figure8MMA;
Print@Plot[((Norm[x1[t]-x2[t],2])^2)((Norm[x1[t]-x3[t],2])^2)((Norm[x2[t]-x3[t],2])^2),{t,0,6.325}];
Table[{t,((Norm[x1[t]-x2[t],2])^2)((Norm[x1[t]-x3[t],2])^2)((Norm[x2[t]-x3[t],2])^2)},{t,0,6.325,(6.325/((2^10)-1))}]
];

NonlinearModelFit[%,a+b Cos[c t]+d Sin[c t],{{a,2.45},{b,.3},{c,6},{d,0}},t]

Show[,Plot[%[t],{t,0,6.325},PlotStyle->{Red,Dashed}]]

Block[{x1,x2,x3},
{x1,x2,x3}=Figure8MMA;
Print@Plot[((Norm[x1[t]-x2[t],2])^2)+((Norm[x1[t]-x3[t],2])^2)+((Norm[x2[t]-x3[t],2])^2),{t,0,6.325}];
Table[{t,((Norm[x1[t]-x2[t],2])^2)+((Norm[x1[t]-x3[t],2])^2)+((Norm[x2[t]-x3[t],2])^2)},{t,0,6.325,(6.325/((2^10)-1))}]
];

NonlinearModelFit[%,a+b Cos[c t]+d Sin[c t],{{a,2.45},{b,.3},{c,6},{d,0}},t]

Show[,Plot[%[t],{t,0,6.325},PlotStyle->{Red,Dashed}]]

Block[{x1,x2,x3},
{x1,x2,x3}=Figure8MMA;
Print@Plot[0.013031263547355615`(((Norm[x1[t]-x2[t],2])^2)((Norm[x1[t]-x3[t],2])^2)((Norm[x2[t]-x3[t],2])^2))-0.1987743855346075`(((Norm[x1[t]-x2[t],2])^2)+((Norm[x1[t]-x3[t],2])^2)+((Norm[x2[t]-x3[t],2])^2)),{t,0,6.325},PlotRange->All];
];

Block[{x1,x2,x3},
{x1,x2,x3}=Figure8MMA;
Print@Plot[x1'[t].x1'[t]+x2'[t].x2'[t]+x3'[t].x3'[t],{t,0,6.325},PlotRange->All];
Table[{t,x1'[t].x1'[t]+x2'[t].x2'[t]+x3'[t].x3'[t]},{t,0,6.325,(6.325/((2^10)-1))}]
];

NonlinearModelFit[%,a+b Cos[c t]+d Sin[c t],{{a,2.45},{b,.3},{c,6},{d,0}},t]

Show[,Plot[%[t],{t,0,6.325},PlotStyle->{Red,Dashed}]]

Block[{x1,x2,x3},
{x1,x2,x3}=Figure8MMA;
Print@Plot[0.1542507160492513`(((Norm[x1[t]-x2[t],2])^2)((Norm[x1[t]-x3[t],2])^2)((Norm[x2[t]-x3[t],2])^2))+0.1987743855346075`(x1'[t].x1'[t]+x2'[t].x2'[t]+x3'[t].x3'[t]),{t,0,6.325}];
];

0.013031263547355615`(3.805301799949746` +0.1987743855346075` Cos[5.960234450409174` t]+0.0004681002925529526` Sin[5.960234450409174` t])-0.1987743855346075`(5.9870653344861955` +0.013031263547355615` Cos[5.960005227008194` t]+0.000021384038472949802` Sin[5.960005227008194` t])/.{5.960005227008194`->5.96,5.960234450409174`->5.96}//Expand

0.1542507160492513`(3.805301799949746` +0.1987743855346075` Cos[5.960234450409174` t]+0.0004681002925529526` Sin[5.960234450409174` t])+0.1987743855346075`(2.574327954867511` -0.1542507160492513` Cos[5.960803788824705` t]-0.0006407232443525828` Sin[5.960803788824705` t])/.{5.960803788824705`->5.96,5.960234450409174`->5.96}//Expand

-0.000055154563884767866` (0.013031263547355615`(3.805301799949746` +0.1987743855346075` Cos[5.960234450409174` t]+0.0004681002925529526` Sin[5.960234450409174` t])-0.1987743855346075`(5.9870653344861955` +0.013031263547355615` Cos[5.960005227008194` t]+0.000021384038472949802` Sin[5.960005227008194` t]))-1.849339171142787`*^-6(0.1542507160492513`(3.805301799949746` +0.1987743855346075` Cos[5.960234450409174` t]+0.0004681002925529526` Sin[5.960234450409174` t])+0.1987743855346075`(2.574327954867511` -0.1542507160492513` Cos[5.960803788824705` t]-0.0006407232443525828` Sin[5.960803788824705` t]))/.{5.960005227008194`->5.96,5.960803788824705`->5.96,5.960234450409174`->5.96}//Expand

-0.000055154563884767866` 0.013031263547355615`-1.849339171142787`*^-6 0.1542507160492513`

0.000055154563884767866` 0.1987743855346075`

-1.849339171142787`*^-6 0.1987743855346075`

Block[{x1,x2,x3},
{x1,x2,x3}=Figure8MMA;
Plot[-0.10039955491885759`(((Norm[x1[t]-x2[t],2])^2)((Norm[x1[t]-x3[t],2])^2)((Norm[x2[t]-x3[t],2])^2))+1.0963314545623987`(((Norm[x1[t]-x2[t],2])^2)+((Norm[x1[t]-x3[t],2])^2)+((Norm[x2[t]-x3[t],2])^2))-3.676012573889878`*^-2(x1'[t].x1'[t]+x2'[t].x2'[t]+x3'[t].x3'[t]),{t,0,6.325},PlotRange->All]
]

Finding periodic orbits

ReturnError[sim_,t0_Real]:=Norm[Flatten@{Through[sim[0]]-Through[sim[t0]],(D[Through[sim[t]],t]/.t->0)-(D[Through[sim[t]],t]/.t->t0)},2]

With[{fig8=MMAIntegrator[3,2,Ones@3,Figure8Init][10.]},
Plot[ReturnError[fig8,t0],{t0,0,10}]
]

MinimumReturnError[sim_,points_Integer]:=
Block[{Ti,Tf,dT,res},
{Ti,Tf}=First@First[sim]["Domain"];
dT=((Tf-Ti)/points);
res=First/@MinimalBy[Table[FindMinimum[{ReturnError[sim,ts],Ti+dT<ts<Tf-dT},{ts,tsi}],{tsi,Ti+dT,Tf-dT,dT}],First]//Min;
(*Print["Return error with p = ",First[sim]'[Ti]," is ",res];*)
res
]

ClearAll@GradientDescent;
GradientDescent[n_Integer,d_Integer,ms_,pot_,{init_,vars_}][T_Real,R_Real][initval_,points_Integer,δ_Real,α_Real,iter_Integer]:=
Block[{state,f,δf,grad,num=Length@vars},
f=MinimumReturnError[NDSIntegrator[n,d,ms,pot,ApplyVariables[{init,vars},initval]][T,R],points];
If[(α/(Sqrt[iter]))<(10^(-5)),Return[{f,First@#->Last@#&/@Transpose@{vars,initval}}]];
(*Print["Current step:\n\tp = ",initval,"\n\tf = ",f,"\n\tα = ",(α/(Sqrt[iter]))];*)
δf=Table[
MinimumReturnError[NDSIntegrator[n,d,ms,pot,ApplyVariables[{init,vars},initval+δ UnitVector[num,i]]][T,R],points]
,{i,num}];
grad= ((δf-f)/δ);
GradientDescent[n,d,ms,pot,{init,vars}][T,R][initval-(α/(Sqrt[iter]))grad,points,δ,α,iter+1]
]

Dynamic@f

GradientDescent[3,2,Ones@3,NewtonianPotential,{SymInitialConditions2D[3][px,py],{px,py}}][10.,2.][{.33,.53},10,(10.^(-3)),(10.^(-4)),1]

Block[{reterror},
reterror[px0_Real,py0_Real]:=MinimumReturnError[NDSIntegrator[3,2,Ones@3,NewtonianPotential,ApplyVariables[{SymInitialConditions2D[3][px,py],{px,py}},{px0,py0}]][10.,2.],10];
FindMinimum[reterror[px0,py0],{{px0,.33},{py0,.53}}]
]

Figures for the problem statement

PlotState[Figure8NSD,5.]

Export["fig8.pdf",%]