From BioUML platform
Jump to: navigation, search

This algorithm is ported for Java from C++ code (last updated: Nov 15, 2002), written by Blake Ashby on the base of code RADAU5 originally written in FORTRAN (version of July 9, 1996, latest small correction: January 18, 2002) by E. Hairer and G. Wanner (Universite de Geneve, Dept. de Mathematiques Ch-1211 Geneve 24, Switzerland) and which is in details described in their book[1].
This code computes the numerical solution of a differential system of first order ordinary differential equations y'=f(x,y).

The method used is an implicit Runge-Kutta method (RADAU IIA) of order 5 with step size control and continuous output. (See Section IV.8 of [1]). The method is designed to solve stiff systems, for non-stiff systems it is rationally to use simpler methods (requiring less computational load), e.g. Dormand-Prince.

User-provided solver parameters Description


itoler = 0: Both rtoler and atoler are scalars. The code keeps, roughly, the local error of y[i] below rtoler*fabs(y[i])+atoler .
itoler = 1: Both rtoler and atoler are vectors. The code keeps the local error of y[i] below rtoler[i]*fabs(y[i])+atoler[i] .


Initial step size guess;
For stiff equations with initial transient, h = 1.0/(norm of f') , usually 1.0e-3 or 1.0e-5, is good. This choice is not very important, the step size is quickly adapted. (if h = 0.0 on input, the code sets h = 1.0e-6 ).


Maximal step size, default (when hmax is set = 0) xend - x .


This is the maximal number of allowed steps. The default value is 100000.


rounding unit, default 1.0e-16.


The safety factor in step size prediction, default 0.9.


Parameters for step size selection the new step size is chosen subject to the restriction 1/facl <= hnew/hold <= 1/facr .
Default values: facl = 5.0 , facr = 1.0/8.0 .


The maximum number of Newton iterations for the solution of the implicit system in each step. The default value is 7.


If startn == 0 the extrapolated collocation solution is taken as starting value for Newton's method. If startn != 0 starting values are used.
The latter is recommended if Newton's method has difficulties with convergence. (This is the case when nstep is larger than naccpt + nrejct ; see output parameters). Default is startn = false .


Switch for step size strategy:

If npred == 1 : mod. predictive controller (Gustafsson)
If npred == 2 : classical step size control.

The default value (for npred = 0 ) is npred = 1 .
The choice npred = 1 seems to produce safer results. For simple problems, the choice npred == 2 produces often slightly faster runs.


If hess != 0 , the code transforms the Jacobian matrix to Hessenberg form. This is particularly advantageous for large systems with full Jacobian. It does not work for banded Jacobian ( mljac < n ) and not for implicit systems ( imas = 1 ).


Stopping criterion for Newton's method, usually chosen < 1. Smaller values of fnewt make the code slower, but safer.
Default min(0.03, sqrt(rtoler)) .


If quot1 < hnew/hold < quot2 , then the step size is not changed. This saves, together with a large thet, lu-decompositions and computing time for large systems.
For small systems one may have quot1 = 1.0 , quot2 = 1.2 , for large full systems quot1 = 0.99 , quot2 = 2.0 might be good. Defaults quot1 = 1.0 , quot2 = 1.2 .


Decides whether the Jacobian should be recomputed. Increase thet, to 0.1 say, when Jacobian evaluations are costly. for small systems thet should be smaller (0.001, say). Negative thet forces the code to compute the Jacobian after every accepted step. Default 0.001.


  1. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Springer Series in Computational Mathematics 14, Springer-Verlag 1991, Second Edition 1996.
Personal tools

BioUML platform
Analysis & Workflows
Collaborative research
Virtual biology