作者schmitt (小密特)
看板Python
标题[问题] debug
时间Fri May 25 21:37:32 2018
我参考网路上面的说明，希望可以画出图来。
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://webptt.com/cn.aspx?n=bbs/Python/M.1527255455.A.E95.html
1F：→ HybridSC: 这code不可能会爆...网格一开始就写死了 05/25 23:16
2F：→ HybridSC: 那三行的语法完全不对，return不是那样用的... 05/25 23:23
3F：→ HybridSC: 啊我看懂前两行了，硬要的话，放def下 05/25 23:26
4F：→ HybridSC: 第三行我就猜不透了，要看boundary condition怎麽设 05/25 23:29