访问 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`

对这个模块有一些经验，我不知道有什么方法可以在不深入研究内部结构的情况下获得步长。

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

Accessing Scipy ode solver's internal steps

python

scipy

ode

`scipy.integrate.odeint`

`scipy.integrate.ode`