Implementation of the Lotke-Volterra system of equations in Python for SAMBa Software Skills Training.
import numpy as np

def evolve(phi, f, Df, t, z0, T, N):
    z = np.zeros(N+1)
    h = T / N
    z[0] = z0
    for idx in range(1, N+1):
        z[idx] = phi(f, Df, t + idx*h, z[idx - 1], h)
    return z
evolve.py
potoole7/lotka_volterra
import numpy as np

def evolve_vec(phi, f, Df, t, z0, T, N):
    z = np.zeros((N+1,2))
    h = T / N
    z[0,:] = z0
    for idx in range(1, N+1):
        z[idx,:] = phi(f, Df, t + idx*h, z[idx - 1,:], h)
    return z
evolve_vec.py
potoole7/lotka_volterra
def fwd_euler(f, Df, t, z, h):
    return z + h*f(t, z)
fwd_euler.py
potoole7/lotka_volterra
import numpy as np
import matplotlib.pyplot as plt
import timestep_functions as ts

# max time
T = 100
# number of time steps
N = 100000
# timestep
# dt=T/N
# time (for plot)
t=np.linspace(0,T,N+1)
# ODE parameters
a=1.7
b=0.9
d=0.2
g=0.3

# arrays to store solutions
sol = np.zeros((N+1, 2))

# initial conditions
sol[0,0] = 10
sol[0,1] = 10

# Forward Euler solution
sol_fe = ts.f_euler(sol, N, T, a, b, g, d)

# Plot
fig, ax = plt.subplots(1,2)
ax[0].plot(t, sol_fe[:,0])
ax[0].plot(t, sol_fe[:,1])
ax[1].plot(sol_fe[:,0], sol_fe[:,1])
plt.show()


lotka.py
potoole7/lotka_volterra
import numpy as np

def lv(t,x, y, alpha, beta, delta, gamma):
    return np.array([
                     alpha*x - beta*x*y,
                     delta*x*y - gamma*y
                    ])

alpha, beta, delta, gamma = 0.1, 0.2, 0.3, 0.4
x0 = 1
y0 = 1
t_span = [0, 5]
lotka_volterra_test.py
potoole7/lotka_volterra
import numpy as np
import matplotlib.pyplot as plt

import fwd_euler
import evolve
import evolve_vec

# solve dz/dt = -z using forward euler method
z = evolve.evolve(phi=fwd_euler.fwd_euler, 
                  f=lambda t, z: -z,
                  Df=lambda t, z: np.array([-1]), 
                  t=0,
                  z0=1,
                  T=1,
                  N=50
                  )

# plot numerical solution of dz/dz=-z against actual e^(-z)
t_span=np.linspace(0,1,51)
plt.plot(t_span,z)
plt.plot(t_span, np.exp(-t_span), "--o")

# lotka-volterra parameters
a=1.1
b=0.5
d=0.4
g=0.3

lotka=evolve_vec.evolve_vec(phi=fwd_euler.fwd_euler,
                    f=lambda t,X: np.array([a*X[0]-b*X[0]*X[1],d*X[0]*X[1]-g*X[1]]),
                    Df=lambda t,x,y: np.array([0,0]),
                    t=0,
                    z0=np.array([10,10], ndmin=2),
                    T=50,
                    N=10000
                    )

fig, ax = plt.subplots()
ax.plot(lotka[:,0], lotka[:,1])

fig_p, ax_p = plt.subplots()
for k in range(5):
    lotka_phase=evolve_vec.evolve_vec(phi=fwd_euler.fwd_euler,
                    f=lambda t,X: np.array([a*X[0]-b*X[0]*X[1],d*X[0]*X[1]-g*X[1]]),
                    Df=lambda t,x,y: np.array([0,0]),
                    t=0,
                    z0=np.array([(k+1)*2,(k+1)*2], ndmin=2),
                    T=50,
                    N=100000
                    )
    ax_p.plot(lotka_phase[:,0], lotka_phase[:,1], label="X0=[{},{}]".format(k*2, k*2))
ax_p.set_xlabel("Prey")
ax_p.set_ylabel("Predator")
ax_p.legend()
plt.show()
main.py
potoole7/lotka_volterra
# Forward Euler Method
# Args: 
# sol: np.array of solutions, where:
#   sol[:, 0] is the prey population (x), and 
#   sol[:, 1] is the predator population (y)
# N: number of time steps
# T: max time
# a, b: prey parameters, max prey per capita growth rate, effect of presence 
# of predators on the prey growth rate.
# g, d: predator parameters, predator's per capita death rate, and the effect of 
# the presence of prey on the predator's growth rate.
def f_euler(sol, N, T, a, b, g, d):

   dt = T / N # time step

   for i in range(N):
     sol[i+1, 0] = sol[i, 0] + a * sol[i, 0]* dt - b * sol[i, 0] * sol[i, 1] * dt
     sol[i+1, 1] = sol[i, 1] - g * sol[i, 1]* dt + d * sol[i, 0] * sol[i, 1] * dt

   return sol

timestep_functions.py
potoole7/lotka_volterra