\% ********** Class DAEs **********
\topic dae_Class {differential-algebraic equations}
The class {\bf differential-algebraic equations} (DAEs) comprises
systems of the differential-algebraic equations
which can by written in the following form:
\par
\qc{$M(\alpha)\dot x\:=\:F(t,x,\alpha),$}
\par
where $t$ is time, $x\in {\bf R}^n$ is a vector of coordinates (phase variables),
$\alpha\in {\bf R}^m$ is a vector of numerical parameters,
$F$ is a vector-function (right-hand sides), and a (possibly singular) 
$n\times n$ matrix $M$ does not depend on $t$ and $x$.
$M$ and $F$ are assumed to be sufficiently smooth but has no other special 
properties.
\par
A system of DAEs is specified or edited
in the \jump kh_RHSWin {System specification dialog box}
which is opened by the {\bf Select|System} command from the
\jump kh_MainWin {Main window.}\par
There are following class-specific fields in the dialog box.\par
{\desc 3
{\bf Coordinates}\>is where you list the names of phase variables.\par
{\bf Parameters}\>is where you list names of parameters your system
depends on.\par
{\bf Time}\>is where you type in the name of time variable
which must be a simple name.\par
}
\par
The vector-function $F$ (right-hand sides, RHS) should be specified
using the following
notation: the value of i-th component of $F$ must be assigned
to $\bf RHS\underscore(i)$ where $i=0,1,\ldots,n-1$ (see example below).
That is the value of the first RHS
must be assigned to $\bf RHS\underscore(0)$ and so on.
\par
To specify the matrix $M=\{m_{i\:j}\}$ (left-hand sides, LHS) use the corresponding
command from the {\bf Specify} menu which puts the cursor to the
appropriate input field. The matrix is given by its non-zero coefficients
as follows. For each coefficient $m_{i\:j}\neq 0$
where $0\leq i,j\leq n-1$
the value of $m_{i\:j}$ must be
assigned to $\bf LHS\underscore(i,j)$. Here $\bf i$ is the number of the equation
and $\bf j$ is the number of j-th coordinate.
For example, if $\bf w$ is third coordinate and the coefficient
of $\bf w\prime$ in the second equation is $\bf p*p$ then
the assignment operator should look like
$\bf LHS\underscore(1,2)\:=\:p*p.$
\par
There are additional commands in the {\bf Specify} menu of the
\jump kh_RHSWin {System specification dialog box} for this class.
Use them if you chose to supply derivatives of RHS and LHS by yourself.
\par
A first order derivative of i-th RHS ($i\geq 0$) with respect to
a coordinate or parameter named $\bf u$ should be assigned to
$\bf der1(i,u)$.
Similary, a second order derivative of i-th RHS ($i\geq 0$)
with respect to names $\bf u$ and $\bf v$  should be assigned to
$\bf der2(i,u,v)$ and so on. {\bf Important note.}
{CONTENT} utilizes the property
$\frac {\partial F}{\partial u \partial v}=\frac {\partial F}{\partial v \partial u}$
to minimize memory used for derivatives. That means, in particulair,
that $\bf u$ and $\bf v$ must be in {\bf non-descendent order} (the order
is given by their appearance in the
{\bf System specification dialog box).}
\par
A first order derivative of non-zero coefficient $m_{i\:j}$
with respect to a parameter $\bf p$ should be assigned to
$\bf der1(i,j,p)$, e.g. $\bf der1(1,2,alpha)=2*alpha$ gives
the derivative with respect to $\bf alpha$
of the coefficient of the third coordinate in the second equation.
\par
The following constants are defined and may be used
in RHS, LHS, its derivatives and local functions:\par
{\desc 2
{\bf COORDINATESDIM}\>the number of components.\par
{\bf PARAMETERSDIM}\>the number of parameters.\par
}
\par
{\bf Example:} For the van der Pol system
\par
\qc{$
\system 2 rcl 1.2 1 
{\dot{\xi}_1}{=}{-\;\xi_2}
{\varepsilon\dot{\xi}_2}{=}{\;\;\xi_1-\;\xi^3_{2}+\;\xi_{2}}
$}
\par
the following specification is valid:
\par
{\li 5
\desc 7
{\bf Coordinates:} \>$x[2]$\par
{\bf Parameters:} \>$eps$\par
{\bf Time:}\> $t$\par
{\bf RHS:} \>$RHS\underscore(0)=-x[1];$\par
{} \>$RHS\underscore(1)=x[0]-x[1]*x[1]*x[1]+x[1];$\par
{\bf LHS:} \>$LHS\underscore(0,0)=1;$\par
{} \>$LHS\underscore(1,1)=eps;$\par
}
\par
If you want to specify the first-order derivatives yourself, you write
for nonzero values
\par
{
\desc 4
{\bf RHS derivatives:} \>der1(0,x[1])=-1;\par
 \>der1(1,x[0])=1;\par
 \>der1(1,x[1])=-3*x[1]*x[1]+1;\par
{\bf LHS derivatives:} \>der1(1,1,eps)=1;\par
}
\end
\% dae_Class

