Verhulst 方程的相图

Posted 2023-02-18

【中文标题】Verhulst 方程的相图

【英文标题】：Phase portrait of Verhulst equation

【发布时间】：2021-12-24 03:58:01

【问题描述】：

我试图举一本书的例子——“Dynamical Systems with Applications using Python”，我被要求绘制 ，然后我看到了这篇文章：How to plot a phase portrait of Verhulst equation with SciPy (or SymPy) and Matplotlib?

我得到的情节与上一篇文章中的用户相同。每当我尝试使用公认的解决方案时，我都会收到“除以零”错误。为什么How to plot a phase portrait of Verhulst equation with SciPy (or SymPy) and Matplotlib? 中接受的解决方案不起作用？

非常感谢您的帮助！

编辑：

使用上一篇文章中的代码和@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()

我们有：

这仍然不同于：

我错过了什么？

【问题讨论】：

该错误是由复制粘贴错误引起的。 OP 和我的本地代码中的顺序是 b,d,c = 1,2,1，现在在链接的答案中更正。通过更正，C1 = 1*1-2*2=-3 不再为零。

是的，为了获得与对角线交叉的适当初始点values = np.linspace(0.05, 0.15, 5)，进行了更多修改，然后从 X0 = [v,0.2-v] 和 X0=6*X0 获得解决方案。

@LutzLehmann 最好用生成绘图的完整代码更新其他答案。最好的问候。

@TrentonMcKinney，我刚刚添加到这篇文章，我不知道你是否想链接到其他帖子。

由于您的问题中有指向该问题的链接，因此它们将在页面上显示为已链接

【参考方案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 的发言）