Modeling an orbital motion

From now on, we will only use velocity Verlet method for solving Newton’s equation of motion. We generalize the method for systems of multiple variables. These systems can have a particle in 2D or 3D, or multiple particles.

# Velocity Verlet method for high dimensional systems (2D, 3D, or multiple particles)
velocity_verlet_generic <- function(f, t0, x0, v0, t.total, dt, ...){
  # f:  2nd derivative function
  # t0: initial time
  # x0: a vector of initial position x
  # v0: a vector of initial velocity v
  # t.total: total simulation time
  # dt: time step size 
  t_all = seq(t0, t.total, by=dt)
  nt_all = length(t_all)
  nx = length(x0)
  x_all = matrix(0, nrow = nt_all, ncol = nx)
  v_all = matrix(0, nrow = nt_all, ncol = nx)
  x_all[1,] = x0
  v_all[1,] = v0 
  for (i in 1:(nt_all-1)) {
    v_half = v_all[i,] + 0.5 * dt * f(t_all[i], x_all[i,],...)
    x_all[i+1,] = x_all[i,] + dt * v_half
    v_all[i+1,] = v_half + 0.5 * dt * f(t_all[i+1], x_all[i+1,],...)
  }
  return(cbind(t_all, x_all, v_all))   # the output is a matrix of t, x, v for all time steps
}

In this new example, a 2D-particle is driven by a gravitational force from the center (e.g., Earth)

\[f(r) = \frac{G}{r} \tag{3}\] \(G\) is a constant. The direction of the force is from the particle to the center. We can write the \(x\) and \(y\) components of the force, \(f_x\) and \(f_y\), as

\[\begin{equation} \begin{cases} f_x(x, y) = - \frac{Gx}{r^2} \tag{4}\\ f_y(x, y) = - \frac{Gy}{r^2} \end{cases} \end{equation}\]

# gravitational force in 2D
f_g2d <- function(t, x, G){
  # x: a vector of size two: x, y positions
  r2 = x[1]**2 + x[2]**2
  return(-x*G/r2)
}

results_g2d = velocity_verlet_generic(f = f_g2d, t0 = 0, x0 = c(4, 0), v0 = c(0, 1), t.total = 100, dt = 0.1, G = 1)
plot(results_g2d[,1], results_g2d[,2], type = "l", col=2,
    xlab="t", ylab="position", xlim=c(0,100), ylim=c(-8,8)) 
lines(results_g2d[,1], results_g2d[,3], col=4)
legend("topleft", inset=0.02, 
       legend = c("x","y"),
       col=c(2,4), lty=1, cex=0.8)

plot(results_g2d[,2], results_g2d[,3], type = "l", col=1,
    xlab="x", ylab="y", xlim=c(-8,8), ylim=c(-8,8)) 

results_g2d2 = velocity_verlet_generic(f = f_g2d, t0 = 0, x0 = c(2, 2), v0 = c(-1, 1), t.total = 1000, dt = 0.1, G = 1)
plot(results_g2d2[,1], results_g2d2[,2], type = "l", col=2,
    xlab="t", ylab="position", xlim=c(0,100), ylim=c(-8,8)) 
lines(results_g2d2[,1], results_g2d2[,3], col=4)
legend("topleft", inset=0.02, 
       legend = c("x","y"),
       col=c(2,4), lty=1, cex=0.8)

plot(results_g2d2[,2], results_g2d2[,3], type = "l", col=1,
    xlab="x", ylab="y", xlim=c(-8,8), ylim=c(-8,8))

Two-body problem

We then consider two particles of equal masses in 2D, a system that can be described by 4 variables.

# two particles in 2D with gravitational forces
f_2body_2d <- function(t, x, G){
  # x: a vector of size six: (x1, y1, x2, y2)
  r12_2 = (x[1] - x[3])**2 + (x[2] - x[4])**2
  
  f1_x = (x[3]-x[1])/r12_2
  f1_y = (x[4]-x[2])/r12_2
  f2_x = (x[1]-x[3])/r12_2
  f2_y = (x[2]-x[4])/r12_2

  return(G*c(f1_x, f1_y, f2_x, f2_y))
}

set.seed(1)  # we use random initial conditions

results_2body = velocity_verlet_generic(f = f_2body_2d, t0 = 0, x0 = runif(4, -5, 5), v0 = runif(4, -0.3, 0.3),
                                        t.total = 1000, dt = 0.1, G = 1)
plot(results_2body[,2], results_2body[,3], type = "l", col=2,
    xlab="x", ylab="y")
lines(results_2body[,4], results_2body[,5], col=4)
legend("topleft", inset=0.02, 
       legend = c("1","2"),
       col=c(2,4), lty=1, cex=0.8)

The time trajectories contain patterns but with a constant shift, because of a non-zero total momentum. To eliminate the constant shift (consider that the system is in another reference with a constant velocity), we can shift the velocities with constant values, so that the system has zero total momentum.

