-
Notifications
You must be signed in to change notification settings - Fork 0
/
lotka-voltera.py
executable file
·127 lines (105 loc) · 3.25 KB
/
lotka-voltera.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#!/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
def odeSolver(t0, y0, dfFunc, h, nSteps, solverStepFunc):
"""This is a general ODE solver that takes the
derivative df/dt (dfFunc) and the algorithm for one time
step (solverStepFunc) as function arguments.
t0 = Initial time
y0 = Initial function values (array)
nSteps = Total number of integration steps
solverStepFunc = Euler Method, Midpoint RK or RK4
"""
yn = y0
tn = t0
tlist = [t0]
ylist = [y0]
for n in range(nSteps):
yn1 = solverStepFunc(tn, yn, h, dfFunc)
tn1 = tn + h
tlist.append(tn1)
ylist.append(yn1)
tn = tn1
yn = yn1
return (np.array(tlist), np.array(ylist))
def eulerStep(tn, yn, h, dfdt):
yn1 = yn + h * dfdt(tn, yn)
return yn1
def MidPointRK2Step(tn, yn, h, dfdt):
y_n1 = yn + h * dfdt(tn + h / 2, yn + (h / 2) * dfdt(tn, yn))
return y_n1
def RK4Step(tn, yn, h, dfdt):
pass
def LotkaVolterra(t, y):
"""
Implements the Lotka Volterra System dm/dt and dfox/dt
where y=(m,fox), dy/dt=(dm/dt, dfox/dt)
"""
K_m = 2
K_mf = 0.02
K_fm = 0.01
K_f = 1.06
m, f = y
dmdt = K_m * m - K_mf * m * f
dfoxdt = -K_f * f + K_fm * f * m
return np.array([dmdt, dfoxdt])
if __name__ == "__main__":
h = 0.001
nSteps = 60000
t0 = 0
y0 = np.array([100, 15]) # y0[0] = mice, y0[1] = foxes
init_cond_phase = np.linspace(
1.0, 100.0, 25
) # initial conditions for mice population (prey)
t, y = odeSolver(t0, y0, LotkaVolterra, h, nSteps, eulerStep)
plt.figure()
plt.grid()
plt.title("Forward Euler method")
plt.plot(t, y[:, 0], "xb", label="Mice") # population versus time
plt.plot(t, y[:, 1], "+r", label="Foxes") # population versus time
plt.xlabel("Time t, [days]")
plt.ylabel("Population")
plt.legend()
plt.savefig("euler_pop.png")
plt.show()
# phase plot
plt.figure()
for mice in init_cond_phase:
X0 = np.array([mice, 15])
t, y = odeSolver(t0, X0, LotkaVolterra, h, nSteps, eulerStep)
plt.plot(y[:, 0], y[:, 1], "-", label="$m_0 =$" + str(X0[0]))
plt.xlabel("Mice")
plt.ylabel("Foxes")
plt.legend()
plt.title("Mice vs Foxes")
plt.savefig("euler_phase.png")
plt.show()
t, y = odeSolver(t0, y0, LotkaVolterra, h, nSteps, MidPointRK2Step)
plt.figure()
plt.grid()
plt.title("Midpoint Runge-Kutta method")
plt.plot(t, y[:, 0], "xb", label="Mice") # population versus time
plt.plot(t, y[:, 1], "+r", label="Foxes") # population versus time
plt.xlabel("Time t, [days]")
plt.ylabel("Population")
plt.legend()
plt.savefig("midpoint_pop.png")
plt.show()
# phase plot
plt.figure()
for mice in init_cond_phase:
X0 = np.array([mice, 15])
t, y = odeSolver(t0, X0, LotkaVolterra, h, nSteps, MidPointRK2Step)
plt.plot(y[:, 0], y[:, 1], "-", label="$m_0 =$" + str(X0[0]))
plt.xlabel("Mice")
plt.ylabel("Foxes")
plt.legend()
plt.title("Mice vs Foxes")
plt.savefig("midpoint_phase.png")
plt.show()
# t, y = odeSolver(t0, y0, LotkaVolterra, h, nSteps, RK4Step)
#
# plt.plot(t, y[:,0], ..)
# ... # phase plot
#
# plt.show()