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

Posted 2023-03-12

【中文标题】访问 Scipy ode 求解器的内部步骤

【英文标题】：Accessing Scipy ode solver's internal steps

【发布时间】：2017-09-11 23:42:10

【问题描述】：

我目前正在为一个涉及求解微分方程的项目从 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。

【问题讨论】：

"然后对那个时间点进行插值并输出" 你确定吗？对我来说听起来有点不寻常。正常的程序是缩短会超调的第一步，因此它会达到要求的点。 - 重新掌握时间步长。如果你真的很绝望，你可以随时检测被积函数并记录调用它的 t。

@PaulPanzer：对于某些可以假设插值足够准确并且缩短时间步长相对乏味的方法来说，插值是很常见的。 scipy.integrate.odeint 插值，从 tcur 数据（请参阅我的回答）以及您建议的日志中可以看出 - 两者都不包含初始时间步长。

@PaulPanzer：如果你真的很绝望，你可以随时检测被积函数并记录调用它的 t。 - 这不会给你想要的结果，因为许多求解器还需要在步骤之间有时评估导数 (odeFun)，例如，explicit midpoint method 举一个简单的例子。

这是因为它们代表了在积分容差范围内完整描述时间序列所需的最少点数。 – 您能否详细说明为什么需要此信息？我感觉到XY problem。

@Wrzlprmft 哇，即使按照我的标准，一条评论中也有很多错误的陈述。感谢您正确设置它们。

【参考方案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

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