Python Matplotlib (Plots, Customisation, Logistic Map, ODEs)

Python Matplotlib

Data visualisation is the graphical representation of information and data. By employing visual elements such as charts, graphs and maps, identify patterns, trends and outliers.

Original page had a screenshot here.

pyplot as plt

import matplotlib.pyplot as plt

matplotlib.pyplot is a collection of 2D plotting functions for visualising data from the simplest plots to highly customisable polished works of art. Comes with the Anaconda distribution.

Normal plots

General formatting

Adding data

x  = [1980, 1990, 2000, 2010, 2020, 2030, 2040, 2050]
y1 = [1.00, 1.18, 1.29, 1.37, 1.44, 1.46, 1.45, 1.40]
y2 = [0.70, 0.87, 1.06, 1.23, 1.38, 1.50, 1.59, 1.64]
y3 = [0.23, 0.25, 0.28, 0.31, 0.33, 0.35, 0.37, 0.38]   # If >1 lines
 
# Generates a numpy array of numbers evenly spaced over a specified interval
x = np.linspace(-2, 5, 50)   # 50 data points including -2 and 5

plot()

plt.plot(x, y)
# x and y are data along the x and y axes respectively
# X VALUES ARE OPTIONAL; if absent, the default will be [0., 1., 2., 3., 4.]
 
# Set colour, marker and line styles for legend
plt.plot(x, y1, "rx-", label="China")
plt.plot(x, y2, "bo--", label="India")
plt.plot(x, y3, "g^:", label="USA")
 
# One-line to plot 2 graphs
plt.plot(x, np.sin(x), "r-", x, np.cos(x), "b-")

plt.xlim()

# Set range for x and y axes
plt.xlim(-2, 20)
plt.ylim(-20, 20)

Different types of graphs

Linear, Quadratic, Cubic

y1 = 2 * x - 7                  # Linear
y2 = x**2 - 5*x - 10            # Quadratic
y3 = x**3 - 3*x**2 - 2*x - 5    # Cubic
 
plt.plot(x, y1, "r-o")
plt.plot(x, y2, "g-s")
plt.plot(x, y3, "b-^")

Trig

x = np.linspace(0, np.pi*4, 100)   # x data is given in radians
 
plt.plot(x, np.sin(x), "r-")   # Plot sine; take note that x is in radians
plt.plot(x, np.cos(x), "b-")
 
# π is unicode for pi symbol

Semi-log

# One axis log, one axis normal
# Year data
x = [1900, 1910, 1920, 1930, 1940, 1950, 1960, 1970, 1980, 1990, 2000, 2010, 2017]
 
# Stock market index data (taken at the end of every decade)
y = [68, 81, 71, 244, 151, 200, 615, 809, 824, 2633, 10787, 11577, 20656]
 
# semilogy() plots y-axis in log (base 10) scale, x-axis linear.
plt.semilogy(x, y)

Multi-line graphs

Single line

import matplotlib.pyplot as plt
 
x = [1600, 1700, 1800, 1900, 2000]
y = [0.2, 0.5, 1.1, 2.2, 7.7]
 
plt.plot(x, y)
 
plt.title("World Population over Time")
plt.xlabel("Year")
plt.ylabel("Population (in billions)")
plt.show()

Multi line

x  = [1980, 1990, 2000, 2010, 2020, 2030, 2040, 2050]
y1 = [1.00, 1.18, 1.29, 1.37, 1.44, 1.46, 1.45, 1.40]
y2 = [0.70, 0.87, 1.06, 1.23, 1.38, 1.50, 1.59, 1.64]
y3 = [0.23, 0.25, 0.28, 0.31, 0.33, 0.35, 0.37, 0.38]
 
plt.figure(figsize=(7, 5))   # Set figure size
 
plt.plot(x, y1, "rx-", label="China")
plt.plot(x, y2, "bo--", label="India")
plt.plot(x, y3, "g^:", label="USA")
 
plt.title("Populations of China, India & USA and their projections",
          fontname="Times New Roman", fontsize=15)
plt.xlabel("Year", fontname="Times New Roman", fontsize=14)
plt.ylabel("Population (in billions)", fontname="Times New Roman", fontsize=14)
plt.legend()
plt.show()

Other plots (Bar, Scatter, Pie, Error)

Bar charts

Customisation

