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.