The first term , ~I is independent of the concentration. That is why it is called a mixed case. With appropriate value of I , C and L this may represent a cell from which a susbtance diffuses out while theres is a constant input I.

Dividing by V and taking the limit as ∆t goes to zero gives the DE

dC(t)/dt = I/V –(L/V) C(t) . (2)

Thre cases are posed.

First the steady state, when dC/dt becomes zero, has the solution ,

C_{steady
state} = I/V (3)

Second suppose the initial concentration C(0)=0. Solve (2) using the Laplace transform. From the previous chapters we anticipate that the answer is an inverted exponential.

From

L { dC(t)/dt = I/V –(L/V) C(t) }

sc(s) –C(0) = I/(s*V) – (L/V) c(s) ,

we find

c(s) = (I/V)* (1/(s*(s+L/V)))

Using Matlab

>> syms I V L s t;

>> cs=(I/V)* (1/(s*(s+L/V)));

>> ilaplace(cs)

ans =

I/L*(1-exp(-L*t/V))

So C(t) = (I/L) { 1 – exp(-Lt/V) }

Let I=.01 moles/minute , L=0.2 liters/minute , C(0)=0.03 moles/liter,

V=10 liters .

I=.01 ; L=0.2 ; c0=.3 ;V=10 ;

tscale=V/L ;tf=3*tscale ;

t=[0:tf/200:tf];

ct=(I/L)*( 1 - exp(-L*t/V) );

>> plot(t,ct),xlabel('t~minutes'), ylabel('c~moles/liter')

Fig 5.1 . An inverted exponential results for C(t).

Lecture 3. The Two-Compartment Problem

Fig 5.2 The two compartment problem.

The meaning and dimensions of the terms are ,

C₁(t) =Q₁/V₁ , C₂(t) = Q₂/V₂ ~moles/volume

I₁ and I₂ ~ moles/time are inputs to the compartments

K(C₁-C₂) ~ diffusion between compartments

~ (volume/time)*(moles/volume)

L₁C₁ , L₂ C₂ ~ concentration dependent loss

~ (volume/time)*(moles/volume)

In an interval ∆t the net changes in Q₁ and Q₂ are ,

∆ Q₁ = ∆t { I₁ - L₁C₁ - K(C₁-C₂) } (4)

∆ Q₂ = ∆t { I₂ – L₂C₂ + K(C₁-C₂) } (5)

Notice that the last term in (5) must have the opposite sign of the same term appearing in(4).

Dividing (4) by ∆t V₁ and (5) by ∆t V₂ give the two DE,

dC₁/dt = (1/V₁)*{ I₁ - L₁C₁ - K(C₁-C₂) } (6)

d C₂/dt = (1/V₂)* { I₂ – L₂C₂ + K(C₁-C₂) } (7)

We solve (6) & (7) using finite differences with the following data and initial conditions (IC). This is not an example from physiology. The values of the parameters are chosen to illustrate a solution.

L ~liters/minute , K~liters/minute , I ~ moles/minute , C~moles/liter,

V~liters

I1=.01;I2=.02 ;L1=.1 ;L2=.2;Kd=.5*L1;

c1(1)=0.20 ; c2(1)=0.1 ;

V1=10 ; V2= 10 ;

Fig 5.3 Concentrations C₁ , C₂ . C₁ is the upper curve.

% chapter5 sec3

I1=.01;I2=.02 ;L1=.1 ;L2=.2;Kd=.5*L1;

c1(1)=0.20 ; c2(1)=0.1 ;

V1=10 ; V2= 10 ;

A=[V1/L1 V1/Kd V2/L2 V2/Kd];

tscale=min(A);tlarge=max(A);

tf=tlarge; dt=tscale/50;

k=0;

for t=[dt:dt:tf];

k=k+1;

g1=(I1-L1*c1(k)-Kd*(c1(k)-c2(k)))/V1;

g2=(I2-L2*c2(k)-Kd*(c2(k)-c1(k)))/V2;

c1(k+1)=c1(k)+dt*g1;

c2(k+1)=c2(k)+dt*g2;

end

t=[0:dt:tf];

plot(t,c1,t,c2),xlabel('t~minutes'),...

