scipy.integrate.odeint 中的时间相关输入?
Time dependent input in scipy.integrate.odeint?
我想使用 odeint 将方波作为输入积分到耦合微分方程。下面的代码似乎可以正常工作,但是会引发错误:
"IndexError: index out of bounds" 第 12 行。
谁能帮我解决这个问题?
提前致谢,
约蒂卡
1 import pylab as pl
2 import numpy as np
3 from scipy.integrate import odeint
4
5
6 def func1(Rates,time,ip):
7 # 0 = x, 1 = y
8 # dx/dt = -x+ ax+by, dy/dt =-y+ cx+dy
9 leak = -0.3
10 a = d = 0.1
11 b = c = -0.2
12 dxbydt = leak*Rates[0]+ a*Rates[0]+b/2.*Rates[1]+ip[time]
13 dybydt = leak*Rates[1]+ c*Rates[0]+d*Rates[1]+ip[time]
14
15 return [dxbydt, dybydt]
16 time = np.arange(0,1000,1)
17 ip = np.zeros((len(time)))
18
19 ip[300:600] = 5.0
20 initR = np.ones((2))*10
21 fR = odeint(func1,initR,time,args=(ip,))
22
23 pl.figure()
24 #pl.plot(time,ip,'k-',label='ip')
25 pl.plot(time,fR[:,0],'b-',label='x')
26 pl.plot(time,fR[:,1],'r-',label='y')
27 pl.legend()
28
29 pl.show()
问题的关键是ip[time]
,因为time
的值在函数内部不断变化,被反复调用。 time
增加,然后达到ip
的最大指标值即1000
。
会有问题,因为时间是一个浮点数,但它被用作 ip
的索引。
这里 time
的目的是什么?
import numpy as np
from scipy.integrate import odeint
def getIp(time):
if time > 300 and time < 600:
x = 5
else:
x = 0
return x
def func1(Rates,time):
# 0 = x, 1 = y
# dx/dt = -x+ ax+by, dy/dt =-y+ cx+dy
leak = -0.3
a = d = 0.1
b = c = -0.2
dxbydt = leak*Rates[0]+ a*Rates[0]+b/2.*Rates[1]+getIp(time)
dybydt = leak*Rates[1]+ c*Rates[0]+d*Rates[1]+getIp(time)
return [dxbydt, dybydt]
timeGrid = np.arange(0,1000,0.01)
ip = np.zeros((len(timeGrid)))
#ip[300:600] = 5.0
initR = np.ones((2))*10
fR = odeint(func1,initR,timeGrid)
pl.figure()
pl.plot(timeGrid,ip,'k-',label='ip')
pl.plot(timeGrid,fR[:,0],'b-',label='x')
pl.plot(timeGrid,fR[:,1],'r-',label='y')
pl.legend()
pl.show()
我想使用 odeint 将方波作为输入积分到耦合微分方程。下面的代码似乎可以正常工作,但是会引发错误: "IndexError: index out of bounds" 第 12 行。
谁能帮我解决这个问题?
提前致谢,
约蒂卡
1 import pylab as pl
2 import numpy as np
3 from scipy.integrate import odeint
4
5
6 def func1(Rates,time,ip):
7 # 0 = x, 1 = y
8 # dx/dt = -x+ ax+by, dy/dt =-y+ cx+dy
9 leak = -0.3
10 a = d = 0.1
11 b = c = -0.2
12 dxbydt = leak*Rates[0]+ a*Rates[0]+b/2.*Rates[1]+ip[time]
13 dybydt = leak*Rates[1]+ c*Rates[0]+d*Rates[1]+ip[time]
14
15 return [dxbydt, dybydt]
16 time = np.arange(0,1000,1)
17 ip = np.zeros((len(time)))
18
19 ip[300:600] = 5.0
20 initR = np.ones((2))*10
21 fR = odeint(func1,initR,time,args=(ip,))
22
23 pl.figure()
24 #pl.plot(time,ip,'k-',label='ip')
25 pl.plot(time,fR[:,0],'b-',label='x')
26 pl.plot(time,fR[:,1],'r-',label='y')
27 pl.legend()
28
29 pl.show()
问题的关键是ip[time]
,因为time
的值在函数内部不断变化,被反复调用。 time
增加,然后达到ip
的最大指标值即1000
。
会有问题,因为时间是一个浮点数,但它被用作 ip
的索引。
这里 time
的目的是什么?
import numpy as np
from scipy.integrate import odeint
def getIp(time):
if time > 300 and time < 600:
x = 5
else:
x = 0
return x
def func1(Rates,time):
# 0 = x, 1 = y
# dx/dt = -x+ ax+by, dy/dt =-y+ cx+dy
leak = -0.3
a = d = 0.1
b = c = -0.2
dxbydt = leak*Rates[0]+ a*Rates[0]+b/2.*Rates[1]+getIp(time)
dybydt = leak*Rates[1]+ c*Rates[0]+d*Rates[1]+getIp(time)
return [dxbydt, dybydt]
timeGrid = np.arange(0,1000,0.01)
ip = np.zeros((len(timeGrid)))
#ip[300:600] = 5.0
initR = np.ones((2))*10
fR = odeint(func1,initR,timeGrid)
pl.figure()
pl.plot(timeGrid,ip,'k-',label='ip')
pl.plot(timeGrid,fR[:,0],'b-',label='x')
pl.plot(timeGrid,fR[:,1],'r-',label='y')
pl.legend()
pl.show()