>    restart;
with(plots):
with(inttrans):
alias(Y(s)=laplace(y(t),t,s)): 

Warning, the name changecoords has been redefined

Cold pills: Periodic impulses and Laplace transforms

About this worksheet

Author and Date

Matthias Kawski
Arizona State University
http://math.asu.edu/~kawski

>   

Original version MAPLE release 8.
June 2007.

>   

Content, Purpose and Use

This worksheet illustrates a naive applications of Laplace transforms to solve
linear differential equations driven by a periodic pulses -- here interpreted as
taking medicine.

>   

A good part of the worksheet experiments with MAPLE's capabilities of

working with infinite series, or simply long series --- but in the end we take
a pragmatic approach, and type in the expected solution by hand when MAPLE's

invlaplace balks.
Suitable for demos, inviting simple modifications ... and some questions such as
how will the final level change when one doubles the daily dose, or takes twice

the dose every other day. The last sections contains a quick analysis of the
quasi-steady state as a function of dosage and schedule.

>   

Updates and log of modifications.

none yet..

Starting from the differential equation

A finite number of days in the treatment

>    de0a:=diff(y(t),t)=sum(Dirac(t-'k'),'k'=1..7):de0a;
AE0a:=subs(y(0)=0,laplace(de0a,t,s)):
AE0a;
Y(s)=solve(AE0a,Y(s));

diff(y(t),t) = Dirac(t-1)+Dirac(t-2)+Dirac(t-3)+Dirac(t-4)+Dirac(t-5)+Dirac(t-6)+Dirac(t-7)

s*Y(s) = exp(-s)+exp(-2*s)+exp(-3*s)+exp(-4*s)+exp(-5*s)+exp(-6*s)+exp(-7*s)

Y(s) = (exp(-s)+exp(-2*s)+exp(-3*s)+exp(-4*s)+exp(-5*s)+exp(-6*s)+exp(-7*s))/s

The infinite treatment

>    de0:=diff(y(t),t)=sum(Dirac(t-'k'),'k'=1..infinity):de0;
AE0:=subs(y(0)=0,laplace(de0,t,s)):
AE0;
Y(s)=solve(AE0,Y(s));

diff(y(t),t) = sum(Dirac(t-k),k = 1 .. infinity)

s*Y(s) = sum(exp(-s*k)-exp(-s*k)*Heaviside(-k)+PIECEWISE([0, k <= 0],[-exp(-s*k), 0 < k]),k = 1 .. infinity)

Y(s) = sum(exp(-s*k)-exp(-s*k)*Heaviside(-k)+PIECEWISE([0, k <= 0],[-exp(-s*k), 0 < k]),k = 1 .. infinity)/s

Adding decay (excretion) to the model

>    de1:=diff(y(t),t)-y(t)=sum(Dirac(t-'k'),'k'=1..7):
de1;
AE1:=subs(y(0)=0,laplace(de1,t,s)):
AE1;
Y(s)=solve(AE1,Y(s));

diff(y(t),t)-y(t) = Dirac(t-1)+Dirac(t-2)+Dirac(t-3)+Dirac(t-4)+Dirac(t-5)+Dirac(t-6)+Dirac(t-7)

s*Y(s)-Y(s) = exp(-s)+exp(-2*s)+exp(-3*s)+exp(-4*s)+exp(-5*s)+exp(-6*s)+exp(-7*s)

Y(s) = (exp(-s)+exp(-2*s)+exp(-3*s)+exp(-4*s)+exp(-5*s)+exp(-6*s)+exp(-7*s))/(s-1)

>   

>   

Restart in the frequency domain (sign error)

This is what the geometric series (almost) correctly sums to ... did not want to bother

to try to get MAPLE do it to above formula (but made a sign mistake)

>    Y1:=1/s/(1+exp(-s));
y1:=invlaplace(Y1,s,t);

Y1 := 1/(s*(1+exp(-s)))

y1 := 1/2-1/2*(-1)^(-1-floor(t))

