Is this a proper implementation of point charge dynamics with ODEs

Question

Since learning about point charges in my physics II class this semester, I want to be able to investigate not only the static force and field distributions but the actual trajectories of movement of electrically charged particles. The first stage in doing this is to build a naive engine for simulating the dynamics of $n$ individual point particles. I've implemented the solution using matrices in python and was hoping someone could comment on whether I've done so correctly. As I don't know what kind of dynamics to expect, I can't tell directly from the videos that my implementation of my equations is correct.

My Particular Problem

In particular, I cannot tell if in my calculation of Force magnitude I am computing the $\frac{1}{r^{3/2}}$ factor correctly. Why? because when I simulate a dipole and use $2/2$ as an exponent the particles start going in an elliptical orbit. which is what I would expect. However, when I use the correct exponent, I get this: Where is my code going wrong? What am I supposed to expect

I'll first write down the equations I'm using:

Given $n$ charges $q_1, q_2, ..., q_n$, with masses $m_1, m_2, ..., m_n$ located at initial positions $\vec{r}_1, \vec{r}_2, ..., \vec{r}_n$, with velocities $\dot{\vec{r}}_1, \dot{\vec{r}}_2, ..., \dot{\vec{r}}_n$ the force induced on $q_i$ by $q_j$ is given by

$$\vec{F}_{j \to i} = k\frac{q_iq_j}{\left|\vec{r}_i-\vec{r}_j\right|^2}\hat{r}_{ij}$$

where

$$\hat{r}_{ij} =\frac{ \vec{r}_i-\vec{r}_j}{\left|r_i-r_j\right|}$$

Now, the net marginal force on particle $q_i$ is given as the sum of the pairwise forces

$$\vec{F}_{N, i} = \sum_{j \ne i}{\vec{F}_{j \to i}}$$

And then the net acceleration of particle $q_i$ just normalizes the force by the mass of the particle:

$$\ddot{\vec{r}}_i = \frac{\vec{F}_{N, i}}{m_i}$$

In total, for $n$ particles, we have an n-th order system of differential equations. We will also need to specify $n$ initial particle velocities and $n$ initial positions.

To implement this in python, I need to be able to compute pairwise point distances and pairwise charge multiples. To do this I tile the $\mathbf{q}$ vector of charges and the $\mathbf{r}$ vector of positions and take, respectively, their product and difference with their transpose.

<!-- language: python -->

def integrator_func(y, t, q, m, n, d, k):
    y = np.copy(y.reshape((n*2,d)))

    # rj across, ri down
    rs_from = np.tile(y[:n], (n,1,1))

    # ri across, rj down
    rs_to = np.transpose(rs_from, axes=(1,0,2))

    # directional distance between each r_i and r_j
    # dr_ij is the force from j onto i, i.e. r_i - r_j
    dr = rs_to - rs_from

    # Used as a mask to ignore divides by zero between r_i and r_i
    nd_identity = np.eye(n).reshape((n,n,1))

    # WHAT I AM UNSURE ABOUT
    drmag = ma.array(
        np.power(
            np.sum(np.power(dr, 2), 2)
        ,3./2)
        ,mask=nd_identity)

    # Pairwise q_i*q_j for force equation
    qsa = np.tile(q, (n,1))
    qsb = np.tile(q, (n,1)).T
    qs = qsa*qsb

    # Directional forces
    Fs = (k*qs/drmag).reshape((n,n,1))

    # Dividing by m to obtain acceleration vectors
    a = np.sum(Fs*dr, 1)

    # Setting velocities
    y[:n] = np.copy(y[n:])

    # Entering the acceleration into the velocity slot
    y[n:] = np.copy(a)

    # Flattening it out for scipy.odeint to work properly
    return np.array(y).reshape(n*2*d)

def sim_particles(t, r, v, q, m, k=1.):
    """
    With n particles in d dimensions:

    t: timepoints to integrate over
    r: n*d matrix. The d-dimensional initial positions of n particles
    v: n*d matrix of initial particle velocities
    q: n*1 matrix of particle charges
    m: n*1 matrix of particle masses
    k: electric constant.
    """

    d = r.shape[-1]
    n = r.shape[0]
    y0 = np.zeros((n*2,d))
    y0[:n] = r
    y0[n:] = v
    y0 = y0.reshape(n*2*d)
    yf = odeint(
        integrator_func,
        y0,
        t,
        args=(q,m,n,d,k)).reshape(t.shape[0],n*2,d)
    return yf

Answer 1

It's an interesting question - how does a dipole actually move? It seems to me you're not entirely sure what to expect, so we need to get a solid test case where we understand what's actually going on.

On a related note, do try to add your imports and initial conditions the next time :)

For now, let's assume both masses equal and opposite charges. We can then test out circular motion (think two point bodies acting on each other by gravity, I'm sure you can see the similarity): $$ F = \frac{mv_{orbital}^2}{r} = \frac{k q^2}{r^2} \implies v_{orbital} = \sqrt{\frac{k q^2}{m r}} $$ So it's easy to see that if we see all the variables on the right equal to 1 (charges equal to $\pm 1$), the orbital velocity should be 1. So if I run these initial conditions:

t = np.linspace(0, 5, 500)
N = 2 #particles
d = 2 #dimensions

q = np.array([-1.,1.])
m = np.ones_like(q)
k = 1.
r = np.array([[0.5,0],[-0.5,0]]) # initial distance equal to 1
distance = np.linalg.norm(r[0] - r[1])
print(distance) # should return 1
v_orbital_value = (k*q**2/m/distance)**0.5
print(v_orbital_value) #should also return 1
v = np.array([[0.,1.],[0.,-1.]]) * v_orbital_value

y = sim_particles(t, r, v, q, m, k)

I get, seemingly disappointingly...

And then you realize that neither particle is stationary, and $v_{orbital}$ should refer to their relative velocity. Upon dividing it by 2, we get

Both particles exhibit circular motion, as it should be if we're to have a stationary center of mass (which you can check easily, I'll just leave this line here):

plt.plot(0.5*(x1+x2), 0.5*(y1+y2), "go")

To sum up, your code is quite fine (and written pretty neatly, I must add)! Although the ellipsoids you get from the $\frac{1}{r}$ force law look pretty, they're not what we were expecting. Always have some preexisting idea of how to test your code's results.