Phase portrait of Verhulst equation

Question

I was trying to an example of the book -"Dynamical Systems with Applications using Python" and I was asked to plot the phase portrait of Verhulst equation, then I came across this post: How to plot a phase portrait of Verhulst equation with SciPy (or SymPy) and Matplotlib?

I'm getting the same plot as the user on the previous post. Whenever, I try to use the accepted solution I get a "division by zero" error. Why doesn't the accepted solution in How to plot a phase portrait of Verhulst equation with SciPy (or SymPy) and Matplotlib? works?

Thank you very much for you help!

Edit:

Using the code from the previous post and the correction given by @Lutz Lehmann

beta, delta, gamma = 1, 2, 1
b,d,c = 1,2,1

C1 = gamma*c-delta*d
C2 = gamma*b-beta*d
C3 = beta*c-delta*b

def verhulst(X, t=0):
    return np.array([beta*X[0] - delta*X[0]**2 -gamma*X[0]*X[1],
                     b*X[1] - d*X[1]**2 -c*X[0]*X[1]])

X_O = np.array([0., 0.])
X_R = np.array([C2/C1, C3/C1])
X_P = np.array([0, b/d])
X_Q = np.array([beta/delta, 0])

def jacobian(X, t=0):
    return np.array([[beta-delta*2*X[0]-gamma*X[1],  -gamma*x[0]],
                     [b-d*2*X[1]-c*X[0],             -c*X[1]]])

values  = np.linspace(0.3, 0.9, 5)                         
vcolors = plt.cm.autumn_r(np.linspace(0.3, 1., len(values)))

f2 = plt.figure(figsize=(4,4))

for v, col in zip(values, vcolors):
    X0 = v * X_R
    X = odeint(verhulst, X0, t)
    plt.plot(X[:,0], X[:,1], color=col, label='X0=(%.f, %.f)' % ( X0[0], X0[1]) )

ymax = plt.ylim(ymin=0)[1] 
xmax = plt.xlim(xmin=0)[1]
nb_points = 20

x = np.linspace(0, xmax, nb_points)
y = np.linspace(0, ymax, nb_points)

X1, Y1  = np.meshgrid(x, y)
DX1, DY1 = verhulst([X1, Y1])  # compute growth rate on the gridt
M = (np.hypot(DX1, DY1))       # Norm of the growth rate
M[M == 0] = 1.                 # Avoid zero division errors
DX1 /= M                       # Normalize each arrows
DY1 /= M

plt.quiver(X1, Y1, DX1, DY1, M, cmap=plt.cm.jet)
plt.xlabel('Number of Species 1')
plt.ylabel('Number of Species 2')
plt.legend()
plt.grid()

We have:

That is still different from:

What am I missing?

It doesn't work, if you check the definition of c1 with the conditions on the accepted answer, we have zero. As R = [c2/c1, c3/c1], we get the error I was talking about. Believe me, I've checked — RFeynman, Nov 11 '21 at 17:10
The error is caused by a copy-paste error. The order in OP and my local code was and is `b,d,c = 1,2,1`, now corrected in the linked answer. With the correction, `C1 = 1*1-2*2=-3` is no longer zero. — Lutz Lehmann, Nov 11 '21 at 17:43
Thank you very much! I'll try now and I'll update you later! — RFeynman, Nov 11 '21 at 17:46
I ran the code again with that correction and doesn't give me the same plot. Do you get the same plot that is in the other post? I'll update my answer to show what I got. — RFeynman, Nov 11 '21 at 17:48
Yes, there were more modifications to get appropriate initial points transversal to the diagonal, `values = np.linspace(0.05, 0.15, 5)`, solutions then from `X0 = [v,0.2-v]` and then again from `X0=6*X0`. — Lutz Lehmann, Nov 11 '21 at 17:58
@LutzLehmann It might be best to update the other answer with the complete code that produces the plot. Best regards. — Trenton McKinney, Nov 11 '21 at 18:08
@TrentonMcKinney, I just added to this post, I don't know if you would like to link to the other post. — RFeynman, Nov 11 '21 at 18:10
Since there's a link to that question in your question, they will show on the page as linked — Trenton McKinney, Nov 11 '21 at 18:12

RFeynman · Accepted Answer · 2021-11-11T18:31:39.150

With the help of @Lutz Lehmann I could rewrite the code to get want I needed.

The solutions is something like this:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(8, 4), dpi= 80, facecolor='whitesmoke', edgecolor='k')
beta, delta, gamma = 1, 2, 1
b,d,c = 1,2,1

t = np.linspace(0, 10, 100)

C1 = gamma*c-delta*d
C2 = gamma*b-beta*d
C3 = beta*c-delta*b

def verhulst(X, t=0):
    return np.array([beta*X[0] - delta*X[0]**2 -gamma*X[0]*X[1],
                     b*X[1] - d*X[1]**2 -c*X[0]*X[1]])

X_O = np.array([0., 0.])
X_R = np.array([C2/C1, C3/C1])
X_P = np.array([0, b/d])
X_Q = np.array([beta/delta, 0])

def jacobian(X, t=0):
    return np.array([[beta-delta*2*X[0]-gamma*X[1],  -gamma*x[0]],
                     [b-d*2*X[1]-c*X[0],             -c*X[1]]])

values  = np.linspace(0.05, 0.15, 5)                      
vcolors = plt.cm.autumn_r(np.linspace(0.3, 1., len(values)))


for v, col in zip(values, vcolors):
    X0 = [v,0.2-v]
    X = odeint(verhulst, X0, t)
    plt.plot(X[:,0], X[:,1], color=col, label='X0=(%.f, %.f)' % ( X0[0], X0[1]) )

for v, col in zip(values, vcolors):
    X0 = [6 * v, 6 *(0.2-v)]
    X = odeint(verhulst, X0, t)
    plt.plot(X[:,0], X[:,1], color=col, label='X0=(%.f, %.f)' % ( X0[0], X0[1]) )

    
ymax = plt.ylim(ymin=0)[1] 
xmax = plt.xlim(xmin=0)[1]
nb_points = 20

x = np.linspace(0, xmax, nb_points)
y = np.linspace(0, ymax, nb_points)

X1, Y1  = np.meshgrid(x, y)
DX1, DY1 = verhulst([X1, Y1])  # compute growth rate on the gridt
M = (np.hypot(DX1, DY1))       # Norm of the growth rate
M[M == 0] = 1.                 # Avoid zero division errors
DX1 /= M                       # Normalize each arrows
DY1 /= M

plt.quiver(X1, Y1, DX1, DY1, M, cmap=plt.cm.jet)
plt.xlabel('Number of Species 1')
plt.ylabel('Number of Species 2')
plt.grid()

We get what we were looking for:

Finally, I would like to thank again @Lutz Lehnmann for the help. I wouldn't have solved without it him.

Edit 1:

I forgot $t = np.linspace(0, 10, 100)$ and if you change figsize = (8,8), we get a nicer shape in the plot. (Thank you @Trenton McKinney for the remarks)

Phase portrait of Verhulst equation

1 Answers1