\topic dae_Types {types of curves}
CONTENT currently supports the continuation of the following curve:
\par
\par
\jump dae_or {Orbit}
\end
\% dae_Types 

\topic dae_or {orbit of DAEs}
Given a \jump dae_Class {system of DAEs}, 
\par
\qc{
$
M\dot{x}=F(x),\: \: x \in {\bf R}^n,
$
}
\qr {(1)}
\par
an {\bf orbit} starting at $x_0 \in {\bf R}^n$ is a curve 
$x=\varphi(t),\; -\infty\; <\; t\; < \; +\infty,\; \varphi(0)=x_0$,
such that 
\par
\qc{
$
M\dot{\varphi}(t) = F(\varphi(t)).
$
}
\par
If the matrix $M$ is singular, the initial data $x_0$ must satisfy the
algebraic part of (1).
\par
Orbits are generated numerically by \jump integr {integrator}.
\end %dae_or

\topic integr {integrator}
An \jump dae_or {orbit} starting at $t=t_0$ from
a point $x_0 \in {\bf R}^n$ of a system of DAEs
is approximated by the sequence of points $x^k,k=0,1,2,\ldots$
corresponding to time moments $t^k,k=0,1,2,\ldots$
using the following method:
\par
$\;\;\;$ \jump method_RD5 {RADAU5}\par
\par
Scalar functions $g_m(x,t)$ can be {\bf monitored} along the orbit. If such a function
changes sign between $(x^{k},t^{k})$ and $(x^{k+1},t^{k+1})$, its zero can be {\bf detected} and 
computed by bi-sections to withing given tolerance.
\par
Parameters of the integrator can be modified via the corresponding
\jump kh_GeneratorWin {generator window}. \jump start_or {Initial data} can be
set in the corresponding \jump kh_StarterWin {starter window}.

\end

\topic method_RD5 {Radau IIA method}
An \jump dae_or {orbit} starting at $t=t_0$ from
a point $x_0 \in {\bf R}^n$ of a system of ODEs
\par
\qc{
$M\dot{x}=F(t,x),\;\; x \in {\bf R}^n,$
}
\par
is approximated in the implicit order 5 {\bf Radau IIA method by Ehle} as a sequence of points $x^k,k=0,1,2,\ldots,$
corresponding to time moments $t^k,k=0,1,2,\ldots$, constructed as follows.
\par
Introduce one matrix and two vectors:
\par 
\qc{$
a=\matrix 3 ccc 2 1 
  {\frac{88-7\sqrt{6}}{360}}{\frac{296-169\sqrt{6}}{1800}}{\frac{-2+3\sqrt{6}}{225}}
  {\frac{296+169\sqrt{6}}{1800}}{\frac{88+7\sqrt{6}}{360}}{\frac{-2-3\sqrt{6}}{225}} 
  {\frac{16-\sqrt{6}}{36}}{\frac{16+\sqrt{6}}{36}}{\frac{1}{9}},\;\;\;
b=\matrix 3 c 2 1 {\frac{16-\sqrt{6}}{36}}{\frac{16+\sqrt{6}}{36}}{\frac{1}{9}},\;\;\;		  
c=\matrix 3 c 2 1 {\frac{4-\sqrt{6}}{10}}{\frac{4+\sqrt{6}}{10}}{1}.		  
$}
\par
Denote by $\omega_{ij}$ the elements of the inverse matrix $a^{-1}$.
\par
{\bf Set}: $x^0=x_0,\; t^0=t_0,\; h^0=H_0$
\par
{\bf Iterate for} $k=0,1,2,\ldots$\par
AGAIN:
{\li 2
\par
{\bf If} $h^k < H_{min}$ {\bf or} $|t^k-t_0|>T_{max}$ {\bf then Stop}\par
{\li 4
  \par
  {\bf Iterate for} $i=1,2,3:\;$
  \par
  {\bf Solve by Newton iterations for $X_i\;:$} 
  $M(X_i-x^k)=h^k\Sum {j=1} {3} {a_{ij} F(t^k+c_{j} h^k,X_j)},$
  \par
  {\bf If} no convergence {\bf then}  $h^{k}:=\frac{h^k}{2}$ {\bf goto} AGAIN 
  \par
  $x^{k+1}=\brackets() { 1-\Sum {i,j=1} {3} {b_i \omega_{ij}} } x^k+ 
  \Sum {i,j=1} {3} {b_i \omega_{ij} X_j},$
}
\par
{\bf Compute local error and new step size} $h^{k+1}$
\par
$t^{k+1}=t^k + h^k$
}
\par
Exact step-control mechanism is described in Section IV.8 of the book
E. Hairer and G. Wanner "Solving Ordinary Differential 
Equations II. Stiff and Differential-Algebraic Problems". 
Springer Series in Computational Mathematics 14, 
Springer-Verlag 1991. 
\par
The \jump integrpar_RD5 {parameters of the integrator} can be modified via the corresponding
\jump kh_GeneratorWin {generator window}. \jump start_or {Initial data} can be
set in the corresponding \jump kh_StarterWin {starter window}.

