-
Notifications
You must be signed in to change notification settings - Fork 3
/
predator_prey.py
111 lines (83 loc) · 2.59 KB
/
predator_prey.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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# Copyright (c) 2018, sumeet92k
# All rights reserved. Please read the "license.txt" for license terms.
#
# Project Title: Materials Modeling Course Tutorials
#
# Developer: Sumeet Khanna
# Contact Info: https://github.com/sumeet92k
#
"""
TUTORIAL 2: predator prey model
Solves the coupled Lodka-Volterra (predator-prey) equations using different
numerical techniques
"""
import numpy as np
import matplotlib.pyplot as plt
alpha = 4
beta = 2
delta = 3
gamma = 3
dt = 0.001
timesteps = 20000
t_end = dt*(timesteps -1)
t_array = np.linspace(0, t_end, timesteps)
def get_dx_dt(x, y):
return alpha*x - beta*x*y
def get_dy_dt(x, y):
return delta*x*y - gamma*y
x_exp = np.zeros(timesteps)
y_exp = np.zeros(timesteps)
x_pc = np.zeros(timesteps)
y_pc = np.zeros(timesteps)
x_rk4 = np.zeros(timesteps)
y_rk4 = np.zeros(timesteps)
x_exp[0] = 10
y_exp[0] = 5
x_pc[0] = 10
y_pc[0] = 5
x_rk4[0] = 10
y_rk4[0] = 5
for t in range(timesteps - 1):
# explicit method
x_exp[t+1] = x_exp[t] + get_dx_dt(x_exp[t], y_exp[t] )*dt
y_exp[t+1] = y_exp[t] + get_dy_dt(x_exp[t], y_exp[t] )*dt
# predictor corrector method
x_tilde = x_pc[t] + get_dx_dt(x_pc[t], y_pc[t])*dt
y_tilde = y_pc[t] + get_dy_dt(x_pc[t], y_pc[t])*dt
x_pc[t+1] = x_pc[t] + 0.5*( get_dx_dt(x_pc[t], y_pc[t]) + get_dx_dt(x_tilde, y_tilde) )*dt
y_pc[t+1] = y_pc[t] + 0.5*( get_dy_dt(x_pc[t], y_pc[t]) + get_dy_dt(x_tilde, y_tilde) )*dt
# RK 4 method
k1_x = dt*get_dx_dt(x_rk4[t], y_rk4[t])
k1_y = dt*get_dy_dt(x_rk4[t], y_rk4[t])
x_tilde = x_rk4[t] + 0.5*k1_x
y_tilde = y_rk4[t] + 0.5*k1_y
k2_x = dt*get_dx_dt(x_tilde, y_tilde)
k2_y = dt*get_dy_dt(x_tilde, y_tilde)
x_tilde = x_rk4[t] + 0.5*k2_x
y_tilde = y_rk4[t] + 0.5*k2_y
k3_x = dt*get_dx_dt(x_tilde, y_tilde)
k3_y = dt*get_dy_dt(x_tilde, y_tilde)
x_tilde = x_rk4[t] + 1*k3_x
y_tilde = y_rk4[t] + 1*k3_y
k4_x = dt*get_dx_dt(x_tilde, y_tilde)
k4_y = dt*get_dy_dt(x_tilde, y_tilde)
x_rk4[t+1] = x_rk4[t] + (k1_x + 2*k2_x + 2*k3_x + k4_x)/6
y_rk4[t+1] = y_rk4[t] + (k1_y + 2*k2_y + 2*k3_y + k4_y)/6
plt.plot(x_exp, y_exp, label='explicit')
plt.plot(x_pc, y_pc, label='predictor-corrector')
plt.plot(x_rk4, y_rk4, label='RK4')
plt.xlabel("deer") # label the axes
plt.ylabel("lion")
plt.legend()
plt.show()
# plotting lion and deer population with time in a different figure
plt.figure()
plt.plot(t_array, x_rk4, label='deer')
plt.plot(t_array, y_rk4, label='lion')
plt.xlabel("time")
plt.ylabel("population")
plt.legend()
plt.show()