Verlet Integration in Python resulting in particles running away

Question

My fallowing code (was supposed to) solve the equation of motion for two bodies but the result is the particles running way and I wasn't able to find where is the error

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('ggplot')

DIM = 2
N = 2
ITER = 1000

def acc(r, v, a):
    for i in range(N - 1):
        for j in range(i+1, N):
            r2 = 0.0
            rij = [0.0, 0.0]
            for k in range(DIM):
                rij[k] = r[i][k] - r[j][k]
                r2 += rij[k] * rij[k]
            fg = -1.0 /np.power(np.sqrt(r2), 3)
            for k in range(DIM):
                a[i][k] += fg * rij[k]
                a[j][k] -= fg * rij[k]
    return a

def verlet(r, v, a, dt):
    for i in range(N):
        for k in range(DIM):
            v[i][k] += 0.5 * a[i][k] * dt 
            r[i][k] += v[i][k] * dt
    a = acc(r, v, a)
    for i in range(N):
        for k in range(DIM):
            v[i][k] += 0.5 * a[i][k] * dt
    return [r,v]


r = np.zeros((N, DIM))
v = np.zeros((N ,DIM))
a = np.zeros((N, DIM))  

r[0] = [0.5,0.0]
v[0] = [0.0,1.0]

r[1] = [-0.5,0.0]
v[1] = [0.0,-1.0]

dt = 0.01
plt.ion()
for i in range(ITER):
    r,v = verlet(r, v, a, dt)
    plt.scatter(r[0][0], r[0][1])
    plt.scatter(r[1][0], r[1][1],color='yellow')
    plt.pause(0.00005)

And I used the algorithm described in velocity Verlet

Welcome to StackOverflow. Please read and follow the posting guidelines in the help documentation. [Minimal, complete, verifiable example](http://stackoverflow.com/help/mcve) applies here. We cannot effectively help you until you post your code and accurately describe the problem. In particular, provide the failing output, your debugging trace(s), and any analysis or references you have on the problem. — Prune, Jun 13 '16 at 23:04
Hi, actually "code error", the error is in the logic (since 2 particles shouldn't accelerate and go away). I was hopping if someone could find were the error. — 0xhfff, Jun 13 '16 at 23:09

score 1 · Accepted Answer · answered Jul 06 '16 at 14:32

Acceleration doesn't accumulate over time, the way velocity and position do: it should be computed from scratch at each step of the process. So,

Remove a from the argument list of both acc and verlet
Initialize a by zeros at the beginning of acc.
Move the call a = acc(r, v) to the beginning of verlet.

You will see the articles rotating around each other, as expected.

correct motions

Unrelated to your issue but important for effective use of NumPy: get rid of loops over vector coordinates, leave it for NumPy to add and subtract vectors. I rewrote acc method in this way:

def acc(r, v):
    a = np.zeros((N, DIM))  
    for i in range(N - 1):
        for j in range(i+1, N):
            rij = r[i, :] - r[j, :]
            r2 = np.dot(rij, rij)
            fg = -1.0 /np.power(np.sqrt(r2), 3)
            a[i, :] += fg * rij
            a[j, :] -= fg * rij
    return a

If you move `a = acc(r, v)` per 3., then you change the Verlet method into a symplectic Euler method, i.e., reduce the error order from 2 to 1. — Lutz Lehmann, Sep 09 '16 at 14:33

Verlet Integration in Python resulting in particles running away

1 Answers1