# velocity is shifted to remove net changes 
center_v <- function(v, ndim){
  # v: a vector of velocities 
  # ndim: dimension of the system
  v_mat = matrix(v, ncol = ndim, byrow = T)
  v_means = colMeans(v_mat)
  return(v - v_means)  # different sizes, the shorter vector iterates
}

set.seed(1)  # we use random initial conditions
v0 = runif(4, -0.3, 0.3)

results_2body = velocity_verlet_generic(f = f_2body_2d, t0 = 0, x0 = runif(4, -5, 5), v0 = center_v(v0,2),
                                        t.total = 1000, dt = 0.1, G = 1)
plot(results_2body[,2], results_2body[,3], type = "l", col=2,
    xlab="x", ylab="y")
lines(results_2body[,4], results_2body[,5], col=4)
legend("topleft", inset=0.02, 
       legend = c("1","2"),
       col=c(2,4), lty=1, cex=0.8)

An R package gganimate allows to generate movies of time trajectories! (Commented as they are slow to execute)

{r,fig.width = 5, fig.height = 5}
library(ggplot2)
library(gganimate)
library(gifski)
library(png)
npoints = nrow(results_2body)
results_2body_dframe = data.frame(t = results_2body[,1], x = c(results_2body[,2],results_2body[,4]),
                                  y = c(results_2body[,3],results_2body[,5]),
                                  particle = c(rep("1", npoints), rep("2", npoints)))

p = ggplot(results_2body_dframe, aes(x = x, y = y, colour = particle)) + 
  geom_point(aes(colour = particle), show.legend = FALSE, size = 5,alpha = 0.7) + labs(x = "x", y = "y") +
  scale_color_manual(values = c("1" = "red", "2" = "blue"))
p
anim = p + transition_time(t, range = c(0, 50)) + labs(title = "t: {frame_time}")
anim = animate(anim, fps = 20, nframes = 200, end_pause = 0, rewind = F, renderer = gifski_renderer())
anim_save("2-body.gif", animation = anim)

Three-body problem

Now we consider three particles of equal masses in 2D, a system that can be described by 6 variables.

# three particles in 2D with gravitational forces
f_3body_2d <- function(t, x, G){
  # x: a vector of size six: (x1, y1, x2, y2, x3, y3)
  r12_2 = (x[1] - x[3])**2 + (x[2] - x[4])**2
  r13_2 = (x[1] - x[5])**2 + (x[2] - x[6])**2
  r23_2 = (x[3] - x[5])**2 + (x[4] - x[6])**2
  
  f1_x = (x[3]-x[1])/r12_2 + (x[5]-x[1])/r13_2 
  f1_y = (x[4]-x[2])/r12_2 + (x[6]-x[2])/r13_2 
  f2_x = (x[1]-x[3])/r12_2 + (x[5]-x[3])/r23_2 
  f2_y = (x[2]-x[4])/r12_2 + (x[6]-x[4])/r23_2 
  f3_x = (x[1]-x[5])/r13_2 + (x[3]-x[5])/r23_2 
  f3_y = (x[2]-x[6])/r13_2 + (x[4]-x[6])/r23_2 
  return(G*c(f1_x, f1_y, f2_x, f2_y, f3_x, f3_y))
}

set.seed(1)  # we use random initial conditions
v0 = runif(6, -0.01, 0.01)

results_3body = velocity_verlet_generic(f = f_3body_2d, t0 = 0, x0 = runif(6, -5, 5), v0 = center_v(v0, 3),
                                        t.total = 1000, dt = 0.1, G = 1)
plot(results_3body[,2], results_3body[,3], type = "l", col="red",
    xlab="x", ylab="y")
lines(results_3body[,4], results_3body[,5], col="orange")
lines(results_3body[,6], results_3body[,7], col="blue")
legend("topright", inset=0.02, 
       legend = c("1","2", "3"),
       col=c("red","orange","blue"), lty=1, cex=0.8)

Animations too (Commented as they are slow to execute)

{r,fig.width = 5, fig.height = 5}

npoints = nrow(results_3body)
results_3body_dframe = data.frame(t = results_3body[,1], x = c(results_3body[,2],results_3body[,4], results_3body[,6]),
                                  y = c(results_3body[,3],results_3body[,5], results_3body[,7]),
                                  particle = c(rep("1", npoints), rep("2", npoints), rep("3", npoints)))

p = ggplot(results_3body_dframe, aes(x = x, y = y, colour = particle)) + 
  geom_point(aes(colour = particle), show.legend = FALSE, size = 5,alpha = 0.7) + labs(x = "x", y = "y") +
  scale_color_manual(values = c("1" = "red", "2" = "orange", "3" = "blue"))
p
anim = p + transition_time(t, range = c(0, 300)) + labs(title = "t: {frame_time}")
anim = animate(anim, fps = 20, nframes = 300, end_pause = 0, rewind = F, renderer = gifski_renderer())
anim_save("3-body.gif", animation = anim)