Introduction
In general, a system of ordinary differential equations can be expressed in the normal form,
The derivatives of the dependent variables x are expressed explicitly in terms of the independent variable, t, and the dependent variables, x. As long as the function f has sufficient continuity, a unique solution can always be found for an initial value problem where 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, it is typical that derivatives of some of the dependent variables will not even appear in the equations at all. The general form of a system of DAEs is
where the Jacobian with respect to x', 
may be singular.
A system of DAEs can be converted to a system of ODEs by differentiating it with respect to the independent variable t. 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.
One issue which is quite different for DAEs than for ODEs is the specification of initial conditions. For ODEs, as mentioned above 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 which satisfy the equations at all. To better understand this issue, consider the following example [AP98]
Here is a system of DAEs with 3 equations, but only one differential term.
The initial conditions are clearly not free; the second equation requires that 
is either 0 or 1.
This solves the system of DAEs starting with a specified initial condition for the derivative of 
.
To get the solution above, NDSolve first searches for initial conditions which 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 which NDSolve found.
This shows a plot of the solution. The solution 
is obscured by the solution 
which has the same constant value of 1.
However, there may not be a solution from all initial conditions which satisfy the equations.
This tries to find a solution with 
starting from steady state with derivative 0.
NDSolve::nderr: "Error test failure at \!\(t\) == \!\(0.`\); unable to continue."
This shows the initial conditions which NDSolve found.
If you look at the equations with 
set to 1, you can see why it is not possible to advance beyond t = 1.
Substitute 
into the equations
The middle equation effectively drops out. If you differentiate the last equation with
, you get the condition that 
, 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 necessary to specify a value for 
also because it is a differential variable.
This shows a plot of the nonzero components of the solution.
In general, it is at the very least necessary to specify initial conditions for the differential variables because it is quite typical that 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 will not always be able to find initial conditions consistent with the equations since, sometimes this is quite a difficult problem to solve, in fact, "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.
NDSolve::icfail: "Unable to find initial conditions which satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions."
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 will use 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 a chemical kinetics problem due to Robertson. Because of the large and small rate constants, the problem is quite stiff.
This solves the Robertson kinetics problem as an ODE by differentiating the balance equation.
The stiffness of the problem can be seen by seeing that 
and 
have their main variation on two completely different time scales.
This shows the solutions 
and 
.
This solves the Robertson kinetics problem as a DAE.
The solutions for a given component will appear to be quite close, but comparing the chemical balance constraint shows a difference between the two solutions.
This loads a package which makes it easier to extract the coordinates at which NDSolve saved data; these will give better information about the actual error in the balance equation.
This shows 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 both t and the magnitude of the error.l
In this case, both solutions satisfied the balance equations well beyond expected tolerances. Take 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 satisfy the constraint better once the range of quick variation is passed.
It may arise that you 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 also NDSolve expects the number of equations to be the same as the number of dependent variables. NDSolve, however, does have a Projection method which will often solve the problem.
A very simple example of such a constrained system is a nonlinear oscillator modelling the motion of a pendulum.
This defines the equation, invariant constraint, and starting condition for s simulation of 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 just use the invariant.
This solves for the motion of a pendulum using just the invariant equation. The SolveDelayed option tells NDSolve not to symbolically solve the quadratic equation for x', but instead to solve the system as a DAE.
However, this solution may not be quite what you are expecting: the invariant equation has the solution x[t] == constant when it starts with x'[t] == 0. In fact it does not have unique solutions from this starting point. This is because if you do actually solve for x', the function does not satisfy the continuity requirements for uniqueness.
This solves for the motion of a pendulum using just the differential equation. The method ExplicitRungeKutta is used because it can be used also as a submethod of the projection method.
This shows the solution plotted over the last several periods.
This shows a plot of the invariant at the ends of the time steps NDSolve took
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.
This shows a plot of the invariant at the ends of the time steps NDSolve took with the projection method.
Methods
IDA
The IDA package is part of the SUNDIALS (SUite of Nonlinear and DIfferential/ALgebraic equation Solvers) developed at the Center for Applied Scientific Computing of Lawrence Livermore National Laboratory. As described in the IDA user guide [HT99], "IDA is a general purpose solver for the initial value problem for systems of differential-algebraic equations (DAEs). The name IDA stands for Implicit Differential-Algebraic solver. IDA is based on DASPK ...". DASPK [BHP94], [BHP98] is a FORTRAN code for solving large scale differential-algebraic systems.
In Mathematica, an interface has been provided to the IDA package so that rather than needing to write a function in C for evaluating the residual and compiling the program, Mathematica generates the function automatically from the equations you input to NDSolve.
IDA solves the system (1) with Backward Differentiation Formula (BDF) methods of orders 1 through 5, implemented using a variable-step form. The BDF of order k is at time 
is given by the formula
The coefficients 
depend on the order k and past step sizes. Applying the BDF to the DAE (1) gives a system of nonlinear equations to solve:
The solution of the system is achieved by Newton-type methods, typically using an approximation to the Jacobian
"Its [IDAs] most notable feature is that, in the solution of the underlying nonlinear system at each time step, it offers a choice of Newton/direct methods or an Inexact Newton/Krylov (iterative) method." [HT99] In Mathematica, you can access these solvers using method options, or use the default Mathematica LinearSolve which switches automatically to direct sparse solvers for large problems.
At each step of the solution, IDA computes an estimate 
of the local truncation error and the step size and order are chosen so that the weighted norm 
is less than 1. The 
component, 
, of 
is given by
The values prec and acc are taken from the NDSolve settings for the PrecisionGoal->prec and AccuracyGoal->acc.
Because IDA provides a great deal of flexibility, particularly in how the nonlinear equations are solved, there are a number of method options which allow you to control how this is done. You can use the method options to IDA by giving NDSolve the option Method->{IDA, ida method options}.
The options for the IDA method are associated with the symbol IDA in the NDSolve` context.
IDA method options.
The remainder of this documentation will focus on suboptions of the two possible settings for the "ImplicitSolver" option, which can be "Newton" or "GMRES". With "Newton" the Jacobian or an approximation to it is computed and the linear system is solved to find the Newton step. On the other hand, "GMRES" is an iterative method and rather than computing the entire Jacobian, a directional derivative is computed for each iterative step.
The "Newton" method has one method option, LinearSolveMethod, which you can use to tell he Mathematica how to solve the linear equations required to find the Newton step. There are several possible values for this option.
Possible settings for LinearSolveMethod option.
The "GMRES" method may be substantially faster, but is typically quite a bit more tricky to use because for it to really be effective, it typically requires a preconditioner, which can be supplied via a method option. There are also some other method options which control the Krylov subspace process. For use of these refers to the IDA user manual [HT99]
GMRES method options.
As an example problem, consider a 2-dimensional Burgers` equation.
This can typically be solved with an ordinary differential equation solver, but in this example, two things are achieved by using the DAE solver. First, boundary conditions are enforced as algebraic conditions. Second, NDSolve is forced to use conservative differencing by using an algebraic term. For comparison, a known exact solution will be used for initial and boundary conditions.
This defines a function which satisfies Burger's equation
This defines initial and boundary conditions for the unit square consistent with the exact solution.
This defines the differential equation.
This sets the diffusion constant ν to a value for which we can find a solution in a reasonable amount of time and shows a plot of the solution at t == 1.
You can see from the plot that with ν = 0.025, there is a fairly steep front. This moves with constant speed.
This solves the problem using the default settings for NDSolve and the IDA method with the exception of the "DifferentiateBoundaryConditions" option for MethodOfLines, which causes NDSolve to treat the boundary conditions as algebraic.
Since there is an exact solution to compare to, the overall solution accuracy can be compared as well.
This defines a function which finds the maximum deviation between the exact and computed solutions at the grid points at all of the time steps.
This computes the maximal error for the computed solution.
In the following, a comparison will be made with different settings for the options of the IDA method. To emphasize the option settings a function will be defined which times the computation of the solution and gives the maximal error.
This defines a function for comparing different option settings of IDA.
No options uses the defaults as above.
This uses the "Band" method
The band method is not very effective because the band with of the Jacobian is relatively large, partly because of the fourth order derivatives used and partly because the one-sided stencils used near the boundary add width at the top and bottom. You can specify the band width explicitly.
This uses the "Band" method with the width set to include the stencil of the differences in only one direction.
While the solution time was smaller, notice that the error is slightly greater and the total number of time steps is a lot greater. If the problem was more stiff, the iterations likely would not have converged because it was missing information from the other direction. Ideally, the band width should not eliminate information form an entire dimension.
This computes the grids used in the X and Y directions and shows the number of points used.
This uses the "Band" method with the width set to include at least part of the stencil in both directions.
With the more appropriate setting of the band width, the solution is still slightly slower than in the default case. The band method can sometimes be effective on 2-dimensional problems, but is usually most effective on 1-dimensional problems.
This computes the solution using the GMRES implicit solver without a preconditioner.
This is incredibly slow! Using the GMRES method without a preconditioner is not recommended for this very reason. However, finding a good preconditioner is not usually trivial. For this example, a diagonal preconditioner will be used.
The setting of the Preconditioner option should be a function f which accepts four arguments which will be given to it by NDSolve such that f[t, x, x', c] returns another function which will apply the preconditioner to the residual vector. (See IDA User Guide [HT99] for details on how the preconditioner is used.) The arguments t, x, x', c are the current time, solution vector, solution derivative vector, and the constant c in formula (2) above. For example, if you can determine a procedure which would generate an appropriate preconditioner matrix P as a function of these arguments, you could use
Preconditioner->Function[{t,x,xp,c}, LinearSolve[P[t,x,xp,c]]]
to produce a LinearSolveFunction object which will effectively invert the preconditioner matrix P. Typically, for each time the preconditioner function is set up, it is applied to the residual vector several times, so using some sort of factorization such as is contained in a LinearSolveFunction is a good idea.
For the diagonal case, the inverse can be effected simply by using the reciprocal. The most difficult part of setting up a diagonal preconditioner is keeping in mind that values on the boundary should not be affected by it.
This finds the diagonal elements of the differentiation matrix for computing the preconditioner
This gets the positions where elements at the boundary, which satisfy a simple algebraic condition are in the flattened solution vector.
This defines the function which sets up the function called to get the effective inverse of the preconditioner. For the diagonal case, the inverse is done simply by taking the reciprocal.
This uses the preconditioned GMRES method to compute the solution.
Thus, even with a crude preconditioner, the GMRES method computes the solution faster than the using the direct sparse solvers.
For PDE discretizations with higher order temporal derivatives or systems of PDEs, you may need to look at the corresponding NDSolve`StateData object to determine how the variables are ordered so that you can get the structural form of the preconditioner correctly.
