Difficulties with higher order Runge-Kutta implementations

Question

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?

Answer 1

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!

And here's how to add $J_2$ as well.

Answer 2

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.