Just for illustration, we make MAPLE expand the periodic delays into a series
-- the sign erro becomes obvious here.

>    Y1=subs(exp(-s)=x,Y1);
series(subs(exp(-s)=x,Y1),x=0,12);
Y1=simplify(subs(x=exp(-1),series(subs(exp(-s)=x,Y1),x=0,7)));

1/(s*(1+exp(-s))) = 1/(s*(1+x))

series(1/s+(-1/s)*x+1/s*x^2+(-1/s)*x^3+1/s*x^4+(-1/s)*x^5+1/s*x^6+(-1/s)*x^7+1/s*x^8+(-1/s)*x^9+1/s*x^10+(-1/s)*x^11+O(x^12),x,12)

1/(s*(1+exp(-s))) = (1-exp(-1)+exp(-2)-exp(-3)+exp(-4)-exp(-5)+exp(-6)+O(exp(-7))*s)/s

>    display([
plot(y1,t=0..12,-3..3,thickness=3,discont=true),
plot(y1,t=0..12,-3..3,thickness=1)
],xtickmarks=2,ytickmarks=3);

[Maple Plot]

Of course the pills always add, never subtract medicine in the body

>   

Restart in the frequency domain with correct sign

This is what the geometric series  correctly sums to ... did not want to bother

to try to get MAPLE do it to above formula

>    Y2:=1/s/(1-exp(-s));
y2:=invlaplace(Y2,s,t);

Y2 := 1/(s*(1-exp(-s)))

y2 := ceil(t)

Just for illustration, we make MAPLE expand the periodic delays into a series

>    Y2=subs(exp(-s)=x,Y2);
series(subs(exp(-s)=x,Y2),x=0,12);
Y2=simplify(subs(x=exp(-s),series(subs(exp(-s)=x,Y2),x=0,7)));

1/(s*(1-exp(-s))) = 1/(s*(1-x))

series(1/s+1/s*x+1/s*x^2+1/s*x^3+1/s*x^4+1/s*x^5+1/s*x^6+1/s*x^7+1/s*x^8+1/s*x^9+1/s*x^10+1/s*x^11+O(x^12),x,12)

1/(s*(1-exp(-s))) = (1+exp(-s)+exp(-2*s)+exp(-3*s)+exp(-4*s)+exp(-5*s)+exp(-6*s)+O(exp(-7*s))*s)/s

>    display([
plot(y2,t=0..12,-2..12,thickness=3,discont=true),
plot(y2,t=0..12,-2..12,thickness=1)
],xtickmarks=2,ytickmarks=3);

[Maple Plot]

This is the correct picture for every day adding another dose, none ever metabolized/excreted

>   

Now add some decay (metabolization / excretion)

>    Y3:=1/(s+1)/(1-exp(-s));
y3:=invlaplace(Y3,s,t);

Y3 := 1/((s+1)*(1-exp(-s)))

y3 := invlaplace(1/((s+1)*(1-exp(-s))),s,t)

>    Y3=subs(exp(-s)=x,Y3);
series(subs(exp(-s)=x,Y3),x=0,8);
Y3a:=convert(
     simplify(subs(x=exp(-s),series(subs(exp(-s)=x,Y3),x=0,11))),
         polynom);

1/((s+1)*(1-exp(-s))) = 1/((s+1)*(1-x))

series(1/(s+1)+1/(s+1)*x+1/(s+1)*x^2+1/(s+1)*x^3+1/(s+1)*x^4+1/(s+1)*x^5+1/(s+1)*x^6+1/(s+1)*x^7+O(x^8),x,8)

Y3a := (1+exp(-s)+exp(-2*s)+exp(-3*s)+exp(-4*s)+exp(-5*s)+exp(-6*s)+exp(-7*s)+exp(-8*s)+exp(-9*s)+exp(-10*s))/(s+1)

>    y3a:=invlaplace(Y3a,s,t);

