%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)) 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) 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) 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] 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) plt.plot(y1s, y2s) plt.ylabel("Foxes") plt.xlabel("Rabbits") 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) 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) 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 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 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) 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) 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) 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) 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]) 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]) 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]) plt.plot(ys[:, 0], ys[:, 2]) from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.plot(ys[:, 0], ys[:, 1], ys[:, 2])