Помощь с симплектическими интеграторами

Я пытаюсь разработать физическое моделирование, и я хочу реализовать метод symplectic integration четвертого порядка. Проблема в том, что я должен ошибаться в математике, поскольку моя симуляция вообще не работает при использовании симплектического интегратора (по сравнению с интегратором Runge-Kutta четвертого порядка, который достаточно хорошо работает для моделирования). Я вел это навсегда, и все, что я могу найти, это научные статьи по этому вопросу. Я попытался адаптировать метод, используемый в статьях, но мне не повезло. Я хочу знать, есть ли у кого-нибудь исходный код для моделирования, который использует симплектические интеграторы, предпочтительно для имитации гравитационного поля, но любой симплектический интегратор будет делать. На каком языке находится источник, это не имеет большого значения, но я бы оценил язык, использующий синтаксис C-стиля. Спасибо!

Ответ 1

Как вы просили исходный код: Из ЗДЕСЬ вы можете скачать код MATLAB и FORTRAN для симплектических методов для гамильтоновых систем и симметричные методы для обратимых проблемы. И множество других методов для решения уравнений диффузии тоже.

И в ЭТО в документе вы можете найти описание алгоритмов.

Edit

Если вы используете Mathematica this, это тоже может помочь.

Ответ 2

Вот исходный код для метода композиции четвертого порядка, основанный на схеме Верле. Линейная регрессия $\ 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);
}

Ответ 3

Я занимаюсь физикой ускорителя (источники синхротронного излучения), а при моделировании электронов, движущихся через магнитные поля, мы регулярно используем симплектические интеграторы. Наша основная рабочая лошадка - симплектический интегратор 4-го порядка. Как отмечалось в вышеприведенных комментариях, к сожалению, эти коды не так хорошо стандартизированы или доступны в свободном доступе.

Один открытый код отслеживания на основе 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 - удар (скажем, от квадрупольных магнитов или других магнитных полей).