Two Body problem in python using RK4

jbay

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

D H

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.

jbay

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

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.

jbay

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

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.

jbay

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

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

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.

jbay

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

rcgldr

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.

D H

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.

Similar discussions for: Two Body problem in python using RK4
Thread	Forum	Replies
Is the energy conserved FOR EACH BODY in a two-body central force problem?	Classical Physics	2
Python ode problem (need help)	Programming & Comp Sci	9
Two Body problem using Verlet Algorithm Python	Engineering, Comp Sci, & Technology Homework	0
Python recursion problem	Programming & Comp Sci	8
Python problem	Programming & Comp Sci	3

rcgldr Recognitions: Homework Help	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.