Can the gravitational force of a many-body system be calculated via matrix operations?

Question

I'm working on a software project at work that calculates the "gravitational attraction" between points in 1, 2, or 3 dimensions. This is an $O(n^2)$ runtime efficiency problem, however, if I can write it as a sequence of matrix operations, then I might be able to leverage the GPU to run my simulation which would be substantially faster. My linear algebra chops aren't that great however. Is it possible to represent the mass and position as a matrix and compute the net gravitational force for each element in the matrix as a vector?

Edit: The problem I'm trying to solve is this: Given an array of positioned-masses, compute the net Gravitational force that each mass experiences from the other elements in the array. I can easily do this using for-loops, but if it can be done using matrices then I can leverage the GPU for faster calculation.

_{Non-OP edit: given masses $m_{1..n}$, positions $\underline{r}_{1..n}$. Calculate the matrix $F$, for which $F_{ij}=m_im_j\frac{\underline{r}_i-\underline{r}_j}{|\underline{r}_i-\underline{r}_j|^3}$. The goal is to do this with "simple" vector/matrix operations.}

Answer 1

The size of the matrix you have to fill with values is $O(n^2)$. Thus, you don't have an option for an algoritmical optimization. If you calculate optimally, you can't go below $O(n^2)$.

But you have various options for linear optimizations, namely:

• $F_{ij} = -F_{ji}$, on Newton's Laws (but also the formulas show the same). Of course, $F_{ii}=0$.
• The most complex part is to calculate $\frac{\underline{r_i}-\underline{r_j}}{|\underline{r_i}-\underline{r_j}|^3}$, compute it only once and then divide with $m_i$ or $m_j$ as you want. Note in a 2D case where gravity decreases with the recipe of the distance and not with its square, this formula will be far simpler. Although it won't be our Universe any more, it might be useful for graph visualisation tasks.
• The denominator, $|\underline{r_i}-\underline{r_j}|^3$ will be coordinate-independent, you have to calculate it only once for all the cordinates. Only the nominator, $\underline{r_i}-\underline{r_j}$ is dependent on the actually used coordinate.

---

If you are hardline in math and like the half page long formulas, you can also make a step further, derive the formulas yet another time, to calculate $\frac{da}{dt}$. This results far more complex formulas, but your simulation will be far more resistant to the time quantization error. It might be particularly useful if the precision is important.

---

If precision is not so important for you, you can also use the Barnes-Hutt approximation. This makes your problem algorithmically simplified to $O(n\log n)$. Essentially, you create a tree of the stars by recursively halving the space where they are.

This is what also the force simulation of the d3.js does:

https://www.youtube.com/watch?v=kASOnVvjeyo

Note, here you will have to count with a lot of anomalies, the worst is that stars switching tree branch result a sudden discontinuity in the force directions.