Solutions to Exercises

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

12.8. Solutions to Exercises#

12.8.1. Logistic Growth and COVID-19#

  1. First, separate the variables and antidifferentiate both sides

∫dyky(1−yM)=∫dt.

The partial fraction decomposition

Aky+B1−yM=1ky(1−yM)

with A=1 and B=1kM

leads to

1kln⁡MyM−y=t+C.

Considering the initial condition y(0)=y0>0, the solution can be obtained using algebra:

y=M1+(My0−1)e−kt

2a)

Hide code cell source 
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

plt.style.use('seaborn-poster')

%matplotlib inline

k=1
y0a=10
y0b=50
y0c=90
y0d=110
M=100

F = lambda t,y: k*y*(1-y/M)

t_eval = np.arange(0, 10, 0.01)
sola = solve_ivp(F, [0, 10], [y0a], t_eval=t_eval)
solb = solve_ivp(F, [0, 10], [y0b], t_eval=t_eval)
solc = solve_ivp(F, [0, 10], [y0c], t_eval=t_eval)
sold = solve_ivp(F, [0, 10], [y0d], t_eval=t_eval)
plt.figure(figsize = (5,5))
plt.plot(sola.t, sola.y[0])
plt.plot(solb.t, solb.y[0])
plt.plot(solc.t, solc.y[0])
plt.plot(sold.t, sold.y[0])


plt.title('Logistic Growth')
plt.xlabel('t (in days)')
plt.ylabel('number infected')
plt.ylim((0,120)) 
plt.xlim((0,11)) 
plt.legend(['y0=10','y0=50','y0=90','y0=110'])
plt.show()

<ipython-input-1-b136c51ef48b>:5: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-poster')
../../_images/4055b95be7e3137f8a6c1772b354cf3fb1cefa4f3415d8201af20ad3dcba0e5b.png

2b)

../../_images/logex.png

12.8.2. The Basic SIR Model#

Exercise 1

a) If the initial population is increased to 3,000,000, the epidemic appears to end sooner (around day 30)

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, R = y      # unpack current values of y
    beta, gamma      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [-beta * S * I,  
             beta * S * I - gamma * I,
             gamma * I]  
    return derivs

# Parameters
beta=.000001
gamma=.2
# Initial values
S0=3000000
I0=1
R0=0

# Bundle parameters for ODE solver
params = [beta,gamma]
# Bundle initial conditions for ODE solver
y0 = [S0,I0,R0]

