Electrostatic Force - Simulate Trajectory of Test Particle using Runge Kutta - Force always Repels

In the center of a 2D-Plane a positive static charge Q is placed with position r_prime. This charge creates a static electrical Field E.

Now i want to place a test particle with charge Q and position vector r in this static E-Field and compute its trajectory using the 4th order Runge-Kutta method.

For the initial conditions

• Q = 1, r_prime=(0,0)
• q = -1, r = (9, 0), v = (0,0)

one would expect, that the negative charged test particle should move towards the positive charge in the center. Instead i get the following result for the time evolution of the test particles x component:

[9.0,
 9.0,
 8.999876557604697,
 8.99839964155741,
 8.992891977287334,
 8.979313669171093,
 8.95243555913327,
 8.906134626441052,
 8.83385018027209,
 8.729257993736123,
 8.587258805422984,
 8.405449606446608,
 8.186368339303788,
 7.940995661159361,
 7.694260386250479,
 7.493501689700884,
 7.420415546859942,
 7.604287312065716,
 8.226733652039988,
 9.498656905483394,
 11.60015461031076,
 14.621662121713964,
 18.56593806599109,
 ....

The results of the first iteration steps show the correct behavior, but then the particle is strangely repelled to infinity. There must be a major flaw in my implementation of the Runge-Kutta Method, but I checked it several times and I cant find any…

Could someone please take a quick look over my implementation and see if they can find a bug.

"""
Computes the trajectory of a test particle with Charge q with position vector r = R[:2] in a 
electrostatic field that is generated by charge Q with fixed position r_prime
"""

import numpy as np
import matplotlib.pyplot as plt

def distance(r, r_prime, n):
    
    """
    computes the euclidean distance to the power n between position x and x_prime
    """
    return np.linalg.norm(r - r_prime)**n


def f(R, r_prime, q, Q):
    
    """ 
    The equations of motion for the particle is given by: 
    
    d^2/dt^2 r(t) = F = constants * q * Q * (r - r_prime)/||r - r_prime||^3   
        
    To apply the Runge-Kutta-Method we transform the above (constants are set to 1) 
    two dimensional second order ODE into four one dimensional ODEs:
        
    d/dt r_x = v_x
    d/dt r_y = v_y
    d/dt v_x = q * Q * (r_x - r_prime_x)/||r - r_prime||^3
    d/dt v_y = q * Q * (r_y - r_prime_y)/||r - r_prime||^3 '''  
    """
    
    r_x, r_y = R[0], R[1]
    v_x, v_y = R[2], R[3]
    
    dist = 1 / distance(np.array([r_x, r_y]), r_prime, 3)

    # Right Hand Side of the 1D Odes
    f1 = v_x
    f2 = v_y
    f3 = q * Q * dist * (r_x - r_prime[0]) 
    f4 = q * Q * dist * (r_y - r_prime[1])
        
    return np.array([f1, f2, f3, f4], float)

# Constants for the Simulation
a = 0.0       # t_0
b = 10.0       # t_end
N = 100    # number of iterations
h = (b-a) / N # time step

tpoints = np.arange(a,b+h,h)

# Create lists to store the computed values
xpoints, ypoints = [], []
vxpoints, vypoints = [], []


# Initial Conditions

Q, r_prime = 1, np.array([0,0], float)  # charge and position of particle that creates the static E-Field

q, R = -1, np.array([9,0,0,0], float)    # charge and its initial position + velocity r=R[:2], v=[2:] 


for dt in tpoints:
    xpoints.append(R[0])
    ypoints.append(R[1])
    vxpoints.append(R[2])
    vypoints.append(R[3])
    
    # Runge-Kutta-4th Order Method
    k1 = dt * f(R, r_prime, q, Q)
    k2 = dt * f(R + 0.5 * k1, r_prime, q, Q)
    k3 = dt * f(R + 0.5 * k2, r_prime, q, Q)
    k4 = dt * f(R + k3, r_prime, q, Q)
    R += (k1 + 2*k2 * 2*k3 + k4)/6
    
plt.plot(tpoints, xpoints) # should converge to 0

Edit 09.02.2021

The equation of motion for the test particle with charge q and position vector r(t) is given by

$$frac{d^2}{dt^2} mathbf{r}(t) = frac{mathbf{F}}{m} = frac{qmathbf{E}}{m} = frac{k}{m} frac{qcdot Q}{|mathbf{r} – mathbf{r}’|^3} cdot (mathbf{r} – mathbf{r}’) qquad (1)$$

Charge Q generates an electrostatic field E and has a constant position r’. For simplicity we set the constant factor k, m equal to one.

Equation (1) is a two dimensional second order ODE. To apply Runge Kutta we transform it into four one dimensional ODEs:

$$begin{align}
frac{d}{dt}r_x &= v_x qquad &(2)\
frac{d}{dt}r_y &= v_y qquad &(3)\
frac{d}{dt}v_x &= frac{qcdot Q}{|mathbf{r} – mathbf{r}’|^3} cdot (r_x – r’_x) qquad &(4) \
frac{d}{dt}v_y &= frac{qcdot Q}{|mathbf{r} – mathbf{r}’|^3} cdot (r_y – r’_y) qquad &(5)
end{align}$$

Except for the 1/6 factor I forgot in the last RK4 step (thanks for pointing it out Lutz), everything looks correct from my point of view.

I mean, the system behaves correctly at the beginning, the negative test charge is attracted by the positive one. But then the charge is repelled for reasons I don’t understand.