Runge Kutta 4 for orbits and Newtonian Mechanics

Question

How can I use the Runge Kutta 4 method to solve orbits of Newtonian Mechanics, with position vector $\mathbf x$, velocity vector $\mathbf v$ and acceleration vector $\mathbf a$? Do I still have to calculate the acceleration by combining the forces?

Answer 1

You need a lot. To answer your question directly, yes you need to calculate the "force" vector that is driving a change of state of the system.

1) write down the second order equations of motion in whatever coordinate you want, Cartesian, spherical, etc.

2) define a new set of variables px = dx/xt, etc.

3) the last step will reduce the order of the ODE system from second to first. Most ODE solvers only work on first order derivative but that is not an issue as we've been applying reduction of order for at least a century. The cost is that your system has twice as many variables to solve for. It is similar to a Hamiltonian formulation of the problem.

4) At this point we cannot help you unless you state whether you are writing your own code or using a canned ODE library like that in MATLAB, or ODEINT, etc. But basically you will have a first order vector ODE of the form dY/dt = F(Y), where Y = [x, y, z, px, py, pz] or whatever, and F(Y) is given by the forces present in the equations, including terms that arise from transforming to spherical coordinates if you decide to do that, e.g. (d(theta)/dt)^2, etc.

5) Now you can use you F(Y) to evaluate the pieces of an RK step and try it.

Word of caution. I don't know why people are so in love with RK when there are many other solvers out there. RK is explicit and implicit solvers can be more stable. If you have never tried this before RK4 is a good starting point to learn. If you have not given any thought to the step size you will need to. RK4 can drift off and you can either (1) use a very very small fixed step size and cross your fingers, or (2) implement a step size control that checks an error estimate then shrinks the step until a good answer is obtained. RK4 will not give you an error estimate. However and embedded RK5(4) will since you are solving a 5th and 4th order RK at the same time. Other possibilities are checking energy conditions, and back propagating with a requirement that you recover the initial conditions to within a given error tolerance.

Numerical methods for ODE are a lot of fun once you get into them. Good luck. If you can elaborate on some of the questions I've posed to you I (and the rest of us) can help you further.

Answer 2

You could do it, but you don't want to because RK4 doesn't conserve energy. You need to look into symplectic integrators, such as Velocity Verlet.

I discuss how such algorithm can be implemented in this post of mine, among others. It would be easiest to store the positions, velocities and accelerations in a vectors/lists/arrays: $$\mathbf x=\left[x,\,y,\,z\right]$$ and then operate on those, depending of course how easy it is to work with arrays in the language you're using.

If you do go the vector route, the acceleration/force function should return the array so that you can use it directly when updating the velocity and position vectors.