martedì 12 agosto 2014

Status of the code: bugfixes and new issues


ODESET and ODEGET


  • odeset and odeget functions have been slightly modified to be compliant with MATLAB. Each MATLAB option is present and all the options are tested. The coding style has been adapted to the GNU-Octave standard.
  • ode_struct_value_check: this function has been introduced by Roberto in addition to odepkg_structue_check. The relation between the two functions has to be clarified: in particular it is necessary to understand if it is really necessary to have two different functions or one is sufficient.


CHANGES TO THE STEPPERS


  • The runge_kutta_78 stepper has been implemented.
  • Two 4th order steppers have been implemented: runge_kutta_45_dopri (Dormand-Prince coefficients) and runge_kutta_45_fehlberg (Fehlberg coefficients).


CHANGES TO THE SOLVERS


  • ode78 solver has been updated to the new structure. It now exploits the runge_kutta_78 stepper.
  • A series of tests has been added to each solver to check all the functionalities and the all options. This has made me possible to detect some bugs that have been corrected. In particular the adaptive timestep evaluation had some issues that lead to the use of too short timesteps. This has been corrected and now the algorithm proposed in [1] is used.
  • Furthermore the current implementation uses linear interpolation to evaluate the solution at the user specified times. This leads to a considerable loss in accuracy and is not consistent with MATLAB (which guarantees the same order of accuracy of the scheme used). In [1] different methods are proposed for the dense output: these will be used as a reference for the implementation of a better solution.
  • In the released version of odepkg some of the solvers perform local extrapolation, that is the higher-order estimate is chosen as the solution. With the new stepper structure, as it is now, this choice effects all the solvers. It have to be decided whether to perform it or not (MATLAB doesn't seem to do it, thus I suggest to avoid it).
  • MATLAB implementation of ode45 uses the Dormand-Prince (DP) coefficients. In the released version of odepkg there exits two solvers: ode45 that uses the Fehlberg coefficients and ode54 that uses the DP coefficients. To be consistent with MATLAB, ode45 now uses the DP method. This makes the runge_kutta_45_fehlberg stepper and the ode54 solver, as it is now, redundant. Either their elimination or a change of the solver might be considered. However one of the advantages of DP coefficients is the possibility to reuse the last function evaluation at a given step as the first evaluation of the subsequent one. This is not easily done with the stepper structure introduced by Roberto.


CHANGES TO THE OPTIONS


  • InitialStep option has been modified to be MATLAB compatible (it must be a positive scalar).
  • RelTol defalut value has been changed to 1e-3 instead of 1e-6 to be MATLAB compatible.
  • MaxStep option has been implemented.
  • NormControl option has been implemented.


TODO

In addition to the general plan, a couple of issues need to be addressed:

  • Clarify the relation between ode_struct_value_check and odepkg_structue_check.
  • Decide if local extrapolation has to be used or not. My opinion (and the current implementation) is to avoid it to be compliant to what MATLAB seems to be doing.
  • Solve the dense output problem in a way that guarantees the consistency with MATLAB.
  • Consider if it's possible to reduce the number of function evaluation for the Dormand-Prince stepper (ode45) and the Bogacki-Shampine stepper (ode23) exploiting the FSAL property (first same as last).
  • Decide if in the future releases of odepkg ode54 has to be removed or maybe changed to become a Fehlberg solver.



[1] E. Hairer, S.P. N{\o}rsett, G. Wanner, Solving Ordinary Differential Equations, 1993, Springer.

Nessun commento:

Posta un commento