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

Question

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{q\mathbf{E}}{m} = \frac{k}{m} \frac{q\cdot 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{q\cdot Q}{|\mathbf{r} - \mathbf{r}'|^3} \cdot (r_x - r'_x) \qquad &(4) \\\\ \frac{d}{dt}v_y &= \frac{q\cdot 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.

Answer 1

The problem is that, first, the vector field has a singularity at the origin. Moving straight towards the singularity results in a catastrophe in finite time. Second and related, that the speed and all stiffness measures (higher derivatives) increase significantly the closer you get to the origin. So while one would expect a Kepler ellipse (after setting the initial velocity to some non-zero value in $y$ direction), the distortion due to the truncation error produces more like a swing-by maneuver with an extreme gain in velocity. Your step size is so large that you do not see this, this happens in the internal stages of the steps.

To get better results, you would need more control on the speed. One variant is to normalize the vector field that you feed to RK4 to something like $\hat v=\frac{v}{1+\|v\|}$. Here you would need to reconstruct the points on the time axis by integrating the scale factor.

The other variant is to use a method with adapting step size such as all the usual library methods.