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

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

背景

这种解释有些离题，但我认为有必要帮助澄清 cmets 中提出的一些内容，因为其中大部分内容有些违反直觉。

这种对引力相互作用的解释取决于point masses 的概念。假设你有两个点质量，它们位于一个孤立的系统中，彼此相隔一定距离，r₁，质量为 m_{1 sub>} 和 m₂ 分别，

m₁创建的gravitational field由

给出

其中G 是universal gravitational constant，r 是距离m₁ 和 r̂是沿m₁和m₂之间连线的单位方向。

这个场施加在m₂上的引力由下式给出

注意   -   重要的是，这适用于 any 距离处的 any 两点质量。¹

引力相互作用的场特性允许我们使用superposition 来计算由于多重相互作用产生的净引力。考虑如果我们在前面的场景中添加另一个质量，m₃，

那么对质量m₂的引力就是每个其他质量所产生的场的引力之和，

与 r_i,j = r_j,i。这适用于任何分离的任何数量的质量。这也意味着，如果您更喜欢这种形式，可以通过 vector sum 聚合由群众集合创建的字段。

现在考虑如果我们有大量的点质量，M，聚集在一起形成一个密度均匀的连续刚体。然后，我们想计算由于总质量 M 而导致的单个空间上不同的点质量 m 上的引力：

然后，我们可以考虑不同大小的质量的面积（或体积），而不是考虑点质量，然后对这些面积（或体积）对点质量的影响进行积分或求和。那么在二维情况下，引力的大小为

其中σ是聚合质量的密度。²这相当于对每个微分质量引起的引力矢量场求和，σdxdy。这种等效性至关重要，因为它意味着对于任何远离质量分布超出的点质量，由于质量分布而产生的引力几乎完全与其相同将是质量为 M 的点质量，位于质量分布的 center of mass。³⁴

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

然而，为了找到由于任何点的质量分布而导致的点质量上的引力，对于任何形状和分离不可知方式的质量分布，我们必须计算引力场此时，通过对质量分布的每个部分的贡献求和。⁶

回到问题

当然，对于任意多边形或多面体，解析解可能非常困难，因此使用求和要简单得多，算法方法也同样使用求和。

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

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

质心 - 方法 1

当从参考点到质量分布的距离相对于分布的角度范围很大时，并且当质量分布（或任何几个分布）没有参考的几何封闭时，这将很有效.

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

其中M是分布的总质量，r_i是几何空间向量i^th 多边形的中心，m_i 是密度乘以包含质量的多边形部分（即 1.00 表示完全填充的多边形，0.00 表示完全空的多边形）。当您增加采样大小（网格点的数量）时，质心的近似值将接近解析解。一旦你有了质心，计算所创建的引力场就很简单了：你只需将一个质量为 M 的点质量放在 R 点并使用equation from above .

为了演示，这里是在 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

在许多情况下，这种方法也过于矫枉过正，尤其是与第一种方法相比，但它会在 any 分布（在经典机制内）下提供最佳近似值。

这里的思路是将每一块质量分布对一个点质量的影响求和，以确定净引力（基于引力场可以独立添加的前提）

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 接近Planck length，但这些情况不是这个问题的主题。

²这在密度不均匀的情况下变得明显更复杂，而在质量分布不能用符号来描述的情况下没有平凡的解析解。

³可能应该指出，这实际上是integral 正在做的事情；找到质心。

⁴对于点质量在质量分布牛顿Shell Theorem，或必须使用场求和。

⁵在天文学中，这被称为barycenter，并且天体总是围绕系统的重心运行，而不是任何给定天体的质心。 p>

⁶在某些情况下，使用 Newton 的 Shell Theorem 就足够了，但这些情况与分布几何无关。