\end

\topic integrpar_RD5 {integrator parameters (RADAU5)}
The \jump method_RD5 {RADAU5 integrator} has the following parameters:
\par
\par
\desc 18
{\bf Integration data}\>\par
{\li 1 Interval} \> Integration interval length $T_{max}$\par
{\li 1 InitStepsize} \>  Initial stepsize $H_0$\par
{\li 1 MinStepsize} \> Minimal stepsize $H_{min}$\par
{\li 1 MaxStepsize} \> Maximal stepsize $H_{max}$\par
{\li 1 Tolerance} \> Absolute step tolerance $\varepsilon_a$\par
{\li 1 RelTolerance} \> Relative step tolerance $\varepsilon_r$\par
{\li 1 MaxNewtonIter} \> Maximum number of Newton iterations $M_{Newt}$\par
{\li 1 Index1} \> First Index1 phase varaiables have index 1\par
{\li 1 Index2} \> Second Index2 phase varaiables have index 2\par
{\li 1 Index3} \> Third Index3 phase varaiables have index 3\par
{\li 1 mljac} \> Lowest diagonal of the RHS Jacobian matrix $F_x$ that is not zero\par
{\li 1 mujac} \> Upper diagonal of the RHS Jacobian matrix $F_x$ that is not zero\par
{\li 1 mlmas} \> Lowest diagonal of the LHS matrix $M$ that is not zero\par
{\li 1 mumas} \> Upper diagonal of the LHS matrix $M$ that is not zero\par
{\bf User function zero}\>\par
{\li 1 UserFuncTolerance} \> Accuracy of zero location of user functions\par
{\li 1 MaxIterFunc} \> Maximal number of bi-sections for zero location\par
\par
The parameters can be modified via the \jump kh_GeneratorWin {generator window}.
\end

\topic start_or {starter parameters for orbit}
The \jump kh_Starter {starter} which sets necessary data to \jump kh_Generator {generate} 
the \jump dae_or {orbit} has the following {\bf parameters}:
\par
{
\desc 8
{\bf Initial point}\>\par
{\li 1 $t$} \> A name of the time variable and its initial value.\par
{\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par
{\li 1 $\alpha$} \> A list of parameter names and their values at the initial point.\par
{\bf User defined functions}\>\par
{\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}).\par
{\bf Set initial point}\> \par
{\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to
 phase variables and parameters. When it is specified and
 activated, it is called by the started before computing the first point.  \par
}
\par
If the LHS matrix $M$ of the \jump dae_Class {DAEs} is singular, the initial 
phase coordinates $x$ must satisfy the algebraic part of system. 
The parameters can be modified via the \jump kh_StarterWin {starter window}.

\end \% start_or

\topic point_or {point}
You can compute 
\par
{\li 2
\desc 1
\jump dae_or {an orbit}\>
}
\par
starting from a {\bf point} in the phase-parameter space. When the LHS matrix $M$
is singular, the point should be rather close to the manifold defined by the 
algebraic part of the \jump dae_Class {DAEs}.

\end