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

Please don't ask for a complete solution. Work though your code on your own. Than ask for a specific problem, if you stuck somewhere. — T. Christiansen, Feb 07 '20 at 12:36
@ I do have my attempt to plot in the code but want one phase portrait for all initial values not 4 as I am getting. — , Feb 07 '20 at 12:55
You should add that the initial points are all equilibria or close to it, only the third is sufficiently unstable to give a non-constant solution, in relation to that the other 3 solutions appear as points in the 3D plot. You should probably add some other initial points to the collection to fill the box. — Lutz Lehmann, Feb 07 '20 at 15:00

Lutz Lehmann · Accepted Answer · 2020-02-07T21:05:20.240

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

I have been trying to use your solution.The two plots are what I want but still I don't know how to use the code provided in the one I posted. My coding skill is limited so it will be helpful to post the full code that gives the two plots. What is the red line in the box by the way? — , Feb 07 '20 at 20:59
The red line (color/line style parameter`'r'`) is the only trace of the first loop over your given initial conditions, from the third one. The other points give constant solutions that are dimensionless in the 3D plot, thus not plotted or erased by the blue (`'b'`) solutions to the random initial conditions. — Lutz Lehmann, Feb 07 '20 at 20:59
sorry for being too demanding . Do you know how to put comma in the initial condition in the title to distinguish the coordinates? And what about plotting the directions arrow . I know that I should use quiver for that but don't know how? — , Feb 07 '20 at 21:32
I'm not aware of any arrows that could be systematically added to a curve. `streamplot` does this somehow, but does not make the method generally available. // `numpy` prints arrays as matrices. The normal python list is printed with commas, `str(list(x0[k]))` works, but gives the fully identifying format with 17 digits precision. You could also build your own such as `"["+", ".join(f"{x:.6g}" for x in x0[k])+"]"`. — Lutz Lehmann, Feb 07 '20 at 21:41
using str(list(x0[k])) with round to not get all digits works. Thanks — , Feb 07 '20 at 21:58
Unfortunately no, not that I'm aware of it. I have not seen the plt.streamplot command used in 3D, this would be the most likely candidate for a systematic solution. — Lutz Lehmann, Apr 05 '21 at 15:49

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

1 Answers1