# Two Body problem in python using RK4

by jbay
Tags: kepler's problem, physics, program, python, rk4
 P: 9 So I am writing a program in python to do RK4 for the two body problem. I want it to display a sphere moving around another. It currently displays one sphere for a split second and then it goes blank. Any suggestions? from __future__ import division from visual import * from visual.graph import * G = 6.67E-11 msat = 1500E21 mearth = 5.98E24 i = 0 h = 100 earth = sphere(pos = (0,0,0), radius = 6378100, color = color.blue) sat = sphere(pos = (384403000,0,0), radius = 3000000, color = color.red) pearth = vector(0,-1000000E21,0) psat = vector(0,1000000E21,0) trail = curve(color = sat.color) satvelocity = psat / msat earthvelocity = pearth / mearth satacceraltion = G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) Ag = arrow(pos=sat.pos, axis = 10000000000 * satacceraltion, color = color.red) vel = arrow(pos=sat.pos, axis = 100000 * (psat / msat), color = color.red) z = [earth.pos, earthvelocity, sat.pos, satvelocity] Force = [earthvelocity, earthacceraltion, satvelocity, satacceraltion] k1 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)] k2 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)] k3 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)] k4 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)] while 1: (100000) psat = psat + msat * Force[3] * h pearth = pearth - msat * Force[3] * h sat.pos = sat.pos + Force[2] * h earth.pos = earth.pos + Force[0] * h trail.append(pos = sat.pos) Ag.pos = sat.pos Ag.axis = 10000000000 * Force[3] vel.pos = sat.pos vel.axis = 100000 * Force[3] while (i < 4): k1[i] = h * Force[i] k2[i] = h *(z[i] + .5 * k1[i]) k3[i] = h *(z[i] + .5 * k2[i]) k4[i] = h *(z[i] + k3[i]) Force[i] = z[i] + h / 6 * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) i = i + 1 i = 0
 Mentor P: 14,432 Your loop and your use of RK4 are a mess. You are only calculating acceleration once, outside the main loop. Acceleration changes constantly. You need to calculate the accelerations with each RK step, and you need to calculate the accelerations of the Earth and satellite simultaneously. (In other words, don't calculate the acceleration of the satellite, move the satellite, calculate the acceleration of the Earth, and finally move the Earth.) What you should be doing:RK4 step 1:Save the position and velocity of the Earth and the satellite in some variable. You'll need this position and velocity state at the start of the integration interval for every subsequent RK step. Calculate the gravitational accelerations of the Earth and the satellite toward one another via Newton's laws. Save the velocity and acceleration of the Earth and the satellite in k1. Advance the position and velocity of the Earth and the satellite by half a time step using the derivatives in k1. RK4 step #2:Calculate the gravitational accelerations of the Earth and the satellite toward one another via Newton's laws. Save the velocity and acceleration of the Earth and the satellite in k2. Advance the position and velocity of the Earth and the satellite at the start of the integration interval by half a time step using the derivatives in k2. RK4 step #3:Calculate the gravitational accelerations of the Earth and the satellite toward one another via Newton's laws. Save the velocity and acceleration of the Earth and the satellite in k3. Advance the position and velocity of the Earth and the satellite at the start of the integration interval by a full time step using the derivatives in k3. RK4 step #4:Calculate the gravitational accelerations of the Earth and the satellite toward one another via Newton's laws. Save the velocity and acceleration of the Earth and the satellite in k4. Advance the position and velocity of the Earth and the satellite at the start of the integration interval by a full time step using the RK4 step #4 weighted average of the derivatives in k1, k2, k3, and k4.
 P: 9 Thank you! I thought I was updating incorrectly. I had a little trouble understanding you description at first but here is the steps coded out instead of a loop, it displays not how i want yet but I assume that is due to scaling and etc. Tell me what you think? while 1: (100000) k1v1 = h * earthvelocity k1v2 = h * satvelocity k1a1 = h * earthacceraltion k1a2 = h * satacceraltion earth.pos = earth.pos + .5 * k1v1 sat.pos = sat.pos + .5 * k1v2 earthvelocity = earthvelocity + .5 * k1a1 satvelocity = satvelocity + .5 * k1a2 satacceraltion = G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) trail.append(pos = sat.pos) Ag.pos = sat.pos Ag.axis = 10000000000 * satacceraltion vel.pos = sat.pos vel.axis = 100000 * satacceraltion k2v1 = h * earthvelocity k2v2 = h * satvelocity k2a1 = h * earthacceraltion k2a2 = h * satacceraltion earth.pos = earth.pos + .5 * k2v1 sat.pos = sat.pos + .5 * k2v2 earthvelocity = earthvelocity + .5 * k2a1 satvelocity = satvelocity + .5 * k2a2 satacceraltion = G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) trail.append(pos = sat.pos) Ag.pos = sat.pos Ag.axis = 10000000000 * satacceraltion vel.pos = sat.pos vel.axis = 100000 * satacceraltion k3v1 = h * earthvelocity k3v2 = h * satvelocity k3a1 = h * earthacceraltion k3a2 = h * satacceraltion earth.pos = earth.pos + k3v1 sat.pos = sat.pos + k3v2 earthvelocity = earthvelocity + k3a1 satvelocity = satvelocity + k3a2 satacceraltion = G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) trail.append(pos = sat.pos) Ag.pos = sat.pos Ag.axis = 10000000000 * satacceraltion vel.pos = sat.pos vel.axis = 100000 * satacceraltion k4v1 = h * earthvelocity k4v2 = h * satvelocity k4a1 = h * earthacceraltion k4a2 = h * satacceraltion earth.pos = earth.pos + h / 6 * (k1v1 + 2 * k2v1 + 2 * k3v1 + k4v1) sat.pos = sat.pos + h / 6 * (k1v2 + 2 * k2v2 + 2 * k3v2 + k4v2) earthvelocity = earthvelocity + h / 6 * (k1a1 + 2 * k2a1 + 2 * k3a1 + k4a1) satvelocity = satvelocity + h / 6 * (k1a2 + 2 * k2a2 + 2 * k3a2 + k4a2) satacceraltion = G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) trail.append(pos = sat.pos) Ag.pos = sat.pos Ag.axis = 10000000000 * satacceraltion vel.pos = sat.pos vel.axis = 100000 * satacceraltion
