python - 访问 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。

Answer 1

scipy.integrate.odeint

odeint 有一个选项 full_output，它允许您获取包含集成信息的字典，包括 tcur，它是:

vector with the value of t reached for each time step. (will always be at least as large as the input times).

(注意第二句:实际的步骤总是和你想要的输出一样精细。如果你想使用最少的必要步骤，你必须要求粗采样。)

现在，这不会为您提供值，但我们可以通过使用以下步骤进行第二次积分来获得这些值:

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

有一些使用此模块的经验，我不知道有什么方法可以在不深入挖掘内部结构的情况下获得步长。