Verhulst方程相图

Question

我正在尝试这本书的一个例子 - “使用 Python 的应用程序的动态系统”，我被要求绘制 Verhulst 方程的相图，然后我遇到了这个 post :

我得到的情节与上一个 post 的用户相同。每当我尝试使用公认的解决方案时，我都会收到“被零除”错误。为什么中接受的解决方案不起作用？

非常感谢您的帮助！

编辑：

使用之前 post 中的代码和@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()

我们有：

这仍然不同于：

我错过了什么？

Answer 1

在@Lutz Lehmann 的帮助下，我可以重写代码以满足我的需要。

解决方案是这样的：

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()

我们得到了我们想要的东西：

最后，再次感谢@Lutz Lehnmann 的帮助。没有他我也解决不了。

编辑 1：

我忘记了 $t = np.linspace(0, 10, 100)$ 并且如果你改变 figsize = (8,8)，我们会在图中得到一个更好的形状。（感谢@Trenton McKinney 的评论）

Verhulst方程相图

Phase portrait of Verhulst equation

python

matplotlib

scipy

ode

differential-equations