# Difficulties with higher order Runge-Kutta implementations

I've been trying to implement a higher order Runge-Kutta numerical integrator in MATLAB for a bit now with ambitions of n-body trajectory fun. Unfortunately, I have not had success in reaping the rewards of such higher order methods.

I have (I believe) a good grasp of how Runge-Kutta methods work and have ample experience with lower order methods like an RK4 or Dormand Prince.

Various answers on the site have suggested RK78 and RK89. I found Fehlberg's Technical Report (NTRS ID: 19680027281) describing both, as well as C++ source code from Trick and GMAT giving the Butcher tableaus (in a C++ structure) for each method. This gives me confidence that my coefficients are correct as I copied them into my MATLAB script (and only "translated" to MATLAB).

My simple test case is an object in a 2D 185 km orbit above "Earth", that is, a point mass with no atmosphere. I used a fixed 30 second time step (just trying to get one of the orders to work!).

For both of my RK78 & RK89 implementations I get nearly identical erroneous outputs of this growing spiral-like trajectory:

Edit: my specific energy plots were not correct, I have fixed those, the problem(s) still remain with the trajectory.

RK78:

RK89:

These results occur regardless of whether I select the 7th or 8th (or 8th or 9th) order output.

Compare these to a much more reasonable result from a Dormand Price integration programmed in the same structure (trajectory plot excluded because it's expectedly benign):

Here is a code snippet of the core iteration scheme:

% k is a m x 4 matrix where m is the number of stages, 16 for RK89 (shown), 13 for RK78.
% multiplying the weights (C, m x 1, transposed) by k (m x 4) gives the output 'slopes' (1 x 4)

while 1
% Runge-Kutta Slopes
for j = 1:16 % k_j (individual k slopes)
input = zi(i,:);
for jj = 1:16 % (component of k slopes, B matrix)
input = input + B(j,jj)*k(jj,:);
end
k(j,:) = dz(input,G,M_E)'; % call physics model for j'th node
end

% Forward local extrapolation
i = i+1; % move time forward

% Break Conditions
...

% Integrate
zi(i,:) = zi(i-1,:) + dt*(C'*k); % next state
end

% force model function 'dz':
function dzdt = dz(z,G,M_E)
aG_E = G*M_E/(z(1)^2 + z(2)^2); % gravitational acceleration (m/s^2)

% State Variable Definitions:
% z1 = x, z2 = y,
% z3 = Vx, z4 = Vy

% 'slopes':
dzdt(1) = z(3);
dzdt(2) = z(4);
dzdt(3) = -aG_E*z(1)/sqrt(z(1)^2 + z(2)^2);
dzdt(4) = -aG_E*z(2)/sqrt(z(1)^2 + z(2)^2);
end


Does anyone have any clues as to what is going (so horribly) wrong?

• BTW, you should use the standard gravitational parameter rather than the mass. Horizons gives 398600.435436 km^3/s^2 for Earth. Jan 27 at 5:58
• I agree 1000% with @PM2Ring. Never use G*M. Use the standard gravitational parameter instead. Using G*M is a way to chop accuracy down to 4 or 5 decimal places. Jan 27 at 13:26
• I've had limited (being nice) success using adaptive RK techniques. There are multiple commercial vendors whose adaptive integrators are <<expletive deleted>> awful. RK techniques address the problem of solving first order differential equations. Gravitation is a second order differential equation. When one piles the position and velocity together the result is oftentimes a pile of <<expletive deleted>>. You are throwing out geometry by doing so (that the underlying equations are second order ODEs is geometry), and the result of flushing geometry down the toilet is oftentimes garbage. Jan 27 at 13:43
• @PM2Ring I also slightly disagree with myself. It is educational to roll your own RK4, or your own low order symplectic integrator. But rolling your own adaptive integrator? That is not a good idea. A much better idea is to try using something like ODEPACK (LSODE, LSODA, ...). This also involves lots of "fun", "fun" in quotes. A lot of the well-developed integration techniques have FORTRAN-style interfaces (aka "fun"). This is educational. It is not fun in any reasonable sense of the word "fun". Jan 27 at 14:34
• However, there are boatloads of numerical techniques that date back to the 1950s / 1960s that inherently have FORTRAN (not just Fortran but FORTRAN) style interfaces. They might be hiding behind f2c (I prefer the more vulgar eff2c) transformations, but they ultimately are FORTRAN-style interfaces. Learning how to deal with that ugliness can be useful. Jan 27 at 14:36

When in doubt, do what my high school physics teacher taught; use dimensional analysis and check your units. dz/dt should be an acceleration in length per time squared. GM is length cubed per time squared.

$$a_G = -GM \frac{\mathbf{r}}{r^3} = -GM \frac{\mathbf{\hat{r}}}{r^2}$$

which usually works out to something like

$$a_x = GM\frac{x}{(x^2+y^2)^{1.5}}$$

$$a_y = GM\frac{y}{(x^2+y^2)^{1.5}}$$

or acc = pos * ((pos**2).sum())**-1.5

Either that, or you've discovered a new force!

• Should that be pos**2?
– Ludo
Jan 27 at 6:59
• @Ludo yes indeed, thanks!
– uhoh
Jan 27 at 8:23
• The units are correct, also note that it is the same physics model producing the nice Dormand Prince output Jan 27 at 12:52

Foiled by my own eyes! (and MATLAB's indexing)

Lackluster, but expected. There was a bug in my code where the indices were off by one because MATLAB indexes from 1 instead of 0 I didn't look hard enough.

The comments are insightful and I offer these insights in return:

• Same test case: 30s fixed step (7th order solution), point mass 2D Earth

Where the spread (min max) of specific energies (w.r.t. Earth's surface) is on the order of 10's of micro-Joules (per kg), though I think this is stretching the limits of double precision. The orbit only drifts inwards by ~16 micrometers.

• All of the specific energy plots are incorrect (including in question), but its just a constant offset (km's used instead of meters for Earth surface potential) so I won't bother fixing it :) Feb 3 at 17:15