Applying Runge-Kutta(RK4) method to planetary velocity

In summary: I'm confused as to how to calculate it. Can someone please help me understand what y' is and how to calculate it? Thanks!Let me go back to your first post.
  • #1
StillmatikX
4
0
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.
 
Astronomy news on Phys.org
  • #2
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.
 
  • #3
D H said:
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.
 
  • #4
Let me go back to your first post.

StillmatikX said:
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.
 
Last edited by a moderator:
  • #5
D H said:
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...
 
Last edited by a moderator:
  • #6
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.
 
  • #7
D H said:
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.

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.
 
  • #8
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.
 
  • #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 by a moderator:
  • #10
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:
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
 

Related to Applying Runge-Kutta(RK4) method to planetary velocity

1. What is the Runge-Kutta method and how is it applied to planetary velocity?

The Runge-Kutta method is a numerical method used to solve ordinary differential equations. In the context of planetary motion, it can be used to approximate the velocity of a planet at a given point in time by taking into account the gravitational forces acting on the planet.

2. Why is the Runge-Kutta method preferred for calculating planetary velocity?

The Runge-Kutta method is preferred because it is a higher order method, meaning it provides more accurate results compared to simpler methods such as Euler's method. It also allows for a larger time step, making it more efficient for calculating planetary motion over longer periods of time.

3. What are the steps involved in applying the Runge-Kutta method to planetary velocity?

The steps involve first setting up the initial conditions (position and velocity) of the planet, then using the equations of motion and the gravitational forces to calculate the acceleration at that point in time. The RK4 method then uses this acceleration to estimate the velocity at the next time step, and this process is repeated to approximate the velocity at subsequent time steps.

4. Can the Runge-Kutta method be used for all planetary systems?

Yes, the Runge-Kutta method can be applied to any planetary system as long as the equations of motion and the gravitational forces are known. It is a general numerical method that can be used to approximate the motion of any object under the influence of a force.

5. Are there any limitations to using the Runge-Kutta method for planetary velocity calculations?

One limitation is that the RK4 method assumes that the acceleration remains constant during each time step, which may not always be the case in complex planetary systems. Additionally, it is only an approximation and may introduce errors over long periods of time. Other numerical methods, such as the adaptive step size RK4 method, can be used to improve accuracy in these situations.

Similar threads

Replies
65
Views
2K
Replies
9
Views
2K
  • Astronomy and Astrophysics
Replies
4
Views
1K
  • Programming and Computer Science
Replies
8
Views
2K
  • Programming and Computer Science
Replies
15
Views
2K
  • Calculus and Beyond Homework Help
Replies
14
Views
2K
  • Astronomy and Astrophysics
Replies
2
Views
2K
  • Differential Equations
Replies
6
Views
3K
  • Differential Equations
Replies
11
Views
3K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
4K
Back
Top