width = 0.25
plt.bar(x, y, color="skyblue", width=0.7)
plt.xticks(x, ['Team A', 'Team B', 'Team C', 'Team D', 'Team E'])

Bar chart

x = [20, 36, 24, 12, 8]
y = ["Cat", "Dog", "Hamster", "Rabbit", "Terrapin"]
 
plt.figure(figsize=(5, 4))
 
plt.bar(x, y)     # vertical display
plt.barh(x, y)    # horizontal display

Clustered bar chart

x  = np.arange(5)
y1 = [45, 56, 32, 89, 67]
y2 = [32, 56, 78, 48, 90]
y3 = [24, 37, 45, 55, 89]
 
plt.figure(figsize=(6, 4))
 
width = 0.25
plt.bar(x - width, y1, width, color="olive")
plt.bar(x,         y2, width, color='darkkhaki')
plt.bar(x + width, y3, width, color='khaki')

Stacked bar chart

x  = ['Jul', 'Aug', 'Sep', 'Oct']
y1 = np.array([36, 28, 26, 40])
y2 = np.array([22, 24, 22, 28])
y3 = np.array([6, 5, 7, 5])
 
plt.figure(figsize=(5, 4))
 
plt.bar(x, y1,            color='olivedrab', width=0.7)
plt.bar(x, y2, bottom=y1, color='khaki',     width=0.7)
plt.bar(x, y3, bottom=y1+y2, color='teal',   width=0.7)

Scatter plot

x = [14.2, 16.4, 11.9, 12.5, 18.9, 22.1, 19.4, 23.1, 25.4, 18.1, 22.6, 17.2]
y = [215.20, 325.00, 185.20, 330.20, 418.60, 520.25, 412.20, 614.60, 544.80, 421.40, 445.50, 408.10]
 
plt.scatter(x, y, c='red', alpha=0.7)   # c — color; alpha — transparency

Pie chart

x = ['Carbon Monoxide', 'Nitrogen oxides', 'Sulfur oxides', 'Volatile organics', 'Particulates']
y = [49.1, 14.8, 16.4, 13.7, 6.0]
 
fig = plt.figure(figsize=(6, 5))
 
plt.pie(data, labels=pollutants, autopct='%.1f%%',
        colors=["lightblue", "green", "khaki", "salmon", "gold"])

Error bars

x = np.arange(11) / 10
y = (x + 0.1)**3
y_error = 0.05
 
plt.plot(x, y)
plt.errorbar(x, y, yerr=y_error, fmt="o")

Multi-plot

plt.subplot

plt.subplot(rows, cols, plot number)
 
plt.subplot(2, 2, 1)
plt.subplot(2, 2, 2)
plt.subplot(2, 2, 3)
plt.subplot(2, 2, 4)

plt.figure(figsize=(7, 3), facecolor='lightblue')
 
plt.subplot(1, 2, 1)
plt.plot(x, np.sin(x), "r-")
plt...
 
plt.subplot(1, 2, 2)
plt.plot(x, np.cos(x), "b-")
plt...
 
plt.tight_layout()
plt.show()

axes.plot()

x = np.linspace(0, np.pi*4, 100)
 
fig, ax = plt.subplots()
 
ax.plot(x, np.sin(x), "r-")
ax.plot(x, np.cos(x), "b-")
 
ax.set_title("Graph of sine and cosine functions")
ax.legend(("sine", "cosine"), loc=0)
ax.set_xlabel("x (in radians)")
ax.set_ylabel("Trigonometry functions")
 
ax.set_xlim(0, 4*np.pi)
ax.set_ylim(-1, 1)
 
ax.set_xticks(np.arange(0, np.pi*5, np.pi), ('0', 'π', '2π', '3π', '4π'))
ax.set_yticks([-1, -0.5, 0, 0.5, 1])
ax.minorticks_on()
ax.grid(which='major', color='black')
ax.grid(which='minor', color='gray', linestyle='--')

Customisation

Markers

Specifier   Line/Marker Style
'-'         # Solid line
'--'        # Dashed line
'-.'        # Dash-dot line
':'         # Dotted line
 
'o'         # Circle marker
'^'         # Triangle marker
's'         # Square marker
'h'         # Hexagon marker
'x'         # Cross marker
 
Colour Abbreviation/Colour
b           # blue
c           # cyan
g           # green
k           # black
m           # magenta
r           # red
w           # white
y           # yellow

plt.plot(x, y1, "r-o")
plt.plot(x, y2, "g-s")
 
