每次更新ODE求解器中的初始条件 [英] Update initial condition in ODE solver each time step

问题描述

我想解决一个ODE系统，在最初的30,000秒内，我希望我的一个状态变量从相同的初始值开始.在这30,000秒之后，我想将该状态变量的初始值更改为其他值，并在其余时间内模拟系统.这是我的代码:

I am wanting to solve a system of ODEs where for the first 30,000 seconds, I want one of my state variables to start from the same initial value. After those 30,000 seconds, I want to change the initial value of that state variable to something different and simulate the system for the rest of time. Here is my code:

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)

因此，换句话说，我希望在前30,000秒内每次调用odeint时都将y0 [1]设置为相同的值(100).我正在有效地尝试让系统平衡一段时间，然后再将信号输入到系统中.我考虑做if t < 30000: y0[1] = 100之类的事情作为我的ode_rhs()函数的第一行，但是我不确定这是否可行.

So, in other words, I want to set y0[1] to the same value (100) each time odeint is called for the first 30,000 seconds. I'm effectively trying to let the system equilibrate for an amount of time before inputing a signal into the system. I thought about doing something like if t < 30000: y0[1] = 100 as the first line of my ode_rhs() function, but I'm not quite sure that works.

推荐答案

听起来您希望y1(t)保持常数(值为100)， 平衡阶段.您可以通过在此期间确保dy1(t)/dt = 0来做到这一点 阶段.有(至少)两种方法可以实现这一目标.首先是 修改ode_rhs中ydot[1]的计算，如下所示:

It sounds like you want y1(t) to remain constant (with the value 100) for the equilibration phase. You can do this by ensuring that dy1(t)/dt = 0 during this phase. There are (at least) two ways you can accomplish that. The first is to modify the calculation of ydot[1] in ode_rhs as follows:

if t < 30000:
    ydot[1] = 0.0
else:
    ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]

，并将初始条件100用于y[1].

and use the intitial condition 100 for y[1].

请注意，这会在系统的右侧引入不连续性， 但是odeint(Fortran代码LSODA)使用的自适应求解器通常足够健壮以应对它.

Note that this introduces a discontinuity in the right-hand side of the system, but the adaptive solver used by odeint (the Fortran code LSODA) is usually robust enough to handle it.

这是一个独立的例子.我已经将p和t1参数设置为ode_rhs. t1是平衡阶段的持续时间.

Here's a self-contained example. I've made p and t1 arguments to ode_rhs. t1 is the duration of the equilibration phase.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def ode_rhs(y, t, p, t1):
    ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
    if t < t1:
        ydot[1] = 0.0
    else:
        ydot[1] = -p[0]*y[0]*y[1] + p[1]*y[2]
    ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
    return ydot


ydot = np.zeros(3)
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t = np.linspace(0, 200, 2001)
t1 = 20.0

sol = odeint(ode_rhs, y0, t, args=(p, t1))


plt.figure(1)
plt.clf()

plt.subplot(3, 1, 1)
plt.plot(t, sol[:, 0])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[0]')


plt.subplot(3, 1, 2)
plt.plot(t, sol[:, 1])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)

plt.subplot(3, 1, 3)
plt.plot(t, sol[:, 2])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')

plt.show()

上述方法的一个细微变化是通过添加一个 参数为0或1.当参数为0时，将求解平衡系统；当参数为1时，将求解整个系统.这样，ydot[1](在我的较小示例中)的代码为

A slight variation of the above method is to modify the system by adding a parameter that is either 0 or 1. When the parameter is 0, the equlibration system is solved, and when the parameter is 1, the full system is solved. The code for ydot[1] (in my smaller example) is then

ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])

其中full是参数.

为处理平衡阶段，系统在0< = t< t< 0时被求解一次.与 full=0.然后将平衡溶液的最终值用作 第二个解决方案的初始条件，使用full=1运行.这种方法的优点是您不必强迫求解器处理不连续性.

To handle the equilibration phase, the system is solved once on 0 <= t < t1 with full=0. Then the final value of the equilibration solution is used as the initial condition to the second solution, run with full=1. The advantage of this method is that you are not forcing the solver to deal with the discontinuity.

这就是它在代码中的样子.

Here's how it looks in code.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def ode_rhs(y, t, p, full):
    ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
    ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
    ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
    return ydot


ydot = np.zeros(3)
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t1 = 20.0  # Equilibration time
tf = 200.0  # Final time

# Solve the equilibration phase.
teq = np.linspace(0, t1, 100)
full = 0
soleq = odeint(ode_rhs, y0, teq, args=(p, full))

# Solve the full system, using the final point of the
# equilibration phase as the initial condition.
y0 = soleq[-1]
# Note: the system is autonomous, so we could just as well start
# at t0=0.  But starting at t1 makes the plots (below) align without
# any additional shifting of the time arrays.
t = np.linspace(t1, tf, 2000)
full = 1
sol = odeint(ode_rhs, y0, t, args=(p, full))

plt.figure(2)
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(teq, soleq[:, 0], t, sol[:, 0])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[0]')

plt.subplot(3, 1, 2)
plt.plot(teq, soleq[:, 1], t, sol[:, 1])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)

plt.subplot(3, 1, 3)
plt.plot(teq, soleq[:, 2], t, sol[:, 2])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')

plt.show()

这是它生成的图(第一个示例中的图是 几乎完全相同):

And here's the plot that it generates (the plot from the first example is almost exactly the same):

这篇关于每次更新ODE求解器中的初始条件的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！