这是用python ODEint正确实现点电荷动态的

自从我在本学期的物理II课中了解了点电荷以后，我希望能够研究静态力和场分布，而不是调查带电粒子运动的实际轨迹。 这样做的第一步是建立一个模拟n个点粒子动力学的朴素引擎。 我已经在Python中使用矩阵实现了解决方案，并希望有人能评论我是否正确地做了这些。 由于我不知道会有什么样的动态，所以我无法直接从视频中得知我的方程实现是正确的。

我的特殊问题

特别是，我无法确定在计算力的大小时，我是否正确地计算了1/r^(3/2)因子。 为什么？ 因为当我模拟偶极子并使用$ 2/2 $作为指数时，粒子就开始进入椭圆轨道。 这是我所期望的。 但是，当我使用正确的指数时，我得到这个： 我的代码在哪里出错？ 我应该期待什么

我将首先写下我正在使用的等式：

鉴于n收费q_1, q_2, ..., q_n ，与群众m_1, m_2, ..., m_n位于初始位置r_1, r_2, ..., r_n ，用速度(d/dt)r_1, (d/dt)r_2, ..., (d/dt)r_n由q_j在q_i引起的力由下式给出

F_(j -> i) = k(q_iq_j)/norm(r_i-r_j)^{3/2} * (r_i - r_j)

现在，粒子$ q_i $的净边际力是作为成对力的总和给出的

F_(N, i) = sum_(j != i)(F_(j -> i))

然后粒子$ q_i $的净加速度只是通过粒子的质量来归一化力：

(d^2/dt^2)r_i = F_(N, i)/m_i

总的来说，对于n粒子，我们有一个微分方程的n阶系统。 我们还需要指定n初始粒子速度和n初始位置。

为了在Python中实现这一点，我需要能够计算成对点距离和成对收费倍数。 要做到这一点，我要平铺电荷的q矢量和位置的r矢量，分别取其产物和它们的转置差异。

def integrator_func(y, t, q, m, n, d, k):
    y = np.copy(y.reshape((n*2,d)))

    # rj across, ri down
    rs_from = np.tile(y[:n], (n,1,1))

    # ri across, rj down
    rs_to = np.transpose(rs_from, axes=(1,0,2))

    # directional distance between each r_i and r_j
    # dr_ij is the force from j onto i, i.e. r_i - r_j
    dr = rs_to - rs_from

    # Used as a mask to ignore divides by zero between r_i and r_i
    nd_identity = np.eye(n).reshape((n,n,1))

    # WHAT I AM UNSURE ABOUT
    drmag = ma.array(
        np.power(
            np.sum(np.power(dr, 2), 2)
        ,3./2)
        ,mask=nd_identity)

    # Pairwise q_i*q_j for force equation
    qsa = np.tile(q, (n,1))
    qsb = np.tile(q, (n,1)).T
    qs = qsa*qsb

    # Directional forces
    Fs = (k*qs/drmag).reshape((n,n,1))

    # Dividing by m to obtain acceleration vectors
    a = np.sum(Fs*dr, 1)

    # Setting velocities
    y[:n] = np.copy(y[n:])

    # Entering the acceleration into the velocity slot
    y[n:] = np.copy(a)

    # Flattening it out for scipy.odeint to work properly
    return np.array(y).reshape(n*2*d)

def sim_particles(t, r, v, q, m, k=1.):
    """
    With n particles in d dimensions:

    t: timepoints to integrate over
    r: n*d matrix. The d-dimensional initial positions of n particles
    v: n*d matrix of initial particle velocities
    q: n*1 matrix of particle charges
    m: n*1 matrix of particle masses
    k: electric constant.
    """

    d = r.shape[-1]
    n = r.shape[0]
    y0 = np.zeros((n*2,d))
    y0[:n] = r
    y0[n:] = v
    y0 = y0.reshape(n*2*d)
    yf = odeint(
        integrator_func,
        y0,
        t,
        args=(q,m,n,d,k)).reshape(t.shape[0],n*2,d)
    return yf

链接地址: http://www.djcxy.com/p/65411.html

上一篇: Is this a proper implementation of point charge dynamics with python ODEint

下一篇: Solar System Simulator Physics Integration Issues (Unreal Engine 4, C++)