plt.plot(x, y, alpha=x)         # 0 < x < 1 -- Opacity
plt.plot(x, y, markersize=x)

# Arrow ----->
plt.arrow(x0, 0, 0, f(x0)/2, head_width=0.02, head_length=0.02,
          fc='g', ec='g', linewidth=0)

Legend position

0   # Best location
1   # Upper right
2   # Upper left
3   # Lower left
4   # Lower right

plt.legend()
plt.legend(loc=0)
plt.legend(("Linear", "Quadratic", "Cubic"))
plt.legend(("Linear", "Quadratic", "Cubic"), loc=1)

Titles and figure

plt.title("World Population over Time")
plt.title("Populations of China, India & USA and their projections",
          fontname="Times New Roman", fontsize=15)
 
plt.figure(figsize=(7, 5))

Axis labels

plt.xlabel("Year")
plt.ylabel("Population (in billions)")
plt.xlabel("Year", fontname="Times New Roman", fontsize=14)
plt.ylabel("Population (in billions)", fontname="Times New Roman", fontsize=14)
plt.xlabel(r'$x_n$')
plt.ylabel("Life Expectancy", fontsize=12); plt.xlabel("Year", fontsize=12)

Ticks

plt.xticks(np.arange(0, np.pi*5, np.pi))
plt.xticks(np.arange(0, np.pi*5, np.pi), ('0', 'π', '2π', '3π', '4π'))
plt.yticks([-1, -0.5, 0, 0.5, 1])
 
plt.minorticks_on()
plt.grid(which='major', color='black')
plt.grid(which='minor', color='gray', linestyle='--')
plt.xticks(rotation=25)

Multiplot helpers

plt.tight_layout()
plt.grid(True)
plt.grid(True, which="both")
 
ax.grid(which='major', color='black')
ax.grid(which='minor', color='gray', linestyle='--')

Others

plt.show()
plt.grid(True)
plt.grid(True, which="both")
plt.tight_layout()

Saving graphs

plt.savefig("Graphs.png", format="png")
plt.savefig("Graphs.jpg", format="jpg")

Handling data sets

Sources: Google Dataset Search, Kaggle, Statista, Our World in Data, Global Health Observatory, World Bank Open Data.

CSV files — most common file format. Stands for Comma Separated Values. Each line is a data record consisting of fields separated by commas (delimiters).

Pandas

import pandas as pd

Every Python program that reads and works with data (like CSV) uses this module. Pandas came from “panel data”, referring to tabular data.

Basic

Analyse books data

df = pd.read_csv('books.csv')
df

Gaining insight

display(df.head())       # Output the first 5 rows; df.head(n): first n rows
display(df.tail())       # Output the last 5 rows; df.tail(n): last n rows
display(df.info())       # Output information about dataframe
display(df.describe())   # Output descriptive statistics
display(len(df))         # Output number of rows
display(df.shape)        # Output rows vs columns
display(df.count())      # Count non-null data for each field or column

Slice using iloc

display(df.iloc[10])         # Output 11th data record
display(df.iloc[3:7])        # Output 4th to 7th records
display(df.iloc[:, 1])       # Output second column
display(df.iloc[3:7, 0:2])

Slice by column name

display(df['Genre'])
display(df.Genre)
display(df[['Title', 'Genre']])

Filter rows

df_history = df[df['Genre'] == 'history']
display(df_history)
df_history.to_csv('history.csv')

Multiple conditions with &

df_penguin_big = df[(df['Publisher'] == 'Penguin') & (df['Height'] > 200)]
df_penguin_big

Multiple conditions with |

df[(df['Publisher'] == 'Penguin') | (df['Height'] > 200)]

More examples

df[(df['Genre'] == 'fiction') & (df['Title'] >= 'T')]
df[(df['Genre'] == 'fiction') & (df['Height'] < 200)]

Example — Life Expectancy data

data = pd.read_csv('LifeExpectancy.csv')
display(data)
display(data.info())
display(data.describe())

df_THA = data[data['Code'] == 'THA']
display(df_THA, type(df_THA))

df_MEX = data[data['Code'] == 'MEX']
display(df_MEX, type(df_MEX))
display(len(df_MEX))

