Choice between DAE or ODE formulation for chemical systems

Question

Consider a simple ODE system describing the evolution of two chemical species undergoing the reaction $A = B$ :

$$ \frac{dn_A}{dt} = - k * n_A $$ $$ \frac{dn_B}{dt} = k * n_A $$

We can discretize that system and solve it using any classical numerical scheme (e.g. implicit Euler).

However, I can also recast that system as :

$$ \frac{dn_A}{dt} = - k * n_A $$ $$ n_A + n_B = n_A^0 + n_B^0 $$

which is a DAE that I can also solve using implicit Euler.

What are the advantages of each formulation ? In particular for more complex chemical systems, is there any interest in solving DAE systems instead of ODE systems ?

Answer 1

The trade-off is generally this:

• ODEs are numerically a lot of easier to solve than DAEs, in particular if the algebraic constraint is nonlinear (yours is linear). So that argues for the ODE formulation.
• But ODE solvers generally do not preserve invariants that may be inherent in the system -- here, conservation of mass $n_A(t)+n_B(t)=n_A(0)+n_B(0)$. So if preservation of this invariant is important to you, then you should use the DAE formulation.

Answer 2

Let's look at the results from 3 different chemical reaction DAE systems:

• Chemical Akzo
• OREGO
• ROBER

These are all implemented using ModelingToolkit.jl in order to generate the various forms. It's 3 forms actually: the singular mass matrix DAE form, the fully implicit DAE form, and the over-determined ODE form (i.e. solving an ODE of only the differential variables by embedding a nonlinear solve of the algebraic equation into the system, somewhat mentioned in the comments). That's a lot of benchmarking to look into, but let me pull out a representative result:

The prob_choices are in the order I described. You can see that generally the version with the embedded nonlinear solves (with tearing) tend to perform the best, followed by the mass matrix DAE form, followed by the fully implicit DAE form. I'm skipping over a lot of details that can happen in DAE simplification to the ODAE form but in general it tends to perform very well in comparison to other methods which is why symbolic systems like ModelingToolkit and Modelica default to this form. The fully implicit DAE solver form that software like IDA uses is pretty clearly outperformed in every case you can avoid it. Giving a bit of structure, i.e. describing the system as an ODE with a mass matrix, is a good intermediate that's not hard to do by hand.

Now how does the mass matrix DAE form relate to the one you described? This tutorial lays out pretty clearly how one formulates it:

Now why is it important? It turns out that you can solve equations in the singular mass matrix form without having extra overhead in the solver (therefore, "ODEs are numerically a lot of easier to solve than DAEs" is not necessarily true). These course notes go into a lot more detail, but the key is that the Newton iterations for this kind of solver simply use a Jacobian of $M - \Delta t J$ instead of the standard $I - \Delta t J$ you would use in a "normal" implicit Euler. That's all you have to do, and now the mass matrix and the hard constraint is encoded in the solving. So it's essentially "free" in the case of a semi-explicit ODE, i.e. where you have $I$ mass matrix plus some constraint equations.

The caveat is not all methods can make use of this formulation. Most SDIRK integrators in particular do not support singular mass matrices, so most of the time you'd have to use a BDF or Rosenbrock method. However, if you can make use of the mass matrix form, you get effectively no cost difference while encoding the hard constraint, so it's just generally a good idea to use that.

In total, the hierarchy of the ways to solve are as follows:

tl;dr

1. If you have a symbolic system with tearing and all of that jazz, use a special torn ODAE formulation if it's stable. There are cases where it won't be stable because this will most strictly enforce the algebraic equations, which can cause an issue if you need branch switching.
2. Use the mass matrix formulation when you can.
3. Fall back to a stiff ODE solver with the constraints implicit in the equation form when you have to.
4. Use a fully implicit DAE solver (i.e. IDA) only when you really really have to. It's generally slower and will have more numerical difficulties doing things like reinitialization.

That said, (4) isn't always true. There are times where the DAE form is truly implicit and moving to an mass matrix ODE form (which always can be done) requires expanding the system and sometimes ruining some of the structure of the sparsity pattern, so there are cases where (4) is faster than (1) and (2). But if you had to guess the numerical difficulty and the performance you'd get on a random problem (where my implicit distribution is questions asked on https://discourse.julialang.org/), you'd be a pretty good gambler if you always bet in this order.

(Note: there is a software specific thing where BDF methods do generally have better scaling to large problems than Rosenbrock methods, but Sundials CVODE doesn't support mass matrices last I checked and so IDA can be preferred if you need to use a BDF, though this is a software limitation and there are ODE BDF implementations like FBDF which support singular mass matrices. But this quirk is one reason why IDA is a strong go-to for large DAEs even in cases where a mass matrix form may be more appropriate, given how ubiquitous SUNDIALS and implementation of its lineage have been)