Interacting Particles (Part 3)

APPH E3200x
Columbia University

Introduction

After study of the motion of one and two interacting particles, this notebook takes a look at many interacting particles. Computationally, the solution to many particles can use the same algorithm used for two particles. But, without careful programming, the compute time becomes very long. illustrate Newton’s laws of motion.

The calculations below illustrate the motion of many particles interacting with a central force. Once we define the equations of motion (EOM), then the solution is computed efficiently from any initial condition.

Where are the particles?

A particle is represented in phase-space by their position vector and their velocity vector. In three dimensions, each particle is represented at each instant of time with six numbers: three numbers for position and three numbers for velocity.

We let the variable “p” represent a list of the particles phase-space locations.

(* regular dynamics *) nBodies = 2;

(* chaotic dynamics *) nBodies = 8;

p = Table[{x[i][t], y[i][t], z[i][t], vx[i][t], vy[i][t],vz[i][t]}, {i, nBodies}] // Flatten;

p (* for two particles, this is 12 functions of time *)

What are the forces?

For a central force, the force is directed along a vector from one point to another and inversely proportional to some power of the distance between particles. For classical gravity and for Coulomb forces, the force scales as the inverse square of the force between particles.

r1 = {x1,y1,z1};
r2 = {x2,y2,z2};

(* distance = *) Sqrt[(r1 - r2).(r1 - r2)]

distance[i_,j_] := Sqrt[({x[j][t],y[j][t],z[j][t]} - {x[i][t],y[i][t],z[i][t]}).({x[j][t],y[j][t],z[j][t]} - {x[i][t],y[i][t],z[i][t]})]

distance[1,2]

(* vector from r1 towards r2 = *) r2 - r1

force[i_,j_] := ({x[j][t],y[j][t],z[j][t]} - {x[i][t],y[i][t],z[i][t]})/distance[i,j]^3

force[1,2] + force[2,1] // Simplify

Why don’t we include the self-force?

Random Initial Conditions (expressed in the COM frame of reference)

For simplicity, we will locate the particles randomly, and then remove the COM motion.

How many initial conditions are required?

v0=Table[RandomReal[{-1.,1.},3],{i,nBodies}]

vcom = MapThread[Plus, v0]/nBodies

r0=Table[RandomReal[{-1.,1.},3],{i,nBodies}];

rcom = MapThread[Plus, r0 ]/nBodies

p0 = Table[{x[i][0] == r0[[i,1]] - rcom[[1]],
y[i][0] == r0[[i,2]] - rcom[[2]],
z[i][0] == r0[[i,3]] - rcom[[3]], vx[i][0] == v0[[i,1]] - vcom[[1]], vy[i][0] == v0[[i,2]] - vcom[[2]],vz[i][0] == v0[[i,3]]- vcom[[3]]}, {i,  nBodies}] // Flatten;

p0

Newton’s Equations of Motion

Complement[Range[nBodies], {1}]
Complement[Range[nBodies], {2}]

