Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Project: pde
Views: 136

Delay Differential Equations: The method of steps.

Example: Water is running at qq liters/minute into the top of a tank intially containing KK liters of brine with pp grams of salt per liter. It runs out at the bottom of the tank at the same rate. Because of the delay in the diffusion of the salt in the water, the concentration of salt at the outflow is the average average concentration of the salt in the tank at time t−t0t - t_0 where t0t_0 is a fixed delay in minutes. Find the average concentration C(t)C(t) of salt in the tank at time t>=0t >=0 .
First, let's speculate on the qualitative structure of the concentration function. We would expect the concentration to be decreasing over time with 00 as a limiting value, depending on the values of the parameters.
Set up: Let S(t)S(t) be the total amount of salt in the tank in grams, C(t)C(t) be the average concentration of salt in the tank in g/l (grams per liter), and Cout(t)C_{out}(t) be the concentration of salt at the outflow. So S(t)=K C(t)S(t)=K\,C(t), C(0)=pC(0)=p, S(0)=K pS(0)=K\,p, and S′(t)=−q C(t−t0)S^\prime(t)=-q\,C(t-t_0), from which we have C′(t)=1K S′(t)=−qK C(t−t0)C^\prime(t)=\dfrac{1}{K}\,S^\prime(t)=-\dfrac{q}{K}\,C(t-t_0), a delay differential equation (DDE).

We need an initial condition: The natural assumption is that C(t)C(t) is constant for t∈[−t0,0]t\in [-t_0,0] at the original concentration pp.
Now the task is to solve the DDE C′(t)=a C(t−t0) C^\prime(t)=a\,C(t-t_0) with initial condition C(t)=pC(t)=p for t∈[−t0,0]t \in [-t_0,0]. We can do this inductively in steps.

Step 1:
We can calculate C(t)C(t) for t∈[0,t0]t \in [0,t_0] by integration: C(t)=C(0)+∫0tC′(s)ds=p+∫0ta C(s−t0)ds=p+∫0ta pds=p (1+a t)C(t)=C(0)+\int_0^t C^\prime(s) ds = p+\int_0^t a\,C(s-t_0) ds = p + \int_0^t a\,p ds = p\,(1+a\,t)
Step n:
We can calulate C(t)C(t) for t∈[n t0,(n+1) t0]t \in [n\,t_0,(n+1)\,t_0] by integration: C(t)=C(n t0)+∫n t0tC′(s)ds=C(n t0)+∫n t0ta C(s−t0)dsC(t)=C(n\,t_0)+\int_{n\,t_0}^t C^\prime(s) ds = C(n\,t_0)+\int_{n\,t_0}^t a\,C(s-t_0) ds .
%md
Here is the method carried out for the first 3 steps.
︠e6bc279a-a8a8-4b9f-8a6a-11c25056e798︠
var('p a t t0 s')
t0=1;a=-.5;p=.5;
C0(t)=p
C1(t)=C0(0)+integral(a*C0(s-t0),(s,0,t));C1
C2(t)=C1(t0)+integral(a*C1(s-t0),(s,t0,t));C2
C3(t)=C2(2*t0)+integral(a*C2(s-t0),(s,2*t0,t));C3
C=piecewise([[(0,t0),C1],[(t0,2*t0),C2],[(2*t0,3*t0),C3]])
C.plot()
(p, a, t, t0, s) t |--> -0.25*t + 0.500000000000000 t |--> 0.0625*t^2 - 0.375*t + 0.5625 t |--> -0.010416666666666666*t^3 + 0.125*t^2 - 0.5*t + 0.6458333333333334

Problem:

In the original salt problem, a=−qK<0a= -\dfrac{q}{K} < 0. We take t0=1t_0=1 minute, a=−.5a=-.5 l/m and p=.5p=.5 g/l.

Note above that C(t)=0C(t)=0 at around t=2.7t=2.7. Discuss. Solve this problem: By adjusting the value of qq only, make C(1)=0C(1)=0.

