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