draw phase portrait of non-linear ODE system

Question

I'm kinda new to python and I need to draw the phase portrait of this system of equations:

x˙ = x(3 − 2x − y)
y˙ = y(2 − x − y)

I solved the system with the odeint function and tried to plot the phase portrait using the 4 fixed points, but end up with a different plot. Here's the code I wrote so far:

# function for the ODE
def f(u, t):
    x, y = u
    dxdt = x * (3 - 2*x - y)
    dydt = y * (2 - x - y)
    return np.array([dxdt, dydt])

#those are my fixed points + 0.01
values = np.array([[0.01,0.01], [0.01,2.01], [1.51, 0.01], [1.01,1.01])

t = np.linspace(0, 15, 1000)

for v in values:
    X0 = [v[0], v[1]]
    X = odeint(f, X0, t)
    plt.plot(X[:,0], X[:,1], color="black")

plt.ylabel("Sheeps")
plt.xlabel("Rabbits")
plt.legend(loc="upper right")
plt.show()

I should get something like this:

but I'm not. Can anyone help?

I count 8 solution curves in the picture. Note that each plot is always only a half-solution. Also explore the `streamplot` command of `pyplot`, see https://stackoverflow.com/questions/50891484/matplotlib-streamplot-with-streamlines-that-dont-break-or-end, https://stackoverflow.com/questions/58701195/plotting-phase-portraits-in-python-using-polar-coordinates and numerous examples on [math.SE](https://math.stackexchange.com/search?q=+plt.streamplot) — Lutz Lehmann, Jun 07 '23 at 15:27

score 1 · Accepted Answer · answered Jun 07 '23 at 15:52

As @Lutz Lehmann says in their comment, the phase portrait you showed has more than 4 starting points (it actually has 10; there are lines on the x and y axes). When you create your plot, you're only exploring 4 possible trajectories. You don't see those 4 trajectories because some are stable fixed points, so the solution will stay there. What you want to do use is use streamplot to explore a larger number of starting points. Following this solution, you can create the streamplot accordingly.

import numpy as np
import matplotlib.pyplot as plt

plt.close("all")

fixed_points = np.array([[0.,0.], [0.,2.], [1.5, 0.], [1.,1.]])

x = np.linspace(0, 3, 10)
y = np.linspace(0, 3, 10)
X, Y = np.meshgrid(x, y)

Xdot = X * (3 - 2*X - Y)
Ydot = Y * (2 - X - Y)

fig, ax = plt.subplots()
ax.streamplot(X, Y, Xdot, Ydot)
ax.scatter(*fixed_points.T, color="r")
fig.show()

As an aside, when solving initial value problems using scipy, you should favor using solve_ivp over the outdated odeint.

Great. If this solution works, you can accept this as the answer. — jared, Jun 07 '23 at 21:48

draw phase portrait of non-linear ODE system

1 Answers1