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

nikolafmf
Messages
112
Reaction score
0

Homework Statement



Solve the Kepler problem by Runge-Kutta method, where initial conditions x(0)=x0, y(0)=y0, vx(0)=vx0, vy(0)=vy0

Homework Equations



dvx/dt = -GMx/(x2+y2)3/2 (1)
dvy/dt = -GMy/(x2+y2)3/2 (2)
dx/dt = vx (3)
dy/dt = vy (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 k2, k3 etc. For example:

k1 = f(t_n, v_n)
k2 = f(t_n + 1/2h, v_n+1/2hk_1)
k3 = f(t_n+ 1/2h, v_n+1/2hk_2)

Now, if I am right, k1 = k2= k3=k4 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 k1 = k2= k3=k4? 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?
 
Last edited:
Physics news on Phys.org
You have a system of 4 first order differential equations with the vector u = (x, y, vx, vy) 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 ki is a 4-element vector and an expression like un+(h/2)k1 means "multiply the k1 vector with the scalar (h/2) and add it to the vector un".

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.
 
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.
 
nikolafmf said:
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 k1 = f(u), to the corresponding set of 4 equations, like

k1[0] = f(un)[0] = (vxn, vyn, -GMx/(xn2+yn2)(3/2), -GMy/(xn2+yn2)(3/2))[0] = vxn,
k1[1] = f(un)[1] = ... = vyn,
k1[2] = f(un)[1] = ... = -GMx/(xn2+yn2)(3/2), and
k1[3] = f(un)[3] = ... = -GMy/(xn2+yn2)(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 k1 to k4 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 )
 
Last edited by a moderator:
  • Like
Likes Simon Jug
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.
 
Hello everyone, I’m considering a point charge q that oscillates harmonically about the origin along the z-axis, e.g. $$z_{q}(t)= A\sin(wt)$$ In a strongly simplified / quasi-instantaneous approximation I ignore retardation and take the electric field at the position ##r=(x,y,z)## simply to be the “Coulomb field at the charge’s instantaneous position”: $$E(r,t)=\frac{q}{4\pi\varepsilon_{0}}\frac{r-r_{q}(t)}{||r-r_{q}(t)||^{3}}$$ with $$r_{q}(t)=(0,0,z_{q}(t))$$ (I’m aware this isn’t...
Hi, I had an exam and I completely messed up a problem. Especially one part which was necessary for the rest of the problem. Basically, I have a wormhole metric: $$(ds)^2 = -(dt)^2 + (dr)^2 + (r^2 + b^2)( (d\theta)^2 + sin^2 \theta (d\phi)^2 )$$ Where ##b=1## with an orbit only in the equatorial plane. We also know from the question that the orbit must satisfy this relationship: $$\varepsilon = \frac{1}{2} (\frac{dr}{d\tau})^2 + V_{eff}(r)$$ Ultimately, I was tasked to find the initial...
Back
Top