dopri5积分器的步长控制

Question

我正在尝试使用scipy.integrate.ode 中的dopri5 积分器解决一个简单的示例。正如文档所述

这是由于 Dormand & Prince 的 (4)5 阶显式 runge-kutta 方法（具有步长控制和密集输出）。

这应该可以。所以这是我的例子：

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def MassSpring_with_force(t, state):
    """ Simple 1DOF dynamics model: m ddx(t) + k x(t) = f(t)"""
    # unpack the state vector
    x = state[0]
    xd = state[1]

    # these are our constants
    k = 2.5 # Newtons per metre
    m = 1.5 # Kilograms

    # force
    f = force(t)

    # compute acceleration xdd
    xdd = ( ( -k*x + f) / m )

    # return the two state derivatives
    return [xd, xdd]


def force(t):
    """ Excitation force """
    f0 = 1  # force amplitude [N]
    freq = 20  # frequency[Hz]
    omega = 2 * np.pi *freq  # angular frequency [rad/s]
    return f0 * np.sin(omega*t)


# Time range
t_start = 0
t_final = 1

# Main program
state_ode_f = ode(MassSpring_with_force)

state_ode_f.set_integrator('dopri5', rtol=1e-6, nsteps=500,
                       first_step=1e-6, max_step=1e-3)

state2 = [0.0, 0.0]  # initial conditions
state_ode_f.set_initial_value(state2, 0)

sol = np.array([[t_start, state2[0], state2[1]]], dtype=float)

print("Time\t\t Timestep\t dx\t\t ddx\t\t state_ode_f.successful()")

while state_ode_f.t < (t_final):
    state_ode_f.integrate(t_final, step=True)
    sol = np.append(sol, [[state_ode_f.t, state_ode_f.y[0], state_ode_f.y[1]]], axis=0)
    print("{0:0.8f}\t {1:0.4e} \t{2:10.3e}\t {3:0.3e}\t {4}".format(
            state_ode_f.t, sol[-1, 0]- sol[-2, 0], state_ode_f.y[0], state_ode_f.y[1], state_ode_f.successful()))

我得到的结果是：

Time         Timestep    dx      ddx         state_ode_f.successful()
0.49763822   4.9764e-01      2.475e-03   -8.258e-04  False
0.99863822   5.0100e-01      3.955e-03   -3.754e-03  False
1.00000000   1.3618e-03      3.950e-03   -3.840e-03  False

带有警告：

c:\python34\lib\site-packages\scipy\integrate_ode.py:1018: UserWarning: dopri5: 需要更大的 nmax self.messages.get(idid, 'Unexpected idid=%s' % idid))

结果不正确。如果我使用vode 积分器运行相同的代码，我会得到预期的结果。

编辑

这里描述了一个类似的问题： Using adaptive step sizes with scipy.integrate.ode

建议的解决方案建议设置nsteps=1，它可以正确求解 ODE，并具有步长控制。但是，集成器将 state_ode_f.successful() 返回为 False。

Answer 1

不，没有错。您告诉集成商对t_final 执行集成步骤，它会执行该步骤。未报告集成商的内部步骤。

---

明智的做法是将所需的采样点作为算法的输入，例如设置dt=0.1并使用

state_ode_f.integrate( min(state_ode_f.t+dt, t_final) )

dopri5 中没有单一的step 方法，只有vode 定义了它，参见源代码https://github.com/scipy/scipy/blob/v0.14.0/scipy/integrate/_ode.py#L376，这可以解释观察到的差异。

正如您在Using adaptive step sizes with scipy.integrate.ode 中发现的那样，可以通过设置迭代界限nsteps=1 来强制执行单步行为。这每次都会产生警告，因此必须抑制这些特定警告才能看到合理的结果。

---

您不应将参数（积分间隔的常数）用于时间相关力。在MassSpring_with_force 内部使用评估f=force(t)。可能您可以将force 的函数句柄作为参数传递。