This software is designed for the numerical solution of DAEs.
Its main feature is to solve high index problems that can't be solved
with usual Runge-Kutta and BDF solvers. Similarly to the general DAEs
integrators developed by Campbell or Kunkel/Mehrman, it is based on
index reduction techniques. However,
some computer algebra is carried out for reducing formally the
DAE to an index one problem, so that the index is
computed and the hidden constraints are explicitely worked out.
This symbolic pre-processing reduces significantly
the amount of work for the computation of
consistent initial conditions and the numerical solution.
For ODE's and semi-explicit index one DAE's, the code coincides
with RADAU5.
What problems does it handle ?
The software handles arbitrary index general DAEs in the form
where may be rank deficient.
For the special case of quasilinear DAE's (linear w.r.t. derivatives) of type
where A is a rank deficient n-dimensional square matrix, and b, y are
n-dimensional vectors, a symbolic pre-processing allows to
compute the index and check the solvability. This also
reduces significantly the amount of work for the
numerical solver. There is no restriction for the entries of A and b,
although theoretical results are only ensured for multivariate polynomials.
The class of quasilinear problems covers most applications.
Computation of the index and constraints
For a quasilinear DAE,
this can be done with the subroutine
which returns
If either the matrix A or the vector b contains some constants,
they must be declared with the optional argument
const={C1=value1, C2=value2 ...}.
Computation of an Index One Form
Before computing consistent initial conditions and a numerical solution,
one must first compute an Index One Form, i.e
a function G describing the constraints variety.
Two subroutines are available for computing this form:
If either the DAE is not quasilinear or no symbolic
pre-processing is desired, the usual derivative array (where the
ith equation of F is differentiated numdiff[i] times) can be
used for describing the constraint variety. This requires the
knowledge of the index apriori. One should use
the subroutine:
For quasilinear DAEs, one can compute a "minimal"
Index One Form ( where G is the set of constraints ) with the
subroutine:
The index need not be known apriori.
In both cases, the output is
where dimalg is the dimension
of the constraints variety, n is the dimension of the initial DAE,
nbeqns/nbvars are the number of equations/variables of G. The function
F is the initial DAE, the functions JacF and JacG
are the jacobians of F and G (all four are optimized procedures).
For quasilinear DAE's, the subroutine MinIndexOneForm is preferable
since it will return smaller dimensional functions G/JacG. For these
two subroutines, the constants must be declared with the optional
argument const (see DAEreduction).
Computation of Consistent Initial Conditions
A set of consistent initial conditions can be computed with the subroutine
A Newton method with optional jacobian reuse and linear research is used.
Some components of the guess ICGuess can be freezed (no more than the
dimension of the algebraic variety).
Numerical solution
The numerical solution over [t0,tf] can be computed with the subroutine
The numerical solver is based on the variable stepsize
RADAU IIA (order 5) method. Tolerances rtol, atol must be
specified by the user. It returns dense output functions for
the stepsize and each component of the solution.