plt.figure(figsize=(7, 5))
plt.plot(df_THA["Year"], df_THA["Life expectancy"], "b-", label="Thailand")
plt.plot(df_MEX.Year,    df_MEX["Life expectancy"], "r-", label="Mexico")
plt.title("Life expectancy of Thailand and Mexico", fontsize=13)
plt.ylabel("Life Expectancy", fontsize=12); plt.xlabel("Year", fontsize=12)
plt.legend()
plt.show()

plt.figure(figsize=(7, 5))
plt.plot(df_THA.iloc[:, 2], df_THA.iloc[:, 3], "b-", label="Thailand")
plt.plot(df_MEX.iloc[:, 2], df_MEX.iloc[:, 3], "r-", label="Mexico")
plt.title("Life expectancy of Thailand and Mexico", fontsize=13)
plt.ylabel("Life Expectancy", fontsize=12); plt.xlabel("Year", fontsize=12)
plt.legend()
plt.show()

Logistic Map — Introduction

The Logistic Map models how a population changes over time, especially with resource limits.

Original page had a screenshot here.

Formula: X_{n+1} = r · X_n · (1 - X_n) where:

• X_n is the population at time n (0..1)
• r is growth rate
• X_{n+1} is the next population

What happens with different r:

• Small r — population can’t grow, shrinks to zero.
• Moderate r — population stabilises.
• Higher r — oscillates between values.
• Very high r — chaotic.

Constraint: r > 0; r < 4 (else explodes).

Example: X0 = 0.5, r = 3.2 → 0.8, 0.512, 0.802, 0.507...

Plotting Logistic Map

Xn against Xn+1

x = np.linspace(0, 1, 100)
r = 4
y = r*x * (1 - x)
 
plt.plot(x, y)
plt.xlabel("$x_n$")
plt.ylabel(r"$x_{n+1}$")
 
# r is raw string formatting
# f is variable formatting
 
r"$\alpha$"        # α
f"Price = {r}"     # Price = 4
rf"$\alpha$ = {r}" # α = 4
 
$X_{n+1}$ # X and subscript n+1
$X_n$     # X and subscript n
 
$\alpha$  # α
$\pi$     # π

Labels and titles

plt.xlabel(r'Time step $n$', fontsize=15)
plt.ylabel(r'$x_n$', fontsize=15)
plt.title(rf'Logistic map at $\alpha$ = {a1} and {a2}', fontsize=15)

Logistic Map — Visual (time)

def logistic(r, y):
    return r * y * (1 - y)
 
a = 1.6              # Try 0.5, 2.5, 3.2, 3.5, 3.55, 3.569, 3.78, 3.83, 3.99
x = [0.2]            # Data starts from 0.2
N = 200              # Number of iterations / points
transients = 150     # Momentary / initial variations
 
# Iterate N times; total N+1 pops including initial pop
for n in range(1, N+1):
    x.append(logistic(a, x[n-1]))
 
# Plot with transients
plt.subplot(1, 2, 1)
plt.plot(x, 'ro', x, 'b')
plt...
 
# Plot without transients
plt.subplot(1, 2, 2)
y = x[transients:]
plt.plot(y, 'ro', y, 'b')
plt...

Behaviour: low r → extinction; moderate r → equilibrium; high r → oscillation; very high r → chaos. Discarding transients lets you focus on the steady state.

Logistic Map — Visual (correlation)

Plot Xn+1 against Xn to distill the nature of oscillations.

a = 1.6
x = [0.2]
N = 200
transients = 150
 
for n in range(1, N+1):
    x.append(logistic(a, x[n-1]))
 
y1 = x[transients : N-2]
print(len(y1))
 
y2 = x[transients+1 : N-1]
print(len(y2))
 
plt.plot(y1, y2, 'ro', alpha=0.2, markersize=5)
 
print('xn\t\tx(n+1)')
print("-" * 30)
for n in range(N-9, N+1):
    print(f"({x[n-1]:.10f}, {x[n]:.10f})")

Sensitivity to initial conditions

x1 = [0.2]
x2 = [0.200001]
N  = 100
a  = 3.7
 
for n in range(0, N):
    x1.append(logistic(a, x1[n]))
    x2.append(logistic(a, x2[n]))
 
plt.plot(x1, '-ro', label=f'$x_0$ = {x1[0]}')
plt.plot(x2, '-bo', label=f'$x_0$ = {x2[0]:.6f}')

Sensitivity to system parameter

x1 = [0.2]
x2 = [0.2]
N  = 100
a1 = 3.9
a2 = 3.900001
 