Solution: If we increase qq the rate of inflow of fresh water, the concentration of salt will decrease faster. Want to solve the equation C(1)=C1(1)=0C(1)=C1(1)=0 for qq. From below, we see that q=Kq=K does the trick. That is, the rate of inflow of fresh water to KK liters/minute. 


︠06e439e9-a2b7-4199-aa80-ffa8ccd2972b︠
var('p a t t0 s K q')
t0=1;p=.5;a=-q/K
C0(t)=p
C1(t)=C0(0)+integral(a*C0(s-t0),(s,0,t));C1
solve(C1(1)-0,q)
(p, a, t, t0, s, K, q) t |--> -0.5*q*t/K + 0.500000000000000 [q == K] 
%md
<b>Solution:</b>
So set $a=-1$.
If we wanted $C(2)=0$, set $a=-2-

Solution: So set a=−1a=-1. If we wanted C(2)=0C(2)=0

C2(t)=C1(t0)+integral(a*C1(s-t0),(s,t0,t));C2
solve(C2(2),a)
t |--> (0.25*a*t^2 + (-0.5*a + 0.5)*t + 0.25*a - 0.5)*a + 0.5*a + 0.500000000000000 [a == -sqrt(2) - 2, a == sqrt(2) - 2] 

︠c5e6844e-4b43-411a-80f2-3bfc342614eci︠
html('Fixing this so that we can plot the concentration out to an arbitrary time and beyond.')
Fixing this so that we can plot the concentration out to an arbitrary time and beyond.
var('p a t t0 s')
t0=1;a=-.5;p=.5
G=Graphics()
C0(t)=p
for i in range(5):
    C1(t)=C0(i*t0)+integral(a*C0(s-t0),(s,i*t0,t));
    G+=plot(C1,(t,i*t0,(i+1)*t0))
    C0=C1
