Plotting the phase portrait for a system of 3 ODE by matplotlib
Estou a tentar traçar a série temporal e o retrato de fase para um sistema dinâmico de 3 ODE. Eu primeiro ploto a série de tempo para 4 condição inicial, mas quero o retrato de fase com setas também. Estou recebendo 4 parcelas 3D que não é o que eu quero. Quero algo como nas últimas páginas deste jornal persistência e estabilidade de duas presas um sistema predador
o código que estou a usar é:
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()
0
Author: T. Christiansen, 2020-02-07
1 answers
Pré-calcular as soluções para os pontos de equilíbrio indicados como
sol=[ odeint(dynamical_model_bacteria,init,t) for init in x0];
Você pode então usá-los em loops separados para construir os dois diagramas. Usando subplot reduz o número de janelas de imagem, os grafos montados então dar
Depois disso, use o mesmo objecto de eixos para desenhar todas as curvas nele, adicione algumas para alguns pontos iniciais aleatórios para preencher o espaço
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)
Então, o enredo salvo é:
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()
0
Author: Lutz Lehmann, 2020-02-07 21:05:20