Solve Differential Algebraic Equations (DAEs)

What is a Differential Algebraic Equation?

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 calledalgebraic, and the presence of algebraic variables means that you cannot write down the equations in the explicit form $y' = f (t, y)$ . Instead, you can solve DAEs with these forms:

Theode15sandode23tsolvers can solve index-1 linearly implicit problems with a singular mass matrix $M (t, y) y' = f (t, y)$ , including semi-explicit DAEs of the form

$\begin{matrix} y' = f (t, y, z) \\ 0 = g (t, y, z) . \end{matrix}$

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.

$M y' = (\begin{matrix} y'_{1} & 0 & \dots & 0 \\ 0 & y'_{2} & 0 & ⋮ \\ ⋮ & 0 & ⋱ & 0 \\ 0 & \dots & 0 & 0 \end{matrix}) .$

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 theMassSingularoption ofodesetto'yes'. With DAEs, you can also provide the solver with a guess of the initial conditions for $y'_{0}$ using theInitialSlopeproperty ofodeset. This is in addition to specifying the usual initial conditions for $y_{0}$ in the call to the solver.
Theode15isolver can solve more general DAEs in the fully implicit form

$f (t, y, y') = 0 .$

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.

$J = \partial f / \partial y' = (\begin{matrix} \frac{\partial f_{1}}{\partial y'_{1}} & \dots & \frac{\partial f_{1}}{\partial y'_{n}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial f_{m}}{\partial y'_{1}} & \dots & \frac{\partial f_{m}}{\partial y'_{n}} \end{matrix})$

Theode15isolver requires that you specify initial conditions for both $y'_{0}$ and $y_{0}$ . Also, unlike the other ODE solvers,ode15irequires the function encoding the equations to accept an extra input:odefun(t,y,yp).

DAEs arise in a wide variety of systems because physical conservation laws often have forms like $x + y + z = 0$ . Ifx,x',y, andy'are defined explicitly in the equations, then this conservation equation is sufficient to solve forz没有一个表达式orz'.

Consistent Initial Conditions

When you are solving a DAE, you can specify initial conditions for both $y'_{0}$ and $y_{0}$ . Theode15isolver requires both initial conditions to be specified as input arguments. For theode15sandode23tsolvers, the initial condition for $y'_{0}$ is optional (but can be specified using theInitialSlopeoption ofodeset). 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 calledinconsistent. The treatment of the initial conditions varies by solver:

ode15sandode23t— If you do not specify an initial condition for $y'_{0}$ , then the solver automatically computes consistent initial conditions based on the initial condition you provide for $y_{0}$ . If you specify an inconsistent initial condition for $y'_{0}$ , then the solver treats the values as guesses, attempts to compute consistent values close to the guesses, and continues on to solve the problem.
ode15i— The initial conditions you supply to the solver must be consistent, andode15idoes not check the supplied values for consistency. The helper functiondeciccomputes consistent initial conditions for this purpose.

Differential Index

DAEs are characterized by theirdifferential 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

$y (t) = k (t) .$

For this equation, you can take a single derivative to obtain the explicit ODE form

$y' = k' (t) .$

An example of an index-2 DAE is

$\begin{array}{l} y'_{1} = y_{2} \\ 0 = k (t) - y_{1} . \end{array}$

These equations require two derivatives to be rewritten in the explicit ODE form

$\begin{array}{l} y'_{1} = k' (t) \\ y'_{2} = k'' (t) . \end{array}$

Theode15sandode23tsolvers 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 seeSolve Differential Algebraic Equations (DAEs)(Symbolic Math Toolbox)for more information.

Imposing Nonnegativity

Most of the options inodesetwork as expected with the DAE solversode15s,ode23t, andode15i. However, one notable exception is with the use of theNonNegativeoption. TheNonNegativeoption does not support implicit solvers (ode15s,ode23t,ode23tb) applied to problems with a mass matrix. Therefore, you cannot use this option to impose nonnegativity constraints on a DAE problem, which necessarily has a singular mass matrix. For more details, see[1].

Solve Robertson Problem as Semi-Explicit Differential Algebraic Equations (DAEs)

Open Script

This example reformulates a system of ODEs as a system of differential algebraic equations (DAEs). The Robertson problem found inhb1ode.mis a classic test problem for programs that solve stiff ODEs. The system of equations is

$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 - 10^4 y_2y_3- (3 \times 10^7)y_2^2\\ y'_3 &= (3 \times 10^7)y_2^2.\end{array}$

hb1ode解决了这个系统的常微分方程与th稳态e initial conditions $y_1 = 1$ , $y_2 = 0$ , and $y_3 = 0$ . But the equations also satisfy a linear conservation law,

$$y'_1 + y'_2 + y'_3 = 0.$$

In terms of the solution and initial conditions, the conservation law is

$$y_1 + y_2 + y_3 = 1.$$

The system of equations can be rewritten as a system of DAEs by using the conservation law to determine the state of $y_3$ . This reformulates the problem as the DAE system

$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 - 10^4 y_2y_3-(3 \times 10^7)y_2^2\\ 0 &= y_1 + y_2 + y_3 - 1.\end{array}$

The differential index of this system is 1, since only a single derivative of $y_3$ is required to make this a system of ODEs. Therefore, no further transformations are required before solving the system.

The functionrobertsdaeencodes this DAE system. Saverobertsdae.min your current folder to run the example.

functionout = 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 inhb1dae.m.

Solve the DAE system usingode15s. Consistent initial conditions fory0are obvious based on the conservation law. Useodesetto set the options:

Use a constant mass matrix to represent the left hand side of the system of equations.

$\left( \begin{array}{c} y'_1\\ y'_2\\ 0 \end{array} \right) = M y' \rightarrow M = \left( \begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 0 \end{array} \right)$

Set the relative error tolerance to1e-4.
Use an absolute tolerance of1e-10第二个解决方案组件,因为规模varies dramatically from the other components.
Leave the'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');

References

[1] Shampine, L.F., S. Thompson, J.A. Kierzenka, and G.D. Byrne. "Non-negative solutions of ODEs."Applied Mathematics and Computation. Vol. 170, 2005, pp. 556-569.