Another way to evaluate the gravitational force from a uniform cube?

Question

Appendix A of Liu, Baoyin, and Ma (2011) Equilibria, periodic orbits around equilibria, and heteroclinic connections in the gravity field of a rotating homogeneous cube shows an analytic expression for the gravitational potential of a uniform cube. I'd like to reproduce the orbit calculation using python, which will require evaluating the gradient of the reduced potential, i.e. acceleration.

Given sufficient coffee I can probably type that into Wolfram alpha, evaluate the gradient, and script as python. Alternately I could use a numerical gradient which would require four potential evaluations and be faster and easier for me to reliably script. I just want to make plots of single orbits, so I do not require very high accuracy.

Is there another way to get a fairly good approximation to the gravitational force from a uniform cube of side $a$ with similar (or less) effort? Fairly good might be say 1E-06 error at distances > $0.1a$ from a face, possibly need to stay farther from the corner for similar error, but it seems that most of the stable orbits tend to do that anyway.

3D direct integration at each time step is extremely slow with scipy's 'triple quadratic' method, not that I'd ever admit trying it. There may be an amazingly clever integration algorithm, and it could compete in speed considering I will be writing python without any of its numerical acceleration options at this point in time.

Answer 1

I am rewriting my comment as an answer.

I think that you have several options. Some of them are:

• You can use SymPy to find the gravitational field from the potential. Then you can generate your Python code from it.

• You could create some kind of mesh and then compute the gradient numerically using finite differences, piecewise polynomials (FEM-like), or Chebyshev polynomials (see pychebfun).

• You can use a Boundary Integral Representation for the gravitational field of your problem.

• You can represent your cube as a set of $n$ discrete particles and consider that the overall field is just the superposition of these "little" contributions. When $n\rightarrow\infty$, the field should approach the real distribution.