I'm trying to visualize the solution to my code by plotting - in python - a graph of the number of predators versus preys, and using a black circle to mark the current state. My code runs okay except I can't get the plotting portion of the graph to print the plots on the graph.
This is my current python code:
import matplotlib.pyplot as plt
import numpy as np
from trajectoryplotsupport import *
def ExactTrajectory(dt, ntimesteps):
tt = np.linspace(0, dt*ntimesteps, ntimesteps+1)
x = np.cos(2*np.pi*tt)
y = np.sin(2*np.pi*tt)
return x,y
def GetVelocity(pos, tn):
u =-2.0*np.pi * pos[1]
v = 2.0*np.pi * pos[0]
vel = np.array([u,v])
return vel
def RK3(p, t):
vel = GetVelocity(p, t) # stage 1 velocity
pt = p + 0.5*dt*vel # stage 1 provisional position
vel = GetVelocity(pt, t+dt) # stage 2 velocity
pt = 0.25 * p + 0.75*( pt + dt*vel ) # stage 2 provisional position
vel = GetVelocity(pt, t+0.5*dt) # stage 3 velocity
p = ( p + 2.0*(pt + dt*vel) )/3.0 # stage 3 provisional position
return p
def PreyPred(p, t):
b = 0.5 # prey birth rate
m = 0.2 # predator mortality
a = 0.04 # effect of predator abundance on prey
r = 0.01 # effect of prey abundance on predator
h, w = p[0], p[1] # initial densities of prey and predator are h(0) = 4.0 and w(0) = 8.0
dht = h*(b - a*w)
dwt = w*(-m + r*h)
return np.array([dht, dwt])
#################MAIN CODE#############
p = np.array([4.0, 8.0]) # initial position
ntimesteps = 120 # number of time steps
finalTime = 480 # final time
dt = 1 # time step
xe,ye = ExactTrajectory(dt, ntimesteps) # call exact trajectory
# Check on computations. Create 2 array to track the positions and initialize
pos = np.zeros(2) # create current position array
pos[0] = xe[0]; pos[1]=ye[0] # initialize current position
fig, lastPosition = InitialPlot(p, xe,ye)
for it in range(finalTime): # time loop
t = it*dt # set current time
p = RK3(p, t)
print(t, p)
UpdatePlot(fig,lastPosition,p,t) # update the plot
ttpoints = np.arange(0, 480, dt) #edited
hpoints, wpoints = [], []
p = np.array([0, 0], float)
for tt in ttpoints:
hpoints.append(p[0])
wpoints.append(p[1])
p += RK3(p, tt)
plt.plot(ttpoints, hpoints)
plt.plot(ttpoints, wpoints)
plt.xlabel("Prey")
plt.ylabel("Predator")
plt.plot(hpoints, wpoints)
plt.show()