我所知道的关于实现的唯一文档ode23t是在文件中，该文件记录了 ode23tb 的​​实现，即 MATLAB 中的 TRBDF2 方法。

正如通常实现的那样，梯形规则并不是严格意义上的一步法，因为截断误差估计利用了两个先前计算的解值。该方法对于非常僵硬的问题效率不高，因为它不是非常稳定的。

它引用了SPICE 书。这里有两件事。首先，梯形规则，也称为 Crank-Nicholson，不是 L 稳定的，通常不适合高刚度方程（但在某些模型中这种缺乏阻尼可能会产生不良结果）。误差估计的不稳定性加剧了这种情况，误差估计来自截断误差估计。基本上，您保存一些先前的点并使用导数近似来匹配领先的截断误差项，然后将其用作标准 P 或 PI 自适应算法的输入。

对于该算法的开源（MIT 许可）实现，您可以查看Trapezoid来自DifferentialEquations.jl的方法。该实现可在此处的源代码中找到，它源自 Shampine 提供的 SPICE 参考书目中的描述。请注意，由于许可问题，尚未直接针对 MATLAB 源代码进行检查，仅直接针对 Shampine 提供的这些源代码进行检查。

回应编辑

继续对话在 SO 中有点不受欢迎，但我会这样做，因为它是有道理的。对于第一部分，链接方法不是第一步中的隐式欧拉。相反，适应性在第一步中不存在，因为它使用两个步骤来获得二阶导数近似。但是，作为 1 步法，即使采用非自适应步长，梯形规则仍会使用。这使它始终保持第二顺序（尽管您可以有一个第一顺序步骤并且仍然是第二顺序）。

但是...如上所述，MatLab 通过对三次插值求微分来估计误差。

哪里说 MATLAB 是这样做的？MATLAB 对 使用三次插值误差dde23，但我从未见过它声明它在ode23t. 在评论中，我展示了从 MATLAB 的文档到 Shampine 的论文和 SPICE 书的论文轨迹。在 SPICE 书中，他们推导出截断误差估计，该估计使用两个步骤来估计作为主要截断误差项的三阶导数项。据我所知，由于我是一名开源开发人员并且不希望出现许可问题，因此没有查看 MATLAB 源代码，这就是根据给定源代码完成的方式。如果有人有时间，他们可以通过查看 MATLAB 源代码来加倍更改并确认或否认这一点。

要回答如何使用插值法得出误差估计的问题，让我向您展示 Hermite 在这里不起作用。我只是指向实时代码，因为我知道它是正确的，而且它比 TeXing 更容易。Hermite 插值在此处定义：

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/v4.16.4/src/dense/generic_dense.jl#L389

( ,k[1]) 和 ( ,k[2]) 分别作为 t 和 t+dt 处的值，在点 t+ dt 处计算。为了得到插值的导数，你用微分并得到$y_0$$y_1$$\theta$$\theta$

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/v4.16.4/src/dense/generic_dense.jl#L418

请注意，当时，由三次 Hermite 多项式近似的导数为 ，这是从本身处的导数估计。因此，它在区间的端点处没有针对导数的自然误差估计，因为它可以完美地达到这些导数值。当使用这个估计时，它被称为残差校正方法，在这里讨论，其中这个导数估计用于在区间中间通过使用特定点，您可以估计最大差异，这在 RK4 上是如何做到的$\theta=1$$k[2]$$t+dt$$f$$f$ddesd作品。RK4 很简单，它的两点最大残差估计的开源实现可以在这里找到。但是，虽然 MATLAB 确实在不同的积分器中使用了它，但 Shampine 写了它，从技术上讲，你可以将它应用于梯形方法，据我所知，这不是在ode23t.