Scipysolve_ivp：在运行过程中可以更改max_step吗？

Question

我正在求解本质上落在弹性表面上的弹性体的运动方程。通过反复试验，我发现，只有当我使用solve_ivp选择0.0001左右的max_step时，我才能得到合理的结果。如果没有这个，即使像 Radau 这样通常好的方法似乎也无法给出结果。 我相信，只有当我的任何身体靠近表面时才需要这样做，并且这些距离在积分过程中是已知的。 我的问题：我可以在集成过程中更改 max_step 的值吗？

感谢您的帮助！

Answer 1

solve_ivp

无法处理动态变化的步长，但

solve_ivp

是执行

ODESolver

类型求解器的包装器。使用

ODESolver

（例如

Radau

），可以更改

max_step

属性。请参阅 scipy 文档，尽管关于如何使用较低级别的基于类的求解器的示例很少。希望下面的例子可以帮助到您。 例如这段代码限制步长不大于0.1。

from scipy.integrate import solve_ivp, Radau
import numpy as np
def fun(t, y):
    return [1]
        
sol = solve_ivp(fun, [0, 1], [0], method="Radau", max_step=0.1)
print(np.diff(sol.t))

产量

[1.00e-04 1.00e-03 1.00e-02 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 8.89e-02]

以下代码是等效的：

rad = Radau(fun, t0=0, y0=[0], t_bound=1, max_step=0.1)
step_size = []
while rad.status == "running":
    rad.step()
    step_size.append(rad.step_size)
print(np.array(step_size))

产生

[1.00e-04 1.00e-03 1.00e-02 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 8.89e-02]

，这是一样的。

但是现在，我们可以在计算过程中修改

ODESolver

对象来更改计算过程中的最大步长。

rad = Radau(fun, t0=0, y0=[0], t_bound=1, max_step=0.1)
step_size = []
while rad.t < 0.5 and rad.status == "running":
    rad.step()
    step_size.append(rad.step_size)
rad.max_step = 0.05  # modify max_step
while rad.status == "running":
    rad.step()
    step_size.append(rad.step_size)
print(np.array(step_size))

产量

[1.00e-04 1.00e-03 1.00e-02 1.00e-01 1.00e-01 1.00e-01 1.00e-01 1.00e-01 5.00e-02 5.00e-02 5.00e-02 5.00e-02 5.00e-02 5.00e-02 5.00e-02 5.00e-02 5.00e-02 3.89e-02]