2
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ You could use SymPy to compute the gradient of the potential energy and ask it the Python code. You could probably use a Boundary Integral method to find the gravitational field. You could also replace your cube for set of particles and superimpose their gravitational field, when the number of particles tends to infinite you should have the right value. $\endgroup$
    – nicoguaro ♦
    Commented Apr 2, 2017 at 5:24
  • 1
    $\begingroup$ I rewrote my comment as an answer $\endgroup$
    – nicoguaro ♦
    Commented Apr 2, 2017 at 14:38
  • $\begingroup$ @nicoguaro indeed you did! :) $\endgroup$
    – uhoh
    Commented Apr 2, 2017 at 14:39

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Those are all excellent suggestions! SymPy has some overhead - I need to learn to use it, but if it generates python and that's very helpful to me. The mesh interpolation is my favorite at the moment because it is so straightforward, and I've done a lot of charged particle beam orbit work in the past. I will look at the boundary integral as a good exercise, related to what I was doing here. Particles sounds fun and also a good exercise if I use a non-uniform density/weighting. This is precisely the kind of answer I needed. Thanks! $\endgroup$
    – uhoh
    Commented Apr 2, 2017 at 14:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.