A vector field approach to calculating gravitational forces

Computing inverse square forces is one of the most fundamental calculations of classical physics. However, despite centuries of work on special cases [1–5] and multipole [6] or Green's function expansions [7] valid only under certain conditions. In the cases where these tricks are not applicable practical calculations are very computationally intensive to perform. The computational complexity often leads to the need to treat extended bodies as point masses after performing tedious case-by-case calculations to show that this assumption has negligible effect on the final result [8, 9]. In this work, taking inspiration from the field of computer graphics [10], we derive a method to reduce the calculation of the force between a point mass and any arbitrary shape to a sum of one-dimensional integrals rather than a three-dimensional integral. This is particularly beneficial as fast, accurate, multi-dimensional numerical integration is not trivial for many shapes and is still an active topic for computer scientists [11, 12]. The solution provided has been formulated in terms of Newtonian gravity but could easily be applied to other fields such as electrostatics.

The force in the x-direction (F_x) from an arbitrary object of constant density, ρ defined by the volume V, on a point mass of mass m (at the origin of our coordinate system) is

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{x}}{{Gm}\rho }={\iiint }_{V}x{({x}^{2}+{y}^{2}+{z}^{2})}^{-3/2}\,{\rm{d}}V.\end{eqnarray} \tag{ 1 }$$

However,

$$\begin{eqnarray}&&x{({x}^{2}+{y}^{2}+{z}^{2})}^{-3/2}={\rm{\nabla }}\cdot \displaystyle \frac{-1}{\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}}\hat{{\boldsymbol{\imath }}},\end{eqnarray} \tag{ 2 }$$

therefore by applying divergence theorem we can reduce the order of the integration to a closed surface integral:

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{x}}{{Gm}\rho }={\large{\unicode{8751}}}_{S}\displaystyle \frac{-1}{\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}}\hat{{\boldsymbol{\imath }}}\cdot \hat{{\boldsymbol{n}}}\,{\rm{d}}S,\end{eqnarray} \tag{ 3 }$$

where S is the bounding surface of the volume V with surface normal $\hat{{\boldsymbol{n}}}$. The surface can then be broken up into n surface elements

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{x}}{{Gm}\rho }=\displaystyle \sum _{n}{\iint }_{{S}_{n}}\displaystyle \frac{-1}{\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}}\hat{{\boldsymbol{\imath }}}\cdot {\hat{{\boldsymbol{n}}}}_{n}\,{\rm{d}}{S}_{n},\end{eqnarray} \tag{ 4 }$$

For efficient calculation we can consider that as the integral includes $\hat{{\boldsymbol{\imath }}}\cdot \hat{{\boldsymbol{n}}}\,{\rm{d}}S$ we can project each surface into the y–z plane. Any surface which obscures part of the same surface when projected must be broken up into multiple surfaces, or the inverse of the projection is not a single valued function. Also each surface with with no projection y–z plane can be ignored. As such we can write the integral over the n projected surfaces (S_n^yz) as

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{x}}{{Gm}\rho }=\displaystyle \sum _{n}{\iint }_{{S}_{n}^{{yz}}}\displaystyle \frac{-1}{\sqrt{{({X}_{n}(y,z))}^{2}+{y}^{2}+{z}^{2}}}\hat{{\boldsymbol{\imath }}}\cdot \hat{{\boldsymbol{\imath }}}\,{\rm{d}}{S}_{n}^{{yz}},\end{eqnarray} \tag{ 5 }$$

where the x-coordinate of the surface S_n^yz can be written as $x={X}_{n}(y,z)$. The full process described above is shown graphically in figure 1.

Zoom In Zoom Out Reset image size

Any arbitrary object's surface can be broken down, either exactly or to a very high precision with triangular meshing, into one of three simple cases. For forces in the x-direction these are: a flat plane of any arbitrary orientation; a surface (flat or singly curved) whose normal meets the requirements ${\hat{{\boldsymbol{n}}}}_{n}\cdot \hat{{\boldsymbol{\jmath}}}=0$ or ${\hat{{\boldsymbol{n}}}}_{n}\cdot \hat{{\boldsymbol{k}}}=0;$ or a surface (flat or singly curved) whose normal meets the requirement ${\hat{{\boldsymbol{n}}}}_{n}\cdot \hat{{\boldsymbol{\imath }}}=0$. Clearly in the third case equation (5) evaluates to zero and the other cases reduce to single integrals as solved below. For doubly curved surfaces we cannot reduce the surface integral to a single integral without triangular meshing.

2.1. Any flat surface

For a flat surface we can write $X(y,z)={ay}+{bz}+c$, and then we note that as

