Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Two Body problem in python using RK4

  1. Nov 11, 2012 #1
    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 (Text):
    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: Nov 12, 2012
  2. jcsd
  3. Nov 12, 2012 #2

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    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.
     
  4. Nov 12, 2012 #3
    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
     
  5. Nov 13, 2012 #4

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    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 (Text):
     # 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.
     
  6. Nov 13, 2012 #5
    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 (Text):
    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))
     
  7. Nov 13, 2012 #6

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    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.

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

    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.
     
  8. Nov 13, 2012 #7
    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 (Text):
    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))
     
  9. Nov 13, 2012 #8

    rcgldr

    User Avatar
    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.
     
  10. Nov 13, 2012 #9

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    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.
     
  11. Nov 13, 2012 #10
    "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: Nov 13, 2012
  12. Nov 13, 2012 #11

    rcgldr

    User Avatar
    Homework Helper

    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.
     
  13. Nov 14, 2012 #12

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Yep. Looks good!

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Two Body problem in python using RK4
  1. Python rk4 (Replies: 0)

  2. Help in Python rk4 (Replies: 1)

  3. N-body problem in Python (Replies: 10)

Loading...