How do i conserve population in a ODE or DAE

I have a set of differential equations dy_1/dt, ..., dy_n/dt in which the variables y_i must remain positive AND the total population y_1+...+y_n is conserved. This conservation is intrinsic to the ODEs, but numerical integration is always problematic.

I've tried 'nonnegative' and 'events' in several ode solvers but none respects conservation of population. I've also written the problem as a DAE, adding the constraint y_1+...+y_n-1=0. This doesn't work either.

Conservation is violated from the time at which the first variable hits 0 and so is likely due to truncation error.

Any suggested means of tackling this problem would be very welcome,

David

[Aykut Satici]

Hello David,

It is a known fact that Runge-Kutta methods (such as ode45) respect the conservation of quantities that can be expressed as linear functions of the solution. For a proof, you can see

Rosenbaum, J.S., "Conservation Properties of Numerical Integration Methods for Systems of Ordinary Differential Equations," Academic Press, 1976

I think when you add the nonnegativity constraint, MATLAB changes the Runge-Kutta method in a way that destroys this property. To impose the constraints and the conservation law, it may be necessary to come up with a custom integration algorithm that is numerically consistent (the solution converges to the real solution in an appropriate function space) with the desired properties.

Could you please provide your vector field (the right hand side of your differential equations) so I can reproduce the issue?

Bests, Aykut