$$\begin{eqnarray}&&\displaystyle \frac{-1}{\sqrt{{({ay}+{bz}+c)}^{2}+{y}^{2}+{z}^{2}}}\hat{{\boldsymbol{\imath }}}\cdot \hat{{\boldsymbol{\imath }}}=-\left({\rm{\nabla }}\times \displaystyle \frac{1}{\sqrt{{b}^{2}+1}}\mathrm{log}(f(y,z)\sqrt{{b}^{2}+1}+{aby}+{b}^{2}z+{bc}+z)\hat{{\boldsymbol{\jmath}}}\right)\cdot \hat{{\boldsymbol{\imath }}}\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&f(y,z)=\sqrt{({a}^{2}+1){y}^{2}+({b}^{2}+1){z}^{2}+2c({ay}+{bz})+2{abyz}+{c}^{2}}.\end{eqnarray} \tag{ 7 }$$

As such, the x-direction force between a point mass and any object composed of n flat faces is simply

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{x}}{{Gm}\rho }=-\displaystyle \sum _{n}{\oint }_{{c}_{n}}\displaystyle \frac{1}{\sqrt{{b}^{2}+1}}\mathrm{log}(f(y,z)\sqrt{{b}^{2}+1}+{aby}+{b}^{2}z+{bc}+z)\hat{{\boldsymbol{\jmath}}}\cdot \,{\rm{d}}{{\boldsymbol{c}}}_{n},\end{eqnarray} \tag{ 8 }$$

where c_n is the contour enclosing the surface S_n^yz.

2.2. A surface with normal perpendicular to the y or z-axis

In the case of surface whose normal is always perpendicular to the y or z-axis we can write the function $X(y,z)$ as a function of either y or z. Considering the case where X is a function of y

$$\begin{eqnarray}&&\displaystyle \frac{-1}{\sqrt{{(X(y))}^{2}+{y}^{2}+{z}^{2}}}\hat{{\boldsymbol{\imath }}}\cdot \hat{{\boldsymbol{\imath }}}=-\left({\rm{\nabla }}\times \mathrm{arsinh}\left(\displaystyle \frac{z}{\sqrt{{(X(y))}^{2}+{y}^{2}}}\right)\mathop{\,\,\,\hat{{\boldsymbol{\jmath}}}}\limits^{}\right)\cdot \hat{{\boldsymbol{\imath }}}.\end{eqnarray} \tag{ 9 }$$

Then by Stokes theorem we can write

$$\begin{eqnarray}&&{\iint }_{{S}_{n}^{{yz}}}\displaystyle \frac{-1}{\sqrt{{({X}_{n}(y))}^{2}+{y}^{2}+{z}^{2}}}\hat{{\boldsymbol{\imath }}}\cdot \hat{{\boldsymbol{\imath }}}\,{\rm{d}}{S}_{n}^{{yz}}=-{\oint }_{{c}_{n}}\mathrm{arsinh}\left(\displaystyle \frac{z}{\sqrt{{(X(y))}^{2}+{y}^{2}}}\right)\mathop{\,\hat{{\boldsymbol{\jmath}}}}\limits^{\,\,\,}\cdot \,{\rm{d}}{{\boldsymbol{c}}}_{n}.\end{eqnarray} \tag{ 10 }$$

For faces where X is a function of z a similar transformation can easily be calculated.

In the case of a right prism with an axis parallel to the z-axis the top and bottom faces integrate to zero and the force is simply calculated as

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{x}}{{Gm}\rho }=-\displaystyle \sum _{n}{\oint }_{{c}_{n}}\mathrm{arsinh}\left(\displaystyle \frac{z}{\sqrt{{(X(y))}^{2}+{y}^{2}}}\right)\mathop{\,\hat{{\boldsymbol{\jmath}}}}\limits^{\,\,\,}\cdot \,{\rm{d}}{{\boldsymbol{c}}}_{n},\end{eqnarray} \tag{ 11 }$$

clearly a similar result exists for right prisms with an axis parallel to the y-axis. A right prism with an axis parallel to the x-axis requires only integrals to be solved for the top and bottom faces, and X is simply a constant

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{x}}{{Gm}\rho }=-\displaystyle \sum _{n=1,2}{\oint }_{{c}_{n}}\mathrm{arsinh}\left(\displaystyle \frac{z}{\sqrt{{X}^{2}+{y}^{2}}}\right)\mathop{\,\hat{{\boldsymbol{\jmath}}}}\limits^{\,\,\,}\cdot \,{\rm{d}}{{\boldsymbol{c}}}_{n},\end{eqnarray} \tag{ 12 }$$

