%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math
def f(k, r, t):
return k * math.e ** (r * t)
k = 4
r = 0.5
x = np.linspace(0, 20, 4000)
plt.plot(x, f(k, r, x))
[<matplotlib.lines.Line2D at 0x10db198d0>]
p = 4
sim_length = 20
delta_t = 1
xs = []
ys = []
growth_rate_per_step = r * delta_t
num_iter = int(sim_length / delta_t) + 1
for i in range(num_iter):
xs.append(i * delta_t)
ys.append(p)
p += growth_rate_per_step * p
plt.plot(x, f(k, r, x))
plt.plot(xs, ys)
[<matplotlib.lines.Line2D at 0x10d693b90>]
as we shrink the delta_t, we get closer and closer to the actual solution. local error (error per step) is proportional to the square of the step size
time = []
pop = []
simLength = 120
population = 4
M = 500.0
growthRate = 0.1
deltaT = 0.025
growthRatePerStep = growthRate * deltaT
numIterations = int(simLength / deltaT) + 1
for i in range(numIterations):
population += growthRatePerStep * population * (1 - population / M)
# population *= 1 + growthRatePerStep
time.append(i * deltaT)
pop.append(population)
plt.plot(time, pop)
[<matplotlib.lines.Line2D at 0x10db79710>]
time = []
temperature = []
simLength = 70
temp = 98.6
growthRate = 0.023079
envTemp = 72
deltaT = 0.0025
growthRatePerStep = growthRate * deltaT
numIterations = int(simLength / deltaT) + 1
for i in range(numIterations):
temp += -growthRatePerStep * (temp - envTemp)
# population *= 1 + growthRatePerStep
time.append(i * deltaT)
temperature.append(temp)
plt.plot(time, temperature)
plt.plot(time, [82] * len(time))
temperature[-1]
77.28713710191359
r1 = 1
p1 = 6
b1 = 0.27
r2 = 1
p2 = 5
b2 = 0.20
sim_length = 50
dt = 0.001
xs = []
y1s = []
y2s = []
num_iter = int(sim_length / dt) + 1
for i in range(num_iter):
xs.append(i * dt)
y1s.append(p1)
y2s.append(p2)
p1old = p1
p2old = p2
p1 += (r1 * p1old - b1 * p1old * p2old) * dt
p2 += (-r2 * p2old + b2 * p1old * p2old) * dt
plt.plot(xs, y1s)
plt.plot(xs, y2s)
[<matplotlib.lines.Line2D at 0x10dcfb690>]
plt.plot(y1s, y2s)
plt.ylabel("Foxes")
plt.xlabel("Rabbits")
<matplotlib.text.Text at 0x10c472190>
a = 1
b = 3
c = 1
d = 5
s = 4
xr = -8.0/5.0
r = 0.001
I = 2.3
p1 = 1
p2 = 1
p3 = 2
sim_length = 500
dt = 0.001
xs = []
y1s = []
y2s = []
y3s = []
num_iter = int(sim_length / dt) + 1
for i in range(num_iter):
xs.append(i * dt)
y1s.append(p1)
y2s.append(p2)
y3s.append(p3)
p1old = p1
p2old = p2
p3old = p3
p1 += (p2old - a * p1old ** 3 + b * p1old ** 2 - p3old + I) * dt
p2 += (c - d * p1old ** 2 - p2old) * dt
p3 += (r * (s * (p1old - xr) - p3old)) * dt
plt.plot(xs, y1s)
plt.plot(xs, y2s)
plt.plot(xs, y3s)
[<matplotlib.lines.Line2D at 0x10c5cecd0>]
p = 0
r = -0.08
dosage = 3
dosage_freq = 6
sim_length = 200
delta_t = 0.01
xs = []
ys = []
num_iter = int(sim_length / delta_t) + 1
for i in range(num_iter):
xs.append(i * delta_t)
ys.append(p)
p += r * delta_t * p
if i % int(dosage_freq / delta_t) == 0:
p += dosage
plt.plot(xs, ys)
[<matplotlib.lines.Line2D at 0x10da3d7d0>]
N = 1000
I = 1
S = N - I
R = 0
beta = 0.0002
gamma = 0.05
sim_length = 120
dt = 0.01
xs = []
y1s = []
y2s = []
y3s = []
num_iter = int(sim_length / dt) + 1
for i in range(num_iter):
xs.append(i * dt)
y1s.append(S)
y2s.append(I)
y3s.append(R)
Sold = S
Iold = I
Rold = R
S += (-beta * Iold * Sold) * dt
I += (beta * Iold * Sold - gamma * Iold) * dt
R += (gamma * Iold) * dt
plt.plot(xs, y1s)
plt.plot(xs, y2s)
plt.plot(xs, y3s)
plt.xlabel("Time")
plt.ylabel("Population")
plt.title("SIR Epidemiology Model")
plt.legend(['Susceptible', 'Infected', 'Recovered'], loc='center right')
print "r0", N * beta / gamma
print "recovery time", 1 / gamma
print "infections per person", beta * N
r0 4.0 recovery time 20.0 infections per person 0.2
N = 1000
I = 1
S = N - I
R = 0
beta = 0.0002
gamma = 0.05
mu = 0.01
sim_length = 120
dt = 0.01
xs = []
y1s = []
y2s = []
y3s = []
num_iter = int(sim_length / dt) + 1
for i in range(num_iter):
xs.append(i * dt)
y1s.append(S)
y2s.append(I)
y3s.append(R)
Sold = S
Iold = I
Rold = R
S += (-beta * Iold * Sold + mu * N - mu * Sold) * dt
I += (beta * Iold * Sold - gamma * Iold - mu * Iold) * dt
R += (gamma * Iold - mu * Rold) * dt
plt.plot(xs, y1s)
plt.plot(xs, y2s)
plt.plot(xs, y3s)
plt.xlabel("Time")
plt.ylabel("Population")
plt.title("SIR Epidemiology Model")
plt.legend(['Susceptible', 'Infected', 'Recovered'])
print "r0", N * beta / gamma
print "recovery time", 1 / gamma
print "time between infections", 1 / beta
r0 4.0 recovery time 20.0 time between infections 5000.0
p = 4
r = 0.1
sim_length = 20
delta_t = 0.1
xs = []
ys = []
growth_rate_per_step = r * delta_t
num_iter = int(sim_length / delta_t) + 1
for i in range(num_iter):
xs.append(i * delta_t)
ys.append(p)
p += growth_rate_per_step * p
plt.plot(xs, ys)
[<matplotlib.lines.Line2D at 0x1144f2b50>]
p = 4
r = 0.1
xs = np.linspace(0, 20, 200)
ys = [p]
for i in range(1, len(xs)):
dt = xs[i] - xs[i - 1]
p += r * p * dt
ys.append(p)
plt.plot(xs, ys)
[<matplotlib.lines.Line2D at 0x113be78d0>]
def growth(u, dt):
return u * 0.1 * dt
p = 4
xs = np.linspace(0, 20, 200)
ys = [p]
for i in range(1, len(xs)):
dt = xs[i] - xs[i - 1]
p += growth(p, dt)
ys.append(p)
plt.plot(xs, ys)
[<matplotlib.lines.Line2D at 0x113e0bb50>]
def euler(func, initial, deltas):
ys = [initial]
for i in range(1, len(xs)):
dt = xs[i] - xs[i - 1]
initial += func(initial, xs[i]) * dt
ys.append(initial)
return ys
def growth(u, t):
return u * 0.1
p = 4
xs = np.linspace(0, 20, 200)
ys = euler(growth, p, xs)
plt.plot(xs, ys)
[<matplotlib.lines.Line2D at 0x113e6bfd0>]
def euler(func, initial, deltas):
ys = [initial]
for i in range(1, len(xs)):
cur = np.array(ys[-1][:])
dt = xs[i] - xs[i - 1]
cur = cur + np.array(func(cur, xs[i])) * dt
ys.append(cur)
return np.array(ys)
beta = 0.0001
gamma = 0.02
def sir(u, t):
s, i, r = u
sd = -beta * s * i
id = beta * s * i - gamma * i
rd = gamma * i
return [sd, id, rd]
init = [999, 1, 0]
xs = np.linspace(0, 200, 2000)
ys = euler(sir, init, xs)
plt.plot(xs, ys[:, 0])
plt.plot(xs, ys[:, 1])
plt.plot(xs, ys[:, 2])
[<matplotlib.lines.Line2D at 0x114094310>]
def growth(u, t):
return u[0] * 0.1
init = [999, 1, 0]
xs = np.linspace(0, 20, 200)
ys = euler(growth, [4], xs)
plt.plot(xs, ys[:, 0])
[<matplotlib.lines.Line2D at 0x1140f5750>]
def lorenz(u, t):
x, y, z = u
xd = (10 * (y - x))
yd = (x * (28 - z) - y)
zd = (x * y - (8.0/3.0) * z)
return [xd, yd, zd]
init = [999, 1, 0]
xs = np.linspace(0, 100, 200000)
ys = euler(lorenz, [6, 6, 6], xs)
plt.plot(xs, ys[:, 0])
plt.plot(xs, ys[:, 1])
plt.plot(xs, ys[:, 2])
[<matplotlib.lines.Line2D at 0x11423f310>]
plt.plot(ys[:, 0], ys[:, 2])
[<matplotlib.lines.Line2D at 0x114d879d0>]
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(ys[:, 0], ys[:, 1], ys[:, 2])
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x114d9f850>]