Which Differential Solver do I need for my system?

[DonaldStallman]

Hello,

I am trying to implement the Anaerobic Digestion Model No.1 (ADM1) in MathCAD. The model has 27 dynamic state variables, for which I have 27 1st order, non-linear differential equations. The model also implements a set of 9 implicit algebraic equations that must be satisfied.

I was able to set up a solve block and solve all of the equations with 'ODESOLVE', but when looking at the MathCAD help, I'm unclear whether this solver is able to handle non-linear systems.

The help says that ODESOLVE is only for linear systems, but then says it dynamically detects which solver to use depending on the system. One of the choices is rkfixed, which the help claims can solve non-linear systems.

The reason I ask is that while most of the variables are returning sensible values, a couple of the values that are calculated by the algebraic constraints are coming out negative. I have been troubleshooting reasons for this for a while now, and one of them is that maybe I am using the wrong solve function for my system.

Any thoughts would be greatly appreciated,

Donnie

[AlanStevens]

Right-click on ODESOLVE and you should be able to choose from a selection of solvers - try them all!

Can you upload the worksheet? - difficult to give more specific advice without it.

Alan

Message was edited by: AlanStevens Correction! If you include algebraic constraints Odesolve automatically chooses the Radau solver - and you can't override this!

[DonaldStallman]

Thanks for your response Alan.

I am seeing that it is forcing me into the Radau solver, does this mean it is capable of nonlinear? When I look into help I'm mostly just seeing how it can solve stiff systems with algebraic constraints, but no mention of its capability to handle nonlinear.
I have attached MC13 and MC15 versions of my file. I apologize as I'm new to MathCAD and so my sheet is a little messy.

The negative values are coming from my charge balance and I end up with a negative proton concentration in solution.

-Donnie

[AlanStevens]

Hell's bells and buckets of blood!!! That's not easy on the eye!

Donald Stallman wrote:

Thanks for your response Alan.

I am seeing that it is forcing me into the Radau solver, does this mean it is capable of nonlinear?

Yes - it is solving your equations - just not giving you the answers you want! (Incidentally, you don't need to specify the number of intervals as 200 - it works without specifying any number here).

Check and double check that your equations are correct (especially the signs of each term -especially SH). I don't have the domain knowledge to help there (nor the patience to trawl through all of them!).

There are some large instantaneous changes from initial values for some parameters - are you sure all your initial conditions are correct and consistent?

You could simplify the equation set somewhat. For example you could replace SHCO3 with Ka_co2.SIc(t)/(Ka_co2+SH(t)) everywhere it occurs inside the solve block, and just define it as a function outside the solve block. You could then reduce the number of variables to solve for. You could do this for some of your other constraint variables as well. Not that this will necessarily solve your problem, but it might make it easier to identify the problem.

I'll have another look when I have a little more time.

Alan