一个常见的测试包括计算一些可以分析检查的整体平均值，并在减少时间步长的同时重复实验多次。您应该观察到，当您根据对数刻度的误差绘制时间步长长度时，可以轻松地将数据点拟合到斜率等于数值方案全局误差阶数的直线上（在Verlet 的情况）。

以下代码使用 Verlet 算法计算不同时间步长的 Lennard-Jones 振荡器的积分曲线，并收集平均总能量（应该是常数）。在下图中，一个简单的线性拟合显示斜率几乎为 2，正如理论所预期的那样。

#include <cmath>
#include <iostream>

using namespace std;

const double total_time = 1e3;

const double sigma = 1.0;
const double sigma6 = pow(sigma, 6.0);
const double epsilon = 1.0;
const double four_epsilon = 4.0 * epsilon;

inline double force(double q, double& potential) {
  const double r2 = q * q;
  const double r6 = r2 * r2 * r2;
  const double factor6  = sigma6 / r6;
  const double factor12 = factor6 * factor6;

  potential = four_epsilon * (factor12 - factor6);
  return -four_epsilon * (6.0 * factor6 - 12.0 * factor12) / r2 * q;
}

int main() {
  const double q0 = 1.5, p0 = 0.1;
  double potential;
  const double f0 = force(q0, potential);
  const double total_energy_exact = p0 * p0 / 2.0 + potential;

  for (double dt = 1e-3; dt <= 5e-2; dt *= 2.0) {
    const long steps = long(total_time / dt);

    double q = q0, p = p0, f = f0;
    double total_energy_average = total_energy_exact;

    for (long step = 1; step <= steps; ++step) {
      p += dt / 2.0 * f;
      q += dt * p;
      f = force(q, potential);
      p += dt / 2.0 * f;

      total_energy_average += p * p / 2.0 + potential;
    }

    total_energy_average /= double(steps);

    const double err = fabs(total_energy_exact - total_energy_average);
    cout << log10(dt) << "\t"
         << log10(err) << endl;
  }

  return 0;
}

一种简单的方法是认识到（比方说）势能在轨迹上的分布是先验已知的，并且以已知的方式随（比方）温度而变化。使用两个积分器观察到分布的相同（和预期）温度依赖性的统计显着观察是证明它们等效性的良好开端。请参阅http://dx.doi.org/10.1021/ct300688p（文章也可在 arxiv.org 上找到；代码在https://simtk.org/home/checkensemble上）。

另见分子动力学轨迹中的再现性讨论

这不是我的领域，但也许你可以用数字检查守恒定律，例如集合平均能量。