solving a system of ODEs using the RKF78 method - Printable Version +- HP Prime Forum (http://forum.hp-prime.de) +-- Forum: HP Prime (/forumdisplay.php?fid=1) +--- Forum: Software (/forumdisplay.php?fid=4) +--- Thread: solving a system of ODEs using the RKF78 method (/showthread.php?tid=102) solving a system of ODEs using the RKF78 method - cdeaglejr - 04-05-2016 11:08 PM 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