Построение фазового портрета для системы из 3 ODE по matplotlib

Question

Я пытаюсь построить временной ряд и фазовый портрет для динамической системы 3 ODE. Сначала я строю временные ряды для 4 начальных условий, но также хочу фазовый портрет со стрелками. Я получаю 4 3D-графики, которые не то, что я хочу. Я хочу что-то похожее на последних страницах этой статьи Стойкость и стабильность системы хищников с двумя жертвами

Код, который я использую:

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


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

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

g12=2
gamma12=22
def dynamical_model_bacteria(X,t):

    x=X[0]
    y1=X[1]
    y2=X[2]




    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))  
x0=np.array([
       [ 0.11111111,  0.88888889,  0.        ],
       [ 0.37237237,  0.        ,  0.62762763],
       [ 0.        ,  0.10813086, -0.22171438],
       [ 0.17247589,  0.35219856,  0.47532555]])
for i in x0:
    t=np.linspace(0,30,30000)
    x_sol=odeint(dynamical_model_bacteria,i,t)
    x=x_sol[:,0]
    y1=x_sol[:,1]
    y2=x_sol[:,2]
    fig=plt.figure(figsize=(15,5))
    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(i))
    plt.plot(t,x_sol[:,0],'r--' ,linewidth=2.5, markersize=1.5)
    plt.plot(t,x_sol[:,1],'b:', linewidth=3, markersize=0.5)
    plt.plot(t,x_sol[:,2],'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)
    fig=plt.figure(figsize=(15,5))
    ax = plt.axes(projection="3d")

    ax.plot3D(x,y1,y2)

    plt.savefig('EP2_fig2.png', dpi=300, bbox_inches='tight')

    plt.show()

Answer 1

Предварительно вычисляя решения для заданных точек равновесия как

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

, вы можете затем использовать их в отдельных циклах для построения двух диаграмм. Использование подзаголовка уменьшает количество изображений windows, тогда собранные графики дают

После этого используйте один и тот же объект осей, чтобы нарисовать все кривые в это, добавить некоторые для некоторых случайных начальных точек, чтобы заполнить пространство

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