Mentor
P: 14,432

## Two Body problem in python using RK4

In the future, please wrap code tags around your code. Use [ code ] (but with no spaces before or after "code") to start a block of code, [ /code ] (but with no spaces again) to end the block.

For example,
[ code ] # This is python code [ /code ]
 # This is python code
There are still a number of problems with your code.
• Your step #4 is incorrect. This is the cause of the immediate problem.
• You aren't saving the initial state, the state at the start of RK4 step #1.
• Your gravitational accelerations are incorrect. Gravitation is attractive, not repulsive.
• You shouldn't be plotting those intermediate states (RK4 steps #1 to 3). They're not quite "real". All that matters is the final state, the result of step #4.
• There's a lot of replicated code here. That's a sign of "code smell".
• Would your code scale up to handle a N-body problem? Even the solar system will be ungainly, let alone a large value of N.
 P: 9 First this is just the loop from my first post, my initial states are calculated above the loop. Code posted below. Why do I need them over and over again? Second, the accelerations, Earths g = - 9.8, so each acceleration with respect to the other should be negative and allows orbit. how are they incorrect? What is wrong with step number 4? Please be specific, i see nothing when comparing it to the general process formula. I do not mind code smell at this point, I will worry about n - bodies when I am done learning RK4, two bodies completely and more python. from __future__ import division from visual import * from visual.graph import * G = 6.67E-11 msat = 1500E21 mearth = 5.98E24 h = 100 earth = sphere(pos = (0,0,0), radius = 6378100, color = color.blue) sat = sphere(pos = (384403000,0,0), radius = 3000000, color = color.red) pearth = vector(0,-1000000E21,0) psat = vector(0,1000000E21,0) trail = curve(color = sat.color) satvelocity = psat / msat earthvelocity = pearth / mearth satacceraltion = -1 * G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = -1 * G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) Ag = arrow(pos=sat.pos, axis = 10000000000 * satacceraltion, color = color.red) vel = arrow(pos=sat.pos, axis = 100000 * (psat / msat), color = color.blue) display1 = gdisplay(x=600, y=0, width=600, height=600, foreground=color.black, background=color.white) while 1: (100000000) k1v1 = h * earthvelocity k1v2 = h * satvelocity k1a1 = h * earthacceraltion k1a2 = h * satacceraltion earth.pos = earth.pos + .5 * k1v1 sat.pos = sat.pos + .5 * k1v2 earthvelocity = earthvelocity + .5 * k1a1 satvelocity = satvelocity + .5 * k1a2 satacceraltion = -1 * G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = -1 * G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) k2v1 = h * earthvelocity k2v2 = h * satvelocity k2a1 = h * earthacceraltion k2a2 = h * satacceraltion earth.pos = earth.pos + .5 * k2v1 sat.pos = sat.pos + .5 * k2v2 earthvelocity = earthvelocity + .5 * k2a1 satvelocity = satvelocity + .5 * k2a2 satacceraltion = -1 * G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = -1 * G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) k3v1 = h * earthvelocity k3v2 = h * satvelocity k3a1 = h * earthacceraltion k3a2 = h * satacceraltion earth.pos = earth.pos + k3v1 sat.pos = sat.pos + k3v2 earthvelocity = earthvelocity + k3a1 satvelocity = satvelocity + k3a2 satacceraltion = -1 * G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = -1 * G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) trail.append(pos = sat.pos) k4v1 = h * earthvelocity k4v2 = h * satvelocity k4a1 = h * earthacceraltion k4a2 = h * satacceraltion earth.pos = earth.pos + 1 / 6 * (k1v1 + 2 * k2v1 + 2 * k3v1 + k4v1) sat.pos = sat.pos + 1 / 6 * (k1v2 + 2 * k2v2 + 2 * k3v2 + k4v2) earthvelocity = earthvelocity + 1 / 6 * (k1a1 + 2 * k2a1 + 2 * k3a1 + k4a1) satvelocity = satvelocity + 1 / 6 * (k1a2 + 2 * k2a2 + 2 * k3a2 + k4a2) satacceraltion = -1 * G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = -1 * G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) trail.append(pos = sat.pos) Ag.pos = sat.pos Ag.axis = 10000000000 * satacceraltion vel.pos = sat.pos vel.axis = 100000 * satvelocity func1 = gdots(gdisplay=display1, color=color.black) func1.plot(pos=(sat.pos.x, sat.pos.y))
