python - matplotlib绘制3 ODE系统的相图

预先计算给定平衡点的解为

sol=[ odeint(dynamical_model_bacteria,init,t) for init in x0];

然后，您可以在单独的循环中使用它们来构造两个图。使用子图可减少图像窗口的数量，然后将组合图给出

之后，使用相同的轴对象在其中绘制所有曲线，为一些随机初始点添加一些以填充空间

for k in range(len(x0)):
    x,y1,y2=sol[k].T
    ax.plot(x,y1,y2,'r', lw=2)

for k in range(80):
    sol = odeint(dynamical_model_bacteria,0.3*np.random.rand(3),t)
    x,y1,y2=sol.T
    ax.plot(x,y1,y2,'b', lw=1)

然后，保存的图是

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import solve_ivp, odeint
s=0  
omega=np.array([2,3.7,4,3])
omega11,omega12,omega21,omega22=omega


beta=np.array([28,24])
beta1,beta2=beta

epsilon=np.array([1.12,3.1])
epsilon1,epsilon2=epsilon

g12=2
gamma12=2
def dynamical_model_bacteria(X,t):
    x, y1, y2 = X
    dxdt=x*(1-x-y1-y2)+s
    dy1dt=y1*(-epsilon[0]*(1+omega11*y1+omega12*y2)-g12*y2+beta1*x)
    dy2dt=y2*(-epsilon[1]*(1+omega21*y1+omega22*y2)+gamma12*y1+beta2*x) 
    return [dxdt,dy1dt,dy2dt]

print( dynamical_model_bacteria([50,50,50],0))  
# compute some solutions in "orderly" fashion
x0=np.array([
       [ 0.11111111,  0.88888889,  0.        ],
       [ 0.37237237,  0.        ,  0.62762763],
       [ 0.        ,  0.10813086, -0.22171438],
       [ 0.17247589,  0.35219856,  0.47532555]])
t=np.linspace(0,30,8000)

sol=[ odeint(dynamical_model_bacteria,init,t) for init in x0];
# first plot
fig=plt.figure(figsize=(15,5))
for k in range(len(x0)):
    print( dynamical_model_bacteria(x0[k],0)) 
    plt.subplot(2,2,k+1);
    x,y1,y2=sol[k].T
    plt.xlabel('Time  (day)')
    plt.ylabel('populations (num/ml)')
    plt.title(r' $\quad\beta_{1}<\epsilon_{1} and\ \beta_{2}<\epsilon_{2}$'+'      '+'$x_{0}=$' +str(x0[k]))
    plt.plot(t,x,'r--' ,linewidth=2.5, markersize=1.5)
    plt.plot(t,y1,'b:', linewidth=3, markersize=0.5)
    plt.plot(t,y2,'g-.', linewidth=3.5, markersize=0.2)
    # set the limits
    #plt.xlim([0, 30])
    plt.ylim([-1, 1])
    plt.legend( ('Bacteria','protozoa','Daphnia'),
           loc='upper center', shadow=True)
plt.tight_layout(); 
plt.savefig('/tmp/EP2_fig1.png', dpi=300, bbox_inches='tight')
#second plot, adding some random solutions with IC in the same size range
fig=plt.figure(figsize=(15,5))
ax = plt.axes(projection="3d")
for k in range(len(x0)):
    x,y1,y2=sol[k].T
    ax.plot(x,y1,y2,'r', lw=2)

for k in range(80):
    sol = odeint(dynamical_model_bacteria,0.3*np.random.rand(3),t)
    x,y1,y2=sol.T
    ax.plot(x,y1,y2,'b', lw=1)

plt.savefig('/tmp/EP2_fig2.png', dpi=300, bbox_inches='tight')
plt.show(); plt.close()