我试图用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

当你计算电荷的时候，你取你的方向向量，它指向另一个原子，然后把它乘以一个负数（因为你不匹配的电荷会得到一个负积），结果就是一个力向量指向另一个原子。你知道吗

您还应该考虑使用这些数字进行建模的风险（float的epsilon大约在1e-16）。如果你打算在埃标度上建模，最好以埃为单位。如果你打算在米秤上做模特，你可能想坚持你现有的。只是小心重新调整你的常数。你知道吗

方向向量就是代码问题，如果你解决了这个问题，你就得到了下一个问题；在t=0，你的例子中的原子的加速度是2e19，在t=0之后的第一个点，它们的速度是2e16（假设你的单位是对的，这比光速快一个数量级）。它们移动得如此之快，以至于它们飞向对方，然后又飞过对方，然后静电的平方距离反作用力在第二个滴答声后在功能上变为0，它们永远不会从超扭曲中减速。你知道吗

有几种方法可以解决这个问题：更短的滴答声（飞秒？）你也可以尝试用连续的曲线而不是离散的点来建模，但是如果你试着把它缩放到多个原子的话，那会下降得太快。。。归根结底，这只是物理建模的一个核心问题。祝你好运！你知道吗