Mentor
P: 14,432
 Quote by jbay First this is just the loop from my first post, my initial states are calculated above the loop.
It's the initial state at the start of each integration cycle (the four steps of RK4) that you are missing. You aren't using the algorithm correctly, and you are getting incorrect results as a consequence.

 Second, the accelerations, Earths g = - 9.8, so each acceleration with respect to the other should be negative and allows orbit. how are they incorrect?
What you had was (post #4) satacceraltion = G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos). This had the wrong sign. You've corrected this sign error in your newest version of your code (post #6).

 What is wrong with step number 4? Please be specific, i see nothing when comparing it to the general process formula.
Once again you appear to have corrected the problem. In post #3 you had sat.pos = sat.pos + h / 6 * (k1v2 + 2 * k2v2 + 2 * k3v2 + k4v2); now you have sat.pos = sat.pos + 1 / 6 * (k1v2 + 2 * k2v2 + 2 * k3v2 + k4v2). Note: You might want to change that 1/6 to 1.0/6.0; it's best not to give a language the chance of screwing up with integer division.
 P: 9 Oh! I believe I have corrected it! But I have a program for Euler that seems to work at extremely high time steps where this RK4 seems to displays the sat. shot straight away form the earth at the same time steps. Also thank you for all your help, I had previously had no experience with RK4, python and this forum lol k1v1 = h * earthvelocity k1v2 = h * satvelocity k1a1 = h * earthacceraltion k1a2 = h * satacceraltion earth.pos1 = earth.pos + .5 * k1v1 sat.pos1 = sat.pos + .5 * k1v2 earthvelocity1 = earthvelocity + .5 * k1a1 satvelocity1 = satvelocity + .5 * k1a2 satacceraltion1 = -1 * G * mearth * norm(sat.pos1 - earth.pos1) / mag2(sat.pos1 - earth.pos1) earthacceraltion1 = -1 * G * msat * norm(earth.pos1 - sat.pos1) / mag2(earth.pos1 - sat.pos1) k2v1 = h * earthvelocity1 k2v2 = h * satvelocity1 k2a1 = h * earthacceraltion1 k2a2 = h * satacceraltion1 earth.pos2 = earth.pos + .5 * k2v1 sat.pos2 = sat.pos + .5 * k2v2 earthvelocity2 = earthvelocity + .5 * k2a1 satvelocity2 = satvelocity + .5 * k2a2 satacceraltion2 = -1 * G * mearth * norm(sat.pos2 - earth.pos2) / mag2(sat.pos2 - earth.pos2) earthacceraltion2 = -1 * G * msat * norm(earth.pos2 - sat.pos2) / mag2(earth.pos2 - sat.pos2) k3v1 = h * earthvelocity2 k3v2 = h * satvelocity2 k3a1 = h * earthacceraltion2 k3a2 = h * satacceraltion2 earth.pos3 = earth.pos + k3v1 sat.pos3 = sat.pos + k3v2 earthvelocity3 = earthvelocity + k3a1 satvelocity3 = satvelocity + k3a2 satacceraltion3 = -1 * G * mearth * norm(sat.pos3 - earth.pos3) / mag2(sat.pos3 - earth.pos3) earthacceraltion3 = -1 * G * msat * norm(earth.pos3 - sat.pos3) / mag2(earth.pos3 - sat.pos3) k4v1 = h * earthvelocity3 k4v2 = h * satvelocity3 k4a1 = h * earthacceraltion3 k4a2 = h * satacceraltion3 earth.pos = earth.pos + 1.0 / 6.0 * (k1v1 + 2 * k2v1 + 2 * k3v1 + k4v1) sat.pos = sat.pos + 1.0 / 6.0 * (k1v2 + 2 * k2v2 + 2 * k3v2 + k4v2) earthvelocity = earthvelocity + 1.0 / 6.0 * (k1a1 + 2 * k2a1 + 2 * k3a1 + k4a1) satvelocity = satvelocity + 1.0 / 6.0 * (k1a2 + 2 * k2a2 + 2 * k3a2 + k4a2) satacceraltion = -1 * G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos) earthacceraltion = -1 * G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos) trail.append(pos = sat.pos) Ag.pos = sat.pos Ag.axis = 10000000000 * satacceraltion vel.pos = sat.pos vel.axis = 100000 * satvelocity func1 = gdots(gdisplay=display1, color=color.black) func1.plot(pos=(sat.pos.x, sat.pos.y))
 HW Helper P: 6,925 I don't know how much this would affect the results, but you calculate velocity2 based on acceleration1. It might be better to calculate accelerations first, then calculate velocites base on the new accelerations, then position based on the new velocities, so that velocity2 would be based on acceleration2, position2 based on velocity2, and so on.
