辛积分器帮助

Question

我正在尝试开发物理模拟，并且想要实现四阶 辛积分 方法。问题是我的数学一定是错误的，因为使用辛积分器时我的模拟根本不起作用（与模拟效果相当好的四阶龙格-库塔积分器相比）。我一直在谷歌上搜索这个问题，我能找到的都是关于这个主题的科学文章。我尝试采用文章中使用的方法，但没有运气。我想知道是否有人有使用辛积分器进行模拟的源代码，最好是模拟引力场，但任何辛积分器都可以。源代码使用什么语言并不重要，但我希望使用 C 风格语法的语言。谢谢！

Answer 1

正如您所要求的源代码：从 此处 您可以下载 MATLAB 和 FORTRAN 代码哈密​​顿系统的辛方法和可逆问题的对称方法。还有许多其他处理微分方程的方法。

在这篇论文中您可以找到算法的描述。

编辑

如果您使用 Mathematica 这也可能有帮助。

Answer 2

这是基于 Verlet 方案的四阶合成方法的源代码。
$\log_{10}(\Delta t)$ 与 $\log_{10}(Error)$ 的线性回归将显示斜率为 4，正如理论所预期的那样（参见下图）。
更多信息请访问此处、此处 或 此处。

#include <cmath>
#include <iostream>

using namespace std;

const double total_time = 5e3;

// Parameters for the potential.
const double sigma = 1.0;
const double sigma6 = pow(sigma, 6.0);
const double epsilon = 1.0;
const double four_epsilon = 4.0 * epsilon;

// Constants used in the composition method.
const double alpha = 1.0 / (2.0 - cbrt(2.0));
const double beta = 1.0 - 2.0 * alpha;


static double force(double q, double& potential);

static void verlet(double dt,
                   double& q, double& p,
                   double& force, double& potential);

static void composition_method(double dt,
                               double& q, double& p,
                               double& f, double& potential);


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-2; dt <= 5e-2; dt *= 1.125) {
    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) {
      composition_method(dt, q, p, f, potential);
      const double total_energy = p * p / 2.0 + potential;
      total_energy_average += total_energy;
    }

    total_energy_average /= double(steps);

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

  return 0;
}

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;
}

void verlet(double dt,
            double& q, double& p,
            double& f, double& potential) {
  p += dt / 2.0 * f;
  q += dt * p;
  f = force(q, potential);
  p += dt / 2.0 * f;
}

void composition_method(double dt,
                        double& q, double& p,
                        double& f, double& potential) {
  verlet(alpha * dt, q, p, f, potential);
  verlet(beta * dt, q, p, f, potential);
  verlet(alpha * dt, q, p, f, potential);
}

Answer 3

我从事加速器物理领域（同步加速器光源），在对穿过磁场的电子进行建模时，我们定期使用辛积分器。我们的基本主力是四阶辛积分器。正如上面的评论所指出的，不幸的是这些代码没有很好地标准化或容易获得。

一种基于 Matlab 的开源跟踪代码称为 Accelerator Toolbox。 Sourceforge 有一个名为 atcollab 的项目。在这里查看凌乱的维基
https://sourceforge.net/apps/mediawiki/atcollab/index.php ?title=Main_Page

对于集成商，您可以查看这里：
https://sourceforge.net/p/atcollab/code- 0/235/tree/trunk/atintegrators/
积分器是用 C 语言编写的，并通过 MEX 链接到 Matlab。
因为电子是相对论性的，所以动能项和势能项看起来与非相对论性情况略有不同，但仍然可以将哈密顿量写为 H=H1+H2，其中 H1 是漂移，H2 是冲击（例如来自四极磁铁）或其他磁场）。