如何监视SciPy.odeint的过程？

Question

SciPy可以通过scipy.integrate.odeint或其他软件包来求解ode方程，但是在函数完全求解后才能给出结果。但是，如果ode函数非常复杂，则该程序将花费大量时间（一到两天）才能得出整个结果。因此，如何控制方程的求解步骤（当方程尚未完全求解时打印结果）？

Answer 1

您可以拆分集成域并集成分段，将上一个的最后一个值作为下一个分段的初始条件。在其间，打印任何所需的内容。如有必要，请使用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框架differentequations.jl，以适当的编译语言直接编码ODE函数。对于较大的步长（因此较小的步长），请使用高阶方法。测试使用隐式或指数/ Rosenbrock方法是否进一步减少了每个固定间隔的步数或ODE函数评估。差异可以是加速的10到100倍。

• 使用上面的python包装器和ODE函数的某些加速友好实现。

[

Answer 2