Re: Mixed Algebraic/differential system
- To: mathgroup@smc.vnet.net
- Subject: [mg11512] Re: [mg11420] Mixed Algebraic/differential system
- From: "C. Woll" <carlw@u.washington.edu>
- Date: Fri, 13 Mar 1998 12:22:02 -0500
Hi Randy, Is it possible to cast your mixed algebraic/differential system into a pure differential system as follows? Let y1, y2, y3 be the unknown functions whose derivatives you know, and z be the unknown function whose derivative you don't know. Thus, you have the system of equations y1' = f(y1,y2,y3,z,t) y2' = g(y1,y2,y3,z,t) y3' = h(y1,y2,y3,z,t) F(y1,y2,y3) = 0 Then taking a total derivative of your constraint function dF = 0 = dF/dy1 f + dF/dy2 g + dF/dy3 h = G(y1,y2,y3,z,t) where G is defined as above. Then, taking another derivative of G gives dG = 0 = dG/dy1 f + dG/dy2 g + dG/dy3 h + dG/dt + dG/dz z' Solving the above equation for z' gives z' = H(y1,y2,y3,z,t) Now you have a pure system of differential equations. The only problem is that you don't have an initial condition for z. Perhaps you already know what this initial condition is, or perhaps you can determine it from your original equations somehow. Once you have a pure system of differential equations, you can use NDSolve or DSolve directly. Note that Mathematica can solve stiff differential equations if that becomes a problem. Perhaps you have already considered this approach and rejected it for some reason, or maybe the approach is somehow flawed. Anyway, good luck! Carl Woll Dept of Physics U of Washington On Thu, 12 Mar 1998, SR Thomas wrote: > 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 > > >