I'm currently writing a program that attempts to simulate particle movement in a gravitational field with more than one object exerting a force on it.

I decided that I'd have the particle move by constantly assigning it with new velocity vectors. This is determined as follows:

  1. I loop through each object inside the scene that isn't the particle. I use the equation: A = GM/R^2 to determine the acceleration that the particle will feel with respect to the current object. Then I multiply the distance vector (difference X, difference Y) between the particle and the object by the Acceleration and add that to a list.

  2. I loop through the final list of all Acceleration Vectors, and add the X's together and the Y's together in order to determine the resultant acceleration vector will be on the particle.

  3. I take this resultant acceleration vector, and then set it as the velocity vector of the particle.

While I think I've gotten the resultant acceleration down okay, I keep finding myself unable to figure out how I would account for the particles current velocity.

I could add it as a vector to the list of acceleration vectors acting on the particle, but it doesn't behave properly.

For example, if I attempt to get a particle to orbit an object, I expect to use the equation:

V = sqrt(G*M/r) to get myself the appropriate velocity needed to orbit the object.

However, even with both the acceleration vector and velocity vector sharing the same magnitude, the particle behaves as follows: enter image description here

In Scene 1, R is the resultant vector of both V and A. This will now be set as the particles new velocity vector.

enter image description here

In Scene 2, you can see that R is now again the resultant between V and A, but the problem is that V is now incorrectly represented and thus the resultant will angle even more towards the surface of the object.

I quickly saw that this was a problem, I couldn't be setting velocity as resultant acceleration. But I can't think of how I would properly implement it either. The current velocity would have to be 'angled' 90 degrees to the acceleration vector, but I cannot envision how I could determine this in the infinite amount of other cases where it won't be a perfect 90 degrees.

I'm beginning to feel that I've made some major mistakes in the entire implementation of this model too, and if anyone could suggest how I could properly account for velocity that would be such a help.

Thanks again.

  • $\begingroup$ You should be using the acceleration to determine the change in velocity for the next timestep, not the other way around. $\endgroup$ – The Photon Jan 13 '15 at 18:49
  • $\begingroup$ Also, you shouldn't add the accelerations vectorially. You should find the net vector sum of forces and then calculate acceleration due to that net force. $\endgroup$ – Gaurav Jan 13 '15 at 18:52
  • $\begingroup$ @Gaurav, how will that give a different result? Is the mass of the particle different for each acting force? $\endgroup$ – The Photon Jan 13 '15 at 18:53
  • $\begingroup$ @ThePhoton It is convention to do so. I can show you numerous cases where direct vector addition of acceleration instead of force gives the incorrect result. $\endgroup$ – Gaurav Jan 13 '15 at 18:56
  • $\begingroup$ @Gaurav I am now calculating resultant Force using GmM/R^2. Then A is acquired at the end by dividing the X and Y components of the force vector by the mass of the current planet. I believe this is what you recommended? $\endgroup$ – Micrified Jan 13 '15 at 19:40

You are trying to model a simple 2-body simulation (a special case of $n$-body simulations). Your basic set up should revolve around two equations: \begin{align} \frac{dx}{dt}&=v\\ \frac{dv}{dt}&=\frac{F}{m} \end{align} which is really breaking up Newton's 2nd law into two steps. Here $$ F=G\frac{mM}{r^2} $$ is the gravitational force and $a=F/m$ your acceleration.

In discrete steps, you should be evolving the position and velocity of the particles as \begin{align} x_{new} &= x_{old}+v_{old}\cdot dt\tag{1}\\ v_{new} &= v_{old}+\frac{F}{m}\cdot dt\tag{2} \end{align} where $x$ and $v$ are 2D vectors.

Your most simple algorithm is the Euler method, described in Equations (1) & (2). A more stable integration scheme, called Verlet integration, goes

  1. Compute the new position step: $$x_{n+1}=x_n+v_{n}dt+\frac12a(x_n)dt^2$$
  2. Compute new velocity: $$v_{n+1}=v_{n}+\frac12\left(a(x_{n+1})+a(x_n)\right)dt$$
  3. Set $x_n=x_{n+1}$ and $v_n=v_{n+1}$, increment time: $t=t+dt$ and repeat until $t=t_{end}$