# Make time array for solution
tStop = 70  # number of days
tInc = .001  
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(f, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(t, psoln[:,0],label='susceptible')
plt.plot(t, psoln[:,1],label='infected')
plt.plot(t, psoln[:,2],label='recovered')
plt.legend()
plt.xlabel('day')
plt.ylabel('number')
plt.savefig("SIR.png")
plt.tight_layout()
plt.show()
../../_images/85c02a95a706c7e97c1fecb850e9670cadc77756cb9ede948ef57ce952af45df.png

If the initial population is reduced to 1,000, there is no epidemic.

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, R = y      # unpack current values of y
    beta, gamma      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [-beta * S * I,  
             beta * S * I - gamma * I,
             gamma * I]  
    return derivs

# Parameters
beta=.000001
gamma=.2
# Initial values
S0=1000
I0=1
R0=0

# Bundle parameters for ODE solver
params = [beta,gamma]
# Bundle initial conditions for ODE solver
y0 = [S0,I0,R0]

# Make time array for solution
tStop = 1000  # number of days
tInc = .001  
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(f, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(t, psoln[:,0],label='susceptible')
plt.plot(t, psoln[:,1],label='infected')
plt.plot(t, psoln[:,2],label='recovered')
plt.legend()
plt.xlabel('day')
plt.ylabel('number')
plt.savefig("SIR.png")
plt.tight_layout()
plt.show()
../../_images/1285466aee23ac521f87f857c5bb4c6866aa230d10bc908e509399d2f509cf72.png

b) If the transmission rate is high (β=.0001), the entire population is quickly infected, and the epidemic ends sooner (around day 20).

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, R = y      # unpack current values of y
    beta, gamma      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [-beta * S * I,  
             beta * S * I - gamma * I,
             gamma * I]  
    return derivs

# Parameters
beta=.0001
gamma=.2
# Initial values
S0=1000000
I0=1
R0=0

# Bundle parameters for ODE solver
params = [beta,gamma]
# Bundle initial conditions for ODE solver
y0 = [S0,I0,R0]

# Make time array for solution
tStop = 70  # number of days
tInc = .001  
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(f, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(t, psoln[:,0],label='susceptible')
plt.plot(t, psoln[:,1],label='infected')
plt.plot(t, psoln[:,2],label='recovered')
plt.legend()
plt.xlabel('day')
plt.ylabel('number')
plt.savefig("SIR.png")
plt.tight_layout()
plt.show()
../../_images/14767bfddad094242a4ec384d70ab18918f39622e69240222bb0d657d76da209.png

If the infection rate is reduced (β=10−7), there is no epidemic.

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, R = y      # unpack current values of y
    beta, gamma      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [-beta * S * I,  
             beta * S * I - gamma * I,
             gamma * I]  
    return derivs

# Parameters
beta=.0000001
gamma=.2
# Initial values
S0=1000000
I0=1
R0=0

# Bundle parameters for ODE solver
params = [beta,gamma]
# Bundle initial conditions for ODE solver
y0 = [S0,I0,R0]

# Make time array for solution
tStop = 70  # number of days
tInc = .001  
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(f, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(t, psoln[:,0],label='susceptible')
plt.plot(t, psoln[:,1],label='infected')
plt.plot(t, psoln[:,2],label='recovered')
plt.legend()
plt.xlabel('day')
plt.ylabel('number')
plt.savefig("SIR.png")
plt.tight_layout()
plt.show()
../../_images/4222cd2c4ff4030ffa27da3f2941bd4e888a5bb207b5653a11e86a5afe104d9f.png

c) If the recovery rate is increased (γ=.3), fewer people get infected.

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, R = y      # unpack current values of y
    beta, gamma      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [-beta * S * I,  
             beta * S * I - gamma * I,
             gamma * I]  
    return derivs

# Parameters
beta=.000001
gamma=.3
# Initial values
S0=1000000
I0=1
R0=0

# Bundle parameters for ODE solver
params = [beta,gamma]
# Bundle initial conditions for ODE solver
y0 = [S0,I0,R0]

# Make time array for solution
tStop = 70  # number of days
tInc = .001  
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(f, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(t, psoln[:,0],label='susceptible')
plt.plot(t, psoln[:,1],label='infected')
plt.plot(t, psoln[:,2],label='recovered')
plt.legend()
plt.xlabel('day')
plt.ylabel('number')
plt.savefig("SIR.png")
plt.tight_layout()
plt.show()
../../_images/215a17fcc5fdc6b18a12d244b1995711e6f4c68b7904b0f9d9620fcff6cf449b.png

If the recovery rate is decreased (γ=.1), more people get infected, and the epidemic persists over a longer number of days.

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, R = y      # unpack current values of y
    beta, gamma      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [-beta * S * I,  
             beta * S * I - gamma * I,
             gamma * I]  
    return derivs

# Parameters
beta=.000001
gamma=.1
# Initial values
S0=1000000
I0=1
R0=0

# Bundle parameters for ODE solver
params = [beta,gamma]
# Bundle initial conditions for ODE solver
y0 = [S0,I0,R0]

# Make time array for solution
tStop = 70  # number of days
tInc = .001  
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(f, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(t, psoln[:,0],label='susceptible')
plt.plot(t, psoln[:,1],label='infected')
plt.plot(t, psoln[:,2],label='recovered')
plt.legend()
plt.xlabel('day')
plt.ylabel('number')
plt.savefig("SIR.png")
plt.tight_layout()
plt.show()
../../_images/57577113d15ec98698ea232143ef16fc60acf65f5d942c527c5dd337ff5c5184.png

d) As we saw in part b), with a high enough transmission rate, the entire population can become infected.

e) As we also saw in part b), with a low enough transmission rate, there is no outbreak.

Exercise 2

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def SIRV(y, t, params):
    S, I, R, V = y      # unpack current values of y
    beta, gamma, v     # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [-beta * S * I - v*V,  
             beta * S * I - gamma * I,
             gamma * I,
             v*S]  
    return derivs

