How to Solve the Kepler Problem Using the Runge-Kutta Method?

[nikolafmf]

Homework Statement

Solve the Kepler problem by Runge-Kutta method, where initial conditions x(0)=x₀, y(0)=y₀, v_x(0)=v_x0, v_y(0)=v_y0

Homework Equations

dv_x/dt = -GMx/(x²+y²)^3/2 (1)
dv_y/dt = -GMy/(x²+y²)^3/2 (2)
dx/dt = v_x (3)
dy/dt = v_y (4)

The Attempt at a Solution

First of all, Runge-Kutta solves equations which have this form:

dv/dt = f(t,v).

But equations (1) and (2) are of the form dv/dt = f(x(t), y(t)). They don't have v on the right side and t is implicit. So, problem arises when I compute k₂, k₃ etc. For example:

k₁ = f(t_n, v_n)
k₂ = f(t_n + 1/2h, v_n+1/2hk_1)
k₃ = f(t_n+ 1/2h, v_n+1/2hk_2)

Now, if I am right, k₁ = k₂= k₃=k₄ in my case, because I have not v on the right side of the equation. If that is true, then Runge-Kutta reduces to ordinary Euler method and brings nothing new. My question is, are really k₁ = k₂= k₃=k₄? And if they are not equal, where is my mistake?

Or, I can ask otherwise, how can I implement the Runge-Kutta method to Kepler problem?

 

[Filip Larsen]

You have a system of 4 first order differential equations with the vector u = (x, y, v_x, v_y) being the independent state, written in short as

du/dt = f(u,t)

and it is this system that your RK method solves. This means that the equations defining the particular RK-method you are using (RK4 in your case) must be considered vector equations, that is, in your case each k_i is a 4-element vector and an expression like u_n+(h/2)k₁ means "multiply the k₁ vector with the scalar (h/2) and add it to the vector u_n".

When you write code for this you may of course choose to write out the vector operations explicitly as scalar operations on four different variables if you do not have convenient access to a vector library.

 

[nikolafmf]

Filip, thank you very much for your response. But I still can't see that du/dt=f(u,t).
If u=(x,y,vx,vy), then du/dt = (vx, vy, -GMx/(x^2+y^2)^(3/2), -GMx/(x^2+y^2)^(3/2)), right? Then what function of u is it? I see that it is a function of the components of u, but not of u itself. And, t doesn't even appear here (at least not explicitely).
Again, thank you for your time.

 

[Filip Larsen]

If u=(x,y,vx,vy), then du/dt = (vx, vy, -GMx/(x^2+y^2)^(3/2), -GMx/(x^2+y^2)^(3/2)), right?

That is correct (well, the last component should be proportional with y and not x, but I assume that is just a typo).

Then what function of u is it? I see that it is a function of the components of u, but not of u itself. And, t doesn't even appear here (at least not explicitely).

The vector u is just a convenient notation for the complete state - a notion used by the RK scheme you mentioned part of. If you are unsure about what a vector is you probably want to refer to a textbook for some background (see [1] for a brief overview). You are free to translate each vector equation, like k₁ = f(u), to the corresponding set of 4 equations, like

k₁[0] = f(u_n)[0] = (vx_n, vy_n, -GMx/(x_n²+y_n²)^(3/2), -GMy/(x_n²+y_n²)^(3/2))[0] = vx_n,
k₁[1] = f(u_n)[1] = ... = vy_n,
k₁[2] = f(u_n)[1] = ... = -GMx/(x_n²+y_n²)^(3/2), and
k₁[3] = f(u_n)[3] = ... = -GMy/(x_n²+y_n²)^(3/2)

where I have used the code-like notation "vector" to mean "the i-th component of vector" for i = 0, 1, 2, 3. So, for each of k₁ to k₄ you can write up such 4 scalar expressions.

You are correct that time t does not explicitly appear in the field, that is, f(u, t) = f(u). This is quite common in dynamical systems and a system with this characteristic is called an autonomous system (see for instance [2]).[1] http://en.wikipedia.org/wiki/Euclidean_vector
[2] http://en.wikipedia.org/wiki/Autonomous_system_(mathematics )

 

[jela]

I must solving a Kepler's problem (for elliptical orbit, for Earth) with Runge-Kutta 4th methods in MATLAB, but I don't know what is my differential equation and how I can solving it, what I must put it in differential equation for initial value.

For examples when I change a steps I must get where will be a Earth on his orbit.

And for last value I must get the correct value for period of Earth, 365.24 days, and position of Earth return in started position.