The author is very grateful to Mr Guillaume Lion at 06130 GRASSE, France who improved all of the MATLAB codes and corrected some bugs in implicit6.m and implicit8.m He also made corrections to the FORTRAN codes. Finally all of the integrators written in F90 and MATLAB were modified so that they could support backward integration. These substantial changes have greatly improved the codes.
Four codes 4irkfinal, 6irkfinal, 8irkfinal and10irkfinal are given for the numerical integration of reversible systems of ODEs using a fixed stepsize of integration and are written in Fortran 77. The codes use a special class of Runge-Kutta-Nystrom formulae and are of orders 4,6, 8 and 10 respectively.
The four problems to be run are given in the code problem.f
It should be clear how to run these codes on reading the driver.
All of these programs are adapted from programs which run Gauss Runge-Kutta codes using a fixed stepsize of integration, which were written by Ernst Hairer. We acknowledge the help received from him particularly in implementing compensated summation. The Hairer codes are given in his recent book with Lubich and Wanner while the Runge-Kutta- Nystrom codes are given here. Comparisons with these two approaches on Kepler's equations and also on the equations of the outer solar system (see book by Hairer at el.) are given here. The codes also run the harmonic equation and pendulum problem. The paper has now been updated and gives algorithms of order up to 10.
This Fortran 90 code allows the integration of reversible ODEs using both fixed and variable stepsizes. The available integration formulae are the order 4, 6 and 8 schemes from above; and an explicit 4th order scheme. The code for this is given here in the order it should be compiled: commonfunctions.f90, userfunctions.f90, explicit4.f90, rkn4.f90, rkn6.f90, rkn8.f90, finalrkn.f90. To enter your own problems and parameters you only need to edit userfunctions.f90 and finalrkn.f90.
To help with the easy compilation of this program we also provide both a Unix make file and a Microsoft visual studio project file.
To solve Kepler"s equations with your own eccentricity, integration period, choice of method and tolerance/stepsize just compile and run the program and you will be asked to input the relevant data.
For larger eccentricities it is much better to use quadruple precision (Real16). To do this just instruct your compiler to assume all reals to be real16, and edit the value of uround at the top of userfunctions.f90 to about 1.e-30.
To hard code you own parameters simply edit mysettings in userfunctions.f90 and perhaps finalrkn.f90 if you wish to set up batch tests.
To input your own problem use myF to input the second differential, myICs to input the initial conditions, and myerror to input you preferred error measurement.
This Matlab code is exactly the same as the Fortran 90 code described above. The files required are correct_h.m, correct_hs.m, explicit4.m, implicit4.m, implicit6.m, implicit8.m, myerror.m, myF.m, myICs.m, mySettings.m, predict_h.m and variablestep.m.
To solve Kepler"s equations with your own eccentricity, integration period, choice of method and tolerance/stepsize put all the above matlab files in your work directory, input your parameters in to mySettings.m and then call the function variablestep().
To input your own problem use myF.m to input the second differential, myICs.m to input the initial conditions, and myerror.m to input you preferred error measurement.
You can use these codes to reproduce the results from our latest paper.
1832 people have visited this page since October 2006.