where c₁ and c₂ are contours around the top and bottom faces respectively. Right prisms of any other orientation can be solved with a coordinate transformation.

The method derived above provides a simple framework for calculation of the gravitational forces between an arbitrary object and a point mass as a one-dimensional integral. There are multiple benefits to the reduction of a multi-dimensional integral to a one-dimensional integral. N-dimensional integrals scale as M^N where M is the number of function evaluations for the integral to converge to the desired precision. Also the limits for N-dimensional integrals, except the most primitive shapes, are functions of $N-1$ variables, thus to calculate such an integral one must invest significant time programming the limits on a case-by-case basis. One can also solve the triple integral by Monte-Carlo integration, but as Monte-Carlo converges asymptotically it is not suitable for high precision calculations [13]. For these reason reduction of N-dimensional integrals to a lower dimensionality is the suggested method where possible [13]. We note that a similar methodology has been applied to the calculation of gravitational potentials of a general polyhedra [14].

For the example of a polyhedron, such as in figure 1, the method presented is simple to compute. From three points defining each face of the polyhedron the function $X(y,z)={ay}+{bz}+c$ is simply calculated as a matrix inversion

$$\begin{eqnarray}&&\left(\begin{array}{c}a\\ b\\ c\end{array}\right)={\left(\begin{array}{ccc}{y}_{1} & {z}_{1} & 1\\ {y}_{2} & {z}_{2} & 1\\ {y}_{3} & {z}_{3} & 1\end{array}\right)}^{-1}\left(\begin{array}{c}{x}_{1}\\ {x}_{2}\\ {x}_{3}\end{array}\right),\end{eqnarray} \tag{ 13 }$$

where $({x}_{i},{y}_{i},{z}_{i})$ are the coordinates of point i. If the matrix is singular the no integrals need to be computed. If the matrix is non-singular then equation (8) can be evaluated with Gaussian quadrature.

To demonstrate calculations for singly curved surfaces we note that the the force from a cylinder can be written down almost by inspection. If its axis parallel to the z-axis and its centre of mass at $(a,0,0)$ equation (11) becomes

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{x}}{{Gm}\rho }=2{\int }_{-R}^{R}\left(\mathrm{arsinh}\left(\displaystyle \frac{H}{2\sqrt{{(a-\sqrt{{R}^{2}-{y}^{2}})}^{2}+{y}^{2}}}\right)-\mathrm{arsinh}\left(\displaystyle \frac{H}{2\sqrt{{(a+\sqrt{{R}^{2}-{y}^{2}})}^{2}+{y}^{2}}}\right)\right)\,{\rm{d}}y,\end{eqnarray} \tag{ 14 }$$

where H and R are the height and radius of the cylinder respectively. We note that this can be shown to evaluate to the same value as calculations performed via more complex procedures such as the method of Chen and Cook [2]. We do note that for the particular case of a complete cylinder Chen and Cook's method is computationally more efficient by an order of magnitude as it reduces the problems to elliptical integrals, however it cannot be adapted to more complex shapes with singly curved surfaces.

The methods presented have been compared in both compiled (C++ using the GNU Scientific Library) and interpreted (GNU Octave) languages. We note that as Octave uses compiled libraries for Gaussian quadrature and matrix inversion the computation time is not significantly different. Computation of simple geometries converges to a relative uncertainty of less than one part in 10¹² in ∼1 ms, two orders of magnitude faster than for 3D integration. To calculate the force between two extended bodies this result itself must be triple integrated over a second body, thus computation of extended bodies can often be reduced from minutes to seconds per calculation. For more complex polyhedra which are non-trivial to compute with triple integrals our method scales linearly with the number of faces, the results are in agreement with Monte-Carlo integration to within the statistical noise of repeated Monte-Carlo integrations. Even given tens of seconds Monte-Carlo integration does not produce results with errors below a part in 10⁵.

Core components of many gravitational experiments are primitive shapes allowing for fast analytical calculations of forces. The forces from many other complex experimental components cannot be solved analytically and thus require slower numerical integration techniques. Yet these components can easily be exported as triangular meshes from most computer aided design or finite element analysis software. We provide a method to calculate gravitational forces given this mesh using one-dimensional rather than three-dimensional integration. The result is simple to code and computationally efficient. While the method provided in this paper is for forces in the x-direction, for a point mass at the origin, any other direction or mass position can be calculated using a simple coordinate transformation. We are, as such, confident that this method will be of use for many practical calculations.

The author would like to thank Stephan Schlamminger, Clive Speake, and Rob Smith for fruitful discussions.