I have created a simple 1 Dimensional two body particle model in C++. In the model, particle 1 starts at position (0,0) and particle 2 starts at position (1,0). The particles are accelerated towards each other at a rate proportional to the distance between them.

The code looks OK to me, but the results look very bazaar and completely nonphysical. I suspect the problem may be related to the discrete time steps taken. I have tried making the rate of acceleration proportional to distance^2 and the results don't look any more realistic.

Here is the code and results:

enter image description here

enter image description here

The above graph is a Position Vs. Time graph.

How can I adjust the code to obtain more reasonable results?

  • 1
    $\begingroup$ I think you should add in your question what you expect to happen. It sounds to me like "accelerated towards each other at a rate proportional to the distance between them" means this should be a harmonic oscillator. Is that correct?? $\endgroup$ – Jared Jun 28 '15 at 5:18
  • 1
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$ – Qmechanic Jun 28 '15 at 5:44
  • 1
    $\begingroup$ @Qmechanic - I agree with Jared: this is about numerical methods for physics. While numerical integration is something that non-physicists study, the techniques needed to tackle this should be of interest to this community. The currently accepted answer doesn't do the question justice... $\endgroup$ – Floris Jun 28 '15 at 6:19
  • 1
    $\begingroup$ I was going to be nasty and write a nasty community wiki answer. Instead I'll just be nasty in a comment and say that this code is awful. $\endgroup$ – David Hammen Jun 28 '15 at 7:06
  • 4
    $\begingroup$ I'm voting to close this question as off-topic because it's about debugging a program and not physics. Perhaps Computational Science might be better suited for this program. $\endgroup$ – Kyle Kanos Jun 29 '15 at 1:16

In situations like this, it is a good idea to adjust the time step based on the gradient of the force - because the whole concept of numerical integration is that "things don't change too much from now until the next time step", and that assumption is violated when you move rapidly through a region with fast-changing force.

This has a risky side-effect: if you make the time step "small enough" to cope with the rapid variation, and the variation becomes "infinitely rapid", your algorithm may get "stuck" - you may run into Xeno's paradox of the Achilles and the tortoise.

To avoid getting stuck, some numerical integrations (like Runge-Kutta) are better at taking into account higher order curvature - allowing a bigger step without losing accuracy.

In your case though, the acceleration for the entire next time step is determined by the current position - so if your particles happen to be close to each other at the start of the step, they get thrown far away, and since their attractive force will then be greatly diminished they will have a hard time getting back.

A few criticisms on your code:

  1. You define particle1[4] but only use two components
  2. You hard code the time step instead of using a variable like dt
  3. You use the most basic integration method... please learn about others
  4. You implicitly define mass = 1 and force constant = 1; consider making those variables (even if you set them to 1)

Regarding the first point, I would use a struct for my particles:

typedef struct{
  double x;
  double v;

and then you can access their properties with

PARTICLE p1, p2;
p1.v = 0.;
p1.x = 0.;
p2.v = 0.;
p2.x = 1.;

Which is immediately more readable... If you want to work in two dimensions, you can either write them explicitly as part of your struct: (double x; double y; double vx; double vy;) or make them arrays.

But really - variable time steps and higher order interpolations would help with the stability of your solution.

  • 1
    $\begingroup$ Thanks, you hit the nail on the head of what is wrong with the program. Thanks for your tips too, I'm just starting out in numerical methods. $\endgroup$ – Kenshin Jun 28 '15 at 6:33
  • 1
    $\begingroup$ Writing a well behaved simulator is hard. But the rewards are awesome as you learn all there is to know. You are in step one and this answer takes you through steps 2, 3 and 4. The more you work at it the more nuances you tackle taking you through steps 5...12. You have earned a PhD at this point. $\endgroup$ – ja72 Jun 28 '15 at 22:37

You're getting particles that get 0 distance from each other, and therefore have infinite acceleration. The computer cannot resolve this, and so it explodes according to how close the step gets to zero.

I'm not sure what "physical" behaviour would be, given that zero distance point masses are difficult to study. But if you added in a y coordinate and gave it an initial velocity you might get better results (basically an elongated ellipse instead of an oscillating line).

You could make the initial y velocity very small, so that you approximate the X movement. Or maybe just add a small offset to your division.

  • $\begingroup$ It's not chaos theory if it's a computer anomaly. $\endgroup$ – Jared Jun 28 '15 at 5:28
  • $\begingroup$ If you want "chaos theory", then you can look at my following animation where this forms a "moon". Each compositie sphere is composed of randomly distributed balls and sometimes you get a moon and sometimes you don't (I just so happened to record one of the times when a moon formed--which, with my fine tuning, is around 1 in 10 runs). $\endgroup$ – Jared Jun 28 '15 at 5:31
  • 1
    $\begingroup$ Pot-a-to, pot-ah-to. The term chaos theory is not well defined anyway. But what I mean is that tiny changes in initial conditions cause large changes in output. $\endgroup$ – Dr Xorile Jun 28 '15 at 5:32
  • 2
    $\begingroup$ So like if the OP changes his time step, the oscillation changes dramatically? ;-) $\endgroup$ – Dr Xorile Jun 28 '15 at 5:35
  • 1
    $\begingroup$ I'm just kidding. You're right, of course. But the point is that the calculation is ill defined because of the divide by zero error that should occur if the time step were just right. $\endgroup$ – Dr Xorile Jun 28 '15 at 5:36

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