# Parameters
beta=.000001
gamma=.2
v=.01
# Initial values
S0=3000000
I0=1
R0=0
V0=0

# Bundle parameters for ODE solver
params = [beta,gamma,v]
# Bundle initial conditions for ODE solver
y0 = [S0,I0,R0,V0]

# Make time array for solution
tStop = 70  # number of days
tInc = .001  
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(SIRV, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(t, psoln[:,0],label='susceptible')
plt.plot(t, psoln[:,1],label='infected')
plt.plot(t, psoln[:,2],label='recovered')
plt.plot(t, psoln[:,3],label='vaccinateded')
plt.legend()
plt.xlabel('day')
plt.ylabel('number')
plt.savefig("SIRV.png")
plt.tight_layout()
plt.show()
../../_images/ddaa314356c3aca23eaaf976dd4536712057e16434793d8589d23295e8178202.png
../../_images/SIRV.png

Exercise 3

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def SIRD(y, t, params):
    S, I, R, D= y      # unpack current values of y
    beta, gamma, mu     # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [-beta * S * I,  
             beta * S * I - gamma * I - mu*I,
             gamma * I,
             mu*I]  
    return derivs

# Parameters
beta=.000001
gamma=.2
mu=.01
# Initial values
S0=3000000
I0=1
R0=0
D0=0

# Bundle parameters for ODE solver
params = [beta,gamma,mu]
# Bundle initial conditions for ODE solver
y0 = [S0,I0,R0,D0]

