Interacting Particles

APPH E3200x
Columbia University

Introduction

Two interacting particles illustrate Newton’s laws of motion. By computing the solution for interacting particles in three dimensions, we also learn about vectors and coordinates.

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

Question: how would you extend these calculations to an arbitrary number of interacting particles?

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.

p = Table[{x[i][t], y[i][t], z[i][t], vx[i][t], vy[i][t],vz[i][t]}, {i, 1, 2}] // 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?

Initial Conditions

For simplicity, we will locate the particles randomly.

How many initial conditions are required?

p0=Flatten[Table[{x[i][0]==RandomReal[{-1,1}],y[i][0]==RandomReal[{-1,1}],
z[i][0]==RandomReal[{-1,1}],
vx[i][0]==RandomReal[{-1,1}],vy[i][0]==RandomReal[{-1,1}],
vz[i][0]==RandomReal[{-1,1}]},{i,1,2}]];

p0

Newton’s Equations of Motion

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

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

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

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

Example Solution

maxTime = 10.0;

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

Visualizing the Solution

Plot[x[1][t] /. sol1, {t, 0, maxTime}]

threeRange = {{-3.0, 3.0}, {-3.0, 3.0}, {-3.0, 3.0}};

objects[tt_] := Graphics3D[{PointSize[0.04],Red, Point[{x[1][t], y[1][t], z[1][t]}], PointSize[0.04],Blue,Point[{x[2][t], y[2][t], z[2][t]}] }/. sol1 /. t -> tt, PlotRange->threeRange, Axes -> True, AxesLabel -> {"x","y","z"}, AxesOrigin-> {0,0,0}, BoxStyle -> LightGray, LabelStyle-> {Bold}]

Show[objects[1.0]]

Animation

Animate[objects[t], {t, 0.0, maxTime, 0.1},AnimationRunning-> False]

Averages

We can find the average radius and the average position.

avgRadius[tt_] = (Plus @@Table[Sqrt[x[i][t]^2 + y[i][t]^2 + z[i][t]^2] /. sol1/. t -> tt, {i, 1, 2}] )/2;

avgRadius[0.2]

Plot[avgRadius[t],{t,0,maxTime},PlotLabel->"Average Radius"]

The average position is a vector...

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

avgPosition[0.2]

ListPointPlot3D[Table[avgPosition[t],{t,0,maxTime,0.1}],PlotStyle->{Blue,PointSize[0.02]},PlotLabel->"Average Position", PlotRange->3 threeRange,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]}/. sol1 /. t -> tt, {i, 1, 2}])/2;

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

ListPointPlot3D[Table[avgVelocity[t],{t,0,maxTime,0.1}],PlotStyle->{Blue,PointSize[0.02]},PlotLabel->"Average Velocity", PlotRange->3 threeRange,Axes -> True, BoxRatios -> Automatic, AxesLabel -> {"x","y","z"}, AxesOrigin-> {0,0,0}]

Removing the Center of Mass

Since the mass of each object is identical, we can easily sum (and average) the initial velocities and positions of the objects in order to find the center of gravity and the center of mass velocity.

Defining the COM

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

vcom = MapThread[Plus, v0]/(2)

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

rcom = MapThread[Plus, r0 ]/(2)

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, 1, 2}] // Flatten;

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

Example Solution

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

Animation

objects2[tt_] := Graphics3D[{PointSize[0.04],Red, Point[{x[1][t], y[1][t], z[1][t]}], PointSize[0.04],Blue,Point[{x[2][t], y[2][t], z[2][t]}] }/. sol2 /. t -> tt, PlotRange->threeRange, Axes -> True, AxesLabel -> {"x","y","z"}, AxesOrigin-> {0,0,0}, BoxStyle -> LightGray, LabelStyle-> {Bold}]

Animate[objects2[t], {t, 0.0, maxTime, 0.1},AnimationRunning-> False]

ParametricPlot3D[{Evaluate[{x[i][t],y[i][t],z[i][t]}/.i->1/.sol2],Evaluate[{x[i][t],y[i][t],z[i][t]}/.i->2/.sol2]},{t,0,maxTime},PlotPoints->200,PlotLabel->"Center of Mass Frame"]

Averages

We can find the average radius and the average position.

avgRadius[tt_] = (Plus @@Table[Sqrt[x[i][t]^2 + y[i][t]^2 + z[i][t]^2] /. sol2/. t -> tt, {i, 1, 2}] )/2;

Plot[avgRadius[t],{t,0,maxTime},PlotLabel->"Average Radius"]

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

ListPointPlot3D[Table[avgPosition[t],{t,0,5,0.1}],PlotStyle->{Blue,PointSize[0.02]},PlotLabel->"Average Position"]

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

avgVelocity[0.0]
avgVelocity[maxTime]

Summary

Newton’s laws of motion completely describe the motion of interacting particles. 

In three dimensions, two interacting particles follow elliptical orbits, as first described by Kepler. 

In general, interacting particles have a relative position with respect to the average position.