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

Question

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

Answer 1

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

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

Редактировать

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

Answer 2

Вот исходный код для метода композиции четвертого порядка на основе схемы врызги. Линейная регрессия $ log_ {10} ( delta t) $ vs. $ log_ {10} (ошибка) $ покажет наклон 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

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

Один с открытым исходным кодом для отслеживания 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 Link на MATLAB. Поскольку электроны являются релятивистскими кинетическими и потенциальными терминами выглядят немного отличаться, чем в нерелятивистском случае, но можно еще писать гамильтониана как H = H1 + H2, где H1 является дрейфом, а H2 - удар (скажем, что из квадрупольных магнитов или другие магнитные поля).