Permalink
Jan 29, 2021
Newer
100644
152 lines (119 sloc)
4 KB
1
#!/usr/bin/env python3
2
# -*- coding: utf-8 -*-
3
"""
4
trirep.py
5
Created on Thu May 9 16:23:30 2019
6
7
@author: nolte
8
Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019)
9
10
Population dynamics
11
"""
12
13
import numpy as np
14
from scipy import integrate
15
from matplotlib import pyplot as plt
16
17
plt.close('all')
18
19
def tripartite(x,y,z):
20
21
sm = x + y + z
22
xp = x/sm
23
yp = y/sm
24
25
f = np.sqrt(3)/2
26
27
y0 = f*xp
28
x0 = -0.5*xp - yp + 1;
29
30
plt.figure(2)
31
lines = plt.plot(x0,y0)
32
plt.setp(lines, linewidth=0.5)
33
plt.plot([0, 1],[0, 0],'k',linewidth=1)
34
plt.plot([0, 0.5],[0, f],'k',linewidth=1)
35
plt.plot([1, 0.5],[0, f],'k',linewidth=1)
36
plt.show()
37
38
39
def solve_flow(y,tspan):
40
def flow_deriv(y, t0):
41
#"""Compute the time-derivative ."""
42
43
f = np.zeros(shape=(N,))
44
for iloop in range(N):
45
ftemp = 0
46
for jloop in range(N):
47
ftemp = ftemp + A[iloop,jloop]*y[jloop]
48
f[iloop] = ftemp
49
50
phitemp = phi0 # Can adjust this from 0 to 1 to stabilize (but Nth population is no longer independent)
51
for loop in range(N):
52
phitemp = phitemp + f[loop]*y[loop]
53
phi = phitemp
54
55
yd = np.zeros(shape=(N,))
56
for loop in range(N-1):
57
yd[loop] = y[loop]*(f[loop] - phi);
58
59
if np.abs(phi0) < 0.01: # average fitness maintained at zero
60
yd[N-1] = y[N-1]*(f[N-1]-phi);
61
else: # non-zero average fitness
62
ydtemp = 0
63
for loop in range(N-1):
64
ydtemp = ydtemp - yd[loop]
65
yd[N-1] = ydtemp
66
67
return yd
68
69
# Solve for the trajectories
70
t = np.linspace(0, tspan, 701)
71
x_t = integrate.odeint(flow_deriv,y,t)
72
return t, x_t
73
74
# model_case 1 = zero diagonal
75
# model_case 2 = zero trace
76
# model_case 3 = asymmetric (zero trace)
77
print(' ')
78
print('trirep.py')
79
print('Case: 1 = antisymm zero diagonal')
80
print('Case: 2 = antisymm zero trace')
81
print('Case: 3 = random')
82
model_case = int(input('Enter the Model Case (1-3)'))
83
84
N = 3
85
asymm = 3 # 1 = zero diag (replicator eqn) 2 = zero trace (autocatylitic model) 3 = random (but zero trace)
86
phi0 = 0.001 # average fitness (positive number) damps oscillations
87
T = 100;
88
89
90
if model_case == 1:
91
Atemp = np.zeros(shape=(N,N))
92
for yloop in range(N):
93
for xloop in range(yloop+1,N):
94
Atemp[yloop,xloop] = 2*(0.5 - np.random.random(1))
95
Atemp[xloop,yloop] = -Atemp[yloop,xloop]
96
97
if model_case == 2:
98
Atemp = np.zeros(shape=(N,N))
99
for yloop in range(N):
100
for xloop in range(yloop+1,N):
101
Atemp[yloop,xloop] = 2*(0.5 - np.random.random(1))
102
Atemp[xloop,yloop] = -Atemp[yloop,xloop]
103
Atemp[yloop,yloop] = 2*(0.5 - np.random.random(1))
104
105
tr = np.trace(Atemp)
106
A = Atemp
107
for yloop in range(N):
108
A[yloop,yloop] = Atemp[yloop,yloop] - tr/N
109
110
else:
111
Atemp = np.zeros(shape=(N,N))
112
for yloop in range(N):
113
for xloop in range(N):
114
Atemp[yloop,xloop] = 2*(0.5 - np.random.random(1))
115
116
tr = np.trace(Atemp)
117
A = Atemp
118
for yloop in range(N):
119
A[yloop,yloop] = Atemp[yloop,yloop] - tr/N
120
121
plt.figure(3)
122
im = plt.matshow(A,3,cmap=plt.cm.get_cmap('seismic')) # hsv, seismic, bwr
123
cbar = im.figure.colorbar(im)
124
125
M = 20
126
delt = 1/M
127
ep = 0.01;
128
129
tempx = np.zeros(shape = (3,))
130
for xloop in range(M):
131
tempx[0] = delt*(xloop)+ep;
132
for yloop in range(M-xloop):
133
tempx[1] = delt*yloop+ep
134
tempx[2] = 1 - tempx[0] - tempx[1]
135
136
x0 = tempx/np.sum(tempx); # initial populations
137
138
tspan = 70
139
t, x_t = solve_flow(x0,tspan)
140
141
y1 = x_t[:,0]
142
y2 = x_t[:,1]
143
y3 = x_t[:,2]
144
145
plt.figure(1)
146
lines = plt.plot(t,y1,t,y2,t,y3)
147
plt.setp(lines, linewidth=0.5)
148
plt.show()
149
plt.ylabel('X Position')
150
plt.xlabel('Time')
151
152
tripartite(y1,y2,y3)