for n in range(0, N):
    x1.append(logistic(a1, x1[n]))
    x2.append(logistic(a2, x2[n]))
 
plt.plot(x1, '-ro', label=rf'$\alpha$ = {a1}')
plt.plot(x2, '-bo', label=rf'$\alpha$ = {a2:.6f}')

Logistic Map — Visual (cobweb)

Plotting lines using 2 points

plt.figure(figsize=(9, 3))
 
plt.subplot(1, 3, 1)
plt.plot([1, 5], [2, 10])
 
plt.subplot(1, 3, 2)
plt.plot([3, 3], [0, 5])
plt.xticks(range(0, 6))
 
plt.subplot(1, 3, 3)
plt.plot([0, 5], [3, 3])
plt.yticks(range(0, 6))

Cobweb diagram

def cobweb(f, x0, iterno, x_min=0, x_max=1):
 
    N = 100
    xx = np.linspace(x_min, x_max, N)
 
    plt.plot(xx, f(xx), 'k')    # Plot logistic function via lambda
    plt.plot(xx, xx, "k--")     # Plot y = x
 
    plt.plot([x0, x0], [0, f(x0)], 'g', linewidth=1)
    plt.arrow(x0, 0, 0, f(x0)/2, head_width=0.02, head_length=0.02,
              fc='g', ec='g', linewidth=0)
 
    x, y = x0, x0
 
    for n in range(iterno):
 
        if 0 < n < (iterno):
 
            plt.plot([x, y], [y, y], 'g', linewidth=1)
            plt.arrow(x, y, (y-x)/2, 0, head_width=0.02, head_length=0.02,
                      fc='g', ec='g', linewidth=0)
 
            plt.plot([y, y], [y, f(y)], 'g', linewidth=1)
            plt.arrow(y, y, 0, (f(y)-y)/2, head_width=0.02, head_length=0.02,
                      fc='g', ec='g', linewidth=0)
 
        x, y = y, f(y)
        plt.scatter(x, f(x), marker='o', s=30, c='b')
 
plt.figure(1, (10, 6))
plt.xlabel("$x_{n}$"); plt.ylabel("$x_{n+1}$")
 
r = 3.99
f = lambda x: r*x*(1 - x)
cobweb(f, 0.2, 7)

Logistic Map — Visual (bifurcation)

Reset variables without prompt

%reset

Bifurcation diagram

rLow, rHigh = 0, 3.99
 
plt.figure(figsize=(15, 10))
 
nTransient = 400
nIteration = 800
nStep      = 400
 
for r in np.linspace(rLow, rHigh, nStep):
 
    x = [np.random.rand()]
 
    for n in range(nIteration):
        x.append(r * x[n] * (1 - x[n]))
 
    y = x[nTransient:]
 
    rvalue = r * np.ones(len(y))
    plt.plot(rvalue, y, "k,")
 
plt.title(f"Bifurcation diagram for $r$ in [{rLow}, {rHigh}] of the Logistic map $f(x) = rx(1-x)$")
plt.xlabel("Intrinsic growth rate $r$", fontsize=12)
plt.ylabel("$x_n$", fontsize=12)

Some r values lead to a few unique Xn (limit cycles). As r increases, solutions bifurcate until they fill the continuum (chaos). Windows of calm at r = 3.83–3.85.

Ordinary Differential Equations (ODE)

odeint()

y = odeint(model, y0, t)
 
# model() returns derivative values at requested y and t values.
# y0 — initial value(s) of y
# t — array of t values at which y is computed
# returns y(t) values at the specific t points.
 
from scipy.integrate import odeint

ODE (easy)

# NUMERICAL SOLUTION
k = 0.3
 
def model(y, t):
    dydt = -k * y
    return dydt
 
y0 = 5.0
t  = np.linspace(0, 20, 50)
 
yns = odeint(model, y0, t)
 
plt.plot(t, yns, "rx")
 
# EXACT SOLUTION
yes = 5.0 * np.exp(-k * t)
plt.plot(t, yes, "b-", label="Exact solutions")

Finding the error

ydiff = yes - yns
plt.plot(t, ydiff, "k-")

Different k values

k = [0, 0.1, 0.2, 0.3, 0.4, 0.5]
 
for x in k:
    yns = odeint(model, y0, t, (x,))
    plt.plot(t, yns, "-", label=f"k = {x:.1f}")

