Introduction

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 transient 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. For example, the following equation

    

does not explicitly contain any derivatives of . Such variables are often referred to as algebraic variables.

The general form of a system of DAEs is

Solving systems of DAEs often involves many steps. The flow chart shown below indicates the general process associated with solving DAEs in NDSolve.

Flow chart of steps involved in solving DAE systems in NDSolve.

Generally, 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 the number of times needed to differentiate the DAEs to get a system of ODEs.

The DAE solver methods built into NDSolve work with index-1 systems, so for higher-index systems an index reduction may be necessary to get a solution. NDSolve can be instructed to perform that index reduction. When a system is found to have an index greater than 1, NDSolve generates a message and number of steps need to be taken in order to solve the DAE.

As a first step for high-index DAEs, the index of the system needs to be reduced. The process of differentiating to reduce the index, referred to as index reduction, can be done in NDSolve with a symbolic process.

The process of index reduction leads to a new equivalent system. In one approach, this new system is restructured by substituting new (dummy) variables in place of the differentiated variables. This leads to an expanded system that then can be uniquely solved. Another approach for restructuring the system involves differentiating the original system number of times until the differentials for all the variables are accounted for. The preceding differentiated systems are treated as invariants.

In order to solve the new index-reduced system, a consistent set of initial conditions must be found. A system of ODEs in normal form, , can always be initialized by giving values for at the starting time. However, for DAEs it may be difficult to find initial conditions that satisfy the residual equation ; this amounts to solving a nonlinear algebraic system where only some variables may be independently specified. This will be discussed in further detail in a a later section. Furthermore, the initialization needs to be consistent. This means that the derivatives of the residual equations also need to be satisfied. Commonly, higher-index systems are harder to initialize. In this case NDSolve cannot, in general, see how the components of interact and is thus not able do an automatic index reduction. NDSolve, however, has a number of different methods, accessible via options, to perform index reduction and find a consistent initialization of DAEs.

The remainder of this document is structured in the following manner: first, some terminology around indices of DAEs is introduced. Following that, the next section treats solving of index-1 DAEs. Next, methods for index reduction of higher-order index DAEs are presented. Consistent initialization of DAEs is treated subsequently. In addition to the examples in the sections mentioned so far, an entire last section with additional examples is given.

Index of a DAE

For a general DAE system , the index of the system refers to the minimum number of differentiations of part or all the equations in the system that would be required to solve for uniquely in terms of and . Note that there are many slightly different concepts of index in the literature; for the purpose of this document the preceding concept will be used.

Consider the following DAE

    

In order to get a differential equation for , the first equation must be differentiated once. This leads to the new system

    

The differentiated system now contains , which must be accounted for. To get an equation for , the second equation must now be differentiated three times. This leads to

    

Since part of the system had to be differentiated three times to get a derivative equation for , the index of the DAE is 3. The final differentiated system is said to be an index-0 system. To see ODEs that are obtained by the differentiations, you can perform appropriate substitutions, in this case subtracting the second equation from the first. The resulting equation can be written as

    

Similarly, the index-1 system can be found by differentiating the second equation twice, resulting in the following system

    

It is typical to reduce the index of the system to 1, since the underlying integration routines can handle them efficiently.

The index of a DAE may not always be obvious. Consider a system that may exhibit index-1 or index-2 behavior, depending on what the initial conditions are.

For this DAE, the initial conditions are not free; the second equation requires that be either 0 or 1.

If , then the two remaining equations make an index-1 system, since differentiation of the third equation gives a system of two ODEs

    

While the initial condition for , the initial condition for can be varied.

On the other hand, if , then the two remaining equations make an index-2 system. Differentiating the third equation gives , and substituting into the first equation gives , which needs to be differentiated to get an ODE

    

There are also no additional degrees of freedom in the initial conditions, since and .

Both the number of initial conditions needed and the index of a DAE may influence the actual solution being found.

Treatment of Differential Equations

In order to solve a differential equation, NDSolve converts the user-specified system into one of three forms. This step is quite important because depending on how the system is constructed, different integration methods are chosen.

You may specify in which form to represent the equations by using the option Method->{"EquationSimplification"->simplification}. The following simplification methods can be specified for the "EquationSimplification" option.

Equation simplification methods.

Consider the system

    

By default, an attempt is made to solve for the derivatives explicitly in terms of the independent and dependent variables and . For efficiency purposes, NDSolve first attempts to try to find the explicit system by using LinearSolve. If that fails, however, then the system is solved symbolically using Solve.

Using the "Solve" method is the same as the default when a symbolic solution is found.

When the system is put in a residual form, the derivatives are represented by a new set of variables that are generated to be unique symbols. You can get the correspondence between these working variables and the specified variables using the "Variables" and "WorkingVariables" properties of NDSolveStateData.

