Commit 6fc6d3701de35538ac327b0533a383931e216970

Some kind of a Python-version of predatorprey.m
demo2/predatorprey.py
(44 / 0)
  
1import numpy as np
2import matplotlib.pyplot as plt
3
4def f(w0, w1, a, b):
5 return np.array([ a*w0*(1-w1), w1*(w0-1)-b*w1**2 ])
6
7def predatorprey(u0, v0, dt, tmax, a, b):
8 print "____"
9 nt = np.floor(tmax/dt)
10 nt = int(nt)
11 print "nt = ", nt
12 ws = np.zeros((nt, 2))
13 for i in range(nt):
14 for j in range(2):
15 ws[i][j] = np.nan
16 w = [u0, v0]
17 print "w = ", w
18 ws[0,:] = w
19 for k in range(nt):
20 #print "k = ", k
21 k1 = f(w[0], w[1], a, b)
22 k2 = f(w[0]+0.5*dt*k1[0], w[1]+0.5*dt*k1[1], a, b)
23 k3 = f(w[0]+0.5*dt*k2[0], w[1]+0.5*dt*k2[1], a, b)
24 k4 = f(w[0]+dt*k3[0], w[1]+dt*k3[1], a, b)
25 w[0] = w[0] + (1.0/6.0)*dt*(k1[0]+2*k2[0]+2*k3[0]+k4[0])
26 #print "w[0] = ", w[0]
27 w[1] = w[1] + (1.0/6.0)*dt*(k1[1]+2*k2[1]+2*k3[1]+k4[1])
28 #print "w[1] = ", w[1]
29 ws[k,:] = w
30 fig = plt.figure()
31 ax = fig.add_subplot(111)
32 #ax.plot(ws[0,0], ws[0,1], 'k0')
33 ax.plot(ws[:,0], ws[:,1], 'k-')
34 #ax.plot(ws[-1,0],ws[-1,1],'k+')
35
36 plt.show()
37
38def main():
39 predatorprey(0.9, 0.9, 0.01, 200, 1, 0)
40 predatorprey(0.9, 0.9, 0.01, 200, 1, +0.1)
41 predatorprey(0.9, 0.9, 0.01, 50, 1, -0.1)
42 predatorprey(0.5, 0.5, 0.01, 50, 1, +0.1)
43
44main()