Register to reply 
Two Body problem in python using RK4 
Share this thread: 
#1
Nov1112, 11:34 PM

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?



#2
Nov1212, 12:09 PM

Mentor
P: 15,066

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:



#3
Nov1212, 04:08 PM

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 


#4
Nov1312, 06:25 AM

Mentor
P: 15,066

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 ]



#5
Nov1312, 08:59 AM

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.



#6
Nov1312, 10:04 AM

Mentor
P: 15,066




#7
Nov1312, 11:48 AM

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



#8
Nov1312, 12:31 PM

HW Helper
P: 7,034

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.



#9
Nov1312, 01:02 PM

Mentor
P: 15,066

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. 


#10
Nov1312, 01:22 PM

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. 


#11
Nov1312, 05:27 PM

HW Helper
P: 7,034




#12
Nov1412, 12:58 PM

Mentor
P: 15,066

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 shortterm 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 shortterm and longterm 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. 


Register to reply 
Related Discussions  
Is the energy conserved FOR EACH BODY in a twobody central force problem?  Classical Physics  2  
Python ode problem (need help)  Programming & Computer Science  9  
Two Body problem using Verlet Algorithm Python  Engineering, Comp Sci, & Technology Homework  0  
Python recursion problem  Programming & Computer Science  8  
Python problem  Programming & Computer Science  3 