State-Space Method for DAEs

Consider a partitioned DAE system of the following form:

If the Jacobian matrices ((∂f)/(∂x')) and ((∂g)/(∂y)) 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 f(t,x, x', y)==x'-(Subscript[f, r])(t,x, y) and ((∂f)/(∂x')) is the identity matrix.

The "StateSpace" time integration method is based on this representation. Given values of x, Newton iterations are used to find y==G(x) and x'==F(t,x,G(x)). The values x' 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 ((∂f)/(∂x')) and ((∂g)/(∂y)) will already have been computed. The right-hand side Jacobian can be computed by differentiating the partitioned system with respect to x

and solving for the right-hand side Jacobian gives

Matrix decompositions for ((∂f)/(∂x')) and ((∂g)/(∂y)) 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 y and x' 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 (RefLink[Method, paclet:ref/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.

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 s.x==b 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 (RefLink[NDSolve, paclet:ref/NDSolve]) examples describing the change in height h and pressure p in tanks connected by pipes:

vars = {h1'[t],h2'[t],h3'[t], p1[t],p2[t],p3[t],pB[t]};eqns = {p1[t]==1.962 h1[t],p2[t]==1.962 h2[t],p3[t]==1.962 h3[t],h1'[t]==-0.19635 (p1[t]-pB[t]),h2'[t]==0.19635 (-p2[t]+pB[t]),h3'[t]==0.0981748 (-p3[t]+pB[t]),0.19635 (p1[t]-pB[t])==0.19635 (-p2[t]+pB[t])+0.0981748 (-p3[t]+pB[t])};

Get the Jacobian with respect to {h',p}:

{b, J} = CoefficientArrays[Apply[Subtract, eqns, {1}], vars];
MatrixForm[J]

Get the BLT ordering:

{br,bc} = SparseArray`BlockTriangularOrdering[J]

Show the reordered matrix:

{r,c} = Map[Flatten,{br,bc}];
MatrixForm[J[[r,c]]]

Show the reordered variables and equations:

TableForm[Transpose[{vars[[c]], eqns[[r]]}]]

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.

An index 1 system of DAEs resulting from index reduction with dummy derivatives of the constrained pendulum problem:

eqns = {x'[t]==x1[t],x1'[t] + x[t] λ[t]==0,yd2[t]+y[t] λ[t]==-1, 2 ((x1[t])^2)+2 x1'[t] x[t]+2 ((yd1[t])^2)+2 y[t] yd2[t]==0,2 x1[t] x[t]+2 y[t] yd1[t]==0,((x[t])^2)+((y[t])^2)==1};

Compute the Jacobian with respect to the highest-order derivatives appearing in the equations:

vars = {x'[t],x1'[t],y[t],yd1[t],yd2[t],λ[t]};
MatrixForm[J = D[Apply[Subtract, eqns, {1}], {vars}]]

Get the BLT ordering:

{br,bc} = SparseArray`BlockTriangularOrdering[J]

Show the reordered matrix:

{r,c} = Map[Flatten,{br,bc}];
MatrixForm[J[[r,c]]]

From this form, it is very clear that the Jacobian is singular if y(t)==0.

Show the reordered variables and equations:

TableForm[Transpose[{vars[[c]], eqns[[r]]}]]

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 (x^2)+(y^2)==1 for y. 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".

Settings for the "BlockLowerTriangular" method option.

An initial condition that will avoid the singularity in the Jacobian at y==0:

init={x[0]==0,x1[0]==5/4,y[0]==-1};

Solve over many periods using the state-space method without using BLT ordering:

AbsoluteTiming[osol = NDSolve[{eqns, init},{x,y},{t,0,10},Method->{"TimeIntegration"->{"StateSpace", "BlockLowerTriangular"->False}}, WorkingPrecision->32]]

Order equations and variables with the BLT ordering:

AbsoluteTiming[sol = NDSolve[{eqns, init},{x,y},{t,0,10},Method->{"TimeIntegration"->{"StateSpace", "BlockLowerTriangular"->"Ordering"}}, WorkingPrecision->32]]

For evaluations, solve by subsystems along the BLT diagonal:

AbsoluteTiming[bsol = NDSolve[{eqns, init},{x,y},{t,0,10},Method->{"TimeIntegration"->{"StateSpace", "BlockLowerTriangular"->"Solve"}}, WorkingPrecision->32]]

Plot the solution:

Plot[Evaluate[{x[t],y[t]} /. bsol], {t,0,10}]

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

Options of the "StateSpace" method.