使用Scipy ode求解器内部步骤的访问

Question

我目前正在从MATLAB切换到Python进行涉及求解微分方程的项目。

在MATLAB中，如果传递的t数组只包含两个元素，则求解器会输出模拟的所有中间步骤。然而，在Python中，你只会得到起点和终点。要获取中间的时间点，你必须明确地指定所需的时间点。

from scipy import integrate as sp_int
import numpy as np

def odeFun(t,y):
    k = np.ones((2))
    dy_dt = np.zeros(y.shape)
    dy_dt[0]= k[1]*y[1]-k[0]*y[0]
    dy_dt[1]=-dy_dt[0]
    return(dy_dt)

t = np.linspace(0,10,1000)
yOut = sp_int.odeint(odeFun,[1,0],t)

我也研究了以下方法：

solver = sp_int.ode(odefun).set_integrator('vode', method='bdf')
solver.set_initial_value([1,0],0)
dt = 0.01
solver.integrate(solver.t+dt)

然而，它仍然需要一个显式的dt。从阅读周围的内容中，我了解到Python的求解器（例如'vode'）计算所请求的dt的中间步骤，然后插值该时间点并输出它。不过，我想直接获得所有这些中间步骤，而不需要插值。这是因为它们代表了在积分容限内完全描述时间序列所需的最小点数。

是否有可用的选项来实现这一点？

我正在使用Python 3。

- sleepyzelaot

---

“然后插值该时间点并输出”这句话你确定吗？听起来有点不寻常。正常的做法是缩短第一步，以避免超过所需点，使其达到请求的点。-- 如果您真的很绝望想获取时间步长，可以在被积函数中加入仪器，并记录调用它的t。 - Paul Panzer

@PaulPanzer：在某些方法中，插值是相当常见的，其中可以假定插值足够准确，并且缩短时间步骤将是相对繁琐的。scipy.integrate.odeint进行插值，可以从tcur数据（请参见我的答案）以及您建议的日志记录中看到 - 两者都不包含初始时间步长。 - Wrzlprmft

@PaulPanzer：如果你真的很绝望，你总是可以对被积函数进行仪器测量，并记录它被调用时的t值。 - 这不会给你想要的结果，因为许多求解器还需要在步骤之间的时间评估导数（odeFun），例如，显式中点法就是一个简单的例子。 - Wrzlprmft

这是因为它们代表了在积分公差内完全描述时间序列所需的最小点数。- 你能详细说明为什么需要这些信息吗？我感觉到了一个 XY 问题。 - Wrzlprmft

@Wrzlprmft 哇，这么多错误的陈述在一个评论中，即使按照我的标准也是如此。感谢您纠正它们。 - Paul Panzer

Answer 1

scipy.integrate.odeint

odeint有一个选项full_output，允许您获取包含有关积分的信息的字典，包括tcur，它是：

每个时间步骤达到的t值的向量。(始终至少与输入时间一样大)。

注意第二句话：实际步骤始终与所需输出一样精细。如果要使用必要的最小步数，则必须要求进行粗略采样。

现在，这并不会给您提供值，但我们可以使用这些步骤再次进行积分来获得这些值：

from scipy.integrate import odeint
import numpy as np

def f(y,t):
    return np.array([y[1]-y[0],y[0]-y[1]])

start,end = 0,10 # time range we want to integrate
y0 = [1,0]       # initial conditions

# Function to add the initial time and the target time if needed:
def ensure_start_and_end(times):
    times = np.insert(times,0,start)
    if times[-1] < end:
        times = np.append(times,end)
    return times

# First run to establish the steps
first_times = np.linspace(start,end,100)
first_run   = odeint(f,y0,first_times,full_output=True)
first_steps = np.unique(first_run[1]["tcur"])

# Second run to obtain the results at the steps
second_times = ensure_start_and_end(first_steps)
second_run   = odeint(f,y0,second_times,full_output=True,h0=second_times[0])
second_steps = np.unique(second_run[1]["tcur"])

# ensuring that the second run actually uses (almost) the same steps.
np.testing.assert_allclose(first_steps,second_steps,rtol=1e-5)

# Your desired output
actual_steps = np.vstack((second_times, second_run[0].T)).T

scipy.integrate.ode

我对这个模块有一些经验，但并不知道如何在不深入研究内部的情况下获取步长。

- Wrzlprmft