这绝不是您问题的答案,而只是一个建议:由于查看 3D 相位肖像也可能有点令人困惑,您是否可以简单地将肖像的 2D 缩小应用到不变的 2D 平面上并查看那?我问这个,因为一般来说,一些 SIR = Susceptible-Infected-Recovered 流行病学模型,就像你提出的那样,具有总人口数是固定的属性,即在你的系统中 x + y + z = c 一些常数 @987654326 @。这是一种低保护类型(保护总人口)。它允许您表达z = c - x - y 并减少您的原始方程组。这是原始系统:
Sucseptible: x, Infected: y, Recovered: z
(the number of susceptible individuals x decreases at a rate which is equal to a portion of all current susceptible individuals x, where this portion is determined by the number of infected individuals y: )
dx/dt = - x*y
(the number of infected individuals y increases at a rate which is equal to the rate at which the susceptible individuals convert into infected, i.e. x*y from the first equation, but decreases by a the number z of formerly infected individuals that have just recovered: )
dy/dt = x*y - z
(the rate of change with which infected individuals recover)
dz/dt = z
注意
dx/dt + dy/dt + dz/dt = - x*y + x*y - z + z = 0
在对t 积分后,得出人口守恒定律:
x + y + z = c
表达z = c - x - y并代入二阶微分方程:
reduced system, showing the conversion dynamics of susceptible to infected individuals:
dx/dt = - x*y
dy/dt = x*y + x + y - c
您甚至可以从这里明确地解决系统问题。但是要快速查看动态:
import numpy as np
import matplotlib.pyplot as plt
X, Y = np.meshgrid(np.linspace(-3.0, 3.0, 30), np.linspace(-3.0, 3.0, 30))
# since x + y + z = c, which is an invariant 2D plane, set up an initial c:
c = 2
# this is simply matrix implementation of your code, since numpy is much faster at
# calculating matrices than running for loops:
u = -X*Y
v = X*Y + X + Y - c
plt.streamplot(X, Y, u, v)
plt.axis('square')#cube???
plt.axis([-3, 3, -3, 3])
plt.show()
并得到以下相图:
如果您想了解易感个体和康复个体之间的动态(或者在感染个体和康复个体之间),表达 y = c - x - z 并将其代入第一个等式并与第三个等式结合起来可能会更好一个。
dx/dt = -c*x + x*x + x*z
dz/dt = z
由于方程dz/dt = z与其他两个解耦,可以求解,然后其他两个方程也可以显式求解。但要从另一个角度快速了解动态:
import numpy as np
import matplotlib.pyplot as plt
X, Z = np.meshgrid(np.linspace(-3.0, 3.0, 30), np.linspace(-3.0, 3.0, 30))
c=2
# x + y + z = c
#y = c - x - z
U = -c*X + X*X + X*Z
V = -Z
plt.streamplot(X, Y, U, V)
plt.axis('square')#cube???
plt.axis([-3, 3, -3, 3])
plt.show()
您也可以对等式 2 和 3 以及第三种观点执行此操作,并对所有变量对的动态有一个不错的了解。