How to simulate motion of 4 particles moving towards each other?

Question

I have a problem that I need to simulate where there are four particles located at the corners of a square room with known lenght a. They are starting to move simultaneously and are supposed to move towards the closest neighbour until one of them reaches some other. They all have constant velocity V. I am supposed to use scipy.integrate.RK45 function.

I am not sure how to describe mathematically the fact that they are all moving towards the closest neightbour. I suppose the function should look something like this (for one particles for now):

def motion_eq(t, y):
        ???

r = np.array([0, 0], 'float')
v = np.array([0, 1], 'float')

dt = 0.1
a = 2

motion_solver = scipy.integrate.RK45(motion_eq, 0, np.hstack([r, v]),
    t_bound = np.inf, first_step = dt, max_step = dt)

particle, *_ = plt.plot(*r.T, 'o')
plt.gca().set_aspect(1)
plt.xlim([-a, a])
plt.ylim([-a, a])
def update():
    motion_solver.step()
    r = motion_solver.y[0:2]
    particle.set_data(*r.T)
    plt.draw()
timer = plt.gcf().canvas.new_timer(interval = 50)
timer.add_callback(update)
timer.start()

plt.show()

Could somebody please give a hint of how to do this? Thanks in advance.

Answer 1

Notes

• We consider a 2D plane for the calculation of the equation of motion.
• We assume four particles, each of which has three neighboring particles, in this case.
• We find the average position (pos) of the four particles (assumed random positions for each particle initially).

Other calculations are in this image:

Code

import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt


def closest_neighbor_dir(r, a):
    all_ptls = np.array([[a, a], [a, -a], [-a, -a], [-a, a]])
    ptls = np.delete(all_ptls, np.where((all_ptls == r).all(axis=1)), axis=0)
    pos = np.mean(ptls, axis=0)
    dir_vec = pos - r
    return dir_vec / np.linalg.norm(dir_vec)


def motion_eq(t, y, a, V):
    r, v = y[:2], y[2:]
    accel = V * closest_neighbor_dir(r, a) - v
    return np.hstack([v, accel])


def update():
    for i in range(4):
        motion_solver[i].step()
        r[i] = motion_solver[i].y[0:2]
    particles.set_data(r[:, 0], r[:, 1])
    plt.draw()


r = np.array([[1.1, 1.7], [-1, 1.3], [-1.2, -1.2], [1.4, -1.7]], dtype=float)
v = np.array([[0, 0], [0, 0], [0, 0], [0, 0]], dtype=float)

dt, a, V = 0.1, 2, 0.2

motion_solver = [scipy.integrate.RK45(lambda t, y: motion_eq(t, y, a, V), 0, np.hstack([r[i], v[i]]),
                                      t_bound=np.inf, first_step=dt, max_step=dt) for i in range(4)]

particles, *_ = plt.plot(*r.T, 'o')
plt.gca().set_aspect(1)
plt.xlim([-a, a])
plt.ylim([-a, a])

timer = plt.gcf().canvas.new_timer(interval=50)
timer.add_callback(update)
timer.start()

plt.show()