ylabel('c1,c2~moles/liter')

We work next a case proposed in page 89 of Reference 1.

Consider the case I₂ =0 , C₁(0)=C₂(0)=0.

MATLAB code

% chapter5 sec3

I1=.01; I2=.0 ;L1=.1 ;L2=.2; Kd=.5*L1;

c1(1)=0.0 ; c2(1)=0.0 ;

V1=10 ; V2= 10 ;

A=[V1/L1 V1/Kd V2/L2 V2/Kd];

tscale=min(A);tlarge=max(A);

tf=tlarge; dt=tscale/50;

k=0;

for t=[dt:dt:tf];

k=k+1;

g1=(I1-L1*c1(k)-Kd*(c1(k)-c2(k)))/V1;

g2=(I2-L2*c2(k)-Kd*(c2(k)-c1(k)))/V2;

c1(k+1)=c1(k)+dt*g1;

c2(k+1)=c2(k)+dt*g2;

end

t=[0:dt:tf];

plot(t,c1,t,c2),xlabel('t~minutes'),...

ylabel('moles/liter')

Fig. 5.4 Concentrations C₁ , C₂ .

We work now a second case proposed in page 91 of Reference 1.

Let K=3*(L₁ + L₂) and L₂ << L₁ . The system will behave as if it were a single compartment.

Data I₁=.1 , I₂ =0 , C₁(0)=C₂(0)=0.

L1=.2 ;L₂ =.05=.2 ; K=3*(L1+L2) , ,V₁ =V₂ =5

The units are again, L ~liters/minute , K~liters/minute , I ~ moles/minute , C~moles/liter, V~liters.

Fig 5.5 Concentration of C₁ , C₂ with a large diffusion from the first compartment to the second.

MATLAB code

I1=.1; I2=.0 ; L1=.2 ; L2=.05 ;

Kd=3*(L1+L2); c1(1)=0.0 ; c2(1)=0.0 ;

V1=5 ; V2= 5 ;

A=[V1/L1 V1/Kd V2/L2 V2/Kd];

tscale=min(A);tlarge=max(A);

tf=tlarge; dt=tscale/50;

k=0;

for t=[dt:dt:tf];

k=k+1;

g1=(I1-L1*c1(k)-Kd*(c1(k)-c2(k)))/V1;

g2=(I2-L2*c2(k)-Kd*(c2(k)-c1(k)))/V2;

c1(k+1)=c1(k)+dt*g1;

c2(k+1)=c2(k)+dt*g2;

end

t=[0:dt:tf];

plot(t,c1,t,c2),xlabel('t~minutes'),...

ylabel('moles/liter')

Section 4. A simple three compartment problem

Fig 5.6. Three compartments with a constant input I₃ at #3 and with outflow

L C₃(t) .There is diffsusion between compartments.

Let C~moles/minute , V ~liters , L~ liters/minute , K~liters/minute ,

and I ~ moles/minute.

The system is described by three coupled first order linear DE.

V₁ dC₁/dt = -K₁₂ ( C₁ – C₂) (8)

liter*(moles/(liter-miniute)) = (liter/minute)*(moles/liter)

Compartment #2 has an incoming diffusion + K₁₂ ( C₁ – C₂) and the

outgoing diffusion -K₂₃ (C₂ – C₃) . The DE is

V₂ dC₂/dt = + K₁₂ ( C₁ – C₂) -K₂₃ (C₂ – C₃) . (9)

The parameters K₁₂ , K₂₃ , need not be the same.

Compartment #3 has a steady input I₃ , an ingoing diffusion

+K₂₃ (C₂ – C₃) and the out flow –LC₃ . The DE is

V₃ dC₃/dt = I₃ + K₂₃ (C₂ – C₃) -L C₃ . (10)

The following data is used for the next plot.

Data

I3=.2; L3=.05 ; ;K12=.01 ;K23=.01;

c1(1)=0.50 ; c2(1)=0.;c3(1)=0. ;

V1=5 ; V2= 5 ;V3=5;

MATLAB

% chapter5 sec4 Three compartment

I3=.2; L3=.05 ; ;K12=.01 ;K23=.01;

