Differential algebraic equations are a type of differential equation where one or more derivatives of dependent variables are not present in the equations. Variables that appear in the equations without their derivative are called algebraic, and the presence of algebraic variables means that you cannot write down the equations in the explicit form . Instead, you can solve DAEs with these forms:
In this form, the presence of algebraic variables leads to a singular mass matrix, since there are one or more zeros on the main diagonal.
By default, solvers automatically test the singularity of the
mass matrix to detect DAE systems. If you know about singularity ahead
of time then you can set the
'yes'. With DAEs,
you can also provide the solver with a guess of the initial conditions
InitialSlope property of
This is in addition to specifying the usual initial conditions for in
the call to the solver.
can solve more general DAEs in the fully implicit form
In the fully implicit form, the presence of algebraic variables leads to a singular Jacobian matrix. This is because at least one of the columns in the matrix is guaranteed to contain all zeros, since the derivative of that variable does not appear in the equations.
ode15i solver requires that you specify
initial conditions for both and .
Also, unlike the other ODE solvers,
the function encoding the equations to accept an extra input:
DAEs arise in a wide variety of systems because physical conservation
laws often have forms like .
y' are defined explicitly in the equations,
then this conservation equation is sufficient to solve for
having an expression for
When you are solving a DAE, you can specify initial conditions
for both and .
ode15i solver requires both initial conditions
to be specified as input arguments. For the
the initial condition for is
optional (but can be specified using the
odeset). In both cases, it is possible that
the initial conditions you specify do not agree with the equations
you are trying to solve. Initial conditions that conflict with one
another are called inconsistent. The treatment
of the initial conditions varies by solver:
If you do not specify an initial condition for ,
then the solver automatically computes consistent initial conditions
based on the initial condition you provide for .
If you specify an inconsistent initial condition for ,
then the solver treats the values as guesses, attempts to compute
consistent values close to the guesses, and continues on to solve
ode15i — The initial conditions
you supply to the solver must be consistent, and
not check the supplied values for consistency. The helper function
decic computes consistent initial conditions
for this purpose.
DAEs are characterized by their differential index, which is a measure of their singularity. By differentiating equations you can eliminate algebraic variables, and if you do this enough times then the equations take the form of a system of explicit ODEs. The differential index of a system of DAEs is the number of derivatives you must take to express the system as an equivalent system of explicit ODEs. Thus, ODEs have a differential index of 0.
An example of an index-1 DAE is
For this equation, you can take a single derivative to obtain the explicit ODE form
An example of an index-2 DAE is
These equations require two derivatives to be rewritten in the explicit ODE form
only solve DAEs of index 1. If the index of your equations is 2 or
higher, then you need to rewrite the equations as an equivalent system
of index-1 DAEs. It is always possible to take derivatives and rewrite
a DAE system as an equivalent system of index-1 DAEs. Be aware that
if you replace algebraic equations with their derivatives, then you
might have removed some constraints. If the equations no longer include
the original constraints, then the numerical solution can drift.
If you have Symbolic Math Toolbox™, then see Step 3. Check and Reduce Differential Index (Symbolic Math Toolbox) for more information.
This example reformulates a system of ODEs as a system of differential algebraic equations (DAEs). The Robertson problem found in hb1ode.m is a classic test problem for programs that solve stiff ODEs. The system of equations is
hb1ode solves this system of ODEs to steady state with the initial conditions , , and . But the equations also satisfy a linear conservation law,
In terms of the solution and initial conditions, the conservation law is
The system of equations can be rewritten as a system of DAEs by using the conservation law to determine the state of . This reformulates the problem as the DAE system
The differential index of this system is 1, since only a single derivative of is required to make this a system of ODEs. Therefore, no further transformations are required before solving the system.
robertsdae encodes this DAE system. Save
robertsdae.m in your current folder to run the example.
function out = robertsdae(t,y) out = [-0.04*y(1) + 1e4*y(2).*y(3) 0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2 y(1) + y(2) + y(3) - 1 ];
The full example code for this formulation of the Robertson problem is available in hb1dae.m.
Solve the DAE system using
ode15s. Consistent initial conditions for
y0 are obvious based on the conservation law. Use
odeset to set the options:
Use a constant mass matrix to represent the left hand side of the system of equations.
Set the relative error tolerance to
Use an absolute tolerance of
1e-10 for the second solution component, since the scale varies dramatically from the other components.
'MassSingular' option at its default value
'maybe' to test the automatic detection of a DAE.
y0 = [1; 0; 0]; tspan = [0 4*logspace(-6,6)]; M = [1 0 0; 0 1 0; 0 0 0]; options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]); [t,y] = ode15s(@robertsdae,tspan,y0,options);
Plot the solution.
y(:,2) = 1e4*y(:,2); semilogx(t,y); ylabel('1e4 * y(:,2)'); title('Robertson DAE problem with a Conservation Law, solved by ODE15S');