Plotting the phase portrait for a system of 3 ODE by matplotlib

Question

I am trying to plot the time series and phase portrait for a dynamical system of 3 ODE. I first plot the time series for 4 initial condition but want the phase portrait with arrows as well. I am getting 4 3D plots which not what I want. I want something like in the last pages of this paper Persistence and stability of a two prey one predator system

The code I am using is:

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

Pre-computing the solutions for the given equilibrium points as

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

you can then use them in separate loops to construct the two diagrams. Using subplot reduces the number of image windows, the assembled graphs then give

After that use the same axes object to draw all the curves in it, add some for some random initial points to fill the space

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)

The saved plot then is

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