# Make time array for solution
tStop = 70  # number of days
tInc = .001  
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(SIRD, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(t, psoln[:,0],label='susceptible')
plt.plot(t, psoln[:,1],label='infected')
plt.plot(t, psoln[:,2],label='recovered')
plt.plot(t, psoln[:,3],label='deceased')
plt.legend()
plt.xlabel('day')
plt.ylabel('number')
plt.savefig("SIRD.png")
plt.tight_layout()
plt.show()
../../_images/a6e5af806aac6129f1ef66bcc53d5c44bc2b0d3b47fc207c1969ece2a2f396ad.png
../../_images/SIRD.png

12.8.3. Cholera in Haiti#

Exercise 1

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, B = y      # unpack current values of y
    H,n,a,K,r,Nb,e      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [n*(H-S)-a*B*S/(K+B),  
             a*B*S/(K+B)-r*I,
             B*Nb+e*I]  
    return derivs

# Parameters
H=10000
n=.0001
K=10**6
r=.2
Nb=-.33
e=10
# Initial values
S0=10000
I0=100
B0=0

exposure_rate=[]
max_infected=[]
    
for a in np.arange(0,100,1):
    exposure_rate.append(a)
    # Bundle parameters for ODE solver
    params = [H,n,a,K,r,Nb,e]
    # Bundle initial conditions for ODE solver
    y0 = [S0,I0,B0]

    # Make time array for solution
    tStop = 10  # number of days
    tInc = .001  
    t = np.arange(0., tStop, tInc)

    # Call the ODE solver
    psoln = odeint(f, y0, t, args=(params,))
    max_infected.append(np.max(psoln[:,1]))
    
# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(exposure_rate, max_infected,label='max infected')
plt.legend()
plt.xlabel('a=exposure rate')
plt.ylabel('max infected population')
plt.savefig("cholerEx1.png")
plt.tight_layout()
plt.show()
../../_images/e4836c961d6d99829f7cc4e2de9b716fc2ef2d35d9ae2951c22f58e2ee52f568.png

The max infected population rises sharply as the exposure rate (a) is increased from 0 to 10. When a=100, the max_infected population has leveled off at approximately 8000 people.

Exercise 2

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, B = y      # unpack current values of y
    H,n,a,K,r,Nb,e      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [n*(H-S)-a*B*S/(K+B),  
             a*B*S/(K+B)-r*I,
             B*Nb+e*I]  
    return derivs

# Parameters
H=10000
n=.0001
K=10**6
r=.2
Nb=-.33
# Initial values
S0=10000
I0=100
B0=0

contam=[]
peak1=[]
peak2=[]
    
for e in np.arange(0,50,1):
    contam.append(e)
    a=40
    # Bundle parameters for ODE solver
    params = [H,n,a,K,r,Nb,e]
    # Bundle initial conditions for ODE solver
    y0 = [S0,I0,B0]

    # Make time array for solution
    tStop = 10  # number of days
    tInc = .001  
    t = np.arange(0., tStop, tInc)

    # Call the ODE solver
    psoln = odeint(f, y0, t, args=(params,))
    peak1.append(np.max(psoln[:,1]))
    
    a=20
        # Bundle parameters for ODE solver
    params = [H,n,a,K,r,Nb,e]
    # Bundle initial conditions for ODE solver
    y0 = [S0,I0,B0]

    # Make time array for solution
    tStop = 10  # number of days
    tInc = .001  
    t = np.arange(0., tStop, tInc)

    # Call the ODE solver
    psoln = odeint(f, y0, t, args=(params,))
    peak2.append(np.max(psoln[:,1]))
    
# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(contam, peak1,label='max infected (a=40)')
plt.plot(contam, peak2,label='max infected (a=20)')
plt.legend(["a=40","a=20"])
plt.xlabel('e=water contamination rate')
plt.ylabel('max infected population')
plt.savefig("cholerEx2.png")
plt.tight_layout()
plt.show()
../../_images/f79cba871952245ccaca1efd30d6f5d104c06051d68d4793ac8427b1e29c6a57.png

Increased exposure rate (a=20 to a=40) increases the max infected population at each level of the water contamination rate (e).

Exercise 3

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, B = y      # unpack current values of y
    H,n,a,K,r,Nb,e      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [n*(H-S)-a*B*S/(K+B),  
             a*B*S/(K+B)-r*I,
             B*Nb+e*I]  
    return derivs

# Parameters
H=10000
n=.0001
K=10**6
e=10
Nb=-.33
# Initial values
S0=10000
I0=100
B0=0

recovery=[]
peak1=[]
peak2=[]
    
for r in np.arange(0,.8,.01):
    recovery.append(r)
    a=40
    # Bundle parameters for ODE solver
    params = [H,n,a,K,r,Nb,e]
    # Bundle initial conditions for ODE solver
    y0 = [S0,I0,B0]

    # Make time array for solution
    tStop = 10  # number of days
    tInc = .001  
    t = np.arange(0., tStop, tInc)

    # Call the ODE solver
    psoln = odeint(f, y0, t, args=(params,))
    peak1.append(np.max(psoln[:,1]))
    
    a=80
        # Bundle parameters for ODE solver
    params = [H,n,a,K,r,Nb,e]
    # Bundle initial conditions for ODE solver
    y0 = [S0,I0,B0]

    # Make time array for solution
    tStop = 10  # number of days
    tInc = .001  
    t = np.arange(0., tStop, tInc)

    # Call the ODE solver
    psoln = odeint(f, y0, t, args=(params,))
    peak2.append(np.max(psoln[:,1]))
    
# Plot results
fig = plt.figure(1, figsize=(5,3))
plt.plot(recovery, peak1,label='max infected (a=40)')
plt.plot(recovery, peak2,label='max infected (a=80)')
plt.xlabel('r=recovery rate')
plt.ylabel('max infected population')
plt.savefig("cholerEx3.png")
plt.tight_layout()
plt.show()
../../_images/7663f5009099b640e0b386299e2d369f687551f31ae1d3ce4af37c953fa79054.png

Peak infection dependency on the cholera recovery rate r for two exposure levels (a=40/day and a=80/day).

Exercise 4

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Equations
def f(y, t, params):
    S, I, B = y      # unpack current values of y
    H,n,a,K,r,Nb,e      # unpack parameters
    #specify dS/dt, dI/dt and dR/dt
    derivs = [n*(H-S)-a*B*S/(K+B),  
             a*B*S/(K+B)-r*I,
             B*Nb+e*I]  
    return derivs

# Parameters
H=10000
n=.0001
K=10**6
r=.2
e=10
Nb=-.33
# Initial values
S0=10000
I0=0


Z=np.eye(101, 201)  
    
for B0 in np.arange(0,100,1):     
    for i in np.arange(0,200,1):  
        a=i/10  # exposure rate
        # Bundle parameters for ODE solver
        params = [H,n,a,K,r,Nb,e]
        # Bundle initial conditions for ODE solver
        y0 = [S0,I0,B0]

        # Make time array for solution
        tStop = 10  # number of days
        tInc = .001  
        t = np.arange(0., tStop, tInc)

        # Call the ODE solver
        psoln = odeint(f, y0, t, args=(params,))
        # Record the peak infection
        Z[B0,i]=np.max(psoln[:,1])
    

# Plot results
# Creating dataset
x=np.arange(0,100,1)
y=np.arange(0,200,1)


# Creating figure
fig = plt.figure()
 
# syntax for 3-D projection
ax = plt.axes(projection ='3d')
 
for i in np.arange(1,100,1):
    for j in np.arange(0,200,1):
        ax.scatter(i, j, Z[i,j],color='lightgray')
ax.set_title('Max Infected')
ax.set_xlabel('a= .1 * i = exposure rate')
ax.set_ylabel('B0=initial bacteria concent')
ax.set_zlabel('Max Infected')
plt.savefig("cholerEx4.png")
plt.show()
../../_images/0e8ee649a8e560d2c574c93b465e594fb0fe02e724153ce5ae49e9721dffa634.png

Exercise 5

a) Though it would be better to revise the model by adding a vaccinated compartment V (t), the basic model might reflect a vaccination program by lowering the initial susceptible population S(0) and raising the value of K, the latter being the bacteria concentration required to give a .5 probability.

