你好

在这里，我们描述某个领域的工作，然后做出一些漂亮的功能（这里的一切都很简单）。

本文将发生什么。

一般情况：

我们描述了基础，即使用向量（对于手头没有麻木的人使用的自行车）
我们描述了相互作用的实质点和领域

特殊情况（基于一般情况）：

让我们可视化电磁场强度的矢量场（第一张和第三张图片）
我们将可视化粒子在电磁场中的运动

切面下见我！

理论背景编程

向量

所有基础的基础都是一个向量（尤其是在我们的案例中）。因此，我们将从向量的描述开始。我们需要什么？向量，距离，模块和一些技术问题上的算术运算。我们将从列表继承的向量。这是其初始化的样子：

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

我们还需要两种方法来更快地生成矢量和实现漂亮的输出。这里没有棘手的问题，因此这里是Vector类的完整代码：

所有向量类代码

 class Vector(list): def __init__(self, *el): for e in el: self.append(e) 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 def __sub__(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 def __mul__(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 def __truediv__(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 def __pow__(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 def __mod__(self, other): return sum((self - other) ** 2) ** 0.5 def mod(self): return self % Vector.emptyvec(len(self)) def dim(self): return len(self) def __str__(self): if len(self) == 0: return "Empty" r = [str(i) for i in self] return "< " + " ".join(r) + " >" def _ipython_display_(self): print(str(self)) @staticmethod def emptyvec(lens=2, n=0): return Vector(*[n for i in range(lens)]) @staticmethod def randvec(dim): return Vector(*[random.random() for i in range(dim)])

物质点

理论上，这里的一切都很简单-该点具有坐标，速度和加速度（为简单起见）。此外，她还有一个质量和一组自定义参数（例如，对于电磁场，电荷不会伤害我们，但没有人打扰您设置至少一个自旋）。

初始化如下：

 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): #        self.acc = self.acc + force / self.mass def clean_acc(self): self.acc = self.acc * 0

做得好，观点本身就做好了。

点代码（输出效果不错）

 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]) 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): self.acc = self.acc + force / self.mass def clean_acc(self): self.acc = self.acc * 0 def __str__(self): r = ["Point {"] for p in self.__params__: r.append(" " + p + " = " + str(getattr(self, p))) r += ["}"] return "\n".join(r) def _ipython_display_(self): print(str(self))

结果：

互动领域

我们将交互场称为包含所有物质点的集合并对其施加作用力的对象。我们将考虑我们奇妙的宇宙的特殊情况，因此我们将进行一个自定义的交互（当然，这很容易扩展）。声明一个构造函数并添加一个点：

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

现在，有趣的部分是声明一个函数，该函数将在此时返回“张力”。尽管此概念指的是电磁相互作用，但在我们的情况下，它是一些抽象矢量，我们将沿着该矢量移动该点。在这种情况下，我们将拥有点q的属性，在特定情况下-点的电荷（通常-我们想要的任何东西，甚至是向量）。那么C点的张力是多少？像这样：

v e c E （ v e c C ） = s u m v e c F_{i} （ v e c C ）

$\ vec E（\ vec C）= \ sum {\ vec F_i（\ vec C）}$

也就是说，此时的张力

v e c C

$\ vec 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) #       (   ,       F,      ) for p in self.points: if coord % p.coords < 10 ** (-10): #      ,        ,   ,     (  ) 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) * pq) p.accelerate(dt) p.move(dt)

这里的一切都很简单。对于每个点，我们确定这些坐标中的张力，然后确定在ETU物料点上的最终力：

v e c F = v e c E * q

$\ vec F = \ vec E * q$

定义缺少的功能。

所有交互字段代码

 class InteractionField: def __init__(self, F): self.points = [] self.F = F def move_all(self, dt): for p in self.points: p.move(dt) def intensity(self, coord): proj = Vector(*[0 for i in range(coord.dim())]) single_point = Point(Vector(), mass=1.0, q=1.0) for p in self.points: if coord % p.coords < 10 ** (-10): continue d = p.coords % coord fmod = self.F(single_point, p, d) * (-1) proj = proj + (coord - p.coords) / d * fmod return proj def step(self, dt): self.clean_acc() for p in self.points: p.accinc(self.intensity(p.coords) * pq) p.accelerate(dt) p.move(dt) def clean_acc(self): for p in self.points: p.clean_acc() def append(self, *args, **kwargs): self.points.append(Point(*args, **kwargs)) def gather_coords(self): return [p.coords for p in self.points]

