Numerical Solution of Differential-Algebraic Equations


In general, a system of ordinary differential equations (ODEs) can be expressed in the normal form,

The derivatives of the dependent variables are expressed explicitly in terms of the independent variable and the dependent variables . As long as the function has sufficient continuity, a unique solution can always be found for an initial value problem where the values of the dependent variables are given at a specific value of the independent variable.

With differential-algebraic equations (DAEs), the derivatives are not, in general, expressed explicitly. In fact, derivatives of some of the dependent variables typically do not appear in the equations. The general form of a system of DAEs is

where the Jacobian with respect to , may be singular.

A system of DAEs can be converted to a system of ODEs by differentiating it with respect to the independent variable . The index of a DAE is effectively the number of times you need to differentiate the DAEs to get a system of ODEs. Even though the differentiation is possible, it is not generally used as a computational technique because properties of the original DAEs are often lost in numerical simulations of the differentiated equations.

Thus, numerical methods for DAEs are designed to work with the general form of a system of DAEs. The methods in NDSolve are designed to generally solve index-1 DAEs, but may work for higher-index problems as well.

This tutorial will show numerous examples that illustrate some of the differences between solving DAEs and ODEs.

This loads packages that will be used in the examples and turns off a message.
Click for copyable input

The specification of initial conditions is quite different for DAEs than for ODEs. For ODEs, as already mentioned, a set of initial conditions uniquely determines a solution. For DAEs, the situation is not nearly so simple; it may even be difficult to find initial conditions that satisfy the equations at all. To better understand this issue, consider the following example [AP98].

Here is a system of DAEs with three equations, but only one differential term.
Click for copyable input

The initial conditions are clearly not free; the second equation requires that be either 0 or 1.

This solves the system of DAEs starting with a specified initial condition for the derivative of .
Click for copyable input

To get this solution, NDSolve first searches for initial conditions that satisfy the equations, using a combination of Solve and a procedure much like FindRoot. Once consistent initial conditions are found, the DAE is solved using the IDA method.

This shows the initial conditions found by NDSolve.
Click for copyable input
This shows a plot of the solution. The solution is obscured by the solution , which has the same constant value of 1.
Click for copyable input

However, there may not be a solution from all initial conditions that satisfies the equations.

This tries to find a solution with starting from steady state with derivative 0.
Click for copyable input
This shows the initial conditions found by NDSolve.
Click for copyable input

If you look at the equations with set to 1, you can see why it is not possible to advance beyond .

Substitute into the equations.
Click for copyable input

The middle equation effectively drops out. If you differentiate the last equation with , you get the condition , but then the first equation is inconsistent with the value of in the initial conditions.

It turns out that the only solution with is , and along this solution, the system has index 2.

The other set of solutions for the problem is when . You can find these by specifying that as an initial condition.

This finds a solution with . It is also necessary to specify a value for because it is a differential variable.
Click for copyable input
This shows a plot of the nonzero components of the solution.
Click for copyable input

In general, you must specify initial conditions for the differential variables because typically there is a parametrized general solution. For this problem with , the general solution is , so it is necessary to give to determine the solution.

NDSolve cannot always find initial conditions consistent with the equations because sometimes this is a difficult problem. "Often the most difficult part of solving a DAE system in applications is to determine a consistent set of initial conditions with which to start the computation" [BCP89].

NDSolve fails to find a consistent initial condition.
Click for copyable input

If NDSolve fails to find consistent initial conditions, you can use FindRoot with a good starting value or some other procedure to obtain consistent initial conditions and supply them. If you know values close to a good starting guess, NDSolve uses these values to start its search, which may help. You may specify values of the dependent variables and their derivatives.

With index-1 systems of DAEs, it is often possible to differentiate and use an ODE solver to get the solution.

Here is the Robertson chemical kinetics problem. Because of the large and small rate constants, the problem is quite stiff.
Click for copyable input
This solves the Robertson kinetics problem as an ODE by differentiating the balance equation.
Click for copyable input

The stiffness of the problem is supported by and having their main variation on two completely different time scales.

This shows the solutions and .
Click for copyable input
This solves the Robertson kinetics problem as a DAE.
Click for copyable input

The solutions for a given component will appear quite close, but comparing the chemical balance constraint shows a difference between them.

Here is a graph of the error in the balance equation for the ODE and DAE solutions, shown in black and blue, respectively. A log-log scale is used because of the large variation in and the magnitude of the error.
Click for copyable input

In this case, both solutions satisfied the balance equations well beyond expected tolerances. Note that even though the error in the balance equation was greater at some points for the DAE solution, over the long term, the DAE solution is brought back to better satisfy the constraint once the range of quick variation is passed.

You may want to solve some DAEs of the form

such that the solution of the differential equation is required to satisfy a particular constraint. NDSolve cannot handle such DAEs directly because the index is too high and NDSolve expects the number of equations to be the same as the number of dependent variables. NDSolve does, however, have a Projection method that will often solve the problem.

A very simple example of such a constrained system is a nonlinear oscillator modeling the motion of a pendulum.

This defines the equation, invariant constraint, and starting condition for a simulation of the motion of a pendulum.
Click for copyable input

Note that the differential equation is effectively the derivative of the invariant, so one way to solve the equation is to use the invariant.

This solves for the motion of a pendulum using the invariant equation. The option tells NDSolve not to symbolically solve the quadratic equation for , but instead to solve the system as a DAE.
Click for copyable input

However, this solution may not be quite what you expect: the invariant equation has the solution constant when it starts with . In fact it does not have unique solutions from this starting point. This is because if you do actually solve for , the function does not satisfy the continuity requirements for uniqueness.

This solves for the motion of a pendulum using only the differential equation. The method is used because it can also be a submethod of the projection method.
Click for copyable input
This shows the solution plotted over the last several periods.
Click for copyable input
This shows a plot of the invariant at the ends of the time steps NDSolve took.
Click for copyable input

The error in the invariant is not large, but it does show a steady and consistent drift. Eventually, it could be large enough to affect the fidelity of the solution.

This solves for the motion of the pendulum, constraining the motion at each step to lie on the invariant.
Click for copyable input
This shows a plot of the invariant at the ends of the time steps NDSolve took with the projection method.
Click for copyable input