帮助符号积分器

Question

我正在尝试开发物理模拟，我想实施四阶 符号整合 方法。问题在于，我必须弄错数学，因为使用Symblectic Integrator（与四阶Runge-Kutta Integrator相比，我的仿真都根本无法正常工作，该集成符可用于模拟相当合理）。我一直在谷歌搜索这个问题，我能找到的只是有关该主题的科学文章。我试图调整文章中使用的方法，但没有运气。我想知道是否有人使用使用符号积分器的仿真源代码，最好是模拟引力场，但是任何符号集成符都可以使用。来源所使用的语言并不重要，但是我会喜欢使用C风格语法的语言。谢谢！

Answer 1

正如您要求提供的源代码：从 这里 您可以为哈密顿系统的符号方法下载MATLAB和FORTRAN代码，以及可逆性问题的对称方法。以及许多其他用于处理DIFF方程的方法。

并在 这 纸您可以找到算法的描述。

编辑

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

Answer 2

这是基于Verlet方案的第四阶组成方法的源代码。 $ log_ {10}（ delta t）$ vs. $ 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

我处于加速器物理（同步子光源）的领域，在对电子通过磁场进行建模时，我们会定期使用符号积分器。我们的基本主力是第四阶符合性集成器。如上所述，不幸的是，这些代码不是很好的标准化或Easilly。

一个基于MATLAB的开源跟踪代码称为Accelerator工具箱。有一个名为AtCollab的SourceForge项目。在这里看到一个混乱的维基https://sourceforge.net/apps/mediawiki/atcollab/index.php?title=main_page

对于集成商，您可以在这里查看：https://sourceforge.net/p/atcollab/code-0/235/tree/trunk/atintegrators/集成器用C与MATLAB的MEX链接编写。因为电子是相对论的，所以动力学和潜在的术语看起来与非权利主义情况有点不同，但是人们仍然可以写入哈密顿量为h = h1+h2，其中h1是h1是漂移，而h2是脚踢或其他磁场）。