U N D E R    D E V E L O P M E N T

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 invertiblility 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 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 method is typically slower than the default 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 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 Mathematica.

SparseArray`BlockTriangularOrdering[s]give an ordering of the rows and columns of the matrix s such that is the k^(th) block along the diagonal of the reordered matrix , 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.

A linear system of DAEs from the NDSolve examples describing the change in height and pressure in tanks connected by pipes.
Click for copyable input
Get the Jacobian with respect to .
Click for copyable input
Get the BLT ordering.
Click for copyable input
Show the reordered matrix.
Click for copyable input
Show the reordered variables and equations.
Click for copyable input

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, are given so getting the derivatives is trivial.

The usefulness of the BLT ordering naturally extends to nonlinear problems.

An index 1 system of DAEs resulting from index reduction with dummy derivatives of the constrained pendulum problem.
Click for copyable input
Compute the Jacobian with respect to the highest-order derivatives appearing in the equations.
Click for copyable input
Get the BLT ordering.
Click for copyable input
Show the reordered matrix.
Click for copyable input

From this form, it is very clear that the Jacobian is singular if .

Show the reordered variables and equations.
Click for copyable input

As before, it is possible to work down the diagonal, solving each subproblem in order. Once , , and 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 , , 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 .

Automaticdetermine how to use the BLT ordering automatically
Falsedo 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 method option.

An initial condition that will avoid the singularity in the Jacobian at .
Click for copyable input
Solve over many periods using the state-space method without using BLT ordering.
Click for copyable input
Order equations and variables with the BLT ordering.
Click for copyable input
For evaluations, solve by subsystems along the BLT diagonal.
Click for copyable input
Plot the solution.
Click for copyable input

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.

Option Summary

option name
default value
"BlockLowerTriangular"Automaticspecify use of the BLT ordering
MethodAutomaticspecify the underlying ODE time integration method

Options of the method.

New to Mathematica? Find your learning path »
Have a question? Ask support »