y3a := exp(-t)+Heaviside(t-1)*exp(-t+1)+Heaviside(t-2)*exp(-t+2)+Heaviside(t-3)*exp(-t+3)+Heaviside(t-4)*exp(-t+4)+Heaviside(t-5)*exp(-t+5)+Heaviside(t-6)*exp(-t+6)+Heaviside(t-7)*exp(-t+7)+Heaviside(t...
y3a := exp(-t)+Heaviside(t-1)*exp(-t+1)+Heaviside(t-2)*exp(-t+2)+Heaviside(t-3)*exp(-t+3)+Heaviside(t-4)*exp(-t+4)+Heaviside(t-5)*exp(-t+5)+Heaviside(t-6)*exp(-t+6)+Heaviside(t-7)*exp(-t+7)+Heaviside(t...
y3a := exp(-t)+Heaviside(t-1)*exp(-t+1)+Heaviside(t-2)*exp(-t+2)+Heaviside(t-3)*exp(-t+3)+Heaviside(t-4)*exp(-t+4)+Heaviside(t-5)*exp(-t+5)+Heaviside(t-6)*exp(-t+6)+Heaviside(t-7)*exp(-t+7)+Heaviside(t...

>    display([
plot(y3a,t=0..10,-2..5,thickness=3,discont=true),
plot(y3a,t=0..10,-2..5,thickness=1)
],xtickmarks=2,ytickmarks=3);

[Maple Plot]

>   

A slight improvement, utilize series, add parameters

Metabolization rate beta, period between taking pills, and amounts taken each time.

>    beta:=-1/5;
p:=1;        
a:=1;   

beta := -1/5

p := 1

a := 1

>    a:='a';
p:='p';
beta:='beta';

a := 'a'

p := 'p'

beta := 'beta'

>    Y40:=a/(s-beta)*sum(exp(-'k'*p*s),'k'=1..infinity);

Y40 := a/(s-beta)/(-1+exp(p*s))

Showing the first few terms for illustration

>    Y4a:=a/(s-beta)*sum(exp(-'k'*p*s),'k'=1..7);
Y4:=a/(s-beta)*sum(exp(-'k'*p*s),'k'=1..25):

Y4a := a/(s-beta)*(exp(-p*s)+exp(-2*p*s)+exp(-3*p*s)+exp(-4*p*s)+exp(-5*p*s)+exp(-6*p*s)+exp(-7*p*s))

>    y4a:=invlaplace(Y4a,s,t);
y4:=invlaplace(Y4,s,t):

y4a := a*(Heaviside(t-p)*exp(beta*(t-p))+Heaviside(t-2*p)*exp(beta*(t-2*p))+Heaviside(t-3*p)*exp(beta*(t-3*p))+Heaviside(t-4*p)*exp(beta*(t-4*p))+Heaviside(t-5*p)*exp(beta*(t-5*p))+Heaviside(t-6*p)*exp...
y4a := a*(Heaviside(t-p)*exp(beta*(t-p))+Heaviside(t-2*p)*exp(beta*(t-2*p))+Heaviside(t-3*p)*exp(beta*(t-3*p))+Heaviside(t-4*p)*exp(beta*(t-4*p))+Heaviside(t-5*p)*exp(beta*(t-5*p))+Heaviside(t-6*p)*exp...
y4a := a*(Heaviside(t-p)*exp(beta*(t-p))+Heaviside(t-2*p)*exp(beta*(t-2*p))+Heaviside(t-3*p)*exp(beta*(t-3*p))+Heaviside(t-4*p)*exp(beta*(t-4*p))+Heaviside(t-5*p)*exp(beta*(t-5*p))+Heaviside(t-6*p)*exp...

My MAPLE does not like these many terms -- thus let us define the solution by hand.

>    y4:=a*sum(Heaviside(t-'k'*p)*exp(beta*(t-'k'*p)),k=1..7);
check=y4-y4a;

y4 := a*(Heaviside(t-p)*exp(beta*(t-p))+Heaviside(t-2*p)*exp(beta*(t-2*p))+Heaviside(t-3*p)*exp(beta*(t-3*p))+Heaviside(t-4*p)*exp(beta*(t-4*p))+Heaviside(t-5*p)*exp(beta*(t-5*p))+Heaviside(t-6*p)*exp(...
y4 := a*(Heaviside(t-p)*exp(beta*(t-p))+Heaviside(t-2*p)*exp(beta*(t-2*p))+Heaviside(t-3*p)*exp(beta*(t-3*p))+Heaviside(t-4*p)*exp(beta*(t-4*p))+Heaviside(t-5*p)*exp(beta*(t-5*p))+Heaviside(t-6*p)*exp(...
y4 := a*(Heaviside(t-p)*exp(beta*(t-p))+Heaviside(t-2*p)*exp(beta*(t-2*p))+Heaviside(t-3*p)*exp(beta*(t-3*p))+Heaviside(t-4*p)*exp(beta*(t-4*p))+Heaviside(t-5*p)*exp(beta*(t-5*p))+Heaviside(t-6*p)*exp(...

check = 0

OK, then let us contionue this for a long time....

>    y4:=a*sum(Heaviside(t-'k'*p)*exp(beta*(t-'k'*p)),k=1..100):

>    y41:=subs(beta=-1/5,p=1,a=1,y4):

>    display([
plot(y41,t=0..10,thickness=3,discont=true),
plot(y41,t=0..10,thickness=1)
],view=-1..12,xtickmarks=2,ytickmarks=3,
title=`One pill daily. How high will it go.`,
titlefont=[TIMES,BOLD,14]);

[Maple Plot]

>    display([
plot(y41,t=0..20,thickness=3,discont=true),
plot(y41,t=0..20,thickness=1)
],view=-1..12,xtickmarks=2,ytickmarks=3,
title=`One pill every day. Metabolization rate 1/5.`,
titlefont=[TIMES,BOLD,14]);

[Maple Plot]

>    y42:=subs(beta=-1/5,p=2,a=2,y4):

>    display([
plot(y41,t=0..20,thickness=3,discont=true),
plot(y41,t=0..20,thickness=1),
plot(y42,t=0..20,thickness=3,color=blue,discont=true),
plot(y42,t=0..20,thickness=1,color=blue)
],view=-1..12,xtickmarks=2,ytickmarks=3,
title=`Two pills every other day. Metabolization rate 1/5.`,
titlefont=[TIMES,BOLD,14]);

[Maple Plot]

>    y43:=subs(beta=-1/5,p=1,a=2,y4):

>    display([
plot([y41,y43],t=0..40,thickness=3,color=[red,blue],
      discont=true,numpoints=600),
plot([5,10,y41,y43],t=0..40,thickness=1,color=[black,black,red,blue])
],view=-1..12,xtickmarks=2,ytickmarks=3,
title=`Double daily dosage -- new level? Metabolization rate 1/5.`,
titlefont=[TIMES,BOLD,14]);

[Maple Plot]

>   

Steady state analysis

A quasi steady-state is characterized by the amount that is metabolized in one time
period is equal to the new dose added. This means that the quasi-steady state ys
must satisfy the following condition:

>    sseq:=ys-ys*exp(-beta*p)=a:
sseq;

ys-ys*exp(-beta*p) = a

>    yssol:=ys=solve(sseq,ys):
yssol;

ys = -a/(-1+exp(-beta*p))

In our case we considered thefixed metabolization rate

>    mrate:=beta=1/5:
mrate;

beta = 1/5

and it is clear that the quasi steady state depends linearly on the dosage,
i.e. we should cosider the ratio

>    ssp:=subs(mrate,yssol/a):
ssp;

1/a*ys = -1/(-1+exp(-1/5*p))

>    display([
   plot([rhs(ssp),rhs(ssp)-1],p=1/5..3,color=[red,blue]),
   seq(plot([[T,0],[T,5/T],[0,5/T]],0..T,color=gray),
       T=[1/2,1,2,3])
   ],view=[0..3,0..30]);

[Maple Plot]

>   

>