Dismiss Notice
Join Physics Forums Today!

Applying Runge-Kutta(RK4) method to planetary velocity

  1. Feb 3, 2008 #1
    Hey guys,
    I'm working on this n-body planetary system simulation and I was just wondering if anyone was familiar enough with the Runge-Kutta(RK4) approximation method that they understood how to apply it to approximating the velocity and position of a planet as it orbits around some central star.

    The problem I'm having with incorporating the method into my simulation is that the acceleration (the derivative of the variable I'm trying to approximate) is not a function of time. Instead, it is being treated as constant (acceleration = total_force / mass).

    Since, I'm not really all that knowledgeable in physics, especially astrophysics, I have absolutely no idea how to define this acceleration as a function of time or even if I really need to in order to use the Runge-Kutta method in my simulation. Any help would be greatly appreciated. Thanks.
  2. jcsd
  3. Feb 3, 2008 #2

    D H

    Staff: Mentor

    The acceleration is not constant in the n-body problem. You instead need to use Newton's law of gravitation to compute the forces. The force acting on one body is the sum of the gravitation forces due to all other bodies.
  4. Feb 3, 2008 #3
    Hey, thanks for the reply. To address your comment, I do understand that the acceleration is dependent on the total forces acting on one body divided by the bodies mass and this is how I am calculating the acceleration for each body in my simulation. I guess I really just don't have a good understanding of the physics required for this. Euler integration is pretty simple and that is what I'm using right now but I need something more precise like RK4. Anyway, can someone help me understand how one would use RK4 in order to approximate the velocity of a planet at a given timestep. Any help would be greatly appreciated. Thanks.
  5. Feb 3, 2008 #4

    D H

    Staff: Mentor

    Let me go back to your first post.

    The phrase that I put in bold is incorrect. The acceleration is not constant. You must recompute at every time step. In fact, you must you must recompute it at every time step. Just because the acceleration does not explicitly depend on time does not mean it doesn't depend on time at all. The gravitational acceleration from object A towards object B depends on the displacement vector from object A to object B, and since this varies with time, so does the acceleration. The acceleration implicitly depends on time.

    The basic idea underlying all Runge Kutta methods is to make multiple Euler steps within a time step, and then take some weighted average of these Euler steps to form the state at the end of the time step. I suggest you read the Wikipedia article on Runge Kutta methods, link here. Note that the first section after the intro covers RK4. If you have any questions after reading that article, fire away.
  6. Feb 4, 2008 #5
    Hey, thanks for the response. To address your comments, I have read the wikipedia article and I do understand that the gravitational acceleration needs to be computed at every time step. I guess my problem before was that I was trying to use the RK4 method to approximate velocity when I should be using it to approximate the position of the planets at a given timestep. The problem I'm having in understanding and applying the Rk4 method, however, is that in the wikipedia article it states that y' = f(t,y). I assume that y' is the derivative of y(the variable one is trying to approximate, in this case, position). Calculating velocity with Euler Integration is not dependent on the position though:

    y = position
    y' = velocity

    y(n+1) = new_position = current_position + current_velocity*time_step + (0.5)*acceleration*(time_step^2)
    y' = f(t,y) = new_velocity = current_velocity + acceleration*time_step

    So my question is: when calculating the k slopes in the RK4 method, what do I plug in for y?

    I can't see the point of plugging in position since the velocity does not depend on it in the above equation.....
  7. Feb 4, 2008 #6

    D H

    Staff: Mentor

    You need to integrate position and velocity. I will call the position and velocity of the ith body in your simulation [itex]{\mathbf x}_i[/itex] and [itex]{\mathbf v}_i[/itex]. These quantities are bolded to indicate that they are vectors. The state, which I will call [itex]{\mathbf y}_i[/itex] to be consistent with the wiki article is an amalgam of the position and velocity:
    [tex]{\mathbf y}_i = \bmatrix {\mathbf x}_i \\ {\mathbf v}_i \endbmatrix[/tex].
    The time derivative of [itex]{\mathbf y}_i[/itex] is then
    [tex]{\mathbf y}^{\prime}_i =
    \bmatrix {\mathbf v}_i \\ \sum_{j\ne i} {\mathbf a}_{i,j} \endbmatrix[/tex]
    where [itex]{\mathbf a}_{i,j}[/itex] is the acceleration due to gravity of the ith body toward the jth body.
  8. Feb 4, 2008 #7
    Hey, again, thanks for the response. Maybe with all this programming that I've been doing over these past few years my applied math skills have become warped and corrupted in some sense so please forgive me. Anyway, I still can't see from what you posted how I am supposed to handle f(t(n), y(n)) and f(t(n) + h/2, y(n) + (h/2)*k1) as they relate to the Euler integration equations I mentioned above. I guess what is really throwing me off is plugging y(n) + (h/2)*k1 into the f function.
  9. Feb 4, 2008 #8

    D H

    Staff: Mentor

    The function "f" is nothing more than the derivative of the state -- i.e., y' -- and it does not depend explicitly on time. In other words, it is just f(y), not f(t,y), in this case.
  10. Feb 4, 2008 #9
    If it helps, here is the C++ code that I use for RK2. It is adapted from the RK4 sample at http://www.gaffer.org/game-physics/integration-basics/:

    // Where the vector_3 kn_acceleration is obtained from Newtonian gravitation() function:
    // ex: A = (GM/r^2) * unit 3-vector pointing from Mercury to the Sun

    inline void proceed_rk2(const double &dt)
    vector_3 k1_velocity = mercury.velocity;
    vector_3 k1_acceleration = gravitation(mercury.position, k1_velocity);
    vector_3 k2_velocity = mercury.velocity + k1_acceleration*dt*0.5;
    vector_3 k2_acceleration = gravitation(mercury.position + k1_velocity*dt*0.5, k2_velocity);

    mercury.velocity += k2_acceleration*dt;
    mercury.position += k2_velocity*dt;
    Last edited: Feb 4, 2008
  11. Feb 4, 2008 #10


    User Avatar
    Science Advisor
    Gold Member

    If you have Visual Basic 6.0, here's some working code. You'll need to add 3 text boxes and 2 buttons to your form. It's only a 2-body problem, with one body massless. I copied this from a book, and modified it for VB6, so don't ask me how it works :) I haven't yet been able to modify it to perform n-bodies.
    Code (Text):

    Option Explicit
    Dim R0 As Double, A As Double, M As Double, EN As Double, nsteps As Double
    Dim EX As Double, X As Double, EY As Double, Y As Double, EVX As Double, VX As Double, EVY As Double, VY As Double
    Dim EAX1 As Double, EAX As Double, EAY1 As Double, EAY As Double
    Dim EVX1 As Double, EVY1 As Double
    Dim EC As Double, SF As Double
    Dim xCenter As Double, yCenter As Double
    Dim t0 As Double, TF As Double, PD As Double, DT As Double, V As Double
    Dim RKxERROR As Double, RKyERROR As Double, RKERROR As Double, RKRELERROR As Double, PI  As Double, j As Integer
    Dim TmStp As Double, EAX2 As Double, EAY2 As Double, EVX2 As Double, EVY2 As Double, EAX3 As Double, EAY3 As Double
    Dim EVX3 As Double, EVY3 As Double, EAX4 As Double, EAY4 As Double, EVX4 As Double, EVY4 As Double
    Dim AVAX As Double, AVAY As Double, AVVX As Double, AVVY As Double, R2 As Double, R As Double, R3 As Double
    Dim delay As Double, oldx As Double, oldy As Double
    Sub RungeKutta4()
    Dim j As Double
    Dim NoOfOrbits As Double
    NoOfOrbits = Val(Text3.Text)
    For j = 1 To nsteps * NoOfOrbits
        EX = X: EY = Y: EVX = VX: EVY = VY 'E = Estimated
        Call CalcEstAccels
        EAX1 = EAX: EAY1 = EAY 'Save the estimates 'A = Acceleration
        EVX1 = EVX: EVY1 = EVY 'Save the estimates 'V = Velocity
        TmStp = DT / 2# 'Half-step estimates, TmStp = Time Step
        Call EstVals
        Call CalcEstAccels
        EAX2 = EAX: EAY2 = EAY 'Save'em
        EVX2 = EVX: EVY2 = EVY 'Save'em
        'Boot-stap up
        Call EstVals
        Call CalcEstAccels
        EAX3 = EAX: EAY3 = EAY 'Save'em
        EVX3 = EVX: EVY3 = EVY 'Save'em
        'Now a full-step & more boot-strapping
        TmStp = DT
        Call EstVals
        Call CalcEstAccels
        EAX4 = EAX: EAY4 = EAY 'Save'em
        EVX4 = EVX: EVY4 = EVY 'Save'em
        'Calculate weighted average accelerations
        AVAX = (EAX1 + 2# * (EAX2 + EAX3) + EAX4) / 6#
        AVAY = (EAY1 + 2# * (EAY2 + EAY3) + EAY4) / 6#
        'Calculate weighted average velocities
        AVVX = (EVX1 + 2# * (EVX2 + EVX3) + EVX4) / 6#
        AVVY = (EVY1 + 2# * (EVY2 + EVY3) + EVY4) / 6#
        'Final integrations
        X = X + DT * AVVX: Y = Y + DT * AVVY
        VX = VX + DT * AVAX: VY = VY + DT * AVAY
        'PSet (oldx, oldy), Me.BackColor
        PSet (X + xCenter, yCenter - Y)
        oldx = X + xCenter: oldy = yCenter - Y
        'For delay = 1 To 1000000: Next delay
    Next j
    End Sub

    Sub CalcEstAccels() 'Inputs EX, EY & M; Outputs EAX & EAY
    R2 = EX * EX + EY * EY 'Faster than EX ^ 2 + EY ^ 2
    R = Sqr(R2): R3 = R * R2 'Faster than R2 ^ 1.5
    EAX = -M * EX / R3
    EAY = -M * EY / R3
    End Sub

    Sub EstVals() 'Inputs X, Y, VX, VY, EAX, EAY, EVX, EVY & TmStp 'Outputs EVX, EVY, EX & EY
    'Must be done in the right order!
    EX = X + TmStp * EVX
    EY = Y + TmStp * EVY
    EVX = VX + TmStp * EAX
    EVY = VY + TmStp * EAY
    End Sub

    Private Sub Command1_Click() ' Do It
    EC = Val(Text1.Text)
    nsteps = Val(Text2.Text)
    Circle (xCenter, yCenter), 3
    M = 200#: R0 = 300# 'M = mass, R0 = apoapsis distance
    Call CalcInitVals
    t0 = 0: TF = PD
    DT = (TF - t0) / nsteps
    X = R0: Y = 100#: VX = 0#: VY = V 'Initial values
    Call RungeKutta4
    RKxERROR = R0 - X: RKyERROR = Y
    RKERROR = Sqr(RKxERROR ^ 2# + RKyERROR ^ 2#)
    End Sub

    Sub CalcInitVals()
    V = Sqr(M * (1 - EC) / R0) 'Initial velocity
    EN = (V * V / 2) - (M / R0) 'Specific mechanical energy
    A = -M / (2 * EN) 'Semi-major axis
    PI = 4# * Atn(1#)
    PD = 2 * PI * A ^ 1.5 / Sqr(M) 'Period
    End Sub

    Private Sub Command2_Click()
    End Sub

    Private Sub Form_Load()
    DrawWidth = 4

    End Sub

    Private Sub Form_Resize()
    xCenter = Me.Width / Screen.TwipsPerPixelX / 2: yCenter = Me.Height / Screen.TwipsPerPixelY / 2

    End Sub

Have something to add?