b) Seasonal behavior of cholera outbreaks which may be related to weather events such as floods or El Nino must be considered. That is, even though there may be a period with no or very few cases, there is a latent possibility for re-emergence of a large number of infected population. Mathematically viable solutions may be overly simplistic in light of the reality of everyday life in Haiti.

From the sensitivity analysis, the contamination rate seems to have the largest effect on the infected population, although exposure rate, recovery rate, and initial concentration of bacteria also have an effect. Within the range of 2.5-30 number of people contaminated per day, the infected population will increase significantly. As the contamination rate increases above 30, the infected population stabilizes and increases more slowly. These observations suggest focusing response and preparation measures on lowering the contamination rate

12.8.4. CWS Model of Alzheimer’s Disease#

Exercise 1

(i=3)0=Eb1b2+Fb4−Eb1b3−Fb3
(i=4)0=Eb1b3+Fb5−Eb1b4−Fb4

Exercise 2

The total exit rate of monomers from the CSF (h2=r21+r23+L2) is the sum of the respective rates monomers move from the CSF to the brain (r21 and from CSF to plasma (r23), these being added to the loss of monomers in the brain by degradation (L2).

Exercises 3

Treatment 2 shows that the plasma and CSF levels of amyloid-beta monomers might increase, but the level in the brain actually decreases.

Exercise 4

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def f(y, t, params):
    b1, c,p = y      # unpack current values of y
    k1,k2,k3, r12,r13,r23,r21,r31,r32,L1,L2,L3,E,F,S = params  # unpack parameters
    derivs = [k1+F*S+r31*p+r21*c-(r13*b1+r12*b1*(1+np.heaviside(t-100,0)-np.heaviside(t-465,0)))-(L1+E*S)*b1-2*E*b1**2,      # list of dy/dt=f functions
             k2+r12*b1*(1+np.heaviside(t-100,0)-np.heaviside(t-465,0))+r32*p-(r21*c+r23*c)-L2*c,
             k3+r13*b1+r23*c-(r31*p+r32*p)-L3*p]
    return derivs

# Parameters
k1 = 7.34*86400         
L1 = 24 
k2=0
L2=0
k3=0
L3=6.7*10**(-3)*86400   
r12=7.6*10**(-6)*86400   
r23=1.7*10**(-3)*86400   
r31=3.7*10**(-5)*86400   
r32=0
r21=0
r13=0
E=0
F=.238
S=0

# Initial values
b10 = 25.7*(10**3)     
c0 = 115
p0 = 29

# Bundle parameters for ODE solver
params = [k1,k2,k3, r12,r13,r23,r21,r31,r32,L1,L2,L3,E,F,S]

# Bundle initial conditions for ODE solver
y0 = [b10,c0,p0]

# Make time array for solution
tStop = 1000
tInc = 0.05
t = np.arange(0., tStop, tInc)

# Call the ODE solver
psoln = odeint(f, y0, t, args=(params,))

# Plot results
fig = plt.figure(1, figsize=(8,8))

# Plot b1 as a function of time
ax1 = fig.add_subplot(311)
ax1.plot(t, psoln[:,0])
ax1.set_xlabel('time')
ax1.set_ylabel('b1')

# Plot c as a function of time
ax2 = fig.add_subplot(312)
ax2.plot(t, psoln[:,1])
ax2.set_xlabel('time')
ax2.set_ylabel('c')

# Plot p as a function of time
ax3 = fig.add_subplot(313)
ax3.plot(t, psoln[:,2])
ax3.set_xlabel('time')
ax3.set_ylabel('p')
plt.savefig("treatment2nonlin.png")
plt.tight_layout()
plt.show()
../../_images/895ba82d400af11e42ba4a8f31b7abf83b03cf1efd66a84c533322cac74192ad.png

The output of the nonlinear is the same as the simplified linear model since the former is equivalent to the latter when E=S=0.

12.8.5. Gravity-Fed Water Delivery#

[2.1] Fgn=Fgcos⁡(θ)=γδVcos⁡(θ)=γδV1−(dzds)2.

[2.2] dpds≈p(s0−δs2)−p(s0)−δs2⇒p(s0−δs2)≈p(s0)−dpdsδs2.

[2.3] Since z has dimension L, we must show the other two terms do as well. Pressure being force per unit area, and force being mass time acceleration, p has dimension MLT2L2=MLT2. Since γ=ρg has dimension ML3LT2=ML2T2, the dimension of pγ is MLT2L2T2M=L. The dimension for u2g is (LT)2LT2=L.

[3.1.1]

[a] Bernoulli’s equation states pγ + u22g + z = c. Since u=0 msec as a result of the tap being closed, zB=0 m as shown in {\bf Figure \ref{fig9}}, and so our equation simplifies to pγ = c. Observing that c = 9 m + 2 m = 11 ft = pkNm29.789kNm3, we obtain p=105.5 kNm2.

[b] From part a we know pγ = c. Observing that c = 30 ft + 6 ft = 36 ft = plbft262.3lbft3, we obtain p=2242.8 lbft2=2242.8lbft2⋅(1ft)2(12in)2=15.6lbin2.

[3.1.2] Referring to Figure (S1), let h be the desired vertical distance from the source reservoir to the faucet. Then h satisfies

1862kNm29.789kNm3=h+3m⇒h≈187m.
../../_images/figS1.png

[3.2.1] [a] At a horizontal distance of 600 m, pressure is greatest since the vertical drop from the source is a maximum (114 m).

[b] The pipe can operate at a pressure head of pγ=917kNm29.789kNm3 = 93.7 m, so 1 intermediate tank is required.

[3.2.2] Pressure at A is equal to pA=[(h1+h2)−(h3+h4)]γ. The pressure at B is equal to pB=(h3+h4)γ. The break-pressure tank effectively reduces the pressure at the faucet by the amount pA.

[3.3.1] .1215 liters/sec.

[3.3.2] [a] A velocity head of 3 m (10 ft) corresponds to a speed u1=[3(2)(9.8)]1/2≈7.7 m/s. By the continuity equation, π(0.025)2u1=π(0.02)2u2 so u2≈11.9 m/s with a corresponding velocity head of roughly 7.3 m (31.6 feet).

[b]

πr12u0=πr22uB⇒u0=(r2r1)22g(h1+h2).

[4.1]

p1γ+u122g+z1=p2γ+u222g+z2+hf⇒0+0+50=p2γ+2.522g+30+25
p2=(50−0.32−25−30)∗9.789kNm2=−52.08kNm2

A negative pressure means that pressure inside the pipe is less than outside. This can cause dirt or air to be sucked into the pipe, resulting in a stoppage of flow.

[4.2] The Reynolds number ρu¯Dμ is dimensionless since all units cancel in the expression ML3LTLMLT.

[4.3.1] Fgs=−Fgsin⁡θ=−mgsin⁡θ=−ρ(πr2L)gsin⁡θ=−γπr2Lsin⁡θ.

[4.3.2] Q=∫02π∫0Ru0(1−(r/R)2)rdrdθ=πR2u02.

[4.3.3] Since the pipe is wider at B than at A, uB is less than uA. Head loss from A to B is represented by hAB.

[4.3.4] $Re=ρu¯Dμ=998.2⋅2⋅D0.001<2100⇒D<0.001m=1mm.$

The flow rate would be

Q=u¯A<u¯π(D2)2=1.6x10−6m3s.

So for a very reasonable u¯ the pipe would have to be unreasonably small. This shows that, while laminar flow is a useful assumption to analyze friction, it is unrealistic for gravity-fed water delivery systems.

[4.4.1] The iron pipe has the greater friction factor.

[4.4.2]

../../_images/figS2.png

[4.5.1] The average velocities associated with the minimum and maximum flow rates are respectively u−=.35(10−3)π(.0125)2=.71(m/s) and u−=1.4(10−3)π(.0125)2=2.85 (m/s). The respective Reynolds numbers are Re=1.23(.71)(.025)1.79(10−5)≈1200, and Re=1.23(2.85)(.025)1.79(10−5)≈4900. This indicates that the flow at the minimum flow rate is laminar, but will be turbulent at the maximum flow rate. Note in {\bf Table 5} the big difference in head loss between these laminar and turbulent flows: 3.45 m per 100 m at the minimum flow rate versus 49.5 m per 100 m at the maximum flow rate.

[4.5.2] Using the equation for f given in {\bf (\ref{Colebrook})} or a Moody chart , we can determine the Reynolds number Re. We may then compute

u−=μReρD,
Q=π(D2)2u−

and

hf=fLDu−22g.

Knowing the head loss, we can easily plot the HGL on the graph shown in Figure (21). The HGL will help us decide whether the pipe may be used for that segment of the system as explained in Section 4.4.

[5.1] There are different possible designs, but a system using a 80 mm diameter is shown here. In 10 years Tesuque Pueblo will have about 1148 people and a required volumetric flow rate of Q=4.27l/s=0.00427m3/s. A 80 mm diameter pipe has a cross-sectional area of 0.005 m2, so our average velocity is 0.85 m/s (low, but within limits). The maximum static pressure will be 2255−1950=335 m at the outlet faucet, which is above the maximum operating pressure head from {\bf Table 5} (146 m). Break pressure tanks at elevations of 2130 m and 2030 m reduce the maximum static pressures to 125 m, 100 m, and 110 m in the first, second, and third sections, respectively (Note that since a 2030 m elevation isn’t included in the survey, we don’t know if there is a suitable building site at this location. In that case a third break pressure tank would be needed). By rounding our flow rate down we get from {\bf Table 5} a head loss of about 0.96 m per 100 m. We will approximate horizontal distance as pipe length, so our HGL is modeled by the equation 2255−0.8522g−0.96x. If we compare this line to the given elevations we find the elevations remain below the HGL line. This means the dynamic pressure will never be negative and the system is feasible.

12.8.6. Earthquake Resistant Construction#

2.1

a) u=c1cos⁡t+c2sin⁡t;