The process of generating an explicit system of ODEs may sometimes become expensive due to the symbolic treatment of the system. For this reason, there is a time constraint of one second put on obtaining the equations. If that time is exceeded, then NDSolve generates a message and attempts to solve the system as a DAE.

The following method options can be specified to the simplification methods to better control the behavior of "EquationSimplification".

As indicated by the message, NDSolve did not solve explicitly for the derivatives because Solve had exceeded the default time constraint of one second. The amount of time that NDSolve spends on obtaining an explicit expression can be controlled using the option "TimeConstraint".

In the preceding example, due to the square terms in the first equation, four possible solutions are encountered. When solving as a DAE, the numerical solution will only follow one of these solution branches, depending on how it is initialized.

To demonstrate the usage of the method option "SimplifySystem", consider the following example:

    

An analytic solution for the above system can be obtained by substituting into the first equation to get a solution for , which in turn is substituted into the second equation to give a solution for . The final solution for the system is

    .

NDSolve cannot solve this system without first performing index reduction of the system and forming a new equivalent system of ODEs. This could prove to be expensive. When the suboption "SimplifySystem" is turned on, NDSolve detects any variable for which an analytic expression/solution can be found and performs repeated substitutions back into the original system. This results in the original system being either simplified or in some cases (like the above) getting back analytic solutions.

If an analytic solution does not exist for certain variables, then NDSolve will return an interpolating function as the solution.

With the suboption turned on, constant parameters can also be directly substituted into the system.

DAE Solution Methods

There are a variety of solution methods built into NDSolve for solving DAEs.

Two methods work with the general residual form of index-1 DAEs, .

• IDA—(Implicit Differential-Algebraic solver from the SUNDIALS package) based on backward differentiation formulas (BDF)

• StateSpace—implicitly solves for derivatives to use an underlying ODE solver

If the system has index higher than 1, both of these solvers typically fail. To accurately solve such systems, index reduction is needed. Note that index-1 systems can be reduced to ODEs, but it is often more efficient to use one of the solvers above.

NDSolve also has some solvers that work with DAEs that can be reduced to special forms.

• MassMatrix—for DAEs of the form

• Projection—for ODEs with invariants

The following subsections describe some aspects of using these methods for DAEs.

ODEs with Invariants

It is common to encounter DAEs of the form

    

where is an invariant that is consistent with the differential equations.

When a system of DAEs is converted to ODEs via index reduction, the equations that were differentiated to get the ODE form are consistent with the ODEs, and it is typically important to make sure those equations are well satisfied as the solution is integrated.

Such systems have more equations than unknowns and are overdetermined unless the constraints are consistent with the ODEs. NDSolve does not handle such DAEs directly because it expects a system with the same number of equations and dependent variables. In order to solve such systems, the "Projection" method built into NDSolve handles the invariants by projecting the computed solution onto after each time step. This ensures that the algebraic equations are satisfied as the solution evolves.

An example of such a constrained system is a nonlinear oscillator modeling the motion of a pendulum.

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

Solving Systems with a Mass Matrix

When the derivatives appear linearly, a differential system can be written in a form

    

where is a matrix often referred to as a mass matrix. Sometimes is also referred to as a damping matrix. If the matrix is nonsingular, then it can be inverted and the system can be solved as an ODE. However, in the presence of a singular matrix, the system can be solved as a DAE. In either case, it may be advantageous to take advantage of the special form.

The system describing the motion of a particle on a cylinder

    

can be expressed in mass matrix form with

    .

Index Reduction for DAEs

The built-in solvers for DAEs in NDSolve currently handle index-1 systems of DAEs fully automatically. For higher-index systems, an index reduction is necessary to get to a solution. This index reduction can be performed by giving appropriate options to NDSolve.

NDSolve uses symbolic techniques to do index reduction. This means that if your DAE system is expressed in the form , where NDSolve cannot see how the components of interact and compute symbolic derivatives, the index reduction cannot be done, and for this reason NDSolve does not do index reduction, by default.

For an example, consider the classic DAE describing the motion of a pendulum:

where there is a mass at the point constrained by a string of length . is a Lagrange multiplier that is effectively the tension in the string. For simplicity in the description of index reduction, take . The figure shows the schematic of the pendulum system.

        Sketch of the forces acting on a pendulum.

If the index of the system is not obvious from the equations, one thing to do is try it in NDSolve, and if the solver is not able to solve it, in many cases it is able to generate a message indicating what the index appears to be for the initial conditions specified.

The message indicates that the index is 3, indicating that index reduction is necessary to solve the system. Note that a complete consistent initialization is used here to avoid possible issues with initialization. This aspect of the solution procedure is treated in "Consistent Initialization of DAEs" with more detail.