(Plus @@ (force[1,#]& /@ Complement[Range[nBodies], {1}]))

(* rate of change in phase-space *) dp = Table[{vx[i][t], vy[i][t], vz[i][t]}~Join~ (Plus @@ (force[i,#]& /@ Complement[Range[nBodies], {i}])), {i,  nBodies}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

(* How many equations? *) Dimensions[eqs]

Example Solution

maxTime = 20.0;

(sol = NDSolve[eqs , p, {t, 0, maxTime}] // Flatten;) // Timing

Visualizing the Solution

x[1][t] /. sol

Plot[{Evaluate[x[1][t] /. sol],Evaluate[x[2][t] /. sol]}, {t, 0, maxTime},PerformanceGoal->"Speed" ]

threeRange = (10.0{{-1,1},{-1,1},{-1,1}});

orbitPlot = ParametricPlot3D[{Evaluate[{x[i][t],y[i][t],z[i][t]}/. i->1/. sol],Evaluate[{x[i][t],y[i][t],z[i][t]}/. i->2/. sol]},{t,0,maxTime},PlotPoints->200, Axes -> True, AxesLabel -> {"x","y","z"}, AxesOrigin-> {0,0,0}, BoxStyle -> LightGray, LabelStyle-> {Bold},PlotLabel->"Center of Mass Frame"]

objects[tt_] := Graphics3D[Flatten[Table[{PointSize[0.02], Point[{x[i][t], y[i][t], z[i][t]}]/. sol /. t -> tt}, {i,1,nBodies}]],PlotRange -> threeRange,  Axes -> True, AxesLabel -> {"x","y","z"}, AxesOrigin-> {0,0,0}, BoxStyle -> LightGray, LabelStyle-> {Bold}]

Show[orbitPlot,objects[1.0]]

Animation

Animate[Show[orbitPlot,objects[t], PlotRange -> (10.0{{-1,1},{-1,1},{-1,1}})], {t, 0.0, maxTime, 0.1},AnimationRunning-> False]

Distance from Origin

r[i_][t_] := Sqrt[({x[i][t],y[i][t],z[i][t]} /. sol).({x[i][t],y[i][t],z[i][t]} /. sol)];

Plot[{r[1][t],r[2][t]}, {t, 0, maxTime}, PlotLabel -> "Distance from COM"]

Plot[Table[r[i][t],{i, nBodies}], {t, 0, maxTime}, PlotLabel -> "Distance from COM", PlotRange -> All]

Conservation Principles

Average Position and Velocity

We can find the average position and the average velocity.

The average position is a vector...

avgPosition[tt_] = (Plus @@ Table[{x[i][t], y[i][t],z[i][t]}/. sol /. t -> tt, {i, nBodies}])/nBodies;

avgPosition[0.2]
avgPosition[1.0]
avgPosition[maxTime]

ListPointPlot3D[Table[avgPosition[t],{t,0,maxTime,0.1}],PlotStyle->{Red,PointSize[0.02]},PlotLabel->"Average Position", PlotRange->(1.0*^-12 {{-1,1},{-1,1},{-1,1}}),Axes -> True, BoxRatios -> Automatic, AxesLabel -> {"x","y","z"}, AxesOrigin-> {0,0,0}]

avgVelocity[tt_] = (Plus @@ Table[{vx[i][t], vy[i][t],vz[i][t]}/. sol /. t -> tt, {i, nBodies}])/nBodies;

avgVelocity[0.0]
avgVelocity[1.0]
avgVelocity[maxTime]

Total Energy

The total energy is the sum of the kinetic energy and the potential energy.

totalKineticEnergy[tt_] := Sum[0.5 (vx[i][t]^2 + vy[i][t]^2 + vz[i][t]^2) /. sol /. t -> tt, {i, nBodies}]

Plot[totalKineticEnergy[t], {t, 0,maxTime},PerformanceGoal->"Speed" ]

potential[i_,j_] := -1/distance[i,j]

(Plus @@ (potential[1,#]& /@ Complement[Range[nBodies], {1}]))

(Plus @@ (potential[2,#]& /@ Complement[Range[nBodies], {2}]))

When summing the potential energy, we need to avoid “double counting”.

totalPotentialEnergy[tt_] := Sum[0.5(Plus @@ (potential[i,#]& /@ Complement[Range[nBodies], {i}]))/. sol /. t -> tt, {i, nBodies}]

Plot[totalPotentialEnergy[t], {t, 0,maxTime},PerformanceGoal->"Speed" ]

initialTotalEnergy = (totalKineticEnergy[0] + totalPotentialEnergy[0])

Plot[100(totalKineticEnergy[t] + totalPotentialEnergy[t]-initialTotalEnergy)/initialTotalEnergy, {t, 0,maxTime}, PlotLabel -> "Percent Variation of Total Energy" , PlotRange -> All,PerformanceGoal->"Speed"]

Angular Momentum

The angular momentum is the vector L = r × v

angularMomentum[tt_] := Sum[{x[i][t],y[i][t],z[i][t]}\[Cross]{vx[i][t],vy[i][t],vz[i][t]} /. sol /. t -> tt, {i, nBodies}]

angularMomentum[0.0]
angularMomentum[1.0]
angularMomentum[10.0]

Sensitivity to Initial Conditions

What happens when we make a single small shift?

v01= v0 + ({{1.0*^-4,0,0}}~Join~Table[{0,0,0},{i,2,nBodies}]);

vcom1 = MapThread[Plus, v01]/nBodies

vcom - vcom1

p01 = Table[{x[i][0] == r0[[i,1]] - rcom[[1]],
y[i][0] == r0[[i,2]] - rcom[[2]],
z[i][0] == r0[[i,3]] - rcom[[3]], vx[i][0] == v01[[i,1]] - vcom1[[1]], vy[i][0] == v01[[i,2]] - vcom1[[2]],vz[i][0] == v01[[i,3]]- vcom1[[3]]}, {i,  nBodies}] // Flatten;

eqs1 = Thread[D[p,t] == dp]~ Join ~ p01;

(sol1 = NDSolve[eqs1, p, {t, 0, maxTime}] // Flatten;) // Timing

Display the difference in motion…

Plot[{Evaluate[x[1][t] /. sol1],Evaluate[x[2][t] /. sol1]}, {t, 0, maxTime},PerformanceGoal->"Speed" ]

Plot[{Evaluate[(x[1][t] /. sol1)-(x[1][t] /. sol)],Evaluate[(x[2][t] /. sol1) - (x[2][t] /. sol)]}, {t, 0, maxTime},PerformanceGoal->"Speed" ]

LogPlot[{Evaluate[Abs[(x[1][t] /. sol1)-(x[1][t] /. sol)]],Evaluate[Abs[(x[2][t] /. sol1) - (x[2][t] /. sol)]]}, {t, 1.0*^-6, 2.0},PerformanceGoal->"Speed" , PlotRange -> {1.0*^-7, 1.0}, PlotLabel -> "Deviation Grows Exponentially"]

δr[t_] := Sqrt[(({x[1][t],y[1][t],z[1][t]} /. sol1)-({x[1][t],y[1][t],z[1][t]} /. sol)).(({x[1][t],y[1][t],z[1][t]} /. sol1)-({x[1][t],y[1][t],z[1][t]} /. sol))];

LogPlot[δr[t], {t, 1.0*^-6, 2.0},PerformanceGoal->"Speed" , PlotRange -> {1.0*^-7, 1.0}, PlotLabel -> "Separation Grows Exponentially"]

Summary

Newton’s laws of motion completely describe the motion of interacting particles. The description of motion is best described in COM coordinates.

In three dimensions, two interacting particles follow elliptical orbits, as first described by Kepler. However, with more than two particles, particle motion is complex. In practice, predicting the motion of many interacting particles is not possible because the dynamics depends very closely on initial conditions. 

Total angular momentum, total linear momentum, and total energy are constant.