c1(1)=0.50 ; c2(1)=0.;c3(1)=0. ;

V1=5 ; V2= 5 ;V3=5;

A=[V1/L3 V1/K12 V1/K23 V2/K23 V2/L3 V3/L3];

tscale=min(A);tlarge=max(A);

tf=tlarge; dt=tf/3000;

k=0;

for t=[dt:dt:tf];

k=k+1;

g1=-K12*(c1(k)-c2(k))/V1;

g2=(K12*(c1(k)-c2(k))-K23*(c2(k)-c3(k)))/V2;

g3=(I3-L3*c3(k)+K23*(c2(k)-c3(k) ) )/V3;

c1(k+1)=c1(k)+dt*g1;

c2(k+1)=c2(k)+dt*g2;

c3(k+1)=c3(k)+dt*g3;

end

t=[0:dt:tf];

plot(t,c1,t,c2,t,c3),xlabel('t~minutes'),...

ylabel('c1,c2,c3~moles/liter')

Fig 5.6 Concentrations C₁ , C₂ , C₃ in the three compartment problem.

Section 5. The delta impulse function

One method for solving inhomogeneous DE is that of the Green’s function.

We explain the method through a concrete example

Let the DE be

d² C/dt² + dC/dt + C = F(t) , t ≥ 0 (11)

with F(t) = cos(2t) and IC C(0)=1 , dC(0)/dt=0.

The complete solution of (11) has two parts , the particular solution C_p(t)

and the complementary or homogeneous solution C_h(t) .

C(t) = C_p(t) + C_h(t) . (12)

The homogenous part satisfies

d² C_h /dt² + dC_h /dt + C_h = 0 and the particular solution satisfies

d² C_p /dt² + dC_p/dt + C_p = F(t) .

The Green function satisfies the DE with a delta function on the right hand side

d² G /dt² + dG/dt + G= δ (t) . (13)

δ (t) is Dirac delta function .Some of its properties are

∫ δ (t) dt = 1 , - ∞ <t < ∞ ,

∫ F(t) δ (t-a) dt = F(a) , - ∞ <t < ∞ .

Once G(t) is obtained , the particular solution is given by the integral

C_p(t) = ∫^t ₀ G(t-t’) F(t’) dt’

The complete solution can be written as

C(t) = ∫^t ₀ G(t-t’) F(t’) dt’ + C_h(t) . (14)

Taking the derivatives upon (14) we verify that

d² C(t) /dt² + dC(t)/dt + C(t) =

∫ {d² G(t-t’) /dt² + dG(t-t’)/dt + G(t-t’)} F(t’) dt’

+ d² C_p(t) /dt² + dC_p(t)/dt + C_h(t)

= ∫ δ(t-t’) F(t’) dt’ + 0 = F(t) . (15)

Example :

We solve (11) in stages. First find the homogeneous solution.

Part 1. Solution to the homogenous equation

d² C/dt² + dC/dt + C =0 with IC , C(0)=1, dC(0)/dt=0 is worked using the

Laplace transform.

L (d² C/dt² + dC/dt + C) = s² c(s) –sC(0)-dC(0)/dt +sc(s)-C(0) +c(s)=0 .

c(s) = (s+1)/(s² +s +1)

MATLAB code

>> syms s t;

>> cs=(1+s)/(s^2+s+1);

>> ilaplace(cs)

gives

C(t) = exp(-1/2*t)*cos(1/2*3^(1/2)*t)+1/3*exp(-1/2*t)*3^(1/2)*sin(1/2*3^(1/2)*t)

t=[0:.05:10];

% Gcomp is the complementary solution with IC g(0)=1

% dG(0)/dt=0

Gcomp=exp(-t./2).*cos(sqrt(3)*t./2) + ...

sqrt(3)*(1/3)*exp(-t./2).*sin(sqrt(3)*t./2);

plot(t,Gcomp),xlabel('t'),ylabel('Gcomplementary')

Fig 5.7 . Complementary solution C _comp. ( or C_homogeneous ) , i.e. the solution of the homogeneous DE.

Part 2. Solution to the delta function impulse

d² G/dt² + dG/dt + G = δ(t) , with G(0)=0 , dG(0)/dt=0

