How to throw a charged particle in a electric vector field?

Question

I am trying to study the dynamics of a charged particle inside potential fields. I tried defining a dummy potential like this:

x = np.linspace(-20, 20, 10)
y = np.linspace(-20, 20, 10)
z = np.linspace(0, 10, 10)
X, Y, Z = np.meshgrid(x, y, z)
V = X**2 + Y**2 + Z**2

and I used gradient function to get Electric field at each point.

Ex, Ey, Ez = np.gradient(V)

This gives me a (10,10,10) 3d array. Now the problem arises when I am studying the motion of a particle with some mass in this vector field.

For example, my force is defined like this:

f_e[i][0] = -Q * Ex_at_r  # Ex but at that r_x value particle is currently
f_e[i][1] = -Q * Ey_at_r
f_e[i][2] = -Q * Ez_at_r

When it comes to forces like gravity you don't need to worry about all this I can just have a 1-D array like this f_g = np.array([0,0,-g]) and call it day, but my f_e is changing as the particle moves. I already tried defining f_e like this, f_e_x = -Q*r[i][0] # X_component of position vector for ith iteration, and this sort of works but this still isn't same as having a vector field and a charge thrown in it. Is there a standard way to go about this?

Answer 1

This is an interesting problem, but I'm not sure you're approach is going to get to a solution.

By considering a 10x10x10 grid, you are forcing the body to only be at one of 1000 discrete positions. What you really want is to figure out where the body is at time t, without restricting its position to only having one of 1000 discrete values.

This is system is, in general, described by a second-order partial differential equation (PDE), where the position of a particle at any time depends upon its acceleration (x''), velocity (x'), and previous position (x, where x = (x, y, z)). In general, there is no exact analytical solution but, by making use of some features in the problem, we can get a numeric solution. (I think there is likely an analytic solution in this case, but I've presented a more general numeric solution here so you can play with changing the V-field.)

As the x, y, and z motions are independent (there are no cross-terms in the V field in the above example), it can effectively be reduced to a series of ordinary differential equations and solved using script.integrate.solve_ivp.

First, a physics recap: the electric potential, V(x,y,z), is a scalar field whose value at any point is proportional to the electrostatic energy a charged particle has at that point. The electrostatic field, E(x,y,z), is a vector field (i.e. it has a size and direction). If you put a particle with charge q in the field, it will experience a force F=Eq, where F points in the direction of E (for positive q).

The first step is to convert the discrete V field into a continuous function. I used the following which I think will re-create the discrete values of V in the original question:

def potential_3d(pos: np.array) -> float:
    x, y, z = pos
    vx = (4 * x - 20) ** 2
    vy = (4 * y - 20) ** 2
    vz = z ** 2
    return vx + vy + vz

Next, we need to calculate the E field. This is the negative of the grad of the field. This can be numerically approximated (and combined with the charge to produce the force) by:

from scipy.optimize import approx_fprime

def e_force_3d(pos: np.array, q: float) -> np.array:
    e_field = -approx_fprime(pos, potential_3d)
    return q * e_field

Then, we need to set up the differential equation: x'' = E(x) q (assuming mass of 1). The initial value problem solver in scipy can only solve first order ODEs, so we need to perform a trick to convert the second order equation into two first order equations. In one dimension, if we define v_x == x' then the differential equations become:

1. x' = v_x
2. v_x' = E(x) q

Generalising to 3 dimensions, we now have six partial differential equations:

1. x' = v_x
2. y' = v_y
3. z' = v_z
4. v_x' = E (x, y, z). i q
5. v_y' = E (x, y, z). j q
6. v_z' = E (x, y, z). k q

where i, j, k are unit vectors along the x, y, and z axis. Coding:

def particle_deriv_3d(t, state, q) -> np.array:
    x, y, z, vx, vy, vz = state
    dpos_dt = np.array([vx, vy, vz])
    dv_dt = e_force_3d(np.array([x, y, z]), q)
    return np.concatenate([dpos_dt, dv_dt])

where state is a 6-tuple containing the 3D position and velocity. This can be solved and plotted using:

from scipy.integrate import solve_ivp


def solver_3d():
    x, y, z = 5, 1, 0
    vx, vy, vz = 3, 0, 0.25
    initial_state = [x, y, z, vx, vy, vz]

    states = solve_ivp(particle_deriv_3d, [0, 10], initial_state, args=(0.5,))
    ax = plt.figure().add_subplot(projection='3d')
    ax.plot(states.y[0], states.y[1], states.y[2])
    plt.show()

Note the first two lines of the function set the initial parameters and should be changed to investigate different solutions. Also, check out the documentation for solve_ivp to get the exact data you require.

The complete code is below, together with a 1-D version that might be useful to get started.

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import approx_fprime
import matplotlib.pyplot as plt


# 1D version
def potential_1d(x: float) -> float:
    # x, y, z = pos
    vx = (4 * x - 20) ** 2
    return vx

def e_force_1d(x: float, q: float) -> np.array:
    e_field = -approx_fprime(x, potential_1d)
    return q * e_field

def particle_deriv_1d(t, state, q) -> np.array:
    x, vx = state
    dx_dt = np.array([vx])
    dvx_dt = e_force_1d(x, q)
    return np.concatenate([dx_dt, dvx_dt])

def solver_1d():
    initial_state = [0, 0]

    states = solve_ivp(particle_deriv_1d, [0, 10], initial_state, args=(0.25,))
    plt.plot(states.t, states.y[0])
    plt.show()


# 3D version
def potential_3d(pos: np.array) -> float:
    x, y, z = pos
    vx = (4 * x - 20) ** 2
    vy = (4 * y - 20) ** 2
    vz = z ** 2
    return vx + vy + vz

def e_force_3d(pos: np.array, q: float) -> np.array:
    e_field = -approx_fprime(pos, potential_3d)
    return q * e_field

def particle_deriv_3d(t, state, q) -> np.array:
    x, y, z, vx, vy, vz = state
    dpos_dt = np.array([vx, vy, vz])
    dv_dt = e_force_3d(np.array([x, y, z]), q)
    return np.concatenate([dpos_dt, dv_dt])

def solver_3d():
    x, y, z = 5, 1, 0
    vx, vy, vz = 3, 0, 0.25
    initial_state = [x, y, z, vx, vy, vz]

    states = solve_ivp(particle_deriv_3d, [0, 10], initial_state, args=(0.5,))
    ax = plt.figure().add_subplot(projection='3d')
    ax.plot(states.y[0], states.y[1], states.y[2])
    plt.show()


if __name__ == '__main__':
    # solver_1d()
    solver_3d()

References: https://pythonnumericalmethods.berkeley.edu/notebooks/chapter22.06-Python-ODE-Solvers.html

https://en.wikipedia.org/wiki/Electric_potential

https://en.wikipedia.org/wiki/Electric_field