solving a system of ODEs using the RKF78 method

[cdeaglejr]

This HP Prime program demonstrates how to interact with a routine called rkf78 which solves a system of ordinary differential equations. The computer program is based on the Runge-Kutta-Felhberg 7(8) numerical method. This algorithm is a variable step size numerical method that allows the user to specify a truncation error tolerance. This tolerance can be used to "calibrate" the algorithm.

This demo program solves the system of Keplerian, unperturbed orbital motion. For this example, a measure of the effectiveness of this algorithm can be determined by how well the orbital closes on itself after one complete orbital period.

Here's the source code for the equations of motion implemented in this example.

EXPORT keqm (t, y)

// first order equations of orbital motion

// Keplerian motion - no perturbations

// input

// t = simulation time (seconds)
// y = state vector

// output

// ydot = integration vector

/////////////////////////////

BEGIN

LOCAL rmag, rcubed;

LOCAL ydot := [0, 0, 0, 0, 0, 0];

rmag := sqrt(y(1) * y(1) + y(2) * y(2) + y(3) * y(3));

rcubed := rmag * rmag * rmag;

// integration vector

ydot(1) := y(4);

ydot(2) := y(5);

ydot(3) := y(6);

ydot(4) := -emu * y(1) / rcubed;

ydot(5) := -emu * y(2) / rcubed;

ydot(6) := -emu * y(3) / rcubed;

return ydot;

END;

The attached PDF document describes other forms of orbital equations of motion.

The following is the output created by this demo program. A typical elliptical orbit is propagated by one orbital period with a tolerance of 1.0e-10. The differences between the initial and final position and velocity vectors is displayed at the end of the printout.

demo_rkf78

initial classical orbital elements

semimajor axis = 12000 km
eccentricity = 0.15
inclination = 51.9434659014 degrees
argument of perigee = 347.329965384 degrees
raan = 212.885279071 degrees
true anomaly = 347.74343648 degrees

initial eci state vector

eci position vector (km)

r_x = −9233.9364482
r_y = −2805.2805302
r_z = −3395.02813242

eci velocity vector (km/sec)

v_x = −0.169315078972
v_y = −4.61622930329
v_z = 4.83421581734

final eci state vector

eci position vector (km)

r_x = −9233.93644828
r_y = −2805.28052594
r_z = −3395.02813705

eci velocity vector (km/sec)

v_x = −0.16931508232
v_y = −4.61622930428
v_z = 4.83421581607

final classical orbital elements

semimajor axis = 12000.0000006 km
eccentricity = 0.150000000001
inclination = 51.9434659015 degrees
argument of perigee = 347.329965374 degrees
raan = 212.885279071 degrees
true anomaly = 347.743436457 degrees

position vector differences (meters)

delta r_x = −0.00008
delta r_y = 0.00426
delta r_z = −0.00463

velocity vector differences (meters/second)

delta v_x = −0.000003348
delta v_y = −0.00000099
delta v_z = −0.00000127