Я пытаюсь разработать физическое моделирование и хочу реализовать метод симплектического интегрирования четвертого порядка. . Проблема в том, что я, должно быть, ошибаюсь в математике, поскольку моя симуляция вообще не работает при использовании симплектического интегратора (по сравнению с интегратором Рунге-Кутты четвертого порядка, который достаточно хорошо работает для симуляции). Я гуглил это вечно, и все, что я могу найти, это научные статьи по этому вопросу. Я пытался адаптировать метод, использованный в статьях, но мне не повезло. Я хочу знать, есть ли у кого-нибудь исходный код для моделирования, в котором используются симплектические интеграторы, предпочтительно для моделирования гравитационного поля, но подойдет любой симплектический интегратор. На каком языке находится источник, не имеет большого значения, но я был бы признателен за язык, использующий синтаксис в стиле C. Спасибо!
Помощь с симплектическими интеграторами
Ответы (3)
Поскольку вы просили исходный код: из ЗДЕСЬ вы можете загрузить код MATLAB и FORTRAN для симплектические методы для гамильтоновых систем и симметричные методы для обратимых задач. А также множество других методов работы с дифференциальными уравнениями.
А в ЭТОМ документе вы можете найти описание алгоритмов.
Изменить
Если вы используете Mathematica, это тоже может помочь.
Вот исходный код метода композиции четвертого порядка, основанного на схеме Верле. Линейная регрессия $\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);
}
Я занимаюсь физикой ускорителей (синхротронные источники света), и при моделировании электронов, движущихся через магнитные поля, мы регулярно используем симплектические интеграторы. Наша основная рабочая лошадка — симплектический интегратор 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 — толчок (скажем, от квадрупольных магнитов). или другие магнитные поля).