python - 如何计算任何2D多边形的重力？

Question

（为简化起见，在2d中工作）我知道由于重力彼此作用在两个球体上的力为
G(m1*m2/r**2)
但是，对于非球形物体，我找不到能够计算相同力的算法或公式。我最初的想法是将圆形填充到对象中，以使重力产生的力等于每个圆形的力之和。例如（伪代码），

def gravity(pos1,shape):
     circles = packCircles(shape.points)
     force = 0
     for each circle in circles:
           distance = distanceTo(pos1,circle.pos)
           force += newtonForce(distance,shape.mass,1) #1 mass of observer
     return force 

这将是一个可行的解决方案吗？如果是这样，我将如何高效，快速地打包圈子？如果没有，是否有更好的解决方案？

编辑：请注意我要如何在特定点上找到物体的力，因此将需要计算圆和观察者之间的角度（并将矢量相加）。它不同于找到施加的总力。

Answer 1

背景

某些解释可能会有些偏离主题，但我认为有必要帮助澄清评论中提到的某些内容，因为其中许多内容都是违反直觉的。

引力相互作用的这种解释取决于point masses的概念。假设您有两个点质量，它们在一个隔离的系统中彼此隔开一定的距离r1，质量分别为m1和m2，



由m1创建的gravitational field由



其中G是universal gravitational constant，r是到m1的距离，r是沿着m1和m2的直线的单位方向。

由该场施加在m2上的重力为



注意-重要的是，这对于任意距离的任何两点质量都是正确的。1

重力相互作用的场性质允许我们采用superposition来计算由于多重相互作用而产生的净重力。考虑是否要在之前的场景中添加另一个质量m3，



那么，作用在质量m2上的引力就是彼此产生的场的引力之和，



与ri，j = rj，i。这适用于任何间距的任何数量的质量。这也意味着，由质量集合创建的字段可以通过vector sum进行汇总，如果您更喜欢这种形式主义的话。

现在考虑是否有大量的点质量M聚集在一起，形成密度均匀的连续刚体。然后，由于合计质量M，我们想计算单个空间上不同的点质量m的重力。



然后，除了考虑点质量外，我们还可以考虑大小不同的质量的区域（或体积），然后对这些面积（或体积）对点质量的影响进行积分或求和。在二维情况下，重力的大小为



其中σ是总质量的密度。2等效于对每个微分质量σdxdy引起的重力矢量场求和。这种等效性至关重要，因为这意味着对于质量分布之外足够远的任何点质量，由于质量分布而产生的重力几乎与位于 质量分布3 4

这就是说，非常近似地，在计算由于任何质量分布引起的引力场时，可以在分布质量中心将质量分布替换为等效质量点质量。这适用于任何数量的空间上不同的质量分布，无论这些分布是否构成刚体。此外，这意味着您甚至可以将分布组聚合到系统质心处的单个点质量中。5只要参考点足够远即可。

但是，为了找到由于任意点的质量分布而在点质量上产生的重力，对于形状和分离不可知的任何质量分布，我们必须通过将每个部分的贡献相加来计算该点处的重力场的质量分布6

回到问题

当然，对于任意多边形或多面体来说，解析解可能会非常困难，因此使用求和要简单得多，并且算法方法将类似地使用求和。

从算法上来讲，这里最简单的方法实际上不是使用 center of mass，而是使用圆形/球形或正方形/立方体。使用压缩并不是不可能，但是从数学上讲，这种方法存在很大的挑战-最好采用依赖更简单数学的方法。一种这样的方法是定义一个包围质量分布空间范围的网格，然后创建简单（多边形/多面体）的多边形或多面体，并以网格点为顶点。这将创建三种多边形或多面体：


不包含任何质量分布的那些
被质量分布完全填充的那些
被质量分布部分填充的那些


geometric packing

质量中心-方法1

当从参考点到质量分布的距离相对于分布的角度范围较大时，以及在没有质量分布（或任何几个分布）的参考几何包围时，这将很好地工作。

然后，通过求和每个多边形的贡献，可以找到分布的重心R



其中M是分布的总质量，ri是第i个多边形的几何中心的空间矢量，mi是密度乘以包含质量的多边形部分的密度（即，对于完全填充的多边形为1.00，对于完全为空的多边形为0.00多边形）。当您增加采样大小（网格点的数量）时，质心的近似值将接近解析解。一旦有了质心，就很容易计算创建的重力场：您只需在R点处放置一个质量为M的点，并使用 即可。

为了演示，这是使用 equation from above库进行多边形操作的Python中二维方法的实现。

import numpy as np
import matplotlib.pyplot as plt
import shapely.geometry as geom

def centerOfMass(r, density = 1.0, n = 100):
    theta = np.linspace(0, np.pi*2, len(r))
    xy = np.stack([np.cos(theta)*r, np.sin(theta)*r], 1)

    mass_dist = geom.Polygon(xy)
    x, y = mass_dist.exterior.xy

    # Create the grid and populate with polygons
    gx, gy  = np.meshgrid(np.linspace(min(x), max(x), n), np.linspace(min(y), max(y), n))
    polygons = [geom.Polygon([[gx[i,j],    gy[i,j]], 
                              [gx[i,j+1],  gy[i,j+1]], 
                              [gx[i+1,j+1],gy[i+1,j+1]], 
                              [gx[i+1,j],  gy[i+1,j]],
                              [gx[i,j],    gy[i,j]]])
                for i in range(gx.shape[0]-1) for j in range(gx.shape[1]-1)]

    # Calculate center of mass
    R = np.zeros(2)
    M = 0
    for p in polygons:
        m = (p.intersection(mass_dist).area / p.area) * density
        M += m
        R += m * np.array([p.centroid.x, p.centroid.y])

    return geom.Point(R / M), M

