如何监控SciPy.odeint进程？

Question

SciPy可以用scipy.integrate.odeint或其他软件包求解常微分方程，但要在函数完全解完后才能给出结果。但是，如果ode函数非常复杂，程序将需要很长时间(一到两天)才能给出完整的结果。那么，我如何控制它求解方程的步骤(当方程还没有完全求解时，打印出结果)？

Answer 1

当我在谷歌上搜索答案时，我找不到一个满意的答案。因此，我使用tqdm项目通过一个概念验证解决方案创建了一个简单的gist。希望这对你有帮助。

编辑:版主要求我对上面的链接给出一个解释。

首先，我使用的是scipy的OOP版本的odeint (solve_ivp)，但您可以将其改编回odeint。假设你想要从time T0集成到T1，并且你想显示每0.1%的进度。您可以修改ode函数以获取两个额外的参数，一个pbar (进度条)和一个state (集成的当前状态)。如下所示：

def fun(t, y, omega, pbar, state):
    # state is a list containing last updated time t:
    # state = [last_t, dt]
    # I used a list because its values can be carried between function
    # calls throughout the ODE integration
    last_t, dt = state
    
    # let's subdivide t_span into 1000 parts
    # call update(n) here where n = (t - last_t) / dt
    time.sleep(0.1)
    n = int((t - last_t)/dt)
    pbar.update(n)
    
    # we need this to take into account that n is a rounded number.
    state[0] = last_t + dt * n
    
    # YOUR CODE HERE
    dydt = 1j * y * omega
    return dydt

这是必要的，因为函数本身必须知道它所在的位置，但是scipy的odeint并没有真正给函数提供这个上下文。然后，您可以使用以下代码集成fun：

T0 = 0
T1 = 1
t_span = (T0, T1)
omega = 20
y0 = np.array([1], dtype=np.complex)
t_eval = np.arange(*t_span, 0.25/omega)

with tqdm(total=1000, unit="‰") as pbar:
    sol = solve_ivp(
        fun,
        t_span,
        y0,
        t_eval=t_eval,
        args=[omega, pbar, [T0, (T1-T0)/1000]],
    )

请注意，args中任何可变的(如列表)都是实例化一次的，并且可以在函数中进行更改。我建议这样做，而不是使用全局变量。

这将显示一个进度条，如下所示：

100%|█████████▉| 999/1000 [00:13<00:00, 71.69‰/s]

Answer 2

您可以拆分集成域并集成这些段，将前一段的最后一个值作为下一段的初始条件。在两者之间，打印出您想要的任何内容。如有必要，请使用numpy.concatenate将这些部件组装起来。

在三体太阳系模拟的标准示例中，替换代码

u0 = solsys.getState0();
t = np.arange(0, 100*365.242*day, 0.5*day);
%timeit u_res = odeint(lambda u,t: solsys.getDerivs(u), u0, t, atol = 1e11*1e-8, rtol = 1e-9)

输出：1 loop, best of 3: 5.53 s per loop

使用进度报告代码

def progressive(t,N):
    nk = [ int(n+0.5) for n in np.linspace(0,len(t),N+1) ]
    u0 = solsys.getState0();
    u_seg = [ np.array([u0]) ];
    for k in range(N):
        u_seg.append( odeint(lambda u,t: solsys.getDerivs(u), u0, t[nk[k]:nk[k+1]], atol = 1e11*1e-8, rtol = 1e-9)[1:] )
        print t[nk[k]]/day
        for b in solsys.bodies: print("%10s %s"%(b.name,b.x))
    return np.concatenate(u_seg)
%timeit u_res = progressive(t,20)

输出：1 loop, best of 3: 5.96 s per loop

显示控制台打印的开销仅为8%。有了更实质性的ODE功能，报告开销的一小部分将显著减少。

也就是说，python，至少在它的标准包中，不是用于工业规模的数字处理的工具。始终使用具有强类型变量的编译版本，以尽可能减少解释开销。

• 使用一些经过大量开发和测试的包，如Sundials或julia-lang框架differentialequations.jl，直接用适当的编译语言编写ODE函数。对于较大的步长，使用高阶方法，因此步长较小。测试使用隐式或指数/Rosenbrock方法是否进一步减少了每个固定间隔的步数或ODE函数评估。不同之处可能是加速比的10到100倍。
• 使用上述的python包装器和你的ODE的一些加速友好的实现function.
• Use准源代码翻译工具JITcode将你的python ODE函数翻译成一个意大利面条状的C指令列表，然后给出一个编译后的函数，可以(几乎)直接从编译的

内核中调用。