RADAU V

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.
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	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] .
hinit	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 ).
hmax	Maximal step size, default (when hmax is set = 0) xend - x .
nmax	This is the maximal number of allowed steps. The default value is 100000.
uround	rounding unit, default 1.0e-16.
safe	The safety factor in step size prediction, default 0.9.
facl facr	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 .
nit	The maximum number of Newton iterations for the solution of the implicit system in each step. The default value is 7.
startn	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 .
npred	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.
hess	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 ).
fnewt	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)) .
quot1 quot2	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 .
thet	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.