show(G)
(p, a, t, t0, s) 
%md
<h3>
Do dde's approach de's when the delay is small?
</h3>
Consider the dde $\left\{\begin{array}{lcl} y'(t) & = & y(t-h)\\ y(t) & = & 1 \mbox{ for }t \in [-h,0] \end{array}\right.$.  Let $y_h$ be the solution on $[0,1]$.  Is $\displaystyle\lim_{h \to 0^+} y_h = exp$?
I think so, but have no idea how to prove it.  Do you?
I can construct functions $y_{1/n}$ for various $n$ and compare the graph with the graph of $exp$.  See below.

Do dde's approach de's when the delay is small?

Consider the dde {y′(t)=y(t−h)y(t)=1 for t∈[−h,0]\left\{\begin{array}{lcl} y'(t) & = & y(t-h)\\ y(t) & = & 1 \mbox{ for }t \in [-h,0] \end{array}\right.. Let yhy_h be the solution on [0,1][0,1]. Is lim⁡h→0+yh=exp\displaystyle\lim_{h \to 0^+} y_h = exp? I think so, but have no idea how to prove it. Do you? I can construct functions y1/ny_{1/n} for various nn and compare the graph with the graph of expexp. See below.
def ddexp(n=10):
    a=1;p=1
    var('s')
    G=Graphics()
    t0=1/n
    C0(t)=p
    for i in range(n+1):
        C1(t)=C0(i*t0)+integral(a*C0(s-t0),(s,i*t0,t));
        G+=plot(C1,(t,i*t0,(i+1)*t0))
        C0=C1
    return (G,C1(1.))


pr1=ddexp(10)
pr2=ddexp(20)
pr3=ddexp(40)
pr4=(plot(exp,(0,1),color='red'),exp(1.))
pr1[1];pr2[1];pr3[1];pr4[1]
show(pr1[0]+pr2[0]+pr3[0]+pr4[0])
2.50100752809690 2.59749527837657 2.65432565421166 2.71828182845905 

︠4d1f19e0-54c1-4d92-b20d-e8726313e776︠
%md
It is pretty convincing to me.

It is pretty convincing to me.

Suppose now that the inflow is not necessarily fresh water, but rather a brine of concentration bb. Our dde now has an added term. C′(t)=b+a C(t−t0) C^{\prime}(t)=b + a\,C(t-t0) We can model this with the Salt procedure below.

var('p a b t t0 s')

def Salt(a=-1,p=1,b=0,t0=1,n=10):
    G=Graphics()
    C0(t)=p
    C=[]
    for i in range(n):
        C1(t)=C0(i*t0)+integral(b+a*C0(s-t0),(s,i*t0,t));
        C+=[[(i*t0,(i+1)*t0),C1]]
        C0=C1
    return piecewise(C)

#Make a call to the default value of Salt()
C=Salt()
(C.plot()+plot(h,(0,4))).show()
(p, a, b, t, t0, s) 
[lambert_w(i,1).n() for i in range(4)]
[0.567143290409784, -1.53391331979357 + 4.37518515306190*I, -2.40158510486800 + 10.7762995161151*I, -2.85358175540904 + 17.1135355394121*I] 
var('t');l=0.567143290409784
h=exp(l*t)
t 

︠dd84c45a-89e1-421a-afa0-d45fcf891be1︠
︠4f9cd073-ee4f-472e-bedf-f35c7b774117i︠
%html
<h4>
    Impossible behavior
</h4>
Note from above that when $t_0=p=1$ and $a=-1$, the graph of $C$ is a damped oscillation around 0, an unrealistic solution   ($C(t)<0$ is impossible) to the salt problem as it stands.  However, suppose the water coming into the tank is not pure but is itself brine of constant concentration $b$ g/l.  Then we might guess that the solution would  oscillate around  $y = b$.

<h5>
    <b>Problem 1:</b>  Let $t_0=1=p$ and $a=-1$.  Find a value for $b$ so that $C(t)>=0$ for $t>=0$
</h5>
<h5>
    <b>Problem 2:</b>  Given that $a<=0$. Find the values for $b>0$ in terms of $a$ so that $C(t)>=0$ for $t>=0$
</h5>

Impossible behavior

Note from above that when t0=p=1t_0=p=1 and a=−1a=-1, the graph of CC is a damped oscillation around 0, an unrealistic solution (C(t)<0C(t)<0 is impossible) to the salt problem as it stands. However, suppose the water coming into the tank is not pure but is itself brine of constant concentration bb g/l. Then we might guess that the solution would oscillate around y=by = b.
Problem 1: Let t0=1=pt_0=1=p and a=−1a=-1. Find a value for bb so that C(t)>=0C(t)>=0 for t>=0t>=0
Problem 2: Given that a[removed]0a[removed]0 in terms of aa so that C(t)>=0C(t)>=0 for t>=0t>=0 
b=var('b')
C=Salt(a=-.5,p=1,b=.019,t0=1,n=10)
C.plot()
for i in range(0,11):
    C(i)
1.0 0.519 0.15825 0.017958333333333062 0.0004218749999998772 0.017707812500001335 0.03250941840277921 0.038583057415668975 0.03950023580886336 0.03885444980735758 0.03824796327777724 
CP=derivative(C)
CP.plot()


︠741ec283-2b01-40d6-a45b-2274325a0df3︠
%md
<h4>
Delay equations with variable delay
</h4>
Let $g$ be any function on $[a,b]$ whose graph lies below the diagonal,  i.e., $g(t)<t$.  The step method extends to the DDE
$f^{\prime}(t)=f(g(t))$.  For example, suppose $g(t)=t^2$ on $[a,b]=[0,1]$.  Now we need to assume that $f(t)$ is given over some initial interval.  Say we assume $f(t)=p$ on $[0,h]$ for some constant $p$ and some small positive number $h$.  Then we can compute $f(t)$ for $t \in [h,\sqrt{h}]$,
$f(t)=f(h)+\displaystyle\int_{h}^t f'(\tau) d\tau = p + \displaystyle\int_{h}^t f(\tau^2) d\tau = p\,(1+t-h)$ .
</br>
For the next step, we have $f(t)=p\,(1+t+\frac{1}{3}t^3 + h^{1/2}k(h)$ on $[h^{1/2},h^{1/4}]$.  We can extend the definition of $f$ in this manner stepping over successive intervals up to 1.  (See below)

If we continue this we note that the polynomial in $t$ generated is the $n^{th}$ term of the power series
$f(t)=\displaystyle 1 + t +\frac{1}{3}t^3 +\frac{1}{3\cdot 7}t^7 +\frac{1}{3\cdot 7\cdot 15}t^{15} + \cdots + \frac{1}{\Pi_{i=1}^{n}2^{i}-1}t^{2^{n}-1} + \cdots $.
The radius of convergence of this series is calculated as $\sqrt{2}$.   The derivative of $f(t)$ is
$\displaystyle f'(t)=1 + t^2 + \frac{1}{3}t^4 + \frac{1}{3\cdot 7}t^6 \cdots = f(t^2)$!  Consequently, we have a solution to the variable dde  initial value problem $f'(t)=f(t^2)$, $f(0)=1$.

Delay equations with variable delay

Let gg be any function on [a,b][a,b] whose graph lies below the diagonal, i.e., g(t)<tg(t)<t. The step method extends to the DDE f′(t)=f(g(t))f^{\prime}(t)=f(g(t)). For example, suppose g(t)=t2g(t)=t^2 on [a,b]=[0,1][a,b]=[0,1]. Now we need to assume that f(t)f(t) is given over some initial interval. Say we assume f(t)=pf(t)=p on [0,h][0,h] for some constant pp and some small positive number hh. Then we can compute f(t)f(t) for t∈[h,h]t \in [h,\sqrt{h}], f(t)=f(h)+∫htf′(τ)dτ=p+∫htf(τ2)dτ=p (1+t−h)f(t)=f(h)+\displaystyle\int_{h}^t f'(\tau) d\tau = p + \displaystyle\int_{h}^t f(\tau^2) d\tau = p\,(1+t-h) .
For the next step, we have f(t)=p (1+t+13t3+h1/2k(h)f(t)=p\,(1+t+\frac{1}{3}t^3 + h^{1/2}k(h) on [h1/2,h1/4][h^{1/2},h^{1/4}]. We can extend the definition of ff in this manner stepping over successive intervals up to 1. (See below)

If we continue this we note that the polynomial in tt generated is the nthn^{th} term of the power series f(t)=1+t+13t3+13⋅7t7+13⋅7⋅15t15+⋯+1Πi=1n2i−1t2n−1+⋯f(t)=\displaystyle 1 + t +\frac{1}{3}t^3 +\frac{1}{3\cdot 7}t^7 +\frac{1}{3\cdot 7\cdot 15}t^{15} + \cdots + \frac{1}{\Pi_{i=1}^{n}2^{i}-1}t^{2^{n}-1} + \cdots . The radius of convergence of this series is calculated as 2\sqrt{2}. The derivative of f(t)f(t) is f′(t)=1+t2+13t4+13⋅7t6⋯=f(t2)\displaystyle f'(t)=1 + t^2 + \frac{1}{3}t^4 + \frac{1}{3\cdot 7}t^6 \cdots = f(t^2)! Consequently, we have a solution to the variable dde initial value problem f′(t)=f(t2)f'(t)=f(t^2), f(0)=1f(0)=1.


︠fc09bba7-7d45-4936-aafa-8adff3182452︠
var('n'); limit((2^(n+1)-1)^(1/(2*n)),n=infinity)
n e^(1/2*log(2)) 
def g(n):
    return
︠15c8352f-d5fa-4ef1-b8ad-5ecb284685f9︠
s,t,h,p=var('s t h p')
C0(t)=p
C=[]
ts=[h]
for i in range(8):
    ts+=[sqrt(ts[-1])]
for i in range(len(ts)-1):
    t0=ts[i]
    t1=ts[i+1]
    C1(t)=C0(t0)+integral(C0(s^2),(s,t0,t));
    C+=[[(t0,t1),C1]]
    C0=C1

C[3][1](t)(h=0)
C[3][1](0)
#C[6][1](t)(h=0)
#C[7][1](t)(h=0)
1/315*p*t^15 + 1/21*p*t^7 + 1/3*p*t^3 + p*t + p 1/315*(64*h - 15)*h^(7/8)*p - 2/21*(4*h - 7)*h^(3/4)*p - 1/3*(2*h - 3)*h^(5/8)*p + 2/3*h^(13/8)*p + 1/3*(2*h - 3)*sqrt(h)*p - (h - sqrt(h))*p + 1/21*h^(7/8)*p - 2/3*h^(3/4)*p - h^(5/8)*p + p 
C[7][1](1.1)(h=0)
4.45675695529338*p 
G=Graphics()
for c in C:
    G+=plot(c[1](p=1,h=.01),(0,1))
G.show()


def vdde_soln(t0=1/4,p=1):
    var('s')
    C0(t)=p
    C=[]
    ts=[t0]
    # construct the sequence t0,t0^.5,t0^.25, .. .999
    while 1-ts[-1] > .001:
        ts+=[sqrt(ts[-1])]
    for i in range(len(ts)-1):
        t0=ts[i]
        t1=ts[i+1]
        C1(t)=C0(t0)+integral(C0(s^2),(s,t0,t));
        C+=[[(t0,t1),C1]]
        C0=C1
    return piecewise(C)


C=vdde_soln()
C.plot()

%md
<h4>
Another Question
</h4>

Another Question 

C=vdde_soln(t0=1/64.)
C2=vdde_soln(t0=1/4.)

C.plot()+C2.plot(color='red')

C(.75)
0.09603041921289346 
var('t');parametric_plot((-t*sin(t),t*cos(t)),(t,-pi/2,pi/2))
t 
var('z')
solve(z+exp(-z)==0,z, to_poly_solve=True)
z [z == lambert_w(---1 
lambert_w(-1).n()
-0.318131505204764 + 1.33723570143069*I) 
lambert_w(-1).n()
-0.318131505204764 + 1.33723570143069*I] 
#help(lambert_w)

[lambert_w(i,1).n() for i in range(4)]
[0.567143290409784, -1.53391331979357 + 4.37518515306190*I, -2.40158510486800 + 10.7762995161151*I, -2.85358175540904 + 17.1135355394121*I] 
var('t');l=0.567143290409784
f=exp(l*t)
t 
f(t=0)
1 
exp(-l)-l
-1.11022302462516e-16 
df=diff(f,t)
g=-f(t=t-1)
g(t=3).n()
df(t=3).n()
-3.10895476357994 3.10895476357994 
[lambert_w(i,-1).n() for i in range(5)]
[-0.318131505204764 + 1.33723570143069*I, -2.06227772959828 + 7.58863117847251*I, -2.65319197403870 + 13.9492083345332*I, -3.02023970816450 + 20.2724576416152*I, -3.28776861154409 + 26.5804714993591*I] 





︠cd5eb1a5-802f-4310-8069-a3350bc1dd47i︠
%md
<p><br /><span style="font-size: large;"><strong>Do dde's approach de's when the delay is small?</strong></span><br /><br />Consider the dde $\left\{\begin{array}{lcl} y'(t) &amp; = &amp; y(t-h)\\ y(t) &amp; = &amp; 1 \mbox{ for }t \in [-h,0] \end{array}\right.$.&nbsp; Let $y_h$ be the solution on $[0,1]$.&nbsp; Is $\displaystyle\lim_{h \to 0^+} y_h = exp$?<br />I think so, but have no idea how to prove it.&nbsp; Do you?<br />I can construct functions $y_{1/n}$ for various $n$ and compare the graph with the graph of $exp$.&nbsp; See below.</p>

<p><span style="font-size: large;"><br /></span></p>