# Numerical Solution of Differential-Algebraic Equations

## 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 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.

In[10]:= |

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].

In[11]:= |

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

In[12]:= |

Out[12]= |

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.

In[13]:= |

Out[13]= |

In[15]:= |

Out[15]= |

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

In[16]:= |

Out[16]= |

In[17]:= |

Out[17]= |

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

In[18]:= |

Out[18]= |

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.

In[19]:= |

Out[19]= |

In[21]:= |

Out[21]= |

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].

In[22]:= |

Out[22]= |

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.

In[23]:= |

In[26]:= |

Out[26]= |

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

In[27]:= |

Out[27]= |

In[33]:= |

Out[33]= |

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

In[34]:= |

Out[37]= |

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.

In[55]:= |

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

In[58]:= |

Out[58]= |

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.

In[59]:= |

Out[59]= |

In[60]:= |

Out[60]= |

In[61]:= |

Out[62]= |

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.

In[63]:= |

Out[63]= |

In[64]:= |

Out[65]= |