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()