Differential equations describe the behavior of a system through rates of change.

Simplest Example:

This can extend to systems with more than one dependent variable.

It is not always possible to find the closed form solution of an ODE problem. This creates a need for solvers that compute a numerical approximation to the solution.

A DAE is a system that is a combination of differential equations and constraints.

One of the simplest examples is the pendulum system, where the m*g (mass times gravitational constant) term has been replaced by Pi^2.

To display derivatives and functions using a simpler notation, use the PDEtools[declare] command.

The typical solution approach for this problem (as presented in textbooks) is to apply a change of variables to  polar co-ordinates, which would proceed as follows.

The last equation implies r=1.

Evaluate and remove lambda(t) as follows, giving an ODE in theta(t).

This is the more commonly known pendulum equation.

As an approximation for small displacements, we can approximate sin(theta) ~ theta giving us the simple harmonic oscillator equation:

This is appropriate for the pendulum equation, but, in general, such a coordinate transform may not exist. Hence, there is the need for numerical solution methods for DAEs.

DAEs have an interesting property.

Consider our constraint:

Differentiate:

Differentiate, eliminate, and isolate:

One more differentiation and elimination yields an ODE for lambda(t).

The above process of successive differentiation until we obtain an ODE for each variable is called index reduction.
The number of required differentiations is called the gear index. It sometimes corresponds to the difficulty in obtaining a numerical solution of the DAE problem. (A greater index corresponds to a more difficult problem.)

Issues:

We obtained 3 algebraic relations for lambda, x, x', y, y': c1,c2,c3 - initial conditions must satisfy these - as a result, we really only have 2 degrees of freedom in the initial conditions for a problem that would normally have 5.
This restricts the initial conditions. They must be checked for consistency.
Note: Constraints are typically nonlinear, and cannot normally be solved in terms of the free initial conditions.

The algebraic constraints must hold along the entire solution.

As previously shown, this involves successive differentiation of system until the DAE becomes an ODE system.

Advantages

Allows us to obtain all constraints directly

Allows use of a standard ODE solver to obtain the solution

Disadvantages

Does not enforce the constraints, other than at the initial point (resulting in constraint drift)

Example:

Find suitable initial conditions:

From the first constraint, we have:

From the second constraint, we have:

And finally from the third constraint, we have:

Now, we can obtain a numerical solution for the index reduced system ignoring the constraints.

(Note: We set maxfun=0 to remove the limit on the number of function evaluations.)

Consider the value at t=1000:

Now consider the constraints at t=1000:

This indicates that the constraints are not satisfied.
This is one of the disadvantages of this approach.

Note: Problems arise when small changes to the problem dictate large changes in the solution at a later time (instability).

The Maple DAE extension solvers (rkf45_dae and rosenbrock_dae) are modifications to ODE solvers, and can operate in this mode:

 This approach is similar to the Method of Index Reduction, but involves the removal of variables that can be solved for.
Advantages and disadvantages are the same as those for the basic index reduction method, but this approach usually provides a better solution.

Examine the two remaining relevant constraints:

We see that the O(1) error for the general index reduction method has been reduced to an O(10^(-2)) error by using the index reduction with removal method.

Note: The Maple DAE extension solvers can operate in this way, but also have additional error control capabilities. They do not remove the index 1 equation for lambda(t), but instead use it in the numerical integration, and monitor it for possible error.
As a result, the solution is more accurate than previously.

The error in the constraints is now O(10^(-4)).

This approach is similar to the Method of Index Reduction, but at each step the solution is projected onto the constraints of the problem.
The magnitude of the projection is monitored as part of the error control strategy.

Advantages

Allows us to obtain all constraints directly

Enforces constraints, removing the constraint drift problem

Disadvantages

Requires a modified solver to add a projection step and error control.

Depending on the stability of system, the projection step may introduce errors in the solution.

The Maple rkf45_dae and rosenbrock_dae solvers operate in this mode by default.

Now the constraints are satisfied well within the tolerance at any solution point.

These methods work directly with the original DAE system, and do not require differentiation of the constraint.
They are typically limited to handling DAEs of index 1, 2, or 3.

Advantages

Does not require symbolic manipulation of the system

Disadvantages

Does not indicate constraints on the initial data, and simply fails if these are not satisfied.

Limited to problems with a maximum index.

The core solver is usually more complex than for the other approaches.

Maple has the mebdfi solver (Modified Extended Backward Difference Formula Implicit).
This solver has the ability to fully solve implicit DAE systems.
For the pendulum example, we get:

Important: The constraints are all satisfied within tolerance even though they are not explicitly used in the solution process.

The pendulum is an excellent example. We have a transformation that converts it to a single nonlinear ODE in theta, which we can solve to high accuracy using our gear method,. This allows us to compare our DAE solutions.

From the last two we get theta(0)=0, and so:

Now obtain a high accuracy solution to this problem using the 'gear' method.

Convert to x, y coordinates, and obtain the lambda value:

Compare with the basic index reduction solution, which is expected to be poor:

Compare with index reduction with removal, which is expected to be moderately poor:

Compare with index reduction with removal and additional error correction (extension methods with projection=false):

Compare with index reduction with projection (extension method in default mode):

Compare with direct method (mebdfi method):

In summary, we see that the Maple default rkf45_dae method and mebdfi method give the best agreement with the solution.

This multi-body system models the behavior of car suspension on a bumpy road.
It consists of 4 second order ODEs and two index-3 DAEs.

Auxiliary variables:

DAE system:

Initial conditions and variables:

Exact solution (solution computed via gear with parameter removal and 1e-22 accuracy):

Solution with default tolerances:

Compare:

Solution with tight tolerances:

Trajectory of x1,y1:

Trajectory of x2,y2:

The double pendulum is a pendulum with an additional weight attached to the first pendulum bob by a string.
It is a system with 4 second order ODEs and two index-3 DAEs.

(x1(t),y1(t)) are the coordinates of the first pendulum bob, and (x2(t),y2(t)) are the coordinates of the second, which is attached to the first.

Set masses, lengths of strings, and g:

DAE system and initial conditions:

Obtain the trajectory for the second mass:

Construct an animation of the system, where tf is the final time and nfr is the number of frames:

Consider the simple pendulum system and the worst rkf45_dae solution (differential=true, projection=false):

Since this is a pendulum, the trajectory (x(t),y(t)) should trace a portion of a circle repeatedly. Constraint drift causes a deviation this behavior.

The default solution mode works much better: