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:
- You define
particle1 but only use two components
- You hard code the time step instead of using a variable like
- You use the most basic integration method... please learn about others
- 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:
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
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.