[問題] ODE 圖形換新的方程式

看板Python作者schmitt (小密特)時間7年前 (2018/07/01 01:34)推噓0(0推 0噓 0→)

留言0則, 0人參與討論串1/1

參考網路來源 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() -- ※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 111.250.133.168 ※ 文章網址: https://www.ptt.cc/bbs/Python/M.1530380051.A.9E6.html