參考網路來源
http://kitchingroup.cheme.cmu.edu/blog/2013/02/21/Phase-portraits-of-a-system-of-ODEs/
我想要把裡面的方程式，換成以下論文裡面的方程式
http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0174946
所以我目前的進度是這樣:
import numpy as np
import matplotlib.pyplot as plt
mu = 0.16
nu= 0.04
D_H = 0.3
D_A=0.02
D_S=0.06
RHO_A = 0.03
RHO_H = 0.0001
epsilon=1.0
gamma=0.02
C=0.002
D=0.008
E=0.1
F=10
C_O=0.02
def laplacian(Z):
Ztop = Z[0:-2, 1:-1]
Zleft = Z[1:-1, 0:-2]
Zbottom = Z[2:, 1:-1]
Zright = Z[1:-1, 2:]
Zcenter = Z[1:-1, 1:-1]
return (Ztop + Zleft + Zbottom + Zright -
4 * Zcenter) / dx**2
deltaA = laplacian(Z)
deltaS = laplacian(Z)
def f(Y, t):
y1, y2,y3,y4 = Y
return [C*A*A* S /H - mu * A + RHO_A * Y+D_A*deltaA,
C * A*S - nu*H+ +RHO_H * Y+D_H*deltaH,
C_O - gamma * S - epsilon * Y * S+D_S*deltaS,
D * A - E * Y + Y * Y / ( 1 + F * Y * Y)
]
y3 = np.linspace(0, 1.0, 20)
y4 = np.linspace(0, 1.0, 20)
Y1, Y2 = np.meshgrid(y3, y4)
t = 0
u, v = np.zeros(Y1.shape), np.zeros(Y2.shape)
NI, NJ = Y1.shape
for i in range(NI):
for j in range(NJ):
x = Y1[i, j]
y = Y2[i, j]
yprime = f([x, y], t)
u[i,j] = yprime[0]
v[i,j] = yprime[1]
v[i,j] = yprime[2]
v[i,j] = yprime[3]
Q = plt.quiver(Y1, Y2, u, v, color='r')
plt.xlabel('$y_1$')
plt.ylabel('$y_2$')
plt.xlim([-2, 8])
plt.ylim([-4, 4])
plt.savefig('e:/phase-portrait.png')
from scipy.integrate import odeint
for y20 in [0, 0.5]:
tspan = np.linspace(0, 10, 200)
y0 = [0.0, y20]
ys = odeint(f, y0, tspan)
plt.plot(ys[:,0], ys[:,1], 'b-') # path
plt.plot([ys[0,0]], [ys[0,1]], 'o') # #start
plt.plot([ys[-1,0]], [ys[-1,1]], 's') # #end
plt.xlim([0, 1])
plt.savefig('e:/phase-portrait-2.png')
plt.show()