N Gravitational Body Simulation

I'm trying to compute an N-Body simulation in Matlab, where I have randomly generated 100 points in 2 dimensions and directly computed the acceleration formula. I believe by doing this I have done it in specifically one time frame, and I want to extend it to multiple time frames but I'm not really sure how...? I have posted the code below and any help of where to go next would appreciated. (Note, in this case I have assumed that the mass M is constant (that of the milky way). The end result it to key in the 127 nearest galaxies to the milky way to see how the galaxies interact but I wanted to see how to code a simple case first)

My code doesn't seem to follow an ODE method like I have seen elsewhere so if my code is completely wrong I'd also be grateful to know (I'm not very experienced in Matlab)!! Thanks.

clear all
N=100; %Number of galaxies
rangex = 100;
rangey = 100;
T=100; %Number of timesteps
G = 6.67408*(10^-11); %Gravitational constant in m^3 kg^-1 s^-2
M = 1.15*1.98855*(10^42); %Mass of MW in kg

for t=1:T
for i=1:N
x(i) = (rand-0.5)*rangex;
y(i) = (rand-0.5)*rangex;
P{i}=[x(i),y(i)];
end

for i=1:N

xs(i) = 0;

for j=1:i-1
dx(i) = P{j}(1,1)-P{i}(1,1);
xs(i)=xs(i)+G*M*((P{j}(1,1)-P{i}(1,1))/(sqrt((P{j}(1,1)-P{i}(1,1)).^2)).^3);
end

for j=i+1:N
dx(i) = P{j}(1,1)-P{i}(1,1);
xs(i)=xs(i)+G*M*((P{j}(1,1)-P{i}(1,1))/(sqrt((P{j}(1,1)-P{i}(1,1)).^2)).^3);
end

end
for i=1:N

ys(i) = 0;

for j=1:i-1
ys(i)=ys(i)+G*M*((P{j}(1,2)-P{i}(1,2))/(sqrt((P{j}(1,2)-P{i}(1,2)).^2)).^3);
end

for j=i+1:N
ys(i)=ys(i)+G*M*((P{j}(1,2)-P{i}(1,2))/(sqrt((P{j}(1,2)-P{i}(1,2)).^2)).^3);
end
end

end

quiver(x,y,xs,ys)
• The code seems harder to read than it should be; are you squaring something, then taking the square root of the result, and then cubing that result? Maybe I lost track of parentheses and read it wrong, or maybe all those parentheses led you to make a mistake. You could probably use vector operations more effectively, especially if you put all your data in a matrix instead of $n$ separate vectors. But those are questions for a MATLAB board. For the purpose of asking a question here, I'd recommend writing out the math formulas you think you want to use. Feb 26 '16 at 20:48
• From trying to reverse-engineer your code, it looks like you're using a simple Euler method, which in my (limited) experience works poorly for gravitational simulations; and besides that, I suspect the equation is wrong. (You can't ignore $y_i - y_j$ while computing the $x$-accelerations, for example.) Feb 26 '16 at 20:51
• Do you really want to re-randomize x and y for each new value of t? Feb 26 '16 at 20:53
• Thanks for your response David. All I was given was the following formula: $-G\sum _{i=1_\neq j}^{N} m_j(\vec x_j - \vec x_i)/(abs(\vec x_j - \vec x_i))^3$ and told to code it (very ambiguous I know!!). Upon researching I am using the Euler method as it is I feel the easiest to use when I haven't been given any direction!! The for t=1:T bit shouldn't be in there yet I accidentally put that in there! Feb 26 '16 at 21:29
• I feel like I'm heading in the complete wrong direction as there is no differential function involved in my code. Any help would be greatly appreciated!!!! Feb 26 '16 at 21:32

Consider the formula you were given, which (based on a comment) was equivalent to this: $$−G \sum_{\substack{1\leq j\leq N \\[.2ex] j\neq i}} m_j \frac{(\vec x_j − \vec x_i)}{(\lVert \vec x_j − \vec x_i\rVert)^3 }$$

Where you had $\mathrm{abs}(\cdots)$, I've written $\lVert \cdots \rVert$ because that's a typical way of denoting the magnitude of a vector. I also noticed that the iteration should be over $j$ for a constant $i$ (so that $m_j$ has to be inside the summation, unlike the constant factor of $G$, which can be moved outside); but I think you got that part right in your implementation.

This formula is essentially one side of the differential equation of motion of the galaxies in your simulation. That's usually all you have in a simulation like this. But the little arrows on top of variables such as $\vec x_i$ means they are vectors, that is, $\vec x_i$ represents both coordinates of one of the galaxies, not just the "$x$" coordinate.

The numerator inside the summation is a simple coordinate-wise sum, so in principle you could compute it for all the coordinates on one axis first, and then for all the coordinates on the other axis. But the denominator cannot be done that way. Suppose you had two galaxies whose coordinates differed by a billion light-years in one dimension, but by only one million light-years in the other dimension. The total force of one galaxy on the other is determined by the billion-light-year distance, but when you look at the second coordinate, you're only using one million light-years as the separation, so you're computing a force that is at least a million ($1000^2$) times as powerful as it should be.

Rather than writing the denominator in the form $(\lVert \vec x_j − \vec x_i\rVert)^3$, try $$\left(\vec r^T \; \vec r\right) ^{1.5},$$ where $\vec r = \vec x_j − \vec x_i.$ The reason this works is that the inner product (aka dot product) of a vector with itself is the square of the vector's magnitude; and as a matrix operation, the inner product of a column vector is written as the product of the transpose of the vector on the left and the vector itself on the right. That is, $\vec r^T \; \vec r = (\lVert \vec x_j − \vec x_i\rVert)^2$. Rather than take the square root of this and then cube it, we can equivalently just take the $3/2$ power.

In vector notation, then, the formula is $$−G \sum_{\substack{1\leq j\leq N \\[.2ex] j\neq i}} m_j \frac{\vec r}{\left(\vec r^T \; \vec r\right) ^{1.5}}, \quad \text{where } \vec r = \vec x_j − \vec x_i.$$

MATLAB gives really good support for matrix and vector operations, so you might as well use them rather than trying to pick apart the components of the vectors.

As for starting with a simple example, a really good place to start would be with just two galaxies. Put them at some distance from each other and give each one a starting velocity so that they will both orbit in the same direction around their common center of mass. The two-body problem is so thoroughly solved that you should be able to find the orbital speed that will make the two orbits circular. So set up those initial conditions, let the simulation run for a few orbits, and see if the bodies maintain circular orbits of constant radius, as they should. (Personally, it's been decades since I tried to do this, but I never had much luck with the obvious Euler-method simulation, which is why I suggest to start with a scenario with a known solution like this before tackling the harder stuff.)