Two Body problem in python using RK4

• Python
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?

Code:
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

Last edited by a moderator:

D H
Staff Emeritus
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.

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

D H
Staff Emeritus
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 ]
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.

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.

Code:
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))

D H
Staff Emeritus
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.

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

Code:
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))

rcgldr
Homework Helper
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.

D H
Staff Emeritus
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.

"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.

Last edited:
rcgldr
Homework Helper
[*]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.

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).

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.

D H
Staff Emeritus