觉得物理太难了？用Python和Matplotlib进行模拟（物理模拟专家）

全文共4206字，预计学习时长8分钟

觉得物理太难了？本文将模拟N维空间中的某些向量场（例如电磁场），以便你更清晰地理解这些概念。

理论基础向量

任何物理场景的基本元素都是向量。我们需要什么？需要向量的算术运算、距离、模块以及一些技术层面的东西。从列表中得到向量。下面是它初始化的样子：

class Vector(list): def __init__(self, *el): for e in el: self.append(e)

现在，创造一个向量：

v = Vector(1, 2, 3)

再进行加法算术运算：

class Vector(list): ... def __add__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self other

结果：

v1 = Vector(1, 2, 3) v2 = Vector(2, 57, 23.2) v1 v2 >>> [3, 59, 26.2]

同样需定义所有的算术运算（向量的完整代码放在最后）。现在需要一个距离函数。创建dist (v1, v2)并不难—但也并不美观，因此重新定义 %运算符。

class Vector(list): ... def __mod__(self, other): return sum((self - other) ** 2) ** 0.5

结果：

v1 = Vector(1, 2, 3) v2 = Vector(2, 57, 23.2) v1 % v2 >>> 58.60068258988115

还需要一些方法来快速生成向量并且完美地输出。这没有什么棘手的问题，因此以上是向量类的整个代码。

粒子

理论上而言，这里的一切都很简单—点具有坐标、速度及加速度。另外，它有大量的自定义参数组（例如，对于电磁场，可以设置电荷）。

下面进行初始化：

class Point: def __init__(self, coords, mass=1.0, q=1.0 speed=None, **properties): self.coords = coords if speed is None: self.speed = Vector(*[0 for i in range(len(coords))]) else: self.speed = speed self.acc = Vector(*[0 for i in range(len(coords))]) self.mass = mass self.__params__ = ["coords", "speed", "acc", "q"] list(properties.keys()) self.q = q for prop in properties: setattr(self, prop, properties[prop])

以下方法可以将点移动、固定及加速：

class Point: ... def move(self, dt): self.coords = self.coords self.speed * dt def accelerate(self, dt): self.speed = self.speed self.acc * dt def accinc(self, force): # Considering exerted force the point obtains acceleration self.acc = self.acc force / self.mass def clean_acc(self): self.acc = self.acc *

很好，这一部分已经完成。

相互作用场

相互作用场可称为物体，它包含来自太空的所有粒子，并且作用于这些粒子。我们将考虑宇宙的一种特殊情况，因此将进行自定义交互（当然，这易于扩展）。表明构造函数并添加点。

class InteractionField: def __init__(self, F): # F - is a custom force, F(p1, p2, r), p1, p2 - points, r - distance inbetween self.points = [] self.F = F def append(self, *args, **kwargs): self.points.append(Point(*args, **kwargs))

现在，有趣的是声明一个函数，该函数在那一点返回“tension”。尽管该概念原指电磁场交互，但是在我们的情况中，它指某种抽象向量，可依据它来移动点。这样将拥有q点的属性，在特殊情况下，我们将拥有点的电荷（通常是我们想要的任何东西，甚至向量）。所以，C点的张力是多少？大概是这样：

C点的张力

C点的电场强度等于作用于某个单位点上所有物质点的力之和。

class InteractionField: ... def intensity(self, coord): proj = Vector(*[0 for i in range(coord.dim())]) single_point = Point(Vector(), mass=1.0, q=1.0) # That's our "Single point" for p in self.points: if coord % p.coords < 10 ** (-10): # Check whether we compare coord with a point P where P.coords = coord continue d = p.coords % coord fmod = self.F(single_point, p, d) * (-1) proj = proj (coord - p.coords) / d * fmod return proj

这时，你已经可以将向量场可视化，但通常最后完成它。下面进行交互。

class InteractionField: ... def step(self, dt): self.clean_acc() for p in self.points: p.accinc(self.intensity(p.coords) * p.q) p.accelerate(dt) p.move(dt)

对于每个点，确定这些坐标的强度，进而确定对该粒子上的最终作用力：

确定对该粒子的最终作用力

质点运动和向量场可视化

最后，来到了最有趣的部分。开始吧！

模拟电磁场中的粒子运动

u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 2 0.1)) for i in range(3): u.append(Vector.randvec(2) * 10, q=random.random() - 0.5)

实际上，系数K应该等于数十亿(9 * 10 ^ (- 9))，但由于在时间t -> 0之前会失超，所以将它们都设为正数会更容易。因此，本物理学中，K=300,000.

接下来，沿着每个轴增加10个点（二维空间），即在坐标轴上从0到10，也给每个点设置-0.25至0.25的电荷。然后，运行一个循环，并依据其坐标（轨迹）绘制点：

X, Y = [], [] for i in range(130): u.step(0.0006) xd, yd = zip(*u.gather_coords()) X.extend(xd) Y.extend(yd) plt.figure(figsize=[8, 8]) plt.scatter(X, Y) plt.scatter(*zip(*u.gather_coords()), color="orange") plt.show(

其结果应该如下：

确定对该粒子的最终作用力

事实上，其绘图是完全随机的，因为目前（2019年）每个点的轨迹是不可预测的。

向量场可视化

需要通过一些步骤来确定坐标，并且在每个坐标上绘制正确的向量。

fig = plt.figure(figsize=[5, 5]) res = [] STEP = 0.3 for x in np.arange(0, 10, STEP): for y in np.arange(0, 10, STEP): inten = u.intensity(Vector(x, y)) F = inten.mod() inten /= inten.mod() * 4 # длина нашей палочки фиксирована res.append(([x - inten[0] / 2, x inten[0] / 2], [y - inten[1] / 2, y inten[1] / 2], F)) for r in res: plt.plot(r[0], r[1], color=(sigm(r[2]), 0.1, 0.8 * (1 - sigm(r[2])))) # Цвет по хитрой формуле чтобы добиться градиента plt.show()

应该已经获得类似这样的：

第一次向量场可视化

你可以加长向量本身，用* 1.5替换* 4：

向量长度增加后的向量可视化程度

变更维度

创建具有200个点的五维空间和一个交互，该交互依赖于第四度而不是距离的平方。

u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 4 0.1)) for i in range(200): u.append(Vector.randvec(5) * 10, q=random.random() - 0.5)

五维已定义所有的坐标、速度等。

模拟如下：

velmod = 0 velocities = [] for i in range(100): u.step(0.0005) velmod = sum([p.speed.mod() for p in u.points]) # Adding sum of modules of all the velocities velocities.append(velmod) plt.plot(velocities) plt.show()