[問題] debug

看板Python作者schmitt (小密特)時間7年前 (2018/05/25 21:37)推噓0(0推 0噓 4→)

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

我參考網路上面的說明，希望可以畫出圖來。 http://kitchingroup.cheme.cmu.edu/blog/2013/02/21/Phase-portraits-of-a-system-of-ODEs/ 但是因為需求，我把圖形的範圍還有原本的方程式做了修改。然後有同學建議可以在解ode時加這個條件，以免它噴掉。然後我就不知道怎麼把預防爆掉的指令放進去。請看下面被 #### 包圍的三行要寫到哪裡去? import numpy as np import matplotlib.pyplot as plt U=0.1 #def f(Y, t): # y1, y2 = Y # return [-y2*(1-y1*y2)+U*y1, y1*(1-y1)+U*y2] def f(Y, t): y1, y2 = Y return [-y2*(1-y1*y2)+U*y1, y1*(1-y1)+U*y2] #######這裡3行，我完全不知道要放哪裡 if abs(Y[0])>2 or abs(Y[1])>2: return[0,0] return[Y[1]-Y[0],Y[0]*Y[1]] ####################### y1 = np.linspace(-3.0, 3.0, 20) y2 = np.linspace(-3.0, 3.0, 20) Y1, Y2 = np.meshgrid(y1, y2) t = 0 u, v = np.zeros(Y1.shape), np.zeros(Y2.shape) NI, NJ = Y1.shape print(NI) 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] Q = plt.quiver(Y1, Y2, u, v, color='r') plt.xlabel('$y_1$') plt.ylabel('$y_2$') plt.xlim([-3, 3]) plt.ylim([-3, 3]) plt.savefig('e:/phase-portrait.png') from scipy.integrate import odeint for y20 in range(-2,2): for y21 in range(-2,2): tspan = np.linspace(0, 4, 100) y0 = [y21, y20] ys = odeint(f, y0, tspan)#, full_output=True) 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([-3, 3]) plt.savefig('e:/phase-portrait-2.png') plt.show() -- ※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 101.8.161.25 ※ 文章網址: https://www.ptt.cc/bbs/Python/M.1527255455.A.E95.html