State-Space Method for DAEs
Consider a partitioned DAE system of the following form:
If the Jacobian matrices 
and 
are both invertible, then by the implicit function theorem, the system can effectively be expressed in the state-space form as
The invertibility requirement is effectively equivalent to the system being of index 1, so a method based on the state-space form is appropriate for index 1 systems of DAEs. The most common form is where 
and 
is the identity matrix.
The "StateSpace" time integration method is based on this representation. Given values of 
, Newton iterations are used to find 
and 
. The values 
are returned for an underlying ODE method that just requires evaluation of the right-hand side function. In between evaluations, previously computed values of the derivatives and algebraic variables are saved to initialize iterations for the next evaluation.
The ODEs resulting from the state-space form are often quite stiff and so require a Jacobian. While this can be computed using finite differences, it is more accurate and efficient to use the problem structure, since the matrices 
and 
will already have been computed. The right-hand side Jacobian can be computed by differentiating the partitioned system with respect to 
and solving for the right-hand side Jacobian gives
Matrix decompositions for 
and 
are typically saved from the Newton iterations so the application of the inverses can be computed efficiently.
The "StateSpace" method is typically slower than the default "IDA" multistep method because it has to solve, possibly iteratively, for 
and 
for each evaluation as opposed to for each time step, and there may be many evaluations per time step. There are some possible advantages that offset the speed in some cases: the iterations local to the function evaluation may be more stable, an A-stable underlying integration method can be used to improve overall stability, the method is effectively a one-step method when the underlying method is one step, and the method is implemented to allow arbitrary-precision solutions.
The Method option of the "StateSpace" method allows you to choose the underlying integration method.
Block Lower Triangular (BLT) Form
Depending on the structure of the system, it may be advantageous to order the system into block lower triangular (BLT) form. Using the BLT form results in the solution of more subproblems of smaller size. If the largest block size can be greatly reduced, using this structure can reduce the computational complexity of the Newton iterations significantly.
The BLT form is used in a number of ways for solving DAEs, including consistent initialization, setting up dummy derivatives for index reduction, and solving sparse linear systems. The description here is focused on its use in the state-space method, but the applicability extends to other areas as well.
In its simplest form, the BLT ordering is a matrix algorithm, and that is the best general way to access it in the Wolfram Language.
| SparseArray`BlockTriangularOrdering[s] | give an ordering {{r1,…,rb},{c1,…,cb}} of the rows and columns of the matrix s such that s[[rk,ck]] is the k th block along the diagonal of the reordered matrix s[[r,c]], where r=Flatten[{r1,…,rb}] and c=Flatten[{c1,…,cb}], and all other elements of the reordered matrix appear below the diagonal |
Getting the block lower triangular order.
When the matrix s is nonsingular, working down the diagonal, the variables in each successive block depend only on the ones in previous blocks. This means that the linear equations 
can be solved by working sequentially down the diagonal, solving each block subsystem in order, possibly saving a huge amount of computational effort.
It is easiest to see how this works in the context of small examples.
For this system, the blocks are all of size 1, which means that it is particularly easy to solve, since you can just pick off the variable values sequentially as can be seen from the listing of the equations in the BLT ordering. For the state space method, h={h1,h2,h3} are given so getting the derivatives h′ is trivial.
The usefulness of the BLT ordering naturally extends to nonlinear problems.
From this form, it is very clear that the Jacobian is singular if 
.
As before, it is possible to work down the diagonal, solving each subproblem in order. Once x', y, and yd1 are found (the nonlinear problems can be done with a scalar Newton solver), then only a 3×3 block needs to be solved to get x1', yd2, and λ.
Note that there are two solutions of the equation 
for 
. Continuity is maintained by using previous values to start the Newton iterations.
The state-space method can be directed to use the block lower triangular ordering in various ways using the method option "BlockLowerTriangular".
| Automatic | determine how to use the BLT ordering automatically |
| False | do not use the BLT ordering |
| "Ordering" | arrange variables and equations using the BLT ordering |
| "Solve" | solve by using the BLT ordering to solve successive subsystems |
Settings for the "BlockLowerTriangular" method option.
For this example, the timings are quite close because the system is so small. For a large system with small blocks in the BLT-ordered Jacobian, the difference can be more substantial.