s² g(s) –sG(0)-dG(0)/dt + sg(s) –G(0) + g(s) = 1 ( 16 )

g(s) = 1/(s² + s + 1) ( 17 )

MATLAB

syms s t;

>> f=1/(s^2+s+1);

>> ilaplace(f)

ans = 2/3*3^(1/2)*exp(-1/2*t)*sin(1/2*3^(1/2)*t)

G(t)=(2/3^1/2)exp(-t/2) sin(3^1/2 t/2)

and

G(t-t’) = (2/3^1/2)exp(-(t-t’)/2) sin(3^1/2 (t-t’)/2)

We plot the solution to the delta function impulse , the Green function for the DE.

MATLAB

t=[0:.05:10];

G=(2/3^(1/2))*exp(-t./2).* sin(3^(1/2)* t./2);

plot(t,G),xlabel('t'),ylabel('G')

Fig 5.8 . Green’s function.

As was said before, the solution to

d² C/dt² + dC/dt + C = F(t) t ≥ 0 is

C(t) = ∫ ^t ₀ G(t-t’) F(t’) dt’ + C(t)_{complementary}

The following FORTRAN code calculates the integral

∫ ^t ₀ G(t-t’) F(t’) dt’ and adds C(t)_{complementary} . The plot is given .

Fig 5.9 Plot of C(t) = ∫ ^t ₀ G(t-t’) F(t’) dt’ + C(t)_{complementary}

FORTRAN code

c Use of Green's function to solve inhomogenous DE 22 nov 2008

data ti ,tf ,nt, nstep/0.,10.,70, 5000/

c g(u=t-t') is the greens'function with L g =delta(t) g(0)=0,dg/dt=0.

g(u)=(2./3.)*sqrt(3.)*exp(-u/2.)*sin(sqrt(3.)*u/2.)

c complentary solutio to L C = 0 , C(0)=1. , dC/dt=0.

comp(t)=exp(-t/2.)*cos(sqrt(3.)*t/2.)+

$ (1./3.)*exp(-t/2.)*sqrt(3.)*sin(sqrt(3.)*t/2.)

c f(t) is the driving force

f(t)=cos(2.*t)

dt=(tf-ti)/float(nt)

do 10 it=0,nt

t=ti+dt*float(it)

dtprime=t/float(nstep)

du=dtprime

sum=0.

do 20 j=1,nstep

tprime=dtprime*float(j)

u=t-tprime

sum=sum+(dtprime/2.)*(g(u)*f(tprime)+g(u-du)*f(tprime-du))

20 continue

ctotal=sum+comp(t)

print 100 ,t,sum,comp(t),ctotal

10 continue

100 format(1x,'t,cinh,comp,c(t)=',4(3x,e10.3))

stop

end

We end this section by solving the problem by the finite difference method.

The generic solution to d² C/dt² + dC/dt + C = F(t) is simply

C (t_n ) = 2 C(t_n-1) – C(t_n-2) + (∆t)² { - (C(t_n-1)-C(t_n))/ ∆t –C(t_n-1) + F(t_n-1) }

The IC is C(0)=1 , dC(0)/dt=0 which is equivalent to C(t₀) =1 , C(t₁) =1 .

The solution is plotted.

Fig 5.10 . C(t) obtained by finite differences. It is identical to fig 5.9.

Fortran code Solution by finite differences

c23 nov 2008 equation d^2c/dt^2+dc/dt+c= cos(2t) IC c(0)=1 ,dc(0)/dt=0.

c solution by finite differences

data ti, tf , nstep/0.,10.,4000 /

data c0,c1/1.,1./

f(t)= cos(2.*t)

a(t)=-(c1-c0)/dt -c1 +f(t)

dt=(tf-ti)/float(nstep)

kp=int(float(nstep)/70.)

kount=kp

print 100,0.,c0

do 10 i=2,nstep

t=ti+dt*float(i)

c2=2.*c1-c0+dt**2*a(t-dt)

if(i.eq.kount)then

print 100,t,c2

kount=kount+kp

endif

c0=c1

c1=c2

10 continue

print 100,t,c2

100 format(1x,'t,C(t)=',2(4x,e10.3))

stop

end