【问题标题】:3D Phase portrait of SIR model in PythonPython中SIR模型的3D相图
【发布时间】:2021-02-11 22:39:12
【问题描述】:

如何使用 Python 创建动态系统的 3D 图? 我在这个question 中找到了我的问题的部分解决方案。但是添加维度似乎很糟糕..我的代码

import numpy as np
import matplotlib.pyplot as plt
X, Y, Z = np.meshgrid(np.linspace(-3.0, 3.0, 30), np.linspace(-3.0, 3.0, 30), np.linspace(-3.0, 3.0, 30))
u, v, w = np.zeros_like(X), np.zeros_like(X), np.zeros_like(X)
NI, NJ, NK = X.shape
for i in range(NI):
for j in range(NJ):
 for k in range(NK):        
  x, y, z = X[i,j,k], Y[i,j,k], Z[i,j,k]
  #LOTKA
       #u[i,j] =  x-x*y
       #v[i,j] = -y+x*y
       #SIR
       u[i,j,k] =  -x*y
       v[i,j,k] =  x*y-z  
       w[i,j,k] =  z         

plt.streamplot(X, Y, Z u, v, w)
plt.axis('square')#cube???
plt.axis([-3, 3, -3, 3,-3,3])
plt.show()

【问题讨论】:

    标签: python numpy math physics ode


    【解决方案1】:

    这绝不是您问题的答案,而只是一个建议:由于查看 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 以及第三种观点执行此操作,并对所有变量对的动态有一个不错的了解。

    【讨论】:

    • 我知道 N=S+I+R 的守恒。感谢您指出了这一点。顺便说一句,系统还有另一个哈密顿:H=R+(a/r)ln(S)
    • @yngabl 我不会称它为 Hamiltonain。这是第一个积分,但那是因为该系统是完全可解的。首先,因为x+y+z=c,然后因为最后一个方程dz/dt = z 与其他两个解耦,并且对于某个常数a 有解z = a*exp(t)。之后,您可以取上面的积分并表示 y = c-x-z = c-x-a*exp(t) 并将其代入第一个方程,这成为完全可解的方程。之后,第二个方程,加上解 x=x(t)z=a*exp(t),变成了 y 的线性方程,因此也可以求解。
    • 它被称为双哈密顿。存在两个辛矩阵,允许从这些守恒量写下 eom..
    猜你喜欢
    • 2021-07-22
    • 2021-04-08
    • 2020-05-09
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2020-05-28
    • 2020-07-27
    相关资源
    最近更新 更多