Difference between revisions of "Figure 4.20: Simulation of the predator-prey system"
Jump to navigation
Jump to search
Line 38: | Line 38: | ||
# Create a nonlinear I/O system | # Create a nonlinear I/O system | ||
− | predprey_sys = ct. | + | predprey_sys = ct.nlsys(predprey_update, states=2) |
# Simulate a trajectory leading to a limit cycle | # Simulate a trajectory leading to a limit cycle | ||
Line 57: | Line 57: | ||
# Generate a phase portrait | # Generate a phase portrait | ||
− | fbs.figure('mlh') | + | fbs.figure('mlh') |
− | def pp_ode(x, t): | + | ct.phaseplot.equilpoints(predprey_sys, [-5, 126, -5, 100]) |
− | return predprey_update(t, x, 0, {}) | + | ct.phaseplot.streamlines( |
− | ct.phase_plot(pp_ode, [0, 60, 7], [0, 50, 6]) | + | predprey_sys, np.array([ |
− | ct.phase_plot(pp_ode, [0, 60, 7], [60, 100, 4]) | + | [0, 100], [1, 0], |
− | ct.phase_plot(pp_ode, [70, 120, 6], [0, 50, 6]) | + | ]), 10, color='b') |
− | ct.phase_plot(pp_ode, [70, 120, 6], [60, 100, 4]) | + | ct.phaseplot.streamlines( |
+ | predprey_sys, np.array([[124, 1]]), np.linspace(0, 10, 500), color='b') | ||
+ | ct.phaseplot.streamlines( | ||
+ | predprey_sys, np.array([[125, 25], [125, 50], [125, 75]]), 3, color='b') | ||
+ | ct.phaseplot.streamlines(predprey_sys, np.array([2, 8]), 6, color='b') | ||
+ | ct.phaseplot.streamlines( | ||
+ | predprey_sys, np.array([[20, 30]]), np.linspace(0, 65, 500), | ||
+ | gridtype='circlegrid', gridspec=[2, 1], arrows=10, color='r') | ||
+ | ct.phaseplot.vectorfield(predprey_sys, [5, 125, 5, 100], gridspec=[20, 20]) | ||
+ | |||
+ | # Add the limit cycle | ||
+ | resp1 = ct.initial_response(predprey_sys, np.linspace(0, 100), [20, 75]) | ||
+ | resp2 = ct.initial_response( | ||
+ | predprey_sys, np.linspace(0, 20, 500), resp1.states[:, -1]) | ||
+ | plt.plot(resp2.states[0], resp2.states[1], color='k') | ||
+ | |||
+ | # Legacy code | ||
+ | # def pp_ode(x, t): | ||
+ | # return predprey_update(t, x, 0, {}) | ||
+ | # ct.phase_plot(pp_ode, [0, 60, 7], [0, 50, 6]) | ||
+ | # ct.phase_plot(pp_ode, [0, 60, 7], [60, 100, 4]) | ||
+ | # ct.phase_plot(pp_ode, [70, 120, 6], [0, 50, 6]) | ||
+ | # ct.phase_plot(pp_ode, [70, 120, 6], [60, 100, 4]) | ||
− | # Plot the limit cycle | + | # # Plot the limit cycle |
− | ct.phase_plot(pp_ode, X0=sim.states[:, -1:].T, T=20) # limit cycle | + | # ct.phase_plot(pp_ode, X0=sim.states[:, -1:].T, T=20) # limit cycle |
− | ct.phase_plot(pp_ode, X0=[[120, 32], [120, 60]], T=20) # outside trajectories | + | # ct.phase_plot(pp_ode, X0=[[120, 32], [120, 60]], T=20) # outside trajectories |
− | ct.phase_plot(pp_ode, X0=[[19, 30]], T=75) # inside trajectories | + | # ct.phase_plot(pp_ode, X0=[[19, 30]], T=75) # inside trajectories |
# Label the plot | # Label the plot | ||
− | |||
plt.xlabel("Hares") | plt.xlabel("Hares") | ||
plt.ylabel("Lynxes") | plt.ylabel("Lynxes") | ||
plt.title("Phase portrait") | plt.title("Phase portrait") | ||
fbs.savefig('figure-4.20-predprey_ctstime-pp.png') # PNG for web | fbs.savefig('figure-4.20-predprey_ctstime-pp.png') # PNG for web | ||
+ | |||
</nowiki> | </nowiki> |
Revision as of 00:48, 7 April 2024
Chapter | Examples |
---|---|
Figure number | 4.20 |
Figure title | Simulation of the predator–prey system |
GitHub URL | https://github.com/murrayrm/fbs2e-python/blob/main/figure-4.20-predprey ctstime.py |
Requires | python-control, |
Figure 4.20: Simulation of the predator–prey system. The figure on the left shows a simulation of the two populations as a function of time. The figure on the right shows the populations plotted against each other, starting from different values of the population. The oscillation seen in both figures is an example of a limit cycle. The parameter values used for the simulations are a = 3.2, b = 0.6, c = 50, d = 0.56, k = 125, and r = 1.6.
# predprey_ctstime.py - Predator-prey model in continuous time # RMM, 28 May 2023 import matplotlib.pyplot as plt import numpy as np import control as ct import fbs # FBS plotting customizations # Define the dynamis for the predator-prey system (no input) def predprey_update(t, x, u, params={}): """Predator prey dynamics""" r = params.get('r', 1.6) d = params.get('d', 0.56) b = params.get('b', 0.6) k = params.get('k', 125) a = params.get('a', 3.2) c = params.get('c', 50) # Dynamics for the system dx0 = r * x[0] * (1 - x[0]/k) - a * x[1] * x[0]/(c + x[0]) dx1 = b * a * x[1] * x[0] / (c + x[0]) - d * x[1] return np.array([dx0, dx1]) # Create a nonlinear I/O system predprey_sys = ct.nlsys(predprey_update, states=2) # Simulate a trajectory leading to a limit cycle timepts = np.linspace(0, 70, 500) sim = ct.input_output_response(predprey_sys, timepts, 0, [25, 20]) # Plot the results fbs.figure('mlh') # FBS conventions plt.plot(sim.time, sim.states[0], 'b-', label="Hare") plt.plot(sim.time, sim.states[1], 'r--', label="Lynx") plt.legend() plt.xlabel("Time $t$ [years]") plt.ylabel("Population") plt.title("Time response") # Save the figure fbs.savefig('figure-4.20-predprey_ctstime-sim.png') # PNG for web # Generate a phase portrait fbs.figure('mlh') ct.phaseplot.equilpoints(predprey_sys, [-5, 126, -5, 100]) ct.phaseplot.streamlines( predprey_sys, np.array([ [0, 100], [1, 0], ]), 10, color='b') ct.phaseplot.streamlines( predprey_sys, np.array([[124, 1]]), np.linspace(0, 10, 500), color='b') ct.phaseplot.streamlines( predprey_sys, np.array([[125, 25], [125, 50], [125, 75]]), 3, color='b') ct.phaseplot.streamlines(predprey_sys, np.array([2, 8]), 6, color='b') ct.phaseplot.streamlines( predprey_sys, np.array([[20, 30]]), np.linspace(0, 65, 500), gridtype='circlegrid', gridspec=[2, 1], arrows=10, color='r') ct.phaseplot.vectorfield(predprey_sys, [5, 125, 5, 100], gridspec=[20, 20]) # Add the limit cycle resp1 = ct.initial_response(predprey_sys, np.linspace(0, 100), [20, 75]) resp2 = ct.initial_response( predprey_sys, np.linspace(0, 20, 500), resp1.states[:, -1]) plt.plot(resp2.states[0], resp2.states[1], color='k') # Legacy code # def pp_ode(x, t): # return predprey_update(t, x, 0, {}) # ct.phase_plot(pp_ode, [0, 60, 7], [0, 50, 6]) # ct.phase_plot(pp_ode, [0, 60, 7], [60, 100, 4]) # ct.phase_plot(pp_ode, [70, 120, 6], [0, 50, 6]) # ct.phase_plot(pp_ode, [70, 120, 6], [60, 100, 4]) # # Plot the limit cycle # ct.phase_plot(pp_ode, X0=sim.states[:, -1:].T, T=20) # limit cycle # ct.phase_plot(pp_ode, X0=[[120, 32], [120, 60]], T=20) # outside trajectories # ct.phase_plot(pp_ode, X0=[[19, 30]], T=75) # inside trajectories # Label the plot plt.xlabel("Hares") plt.ylabel("Lynxes") plt.title("Phase portrait") fbs.savefig('figure-4.20-predprey_ctstime-pp.png') # PNG for web