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

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的概念。假设有两个点质量，它们在一个孤立系统中，彼此相隔一定距离，r₁，质量分别为m₁和m₂

m₁创建的gravitational field由

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

该场施加在m₂上的重力由下式给出：

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

引力相互作用的场性质允许我们利用superposition计算多重相互作用产生的净引力。考虑到如果添加另一个质量， M EEM>到以前的场景，

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

随着r_i，j=r_j，i。这适用于任何分离处的任何数量的质量。它还意味着，如果您喜欢这种形式主义，那么由质量集合创建的字段可以通过vector sum聚合

现在考虑，如果我们有大量的点质量，MeEEM >，在一个连续的刚性密度的统一体中聚集在一起。然后我们想计算单个空间上不同的点质量，m，由于聚合质量，m，上的引力：

而不是考虑点质量，我们可以考虑不同尺寸的质量的区域（或体积），或者将这些区域（或体积）对积分质量的影响进行积分或求和。在二维情况下，引力的大小是

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

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

然而，为了找到由于任何点的质量分布而对点质量产生的引力，对于形状和分离不可知的任何质量分布，我们必须通过将质量分布各部分的贡献相加来计算该点的重力场。⁶

回到问题上来

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

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

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

重心-方法1

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

然后，通过将每个多边形的贡献相加，可以找到分布的质心，R

其中，M是分布的总质量，r_{i是指向i^{th多边形几何中心的空间向量，Mi是密度乘以包含质量的多边形部分（即，对于完全填充的多边形，1.00，对于完全空多边形，0.00）。随着采样大小的增加（网格点的数量）质心的近似值将接近解析解。一旦你有了质心，计算所产生的引力场就很简单了：你只需在点R处放置一个质量M的点，然后使用equation from above}}

为了演示，下面是使用多边形操作的shapely库在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()

这种方法有点过分：在大多数情况下，只要找到多边形的质心和面积乘以质量中心和总质量的密度就足够了。然而，它甚至适用于非均匀质量分布——这就是为什么我将其用于演示的原因

场总和-方法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接近Planck length时，但这些情况不是这个问题的主题

²在密度不均匀的情况下，这会变得非常复杂，在质量分布无法用符号描述的情况下，没有简单的解析解

³可能应该注意到，这就是integral正在做的事情；找到m的中心屁股

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

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

⁶在某些情况下，使用牛顿的Shell Theorem就足够了，但是这些情况不是分布几何不可知的