Установите x
как интеграл v
, чтобы x'=v
, x''=v'
, аналогично x1
для v1
, чтобы ваше уравнение читалось как дифференциальное уравнение второго порядка
a*x''+(b+k1)*x'+c*x=k1*v1+k2*x1
Учитывая v1
в качестве входных данных, вектор состояния должен содержать три интегрированные переменные x, v, x1
, которые дают функцию ODE
def odesys(y,t):
x, v, x1 = y
v1 = eval_v1(t)
return [ v, (k1*v1+k2*x1 - (b+k1)*v-c*x )/a, v1 ]
Чтобы использовать его с odeint, вы можете, например, сделать
t = np.linspace(0,T,2001); # define the end time T before
y0 = [ 0, 0, 0 ] # standard convention is that everything is zero for negative times
y = odeint(odesys, y0, t) # add arguments for higher accuracy if needed
x, v, x1 = y.T # transpose of a list of tuples is a tuple of lists
plt.plot(t,x); plt.show() # as example that should work