Mentor
P: 14,432
 Quote by rcgldr I don't know how much this would affect the results, but you calculate velocity2 based on acceleration1. It might be better to calculate accelerations first, then calculate velocites base on the new accelerations, then position based on the new velocities, so that velocity2 would be based on acceleration2, position2 based on velocity2, and so on.
That's wrong. What you are describing is no longer RK4, and it turns out to be a low accuracy technique.

RK4 is a first order ODE technique. To apply it, or any other first order ODE technique, to a second order ODE, you essentially have to make the state that is being integrated have 2N elements (6 elements in this case).

What jbay is doing is getting close to being correct. The problem is that the algorithm he/she is using also isn't RK4 (yet). jbay needs to do just what I said, which is to save the state at RK4 step #1 and use that as the basis for propagating the state for steps 2 to 4.
 P: 9 "What jbay is doing is getting close to being correct. The problem is that the algorithm he/she is using also isn't RK4 (yet). jbay needs to do just what I said, which is to save the state at RK4 step #1 and use that as the basis for propagating the state for steps 2 to 4. " I am confused to what you are referring. If you can can you pick out what I am doing that is incorrect from the code then this would be helpful, because I believe I use that step for the rest of the steps.
HW Helper
P: 6,925
 Quote by D H [*]RK4 step 1:Save the position and velocity of the Earth and the satellite in some variable. You'll need this position and velocity state at the start of the integration interval for every subsequent RK step. Calculate the gravitational accelerations of the Earth and the satellite toward one another via Newton's laws. Save the velocity and acceleration of the Earth and the satellite in k1. Advance the position and velocity of the Earth and the satellite by half a time step using the derivatives in k1.
 Quote by rcgldr It might be better to calculate accelerations first, then calculate velocites base on the new accelerations, then position based on the new velocities, so that velocity2 would be based on acceleration2, position2 based on velocity2, and so on.
 Quote by D H That's wrong. What you are describing is no longer RK4, and it turns out to be a low accuracy technique. RK4 is a first order ODE technique. To apply it, or any other first order ODE technique, to a second order ODE, you essentially have to make the state that is being integrated have 2N elements (6 elements in this case).
