如何计算任意二维多边形的重力？

Question

如何计算任意二维多边形的重力？

4

（为简单起见，以下内容仅讨论二维情况）我知道两个球体之间由于引力产生的力是 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

这个方案可行吗？如果是，如何高效快速地打包圆形？如果不可行，是否有更好的方案？

编辑：请注意，我想找到物体在特定点的受力情况，因此需要计算圆形和观察者之间的角度（并将向量相加）。这与寻找受到的总力不同。

- Max

1

这个公式难道不适用于每个人吗？ - Clément

6

似乎更像是物理问题而不是编程问题。相关链接。 - Ahmet

2

本质上，它是两个物体体积之间的积分（每个对象的“点”吸引另一个对象的“点”）。 - Willem Van Onsem

3

质心和 G(m1*m2/r**2) 只适用于足够远的物体......对于接近的物体，您需要积分（面积或体积）......或者将您的物体转换为一组质点（精度降低但计算速度更快且更简单）。 - Spektre

1

这是否是可行的解决方案？是的，这本质上是使用叠加原理来找到两个物体之间的净重力相互作用。然而，这是极度浪费的——对于任何形状（2D或3D）的均匀密度物体，你只需要找到它们的质心即可。引力作用将始终沿着从一个质心到另一个质心的射线方向作用。这对于任何数量的物体系统中的任何两个物体都成立，并且它们的相互作用是相互独立的。 - William Miller

显示剩余20条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- William Miller · Accepted Answer

背景

这个解释有些偏题，但我认为它有必要帮助澄清评论中提出的一些问题，并且因为其中许多内容有些反直觉。

这种引力相互作用的解释依赖于质点的概念。假设你有两个质点，在一个孤立系统中彼此分开一定距离r₁，分别具有质量m₁和m₂。

引力场是由m₁产生的，其表达式为

其中G是普适引力常数，r是从m₁到m₂的距离，r̂是沿着它们之间连线的单位方向。

这个场对m₂施加的引力为：

注意 - 重要的是，这对于任何两个质点在任何距离上都是成立的。¹

引力相互作用的场性质使我们能够应用叠加原理来计算多重相互作用导致的净引力。考虑如果我们向之前的情境中添加另一个质量为m₃的物体，

那么对于质量为m₂的物体，其所受到的引力是由每个其他物体产生的引力场相加而成的。

对于任意数量和分离程度的物体，都有r_i,j = r_j,i。这也意味着由一组物体创建的场可以通过矢量求和进行聚合，如果您更喜欢这种形式。

现在考虑如果我们有一个非常大的点质量数M，聚集在一个连续的、均匀密度的刚体中。然后我们想要计算单个空间不同点质量m由于聚合质量M而产生的引力。

那么，我们可以考虑质量差异大小的面积（或体积），而不是考虑点质量，并且将这些区域（或体积）对点质量的作用进行积分或求和。在二维情况下，引力的大小为

其中σ是聚合物质的密度。²这相当于对每个微分质量的重力向量场进行求和，即σdxdy。这种等价性非常重要，因为它意味着对于任何足够远离质量分布的点质量来说，由于质量分布产生的引力力量几乎与位于质量分布center of mass处的质量M的点质量相同。³⁴

这意味着，在计算任何质量分布所产生的引力场时，可以用质量分布在重心处的等效质点代替质量分布，从而得到非常好的近似结果。这适用于任何数量的空间不同的质量分布，无论这些分布是否构成刚体。此外，这意味着您甚至可以将分布组聚集到系统重心处的单个质点中。⁵只要参考点足够远。

然而，为了找到质量分布在任意点上对点质量的引力，无论质量分布的形状和分离方式如何，我们都必须通过将来自质量分布每个部分的贡献相加来计算该点处的引力场。⁶

回到问题

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

从算法角度来看，这里最简单的方法实际上不是几何包装（使用圆/球或正方形/立方体）。使用包装并非不可能，但从数学上讲，这种方法存在重大挑战-最好采用依赖于更简单数学的方法。其中一种方法是定义一个覆盖质量分布空间范围的网格，然后创建以网格点为顶点的简单（正方形/立方体或矩形/长方体）多边形或多面体。这将创建三种类型的多边形或多面体：

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

质心 - 方法1

当参考点到质量分布的距离相对于分布的角度范围很大，并且没有任何几个分布将参考点围住时，这种方法会很有效。

然后，您可以通过对每个多边形的贡献求和来找到分布的质心R。

其中M是分布的总质量，r_i是到第i个多边形几何中心的空间向量，m_i是密度乘以包含质量部分的多边形（即对于完全填充的多边形为1.00，对于完全空的多边形为0.00）。随着采样大小（网格点数）的增加，重心的近似值将接近解析解。一旦你有了重心，就可以轻松地计算所创建的引力场：只需在点R处放置质量为M的点质量，并使用上面的方程。

为了演示，在Python中使用shapely库进行多边形操作的二维实现如下：

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

这种方法有些过于复杂：在大多数情况下，只需找到多边形的质心和面积乘以密度即可得到重心和总质量。然而，它也适用于非均匀分布的质量 - 这就是为什么我用它进行演示的原因。

场总和 - 方法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()

相对简单的测试案例表明，这种方法得到的结果非常接近第一种方法：

但是，虽然第一种方法无法正确处理某些情况，但在传统的情况下，第二种方法没有失败的案例，因此建议优先考虑这种方法。

¹在极端情况下，例如黑洞事件视界之外或当r趋近于普朗克长度时，它会发生瓦解，但这些情况不是本问题的主题。

²当密度不均匀且无法用符号描述质量分布时，情况变得更加复杂，没有显然的解析解。

³应该注意，这实际上就是积分正在做的事情；寻找质心。

⁴对于一个质点位于质量分布内部时，必须使用牛顿的壳定理或场的总和。

在天文学中，这被称为重心，天体总是绕系统的重心而不是任何给定天体的质心运动。

⁶在某些情况下，使用牛顿的外壳定理就足够了，但这些情况并不是分布几何无关的。