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)
plt.plot(tpoints, xpoints) # should converge to 0