I wasn't sure (I've seen "RK4" done both ways), and I was hoping you would clarify this. In your "step 1", I wasn't sure if "save the velocity and acceleration in k1" ... "advance position and velocity ... by half time step using derivatives in k1", meant to advance the velocites based on the initial accelerations or based on the accelerations calculated for k1.
Mentor
P: 14,432
 Quote by jbay Oh! I believe I have corrected it!
Yep. Looks good!

 But I have a program for Euler that seems to work at extremely high time steps where this RK4 seems to displays the sat. shot straight away form the earth at the same time steps.
By "Euler" I assume you mean symplectic Euler, aka Euler-Cromer, aka a whole lot of other names. Simple Euler, where one updates position and velocity simultaneously, is almost always bad. About the only place it isn't bad is uniform linear motion. Symplectic Euler simply changes the update procedure so that velocity is updated first and position is updated using the updated velocity. This simple change has some remarkable qualities. It comes very close to conserving energy and angular momentum. RK4 does not do so.

RK4 is going to do a very good job for quite some time if you use a "reasonable" time step (between 100 to 300 seconds with your problem). You won't see that it isn't conserving energy and angular momentum until after thousands of orbits have passed. Symplectic Euler is going to do a lousy job with that time step in the sense that exactly where the objects are in their orbits quickly becomes very bad. It will however keep the objects in more or less the right orbits.

This is an inherent problem with numerical integration techniques. There is often a tradeoff between stability and accuracy. RK4 does a good job on short-term accuracy, but at the cost of long term stability (and hence long term accuracy). Symplectic Euler does a really lousy job with regard to short term accuracy, but does a much better job with regard to stability.

There are a number of techniques that are considerably better than symplectic Euler in terms of both short-term and long-term accuracy. For example, leap frog and its variants (position verlet and velocity verlet are essentially leap frog) are much, much better than symplectic Euler. The problem is that these techniques are still second order in error (RK4 is forth order). Going beyond second order in error while maintaining symplecticity is very, very hard.

 Related Discussions Classical Physics 2 Programming & Computer Science 9 Engineering, Comp Sci, & Technology Homework 0 Programming & Computer Science 8 Programming & Computer Science 3