b) u(t)=e−12t[c1cos⁡(32t)+c2sin⁡(32t)];

c)u(t)=e−t(c1+c2t);

d) u(t)=c1e(−2+3)t+c2e(−2−3)t.

2.2 For example, to check the solution for c=4, the initial conditions u(0)=1, u.(0)=0 imply that the constants c1 and c2 in the general solution u(t)=c1e(−2+3)t+c2e(−2−3)t must satisfy

c1+c2=1
c1(−2+3)+c2(−2−3)=0.

It follows easily that c1=2+323 and c2=−2+323.

2.3 For 0<c<4, u=hh+up where uh satisfies

uh..+cuh.+4uh=0,

and up=Acos⁡(2t)+Bsin⁡(2t) is a particular solution to the non-homogeneous equation. Solving, we have u=12csin⁡(2t)−2c16−c2e−ct2sin⁡(16−c22t).

../../_images/figE1.png

At the critical value c=4, the solution becomes u=−14te−2t+18sin⁡(2t)

../../_images/figE2.png

Note that in general, increasing the damping coefficient c will reduce the amplitude of the periodic term 12csin⁡(2t) in the solution and cause a more rapid decay of the transient term 2c16−c2e−ct2sin⁡(16−c22t) in the solution.

2.4 Since Rmax=2=kumaxFmax with k=1010kgsec and umax=5x10−3, we obtain the maximum allowable blast force Fmax=2.5(107)N.

