换个姿势学物理!用Python和Matplotlib进行模拟

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

换个姿势学物理!用Python和Matplotlib进行模拟

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

换个姿势学物理!用Python和Matplotlib进行模拟

理论基础

向量

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

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点的张力是多少?大概是这样:

换个姿势学物理!用Python和Matplotlib进行模拟

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)

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

换个姿势学物理!用Python和Matplotlib进行模拟

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

质点运动和向量场可视化

最后,来到了最有趣的部分。开始吧!

模拟电磁场中的粒子运动

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(

其结果应该如下:

换个姿势学物理!用Python和Matplotlib进行模拟

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

事实上,其绘图是完全随机的,因为目前(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()

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

换个姿势学物理!用Python和Matplotlib进行模拟

第一次向量场可视化

你可以加长向量本身,用* 1.5替换* 4:

换个姿势学物理!用Python和Matplotlib进行模拟

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

变更维度

创建具有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()

换个姿势学物理!用Python和Matplotlib进行模拟

给定时间的速度总和

这是在任意给定时间中的速度总和。正如你所看到的,随着时间的推移,它们在缓慢加速。

以上是关于简单模拟基础物理的简短操作说明。

换个姿势学物理!用Python和Matplotlib进行模拟

留言 点赞 关注

我们一起分享AI学习与发展的干货

如需转载,请后台留言,遵守转载规范

相关推荐