Mixed Algebraic/differential system
- To: mathgroup@smc.vnet.net
- Subject: [mg11420] Mixed Algebraic/differential system
- From: SR Thomas <srthomas@ceylan.necker.fr>
- Date: Thu, 12 Mar 1998 01:32:38 -0500
- Organization: CITI2 - Universite Rene Descartes, Paris
Dear MathGroup,
I want to work out an algorithm for a mixed algebraic/ode system such
as:
y'[t]=f(y,t)
is system of odes and
g(y)=0
is a set of algebraic constraints that must be satisfied at all times.
For example, in a compartment containing CO2 and assuming
instantaneous buffer equilibrium between co2, hco3 and protons, plus a
time-dependent exchange of some or all of these with the environment,
we have the odes:
co2'[t] == jco2IN[t] - jco2OUT[t] - jr[t]
(1)
hco3'[t] == jhco3IN[t] -jhco3OUT[t]+ jr[t]
(2)
h'[t] == jhIN[t] - jhOUT[t]+
jr[t]
(3)
taking the reaction rate jr to be positive in the direction of hco3 and
protons.The jiIN and jiOUT are input and output flows that depend on
time and possibly on co2, hco3, and h concentrations. They are given
explicitly, e.g., if co2 input is constant and output is by diffusion
exchange with an outside bath, then jco2IN = 5e-3 mmoles/min, and
jco2OUT=Pco2 (co2[t]-co2BATH), where Pco2 is the permeability to co2.
jr[t] is unknown, but we can't solve for it within the system of
odes because we don't have enough equations. The 4th equation we have
is the algebraic constraint for (effective) reaction equilibrium at all
times (Henderson-Hasselbalch).
Keq == hco3[t] h[t]/co2[t]
or
Keq - (hco3[t] h[t]/co2[t]) ==
0
(4)
which must be valid at all times. Keq is the reaction's equilibrium
constant.
Strictly speaking, reaction equilibrium would imply zero reaction
rate, but in reality it means the rx rate constants are so high that we
get significant reaction flow on the scale of the other system inputs
and outputs without displacing the reactant concentrations measureably
away from their equilibrium distribution. (We could of course treat the
reaction explicitly, giving its rate equation in terms of the
concentrations, but that makes the system very stiff and renders the
numerical analysis more difficult.)
I know we can do this (I've done it) by writing the odes as finite
differences and solving the resulting algebraic system for four
variables: i.e., jr plus the three concentrations at t+Dt. However,
this means choosing a constant delta t or writing some adjustable
step-size algorithm. It's easier to use a canned fortran routine like
lsodi.f.
I would prefer doing it neatly in Mathematica, if there's a way.
So here's my proposed scheme:
To solve for the three concentrations co2, hco3 and h from time t1 to
t2:
0) set or initialize: epsilonG, concentrations, permeabilities, input
rates, etc.
1) set t=t1
2) choose a small time increment, delta t
3) give a guess for jr (after the first trip through the loop, start
with the last successful value)
4) solve the odes (1) to (3) for the three concentrations from t to
t+delta t, (use NDSolve?)
5) calculate deviation from g(y)=0 (equation (4)), i.e., how well is
Hend-Hass. obeyed?
6) if deviation is less than epsilonG, move to the next time step (go
to step (6));
7) if deviation > epsilonG, adjust the guess for jr, and go back to 2)
8) if t=t2 we're done, else
9) increment t (but not past t2) and goto 1
The mysterious part is step 7, adjusting the jr guesses.
So my questions are:
- is there a better choice than NDSolve for the incremental ode
solutions?, and MORE IMPORTANTLY,
- how best to adjust jr, and is there a better approach to this problem
using Mathematica?
Thanks for any help,
Randy Thomas, Ph.D.
INSERM U.467, Necker Medical Faculty Paris France
srthomas@necker.fr