Atualizar a condição inicial no solucionador de ODEs a cada etapa de tempo
Eu estou querendo resolver um sistema de ODEs onde, nos primeiros 30.000 segundos, eu quero que uma das minhas variáveis de estado comece com o mesmo valor inicial. Depois desses 30.000 segundos, quero alterar o valor inicial dessa variável de estado para algo diferente e simular o sistema pelo resto do tempo. Aqui está o meu código:
def ode_rhs(y, t):
ydot[0] = -p[7]*y[0]*y[1] + p[8]*y[8] + p[9]*y[8]
ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]
ydot[2] = -p[10]*y[2]*y[3] + p[11]*y[9] + p[12]*y[9]
ydot[3] = -p[13]*y[3]*y[6] + p[14]*y[10] + p[15]*y[10] - p[10]*y[2]*y[3] + p[11]*y[9] + p[9]*y[8] - p[21]*y[3]
ydot[4] = -p[19]*y[4]*y[5] - p[16]*y[4]*y[5] + p[17]*y[11] - p[23]*y[4] + y[7]*p[20]
ydot[5] = -p[19]*y[4]*y[5] + p[15]*y[10] - p[16]*y[4]*y[5] + p[17]*y[11] + p[18]*y[11] + p[12]*y[9] - p[22]*y[5]
ydot[6] = -p[13]*y[3]*y[6] + p[14]*y[10] - p[22]*y[6] - p[25]*y[6] - p[23]*y[6]
ydot[7] = 0
ydot[8] = p[7]*y[0]*y[1] - p[8]*y[8] - p[9]*y[8]
ydot[9] = p[10]*y[2]*y[3] - p[11]*y[9] - p[12]*y[9] - p[21]*y[9]
ydot[10] = p[13]*y[3]*y[6] - p[14]*y[10] - p[15]*y[10] - p[22]*y[10] - p[21]*y[10] - p[23]*y[10]
ydot[11] = p[19]*y[4]*y[5] + p[16]*y[4]*y[5] - p[17]*y[11] - p[18]*y[11] - p[22]*y[11] - p[23]*y[11]
ydot[12] = p[22]*y[10] + p[22]*y[11] + p[22]*y[5] + p[22]*y[6] + p[21]*y[10] + p[21]*y[3] + p[21]*y[9] + p[24]*y[13] + p[25]*y[6] + p[23]*y[10] + p[23]*y[11] + p[23]*y[4] + p[23]*y[6]
ydot[13] = p[15]*y[10] + p[18]*y[11] - p[24]*y[13]
return ydot
pysb.bng.generate_equations(model)
alias_model_components()
p = np.array([k.value for k in model.parameters])
ydot = np.zeros(len(model.odes))
y0 = np.zeros(len(model.odes))
y0[0:7] = p[0:7]
t = np.linspace(0.0,1000000.0,100000)
r = odeint(ode_rhs,y0,t)
Então, em outras palavras, eu quero definir y0 [1] para o mesmo valor (100) de cada vezodeint
é chamado pelos primeiros 30.000 segundos. Estou efetivamente tentando deixar o sistema equilibrar por um período de tempo antes de inserir um sinal no sistema. Eu pensei em fazer algo comoif t < 30000: y0[1] = 100
como a primeira linha do meuode_rhs()
função, mas não tenho certeza de que funciona.