Solving a system of DAEs versus ODEs (which is preferable)

Assume I have the following system of equations. I am trying to figure out if it's better for me to be solving this as a system of ODEs or as a system of DAEs. The real code can have up to a dozen or so equations, and the system is stiff (I am dealing with chemical kinetics, where the rate constants can vary by many orders of magnitude).

The governing ODEs are: $$\frac{d\theta_{A}}{dt}=2R_{1}(\theta_{A},\theta_{*})-R_{2}(\theta_{A},\theta_{*})$$

$$\frac{d\theta{*}}{dt}=-2R_{1}(\theta_{A},\theta_{*})+R_{1}(\theta_{A},\theta_{*})$$

I can solve for $\theta_{A}$ and $\theta_{*}$ by solving this system. However, it also happens that there is a conservation law such that

$$\theta_{A}+\theta_{*}=1$$

This means I can solve the system of equations consisting of, for example, the first ODE and the algebraic conservation law to get $\theta_{A}$ and $\theta_{*}$ as well.

Computationally, which is easier to solve numerically? I am more concerned with accuracy than speed but am curious either way. Is it easier to solve the system of two ODEs or one of the ODEs and the algebraic conservation law?

In case it's relevant, I'm currently thinking about using ODE15s in MATLAB, which can solve a stiff system of ODEs or DAEs.

• Generally ODEs are easier to compute. But... from that conservation law, isn't there only a single ODE? By substitution the ODEs decouple – Chris Rackauckas Apr 2 '17 at 8:20
• I believe that if you solve your two ODE, ignoring the constraint equation, you will find that the constraint will be only approximately satisfied and that error will increase with time. – Bill Greene Apr 2 '17 at 12:33
• @ChrisRackauckas Yes, that is correct. However, I could solve either the two ODEs (ignoring the constraint equation) or I could solve one ODE and the constraint equation (since the constraint equation decouples one of the ODEs). I'm just not sure which is preferable for a more complicated system. – Argon Apr 2 '17 at 16:31
• @BillGreene This was my gut feeling as well. The constraint equation is a physical law that must hold, and I I'd prefer that to be accurate. – Argon Apr 2 '17 at 16:32
• If you care about having conservation law perfectly (to numerical precision) satisfied, then use the DAE. If you're fine with some (small) error, then use the ODE since it will be easier to find solvers or to even write your own. – spektr Apr 3 '17 at 17:20

Depends on the level from where you consider your problem.

As far as the mathematical equations are concerned, the ODE formulation is preferable, since there are a lot of fine-tuned codes out there, that solve ODEs.

If you have concerns about the validity of your solutions, you should step back and check where the equations came from. As you said, often there are (algebraic) conservation laws in the background. In theory, these laws are implicitly satisfied. In practice, however, there are deviations. From this point of view, it is preferable to have the constraints explicitly in the equations that are solved.

One more general remark. It is valid to say that under certain continuity assumptions, a sufficiently accurate approximation of the ODEs will satisfy the implicit constraints with arbitrary accuracy. However,

• in applications, the constraints (think of the train that should stay on the rails) need to be fulfilled exactly
• in theory, sometimes the ODEs are derived with the assumption that the conservation law is exactly fulfilled (e.g. the Navier--Stokes equations in the standard formulation of the convection term). Then an approximate solution may render the model itself invalid.