2.4.1

u(t)=u2(t)=1mω0∫0tF(τ)sin⁡(ω0(t−τ))dτ
=1mω0∫0TF(τ)sin⁡(ω0(t−τ))dτ+1mω0∫TtF(τ)sin⁡(ω0(t−τ))dτ
=Fmaxmω0∫0T(−τT+1)sin⁡(ω0(t−τ))dτ
=Fmaxmω0[(−τT+1)cos⁡(ω0(t−τ))ω0∣0T−∫0Tcos⁡(ω0(t−τ))ω0(−1T)dτ]
=Fmaxmω02[(−τT+1)cos⁡(ω0(t−τ))∣0T−(sin⁡(ω0(t−τ)ω0T)∣0T]
Fmaxk[(−TT+1)(cos⁡(ω0(t−T))−(−0T+1)(cos⁡(ω0(t))
=−1Tsin⁡(ω0(t−T))+1Tsin⁡(ω0(t))
=Fmaxk[0−cos⁡(ω0t)−1Tsin⁡(ω0(t−T))+1Tsin⁡(ω0(t))]
=Fmaxk[−cos⁡(ωt)]

2.4.2 a)

s2L[u]+ω02L[u]=Fmaxm{L[−tT+1]−L[(−tT+1)UT(t)]}
(s2+ω02)L[u]=Fmaxm{−1T(1s2)+1s−eTs(−1Ts2)}
L[u]=Fmaxm{−1(Ts2)(s2+ω02)+1s(s2+ω02)+e−TsTs2(s2+ω02)}

b) $L[u]=Fmaxm{−1T+ss2(s2+ω02)+1T+e−Tss2(s2+ω02)}$

L[u]=FmaxmTω02{Ts−1s2−Tss2+ω02+1s2+ω02+e−sT(1s2−1s2+ω02)}
L[u]=FmaxmTω02{L[T−t−Tcos⁡(ω0t)+1ω0sin⁡(ω0t)]+e−sTL[t−1ω0sin⁡(ω0t)]}
L[u]=FmaxmTω02{L[T−t−Tcos⁡(ω0t)+1ω0sin⁡(ω0t)+UT(t)[t−T−1ω0sin⁡(ω0(t−T))]]}

c)

u=FmaxmTω02{T−t−Tcos⁡(ω0t)+1ω0sin⁡(ω0t)+UT(t)[t−T−1ω0sin⁡(ω0t)]}
u=Fmaxk{1−tT−cos⁡(ω0t)+1ω0Tsin⁡(ω0t)+UT(t)[tT−1−1ω0Tsin⁡(ω0t)]}

2.5 The ground motion z(t) effectively creates a reversed inertial force −mz...

3.1 W=(w1w2w3) represents the relative displacements, M=(m1000m2000m3) represents the mass of the floors, K=(2(k1+k2)−2k20−2k22(k2+k3)−2k30−2k32k3) represents the spring constants, and F=(F1−m1z..F2−m2z..F3−m3z..) the represents wind force and the inertial force caused by the ground motion.

3.3 x1(t)=6+217517+17+c11cos⁡5+174t+c12sin⁡5+174t and x2(t)=−6+217517−17+c21cos⁡5−174t+c22sin⁡5−174t. The initial conditions w1(0)=w2(0)=w1.(0)=w2.(0)=0 imply that x1(0)=x2(0)=x1.=(0)=x2.(0)=0, giving c11=−6−217517+17,c12=0,c22=6−217517−17,c21=0. Finally, using W=ΦX gives w1(t)=x1(t)+x2(t), and w2(t)=(3−174)x1(t)+(3+174)x2(t).

3.4 U=[u(t)θ(t)], M=[M+mMaMaIG+Ma2] , C=[c000], K=[k00K−Mga] and F(t)=[cz.+kz0]

Contents