我正在尝试在 Python 中模拟电力，但电力似乎总是排斥

Question

编辑：抱歉，我没有考虑编写测试。我会这样做，看看我是否无法找出我做错了什么。感谢建议我编写测试的人！

我正在尝试用 Python 编写一个计算机模拟程序来模拟电力以及原子如何与电力相互作用。对于那些不知道的人，本质上，具有相反（正和负）电荷的东西会吸引，而相同的电荷会排斥，力的大小下降为 1 /（距离的平方）。我尝试将一个带负电的粒子（氧离子）和一个带正电的粒子（氢离子）放入一个坐标系中，并期望它们相互吸引并靠得更近，但实际上它们却是排斥的！当然，我认为这一定是一个错字，所以我在我的电力模型中添加了一个负号，但它们仍然排斥！我不知道发生了什么，并希望你们中的一些人可能对这里出了什么问题有所了解。下面是代码。我已经包含了所有内容，以便您可以从自己的终端运行它，以亲自查看会发生什么。

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import math

Å = 10 ** (-10)
u = 1.660539040 * (10 ** (-27))
k = 8.987551787368 * (10 ** 9)
e = 1.6021766208 * (10 ** (-19))

hydrogen1 = dict(
    name='Hydrogen 1',
    charge=1 * e,
    mass=1.00793 * u,
    position=[0.5 * Å, 0.5 * Å, 0.5 * Å],
    velocity=[0, 0, 0],
    acceleration=[0, 0, 0],
    force=[0, 0, 0]
)
oxygen1 = dict(
    name='Oxygen 1',
    charge=-2 * e,
    mass=15.9994 * u,
    position=[0, 0, 0],
    velocity=[0, 0, 0],
    acceleration=[0, 0, 0],
    force=[0, 0, 0]
)

atoms = [hydrogen1, oxygen1]


def magnitude(vector):
    magnitude = 0
    for coordinate in vector:
        magnitude += (coordinate ** 2)
    return math.sqrt(magnitude)


def scale_vector(vector, scalefactor):
    scaled_vector = vector
    i = 0
    while i < len(vector):
        vector[i] *= scalefactor
        i += 1
    return scaled_vector


def sum_vectors(vectors):
    resultant_vector = [0, 0, 0]
    for vector in vectors:
        i = 0
        while i < len(vector):
            resultant_vector[i] += vector[i]
            i += 1
    return resultant_vector


def distance_vector(point1, point2):
    if type(point1) is list and type(point2) is list:
        pos1 = point1
        pos2 = point2
    elif type(point1) is dict and type(point2) is dict:
        pos1 = point1['position']
        pos2 = point2['position']
    vector = []
    i = 0
    while i < len(pos1):
        vector.append(pos2[i] - pos1[i])
        i += 1
    return vector


def distance(point1, point2):
    return magnitude(distance_vector(point1, point2))


def direction_vector(point1, point2):
    vector = distance_vector(point1, point2)
    length = magnitude(vector)
    return scale_vector(vector, 1 / length)


def eletric_force(obj1, obj2):
    length = k * obj1['charge'] * \
        obj2['charge'] / ((distance(obj1, obj2)) ** 2)
    force_vector = scale_vector(direction_vector(obj1, obj2), length)
    return force_vector


def force_to_acceleration(force, mass):
    scalefactor = 1 / (mass)
    return scale_vector(force, scalefactor)


time = 10

t = 0
period = 1 / 1000

while t < time:
    i = 0
    while i < len(atoms):
        atom = atoms[i]
        position = atom['position']
        velocity = atom['velocity']
        acceleration = atom['acceleration']

        # Moving the atom
        atom['position'] = sum_vectors(
            [position, scale_vector(velocity, period)])

        # Accelerating the atom using its current acceleration vector
        atom['velocity'] = sum_vectors([
            velocity, scale_vector(acceleration, period)])

        # Calculating the net force on the atom
        force = [0, 0, 0]
        j = 0
        while j < len(atoms):
            if j != i:
                force = sum_vectors([force, eletric_force(atoms[i], atoms[j])])
            j += 1

        # Updating the force and acceleration on the atom
        atoms[i]['force'] = [force[0], force[1], force[2]]
        atom['acceleration'] = force_to_acceleration(
            [force[0], force[1], force[2]], atom['mass'])

        i += 1

    t += period

np.random.seed(19680801)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

for atom in atoms:
    name = atom['name']
    position = atom['position']
    X = position[0]
    Y = position[1]
    Z = position[2]
    print(
        f'Position of {name}: [{X}, {Y}, {Z}]')
    color = 'green'
    if 'Oxygen' in atom['name']:
        color = 'red'
    ax.scatter(X, Y, Z, color=color)

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()

以下是测试：

from functions import *
from constants import *
from atoms import hydrogen1, hydrogen2, oxygen1
import math


def test_magnitude():
    assert magnitude({1, 3, -5}) == math.sqrt(35)


def test_sum_vectors():
    assert sum_vectors([[1, 2, 3], [0, -4, 8]]) == [1, -2, 11]


def test_scale_vector():
    assert scale_vector([1, 4, -3], -2) == [-2, -8, 6]


def test_distance_vector():
    assert distance_vector([1, 4, 3], [0, 1, 1]) == [-1, -3, -2]
    assert distance_vector(hydrogen2, hydrogen1) == [Å, Å, Å]


def test_distance():
    assert distance([1, 2, 3], [3, 2, 1]) == math.sqrt(8)
    assert distance(hydrogen1, oxygen1) == math.sqrt(0.75) * Å
    assert distance(hydrogen1, hydrogen2) == Å * math.sqrt(3)


def test_direction_vector():
    assert direction_vector([1, 1, 1], [7, 5, -3]) == [6 /
                                                       math.sqrt(68), 4 / math.sqrt(68), -4 / math.sqrt(68)]
    m = 1 / math.sqrt(3)
    for component in direction_vector(hydrogen2, hydrogen1):
        assert abs(component - m) < 10 ** (-12)


def test_electric_force():
    m = 4.439972744 * 10 ** (-9)
    for component in electric_force(hydrogen1, hydrogen2):
        assert abs(component - m) < 10 ** (-12)


def test_force_to_acceleration():
    assert force_to_acceleration(
        [4, 3, -1], 5.43) == [4 / 5.43, 3 / 5.43, -1 / 5.43]

Answer 1

当您计算电荷时，您将指向伙伴原子的方向矢量乘以一个负数（因为您的不匹配电荷会产生负积），这会导致力矢量指向远离你的伙伴原子。

您还应该考虑使用这些数字进行此建模的风险（浮点数的 epsilon 大约在 1e-16 左右）。如果您打算在埃尺度上建模，最好以埃为单位。如果您打算在米尺上建模，您可能希望坚持使用现有的。请小心重新调整常量。

方向向量是代码问题，如果你解决了这个问题，你就会得到下一个问题；在t=0，您的示例中的原子的加速度为2e19，在t=0 之后的第一点，它们的速度为2e16（假设我的单位正确，是一堆比光速快几个数量级）。它们移动得如此之快，以至于它们飞快地飞向对方，然后相互掠过，然后在第二个滴答声之后，静电的平方反比力在函数上变为 0，它们永远不会从超翘曲中减速。

有一些选项可以解决这个问题；更短的滴答声（飞秒？），更改为相对论速度计算等。您还可以尝试使用连续曲线而不是离散点进行建模，但是如果您尝试将其缩放到多个原子，那么它会下降得如此之快。 ..最终这只是物理建模的核心问题。祝你好运！