How to solve DAE system with non-negativity constraint

Hi

I am encountering problems when trying to solve an index one DAE system in the form of a reactive transport model using the ode15s solver. Issues arise because the solution has to be non-negative to make sense physically, while ode15s allow the 'NonNegative' option only for those problems for which there is no mass matrix. Without the 'NonNegative' option, components of the solution becomes negative and the program ultimately fails. I have read a number of articles concerning strategies for solving index on DAE with non-negative constraints, mostly concerned with so-called "clipping methods", "Newton-damping methods" and "CODOSOL method" to alter iterations outside the feasible region.

And so my question is: Does any Matlab ode-solver allow nonnegative constraints for index one DAE systems for which a mass matrix is supplied? Alternatively, are you aware of user-written programs which could solve my problem (ideally my program should be able to run using matlab solvers)? Finally, could the use of event functions solve my problem by resetting negative components of the solution? (I imagine this will have a huge impact on computation time).

Best regards

Anders

Answer 1

[Torsten]

The usual strategies are:

1. Choose smaller error tolerances for RelTol and/or AbsTol

2. Make a copy of the solution variables in the subroutine where you supply the time derivatives. Then check the copy for negative values and reset them to small positive values within the forthcoming calculations, e.g.

conc_A = y(1);

if conc_A < 0

   conc_A=1e-4;

end

Best wishes

Torsten.