Re: Problem with a system of equations describing an exposure to lead...

*To*: mathgroup at smc.vnet.net*Subject*: [mg50386] Re: [mg50375] Problem with a system of equations describing an exposure to lead...*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Wed, 1 Sep 2004 01:49:26 -0400 (EDT)*References*: <200408311028.GAA18673@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Richard Palmer wrote: > I have a system of equations > > Eqns={u[t]==a bl[t-1], > bo[t]==b bl[t-1]-c bo[t-1], > bl[t]==(1-a-b) bl[t-1]-c bo[t-1], > u[0]==0, > bo[0]==0, > bl[0]==pb} > > They are supposed to represent a system of what happens when a person > suffers a one-time exposure to lead (pb). At any given time, the person > will pass some portion of the lead that is in the blood and soft tissue > (bl[t]) through the urine (u[t]) while some portion moves to the bone > (bo[t]) and a small portion moves from the bone back into the blood. This > is easily solved with RSolve. However, I am concerned the solution may not > be correct: I believe that for any time t, the amount of lead in the bone > and blood less all lead passed to date through the urine should equal the > level of exposure (pb). Taking the Rsolve solutions and making a table of > u[t], bl[t], and bo[t] for times 0-5 will show the problem (substitute > random fractions for a,b, and c and take pb to be 1). Is the mistake in the > solution or in my formulation of the model? I am using Mathematica 5. > > Thanks in advance. There are a few problems. First, you want +c, not -c, in the equation for bl[t]. Also you want (1-c)*bo[t-1] rather than -c*bo[t-1] in the equation for bo[t]. Last, your sanity check needs to work with total passed through urine up to time t rather than an "instantaneous" (or rather, discrete time approximation) rate, hence it needs to integrate u[tau] to time t. Here is what I did to test it. a = 1/5; b = 1/7; c = 1/11; pb = 1; eqns = {u[t]==a*bl[t-1], bo[t]==b*bl[t-1]+(1-c)*bo[t-1], bl[t]==(1-a-b)*bl[t-1]+c*bo[t-1], u[0]==0, bo[0]==0, bl[0]==pb}; InputForm[soln = First[RSolve[eqns, {bo[t],bl[t],u[t]}, t]]] Out[85]//InputForm= {u[t] -> (3^(-2 + t)*5^(-1 - t)*7^(1 - t)* (-165*Sqrt[1901]*(201 - Sqrt[1901])^t + 165*Sqrt[1901]* (201 + Sqrt[1901])^t + 1901*(201 - Sqrt[1901])^t*UnitStep[-1 + t] + 201*Sqrt[1901]*(201 - Sqrt[1901])^t*UnitStep[-1 + t] + 1901*(201 + Sqrt[1901])^t*UnitStep[-1 + t] - 201*Sqrt[1901]*(201 + Sqrt[1901])^t*UnitStep[-1 + t]))/(1901*22^t), bo[t] -> -(((55/3)^(1 - t)*((201 - Sqrt[1901])^t - (201 + Sqrt[1901])^t))/ (14^t*Sqrt[1901])), bl[t] -> (2^(-1 - 2*t)*(-603*(70*(603 - 3*Sqrt[1901]))^t + 3*Sqrt[1901]*(70*(603 - 3*Sqrt[1901]))^t + 7^(1 + t)*10^(2 + t)* (3*(201 - Sqrt[1901]))^t - 7^(1 + t)*10^(2 + t)* (3*(201 + Sqrt[1901]))^t + 67*3^(2 + t)*(70*(201 + Sqrt[1901]))^t + 3^(1 + t)*Sqrt[1901]*(70*(201 + Sqrt[1901]))^t))/(3*Sqrt[1901]*13475^t)} uu[t_] := Evaluate[u[t]/.soln] bbo[t_] := Evaluate[bo[t]/.soln] bbl[t_] := Evaluate[bl[t]/.soln] utot[t_] := NIntegrate[uu[tau], {tau,0,t}] In[90]:= Table[N[{bbo[t],bbl[t],utot[t],bbo[t]+bbl[t]+utot[t]}], {t,0,10}] Out[90]= {{0., 1., 0., 1.}, {0.142857, 0.657143, 0.10913, 0.90913}, {0.223748, 0.444824, 0.272105, 0.940677}, {0.266953, 0.312653, 0.380621, 0.960227}, {0.28735, 0.229726, 0.455336, 0.972412}, {0.294045, 0.177086, 0.508939, 0.98007}, {0.292612, 0.143102, 0.54923, 0.984943}, {0.286454, 0.12064, 0.581008, 0.988101}, {0.277647, 0.105319, 0.607233, 0.990198}, {0.267452, 0.09445, 0.629736, 0.991637}, {0.256631, 0.0863809, 0.649655, 0.992666}} It is encouraging that totals seem to head back toward 1, but probably there is still something wrong insofar as they are not approximately constant. Offhand I do not know whether this is a problem in the formulation or in RSolve. Daniel Lichtblau Wolfram Research