如何在电场中投掷带电粒子？

Question

我正在研究带电粒子在潜在场中的动态行为。我尝试定义一个简单的潜在场，像这样：

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

然后我使用梯度函数来计算每个点的电场。

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

这样我得到了一个 (10,10,10) 的三维数组。现在的问题是，当我研究一个有质量的粒子在这个向量场中的运动时，出现了一些困难。

比如，我的力是这样定义的：

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

对于像重力这样的力，你不需要担心太多，我只需要一个一维数组，比如 f_g = np.array([0,0,-g])，这样就可以了。但是我的电场力 f_e 随着粒子的移动而变化。我已经尝试像这样定义 f_e：f_e_x = -Q*r[i][0] # 第i次迭代的位置向量的X分量，这样做有点效果，但这仍然和在一个向量场中投放一个电荷不太一样。有没有什么标准的方法来处理这个问题呢？

Answer 1

这个问题很有意思，但我不太确定你的方法能否找到解决方案。

通过考虑一个10x10x10的网格，你实际上是把物体的位置限制在1000个固定的位置上。但你真正想要的是在任意时刻t，找出物体的具体位置，而不是仅仅局限于这1000个位置。

这个系统通常是通过一个二阶偏微分方程来描述的，粒子在任何时刻的位置取决于它的加速度（x''）、速度（x'）和之前的位置（x，其中x = (x, y, z)). 通常情况下，没有精确的解析解，但通过利用问题中的一些特性，我们可以得到一个数值解。（我觉得在这个情况下可能有解析解，但我这里提供了一个更通用的数值解，这样你可以尝试改变电场的情况。）

由于x、y和z的运动是独立的（在上面的例子中，电场没有交叉项），所以可以有效地将其简化为一系列常微分方程，并使用script.integrate.solve_ivp来求解。

首先，简单回顾一下物理知识：电势V(x,y,z)是一个标量场，在任何一点的值与带电粒子在该点的静电能量成正比。静电场E(x,y,z)是一个矢量场（也就是说，它有大小和方向）。如果你把一个带电量为q的粒子放在这个场中，它会受到一个力F=Eq的作用，其中F的方向与E相同（对于正的q）。

第一步是将离散的电势V场转换为一个连续的函数。我使用了以下方法，我认为这可以重现原问题中V的离散值：

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

接下来，我们需要计算E场。这是场的负梯度。可以通过数值方法近似计算（并结合电荷来产生力），方法如下：

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

然后，我们需要建立微分方程：x'' = E(x) q（假设质量为1）。在scipy中，初值问题求解器只能解决一阶常微分方程，所以我们需要做一个小技巧，把二阶方程转换为两个一阶方程。在一维情况下，如果我们定义v_x == x'，那么微分方程变为：

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

推广到三维，我们现在有六个偏微分方程：

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

其中i, j, k是沿着x、y和z轴的单位向量。代码如下：

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])

这里的state是一个包含三维位置和速度的六元组。可以使用以下方法求解并绘制结果：

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()

注意，函数的前两行设置了初始参数，应该根据需要进行更改，以便探索不同的解决方案。同时，查看一下solve_ivp的文档，以获取你所需的确切数据。

完整的代码在下面，还有一个一维版本，可能对你入门有帮助。

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()

参考资料： 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