pythonsimulationrunge-kutta

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

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.

Solution

  • Notes

    Other calculations are in this image:

    enter image description here

    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()