ODE (more complicated)

def model(y, t):
    u = 0.0
    if (t >= 10.0):
        u = 2.0
    dydt = (u - y) / 5.0
    return dydt
 
y0 = 1.0
t  = np.linspace(0, 40, 100)
 
yns = odeint(model, y0, t)
 
plt.plot(t, yns, "b--", label="$y(t)$")

ODE — simultaneous differential equations

def model(z, t):
    x, y = z[0], z[1]
    dxdt = 3 * np.exp(-t)
    dydt = 3 - y
    dzdt = [dxdt, dydt]
    return dzdt
 
x0 = 0.0
y0 = 5.0
 
z0 = [x0, y0]
t  = np.linspace(0, 40, 50)
 
zns = odeint(model, z0, t)

ODE — simultaneous w/ piecewise

def model(z, t):
    x, y = z[0], z[1]
 
    if (t >= 5.0):
        u = 2.0
    else:
        u = 0.0
 
    dxdt = 0.5 * (-x*y + u)
    dydt = 0.2 * (-y + x)
    dzdt = [dxdt, dydt]
    return dzdt
 
 
z0 = [0.0, 1.0]
t  = np.linspace(0, 40, 100)
 
zns = odeint(model, z0, t)
 
plt.plot(t, zns[:, 0], "r--", label="$x(t)$")
plt.plot(t, zns[:, 1], "b--", label="$y(t)$")

Predator-Prey modelling (Lotka-Volterra)

Variables: x = rabbit pop, y = fox pop, dxdt/dydt = growth rates, a/b/c/d = rate parameters.

Assumptions: prey finds food always; predator food depends only on prey; rate proportional to size; no environmental change; limitless predator appetite.

print(f"Equilibrium population: prey = {c/d: .2f}, predator = {a/b:.2f}")

Lotka-Volterra

def LV(z, t):
    x, y = z[0], z[1]
    dxdt = x * (a - b*y)
    dydt = y * (d*x - c)
    dzdt = [dxdt, dydt]
    return dzdt
 
    # return [a*z[0] - b*z[0]*z[1], -c*z[1] + d*z[0]*z[1]]

tmax  = 10
ticks = 20 * tmax
ts    = np.linspace(0, tmax, ticks)
 
a, b, c, d = 0.4807, 0.0, 0.9272, 0.0   # b and d are 0 → no interaction
x0 = 30
y0 = 30
z0 = [x0, y0]
 
zs = odeint(LV, z0, ts)
 
plt.title("Lotka-Volterra equations,\n" +
          rf"$\alpha={a: .2f}, \beta={b}, \gamma={c: .2f}, \delta= {d: .2f}, (x_0,y_0)=({x0: .2f},{y0: .2f})$",
          fontsize=10)
plt.plot(ts, zs[:, 0], label='prey')
plt.plot(ts, zs[:, 1], label='predator')
 
plt.plot(zs[:, 0], zs[:, 1])

Lengthening integration tmax, adding interactions

tmax  = 25
ticks = 20 * tmax
ts    = np.linspace(0, tmax, ticks)
 
a, b, c, d = 0.480, 0.025, 0.930, 0.027
x0 = 20
y0 = 40
z0 = [x0, y0]
 
zs = odeint(LV, z0, ts)
 
plt.plot(ts, zs[:, 0], label='prey')
plt.plot(ts, zs[:, 1], label='predator')
plt.plot(zs[:, 0], zs[:, 1])

Logistic Lotka-Volterra

def Logistic_LV(z, t):
    return [a*z[0]*(1 - z[0]/K) - b*z[0]*z[1], -c*z[1] + d*z[0]*z[1]]
 
tmax  = 300
ticks = 20 * tmax
ts    = np.linspace(0, tmax, ticks)
 
a, b, c, d = 0.480, 0.025, 0.930, 0.027
K = 50
 
x0 = 20
y0 = 10
z0 = [x0, y0]
 
zs = odeint(Logistic_LV, z0, ts)
 
prey     = zs[:, 0]
predator = zs[:, 1]
 
plt.plot(ts, prey, label='prey')
plt.plot(ts, predator, label='predator')
 
plt.plot(prey, predator)
plt.scatter([c/d], [a/b * (1 - c/(K*d))], color='red', s=40)

See next

• Python-NumPy — arrays underlying every plot
• Python-Statistics — basic stats outside Matplotlib