我一直在努力创建一些自制的数值方法，并使用它们来可视化我的 Strogatz 动力学教科书中的教科书问题。这感觉像是同时学习数值方法和动力学的好方法。

书中最近的一个问题涉及“严重减速”，并通过比较来举例说明$\dot{x} = -x^3$和$\dot{x} = -x$. 我不能假装理解任何事情，除了一个比另一个衰减得更快。

我决定通过我自制的小型 RK4 算法运行这两个算法并绘制它们。我发现我会因为步长太大而溢出，我想知道为什么。从概念上讲，是什么导致了这种行为？

我知道应该根据解决方案变化的速度来选择步长，否则误差会非常高。我没有得到的是我得到的溢出错误。

代码（ODE 在底部定义和求解。要查看溢出，请将步长参数中的 .001 替换为 .1。：

def rk4(dt, t, field, y_0):

    y = np.asarray(len(t) * [y_0]) 

    for i in np.arange(len(t) - 1):
        k1 = dt * field(t[i], y[i])
        k2 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k1)
        k3 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k2)
        k4 = dt * field(t[i] + dt, y[i] + k3)
        y[i + 1] = y[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6

    return y

def vf_grapher(fn, t_0, t_n, dt, y_0, lintype='-r', sup_title=None,
               title=None, xlab=None, ylab=None):

    t = np.arange(t_0, t_n, dt)
    y_min = .0
    y_max = .0

    fig, axs = plt.subplots()
    fig.suptitle(sup_title)

    axs.set_title(title)
    axs.set_ylabel(ylab)
    axs.set_xlabel(xlab)

    for iv in np.array(y_0, ndmin=1, copy=False):
        soln = rk4(dt, t, fn, iv)
        plt.plot(t, soln, lintype)
        if y_min > np.min(soln):
            y_min = np.min(soln)
        if y_max < np.max(soln):
            y_max = np.max(soln)

    x = np.linspace(t_0, t_n + dt, 11)
    y = np.linspace(y_min, y_max, 11)

    X, Y = np.meshgrid(x, y)

    theta = np.arctan(fn(X, Y))

    U = np.cos(theta)
    V = np.sin(theta)

    plt.quiver(X, Y, U, V, angles='xy')
    plt.xlim((t_0, t_n - dt))
    plt.ylim((y_min - .05 * y_min, y_max + .05 * y_max))
    plt.show()


if __name__ == '__main__':

    def crit_slow(t, x):
        return -x**3


    def exp_decay(t, y):
        return -y

    # graph an ODE modelling critical slow down (has discontinuities in it's 
    second derivative)
    vf_grapher(crit_slow, 0., 10., .001, 10.)

    # graph an ODE modelling exponential decay (compare to critical slowdown)
    vf.vf_grapher(exp_decay, 0., 10., .001, 10.)