density = 1.0     # kg/m^2
G = 6.67408e-11   # m^3/kgs^2
theta = np.linspace(0, np.pi*2, 100)
r = np.cos(theta*2+np.pi)+5+np.sin(theta)+np.cos(theta*3+np.pi/6)

R, M = centerOfMass(r, density)
m = geom.Point(20, 0)
r_1 = m.distance(R)
m_1 = 5.0         # kg
F = G * (m_1 * M) / r_1**2
rhat = np.array([R.x - m.x, R.y - m.y])
rhat /= (rhat[0]**2 + rhat[1]**2)**0.5

# Draw the mass distribution and force vector, etc
plt.figure(figsize=(12, 6))
plt.axis('off')
plt.plot(np.cos(theta)*r, np.sin(theta)*r, color='k', lw=0.5, linestyle='-')
plt.scatter(m.x, m.y, s=20, color='k')
plt.text(m.x, m.y-1, r'$m$', ha='center')
plt.text(1, -1, r'$M$', ha='center')
plt.quiver([m.x], [m.y], [rhat[0]], [rhat[1]], width=0.004, scale=0.25, scale_units='xy')
plt.text(m.x - 5, m.y + 1, r'$F = {:.5e}$'.format(F))
plt.scatter(R.x, R.y, color='k')
plt.text(R.x, R.y+0.5, 'Center of Mass', va='bottom', ha='center')
plt.gca().set_aspect('equal')
plt.show()

shapely

这种方法有点过大：在大多数情况下，找到质心和多边形面积乘以质心和总质量的密度就足够了。但是，它甚至可以用于不均匀的质量分布-这就是为什么我将其用于演示。

现场求和-方法2

在许多情况下，这种方法也是过大的，尤其是与第一种方法相比，但是，它将在任何分布下（经典体制内）提供最佳近似。

这里的想法是总结质量分布的每个块对点质量的影响，以确定净重力（基于可以独立添加重力场的前提）

class pointMass:
    def __init__(self, mass, x, y):
        self.mass = mass
        self.x = x
        self.y = y

density = 1.0     # kg/m^2
G = 6.67408e-11   # m^3/kgs^2

def netForce(r, m1, density = 1.0, n = 100):
    theta = np.linspace(0, np.pi*2, len(r))
    xy = np.stack([np.cos(theta)*r, np.sin(theta)*r], 1)

    # Create a shapely polygon for the mass distribution
    mass_dist = geom.Polygon(xy)
    x, y = mass_dist.exterior.xy

    # Create the grid and populate with polygons
    gx, gy  = np.meshgrid(np.linspace(min(x), max(x), n), np.linspace(min(y), max(y), n))
    polygons = [geom.Polygon([[gx[i,j],    gy[i,j]], 
                              [gx[i,j+1],  gy[i,j+1]], 
                              [gx[i+1,j+1],gy[i+1,j+1]], 
                              [gx[i+1,j],  gy[i+1,j]],
                              [gx[i,j],    gy[i,j]]])
                for i in range(gx.shape[0]-1) for j in range(gx.shape[1]-1)]

    g = np.zeros(2)
    for p in polygons:
        m2 = (p.intersection(mass_dist).area / p.area) * density
        rhat = np.array([p.centroid.x - m1.x, p.centroid.y - m1.y]) 
        rhat /= (rhat[0]**2 + rhat[1]**2)**0.5
        g += m1.mass * m2 / p.centroid.distance(geom.Point(m1.x, m1.y))**2 * rhat
    g *= G

    return g

theta = np.linspace(0, np.pi*2, 100)
r = np.cos(theta*2+np.pi)+5+np.sin(theta)+np.cos(theta*3+np.pi/6)
m = pointMass(5.0, 20.0, 0.0)
g = netForce(r, m)

plt.figure(figsize=(12, 6))
plt.axis('off')
plt.plot(np.cos(theta)*r, np.sin(theta)*r, color='k', lw=0.5, linestyle='-')
plt.scatter(m.x, m.y, s=20, color='k')
plt.text(m.x, m.y-1, r'$m$', ha='center')
plt.text(1, -1, r'$M$', ha='center')
ghat = g / (g[0]**2 + g[1]**2)**0.5
plt.quiver([m.x], [m.y], [ghat[0]], [ghat[1]], width=0.004, scale=0.25, scale_units='xy')
plt.text(m.x - 5, m.y + 1, r'$F = ({:0.3e}, {:0.3e})$'.format(g[0], g[1]))
plt.gca().set_aspect('equal')
plt.show()

对于所使用的相对简单的测试用例，得出的结果与第一种方法非常接近：



但是，尽管在某些情况下第一种方法无法正常工作，但在第二种方法中也没有第二种方法会失败的情况（在传统制度下），因此建议您使用这种方法。



1这确实在极端情况下分解，例如越过黑洞的事件视界或r接近 时，但这些情况不是此问题的主题。

2在密度不均匀的情况下，这变得非常复杂，在无法用符号表示质量分布的情况下，则没有简单的解析解。

3应该注意的是，这实际上是 Planck length所做的；寻找质心。

4对于质量分布内牛顿 integral内的点质量，必须使用场求和。

5在天文学中，这称为 Shell Theorem，并且物体始终绕系统的重心运行-而不是任何给定物体的质心。

6在某些情况下，使用牛顿 barycenter就足够了，但是这些情况与分布几何无关。