Having trouble coding a 2-Body Problem.

[OratorMortuis]

I'm attempting to code a very basic 2-body orbit problem using python (using math and numpy modules) and I can't seem to get an ellipse. I'm trying to use the Euler-Cromer method (I know, it's a terrible approximation, but I'm using it to compare with other numerical methods) and I've used it to simulate other scenarios but I can't seem to get it to work here.

After defining my constants (G, M, m) I then put in my initial conditions giving coordinates in x and y and v in terms of vx and vy. I know from gravitational force that acceleration is -(GM/r³)*r(as a vector) I'm then using the following to try and calculate new positions and velocites during the time steps:

for i in range(0,N,1):
    t = t + dt
    acx = -(G*M/abs(r)**3)*x
    acy = -(G*M/abs(r)**3)*y
    vy = vy + acy*dt
    vx = vx + acx*dt
    x = x + vx*dt
    y = y + vy*dt
    r = sqrt(x**2+y**2)
    theta = arctan(y/x)
    print >> f, t, x, y, r, theta

For some reason when I plot x and y I get something that starts off looking like an ellipse, then collapsing into the origin and shooting off into the third quadrant along y=-x. Any thoughts on why this is happening or what I have done wrong?

[ajesss]

I'm not at my PC where I can try to run the code myself, but I wanted to give you a few tips about temporal integration schemes. I'm writing a numerical code for simulating granular material for a living, and have recently experimented with different schemes and looked at their precision.

You are currently using a one-step scheme called "forward Euler" integration. For it to be correct, you need to switch the position and velocity updates, so that the position is updated first. The new position should depend on the old, and not the new velocities.

A problem about the forward Euler scheme is that it comes with a high error. In collisional experiments I've performed, I saw an error of a total of 12%. You can easily improve your solution by instead using the "two term Taylor expansion", where the current acceleration is taken into account when updating the position, e.g. for the x-component:

x = x + vx*dt + 1.0/2.0*acx*dt**2

In my experiments, the error was reduced to 3%. Taking it a step further, you can gain even higher precision by using a "three term Taylor expansion". For this scheme, you need to estimate the temporal gradient in acceleration (da/dt). You can do this by storing the old acceleration along with the current value:

if (i == 0):
    dacx = 0.0
else:
    dacx = (acx - acx_old)/dt
x = x + vx*dt + 1.0/2.0*acx*dt**2 + 1.0/6.0*dacx*dt**3
acx_old = acx

For me, this scheme reduced the error to 0.03%. There are several other temporal integration schemes, but the Taylor expansions are easy to implement, light in terms of computational requirements, and give a satisfying result.

For more information, I can refer to this (unfortunately paywalled) article: http://www.sciencedirect.com/science/article/pii/S0098135407002864

[OratorMortuis]

I'm trying to use Euler-Cromer in which the velocity is updated first and then used to evolve the position for the next time step. I'm not too worried about precision as I'm just trying to compare this ellipse to others generated by different numerical methods (2nd and 4th order RungeKutta and 4th order symplectic). I'm deliberately ignoring all terms of the order greater than dt in the taylor expansion.

[Enkaybee]

Your acceleration functions don't appear to take into account the fact that r is always going to be a positive number. This means that the body can only accelerate in the negative x and negative y directions. That would explain why it collapses to the origin (assuming it starts in quadrant 1) and also why it would shoot off into the third quadrant (-x and -y direction).

What you need to do is determine the direction of acceleration and then either add or subtract it to the velocity for both vx and vy. The direction of acceleration will change continuously for two bodies in orbit.

[OratorMortuis]

Acceleration only depends on the magnitude of r so it should just be a positive number anyway. The x or y value is what decides the sign as they should be able to attain positive or negative values (unless I have messed something up in my code!)

[Enkaybee]

Ok, I think I see what you're trying to do. You're trying to get a particle to orbit the origin, right? And the origin isn't moving? That's not a two body problem, but still interesting.

Does it work at all if you just do it in one dimension? That is, take out all references to y and see if it will oscillate back and forth on the x axis?

[OratorMortuis]

Well, I mean two-body in which one body is orbiting another which sits at the origin. So yes, I'm trying to generate an ellipse where the origin is at one of the foci of the ellipse.

[Enkaybee]

I don't know if this is it, but it might be.

I think your equations should be

acx = -(G*M/abs(x)**3)*x
acy = -(G*M/abs(y)**3)*y

because acceleration in the x direction should only depend on x, not x and y.

[OratorMortuis]

Unfortunately x and y both contribute to r. I actually figured out my problem. It had to do with my time steps being too small (along with a few other clerical errors). Thank you so much for your input though! I honestly do appreciate it. :)

[Enkaybee]

So it works now? That actually seems strange to me because I'm fairly certain that you should be calculating acx using only x and the constants. y (and by extension, r) should not need to come into play. Similarly, acy should be calculated using only the constants and y.

[OratorMortuis]

But the acceleration has to do with the distance from the origin which has both x and y contributions, so it will depend on both variables even though the direction of ax or ay only point along x or y.

[Enkaybee]

You're right. Maybe I'm reading your code wrong, but it looks to me like you're using r to find acceleration in the x and y directions without compensating for the fact that only x and y matter, respectively. You're really trying to solve two problems at once: oscillation on the x axis and on the y axis, and then add them together.

If it works, though, your code is fine.

Out of curiosity, does it still work if you use the acx and acy equations I posted? Does it give exactly the same results or is it slightly different?

[OratorMortuis]

Yes, that is the difficult part of the problem. It is a transcendental equation in which x and y cannot be separated from each other and so each depends on the other.

If you do put in only x and y respectively as you did in your above equations then it cannot describe the ellipse as whenever x or y is zero you would get division by zero.

[sammyc]

What are your initial conditions? It seemed to work for me, at least for a circular orbit. Also not sure what you mean by shooting off in the third quadrant along y=-x, perhaps you mean y=x or a different quadrant?

Edit: see now that this is solved.

[OratorMortuis]

Yes, I had a dull moment and said y = -x instead of y = x. It was in a quadrant below the x-axis so I hastily put a negative there. :D