And there are other integration methods available to you, but either of these two are simple enough to implement and should work for your purposes (provided $dt$ is small enough).

  • $\begingroup$ How should I implement this? I currently set particle velocities, and the underlying library takes care of moving it. Should I disable this and move them myself by changing particle positions? $\endgroup$ – Micrified Jan 14 '15 at 0:18
  • $\begingroup$ Well if your library only needs velocities, then give it $v_{n+1}$ from Step 2. You keep $x_{n+1}$ local and use it to compute the accelerations, $a(x_{n+1})$ & $a(x_n)$, needed to compute $x_{n+1}$ & $v_{n+1}$. $\endgroup$ – Kyle Kanos Jan 14 '15 at 2:25
  • $\begingroup$ Okay, well. I've implemented the Euler method for determining my new velocities, and while I've made a step in improving the sim, the problem outlined in the images in my OP are still relevant. Notice how the direction of velocity slants more and more towards the center of the planet. It doesn't maintain the 90 degree angle between it and the centripetal force. I'm going to choose your answer and allow this problem to be marked as solved, so that it's not a headache for you. I may repost this in another question some other time.Here is the current result: imgur.com/a/nluDb $\endgroup$ – Micrified Jan 14 '15 at 4:21
  • $\begingroup$ @Owatch how you describe the way you calculate the acceleration it seems like you do this, $\vec{a}=\frac{GM}{r^2}\vec{r}$, but the magnitude of this would be equal to $\frac{GM}{r}$. To correct this you can use $\vec{a}=\frac{GM}{r^3}\vec{r}$. $\endgroup$ – fibonatic Jan 14 '15 at 13:02
  • $\begingroup$ @Owatch: I've uploaded my own two_body C code to github. Take a look to see how I implement the solver and see what you can use for your code. $\endgroup$ – Kyle Kanos Jan 14 '15 at 16:42

Just to be clear, this part

I take this resultant acceleration vector, and then set it as the velocity vector of the particle.

is where your simulation diverges from reality. It sounds like you know that. I've described how to fix your problem below - other than that, it sounds like you have a fine simulation design.

The solution to your problem can be found in the definition of acceleration.

Acceleration is the rate of change of velocity. In mathematical symbols,

$$ a = \frac{\Delta v}{\Delta t} $$

Where $\Delta v$ is the change in velocity for some small change in time $\Delta t$. We can expand this a little by writing out what $\Delta v$ means:

$$ a = \frac{v_1 - v_0}{\Delta t} $$

Since you want to find the velocity at the end of each time step in your simulation, let's solve for $v_1$:

$$ v_1 = v_0 + a \Delta t $$

This is your answer, but let's interpret it carefully. $v_1$ is the velocity at the end of each of your simulation's time steps. $v_0$ is the velocity at the end of the last time step, and at the beginning of the current time step. $\Delta t$ is the duration of your simulation's time step - this is something you have to set for your simulation - how long each discrete time step will be.

So, what this equation tells you is that to "update" the velocity (get the current velocity $v_1$), you take the previous velocity of your particle $v_0$ and add the amount $a \Delta t$.

Now, it sounds like you're doing this in two dimensions, so these are all really vectors:

$$ \vec{v}_1 = \vec{v_0} + \vec{a} \Delta t $$

which means we have two equations, one for each dimension:

$$ v_{1,x} = v_{0,x} + a_x \Delta t $$


$$ v_{1,y} = v_{0,y} + a_y \Delta t $$

  • 1
    $\begingroup$ OP should be aware that this is essentially Euler's method of integration (you are integrating velocity to get position) and it careful attention to the timestep size is needed to get accurate results. $\endgroup$ – The Photon Jan 13 '15 at 18:57
  • $\begingroup$ True enough - every simulation method diverges from reality in its own special way. $\endgroup$ – Brionius Jan 13 '15 at 19:01
  • $\begingroup$ @ThePhoton I am using a constant for a timestep. Is the size of the time step what determines it's accuracy? I have a pretty small timestep. $\endgroup$ – Micrified Jan 13 '15 at 19:43
  • $\begingroup$ Yes, smaller timestep gives better accuracy. If you use more complex integration techniques you can use a larger timestep and get the same accuracy (and also probably an estimate of the error). If Euler isn't good enough, the typical next step is called "4th-order Runge-Kutta". If that isn't good enough, there's tons of literature out there but it's beyond my knowledge. $\endgroup$ – The Photon Jan 13 '15 at 19:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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