特例。粒子运动和矢量场可视化

因此，我们来到了最有趣的地方。让我们从...开始

模拟电磁场中粒子的运动

 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。关于其他所有内容，我认为这很清楚。

r ** 2 + 0.1

-这是避免被零除的问题。我们当然可能会感到困惑，解决了一个很大的扩散系统，但是首先，对于两个以上的物体，没有运动方程，其次，“初学者的文章”概念中显然没有包含这个方程。

接下来，我们沿每个轴添加十个点（二维空间），坐标从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()

应该发生什么：

实际上，绘图是完全随机的，因为在力学发展的时刻，每个点的轨迹都是不可预测的。

矢量场可视化

这里的一切都很简单。我们需要经过一些步骤，并在正确的方向上在每个坐标上绘制一个向量。

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

大概应该得出这个结论。

您可以加长向量本身，将* 4替换为* 1.5：

我们玩的是立体感和造型

创建一个具有200个点的5维空间，并且该交互不取决于距离的平方，而取决于4度。

 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]) #       velocities.append(velmod) plt.plot(velocities) plt.show()

这是在任何给定时间所有速度总和的图表。如您所见，随着时间的推移，它们正在缓慢加速。

好吧，这是关于如何使这种简单的事情的简短说明。但是，如果您玩花会发生什么：

所有代码与演示

 import random class Vector(list): def __init__(self, *el): for e in el: self.append(e) 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 def __sub__(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 def __mul__(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 def __truediv__(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 def __pow__(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 def __mod__(self, other): return sum((self - other) ** 2) ** 0.5 def mod(self): return self % Vector.emptyvec(len(self)) def dim(self): return len(self) def __str__(self): if len(self) == 0: return "Empty" r = [str(i) for i in self] return "< " + " ".join(r) + " >" def _ipython_display_(self): print(str(self)) @staticmethod def emptyvec(lens=2, n=0): return Vector(*[n for i in range(lens)]) @staticmethod def randvec(dim): return Vector(*[random.random() for i in range(dim)]) 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]) 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): self.acc = self.acc + force / self.mass def clean_acc(self): self.acc = self.acc * 0 def __str__(self): r = ["Point {"] for p in self.__params__: r.append(" " + p + " = " + str(getattr(self, p))) r += ["}"] return "\n".join(r) def _ipython_display_(self): print(str(self)) class InteractionField: def __init__(self, F): self.points = [] self.F = F def move_all(self, dt): for p in self.points: p.move(dt) def intensity(self, coord): proj = Vector(*[0 for i in range(coord.dim())]) single_point = Point(Vector(), mass=1.0, q=1.0) for p in self.points: if coord % p.coords < 10 ** (-10): continue d = p.coords % coord fmod = self.F(single_point, p, d) * (-1) proj = proj + (coord - p.coords) / d * fmod return proj def step(self, dt): self.clean_acc() for p in self.points: p.accinc(self.intensity(p.coords) * pq) p.accelerate(dt) p.move(dt) def clean_acc(self): for p in self.points: p.clean_acc() def append(self, *args, **kwargs): self.points.append(Point(*args, **kwargs)) def gather_coords(self): return [p.coords for p in self.points] #  import matplotlib.pyplot as plt import numpy as np import time #     if False: u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 2 + 0.1)) for i in range(10): u.append(Vector.randvec(2) * 10, q=(random.random() - 0.5) / 2) 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() def sigm(x): return 1 / (1 + 1.10 ** (-x/1000)) #    if False: 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) 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() * 1.5 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() #    5-  if False: 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]) velocities.append(velmod) plt.plot(velocities) plt.show()

下一篇文章可能会涉及更复杂的建模，以及流体和Navier-Stokes方程。
UPD：我的同事在这里写的一篇文章

感谢MomoDev帮助渲染视频。

python + matplotlib中的自定义物理交互的基本模拟

理论背景编程

向量

物质点

互动领域

特例。粒子运动和矢量场可视化

模拟电磁场中粒子的运动

矢量场可视化

我们玩的是立体感和造型

More articles:

python + matplotlib中的自定义物理交互的基本模拟

理论背景编程

向量

物质点

互动领域

特例。 粒子运动和矢量场可视化

模拟电磁场中粒子的运动

矢量场可视化

我们玩的是立体感和造型

More articles:

特例。粒子运动和矢量场可视化