# Runge Kutta 4 for orbits and Newtonian Mechanics

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?

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.

• I am using Python and i have already used symplectic methods like the velocity verlet method, but I am currently creating a computer experiment to look at the differences between a few methods. This is why I want to add in the RK4. Also when looking at the RK4 method you have e.g. k2 = h*f(t+h/2, y + k/2) and how I should use this to input the velocity, because I then don't know what to do with the t + h/2 and the y + k/2 how would I input those? Also, I do not understand Hamiltonian Mechanics so any implementation containing this, I would not be able to understand (and how to code it). – Ian Ronk Dec 16 '18 at 19:07
• The symplectic method is Hamiltonian so if you've done that you have seen the Hamiltonian approach. I am familiar with all these methods. The RK is taking a weighted average of Euler steps. You need to evaluate the acceleration (the force) at several points within your step, this requires serial iteration. You need the i-th evaluation to get the (i+1)-th estimate and you need 4 estimates for the RK4. I would recommend looking at Numerical Recipes for the algorithm (it's spelled out step by step). Then you can code it in any language. – ggcg Dec 16 '18 at 19:12
• I have done similar experiments with ray tracing and geodesic solvers, which are nothing more than trajectories. Each method has unique idiosyncrasies. If you get your hands on the actual algorithm (a complete representation) you should have no problem implementing it exactly as written. What exactly do you hope to demonstrate by comparing them all? – ggcg Dec 16 '18 at 19:15
• I have used the Velocity Verlet integration without hamiltonian mechanics, I could show you how I did it, because maybe it is wrong. I have used the algorithm explained on wikipedia and this yielded good results, without the hamiltonian – Ian Ronk Dec 16 '18 at 19:16
• I have a two body problem to solve numerically and I know where they should end up so that I can show which method to use – Ian Ronk Dec 16 '18 at 19:17

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.

• I have already made an approximation with the velocity verlet but wanted to add runge kutta 4 for completion as I have already approximated a differential equation with it. But do you think that it is not worth it to implement it? Also I found this paper: [link] (spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf) would this be the right implementation if I decide to do create the algorithm? – Ian Ronk Dec 16 '18 at 19:01
• @IanRonk I'm not sure why you'd want to mix different methods like that. Just go with Verlet and ignore anyone telling you to use RK$n$ for such orbital simulations due to the lack of energy conservation. – Kyle Kanos Dec 16 '18 at 19:15
• And that is possible without using hamiltonian mechanics? – Ian Ronk Dec 16 '18 at 19:17
• Because I did, but I know that symplectic methods are only for hamiltonian mechanics, so is it possible to use it only for Newtonian Mechanics – Ian Ronk Dec 16 '18 at 19:18
• @IanRonk I don't know what you mean. This method works so long as you consider $\dot x=v$ and $F=m\dot v$. – Kyle Kanos Dec 16 '18 at 19:18