Applying Runge-Kutta(RK4) method to planetary velocity

[StillmatikX]

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.

 

[D H]

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.

 

[StillmatikX]

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.

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.

 

[D H]

Let me go back to your first post.

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

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 http://en.wikipedia.org/wiki/Runge-Kutta_methods" . Note that the first section after the intro covers RK4. If you have any questions after reading that article, fire away.

 

[StillmatikX]

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 http://en.wikipedia.org/wiki/Runge-Kutta_methods" . Note that the first section after the intro covers RK4. If you have any questions after reading that article, fire away.

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

 

[D H]

You need to integrate position and velocity. I will call the position and velocity of the i^th body in your simulation ${\mathbf x}_i$ and ${\mathbf v}_i$. These quantities are bolded to indicate that they are vectors. The state, which I will call ${\mathbf y}_i$ to be consistent with the wiki article is an amalgam of the position and velocity:
$${\mathbf y}_i = \bmatrix {\mathbf x}_i \\ {\mathbf v}_i \endbmatrix$$.
The time derivative of ${\mathbf y}_i$ is then
$${\mathbf y}^{\prime}_i = \bmatrix {\mathbf v}_i \\ \sum_{j\ne i} {\mathbf a}_{i,j} \endbmatrix$$
where ${\mathbf a}_{i,j}$ is the acceleration due to gravity of the i^th body toward the j^th body.

 

[StillmatikX]

You need to integrate position and velocity. I will call the position and velocity of the i^th body in your simulation ${\mathbf x}_i$ and ${\mathbf v}_i$. These quantities are bolded to indicate that they are vectors. The state, which I will call ${\mathbf y}_i$ to be consistent with the wiki article is an amalgam of the position and velocity:
$${\mathbf y}_i = \bmatrix {\mathbf x}_i \\ {\mathbf v}_i \endbmatrix$$.
The time derivative of ${\mathbf y}_i$ is then
$${\mathbf y}^{\prime}_i = \bmatrix {\mathbf v}_i \\ \sum_{j\ne i} {\mathbf a}_{i,j} \endbmatrix$$
where ${\mathbf a}_{i,j}$ is the acceleration due to gravity of the i^th body toward the j^th body.

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.

 

[D H]

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.

 

[shalayka]

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;
}

 

[tony873004]

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.

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
    DoEvents
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#)
RKRELERROR = RKERROR / A
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()
Me.Cls
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