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

Question

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

Answer 1

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

И в статье THIS вы можете найти описание алгоритмов.

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

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

Answer 2

Я занимаюсь физикой акселератора (источники синхротронного излучения), а при моделировании электронов, движущихся через магнитные поля, мы регулярно используем симплектические интеграторы. Наша основная рабочая лошадка - симплектический интегратор 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/ Интеграторов написаны на языке С с MEX ссылкой на Matlab. Поскольку электроны релятивистски, кинетический и потенциальный члены выглядят несколько иначе, чем в нерелятивистском случае, но все же можно записать гамильтониан как H = H1 + H2, где H1 - дрейф, а H2 - удар (скажем, из квадрупольного магниты или другие магнитные поля).

Answer 3

Вот исходный код метода композиции четвертого порядка, основанный на схеме Верле. Линейная регрессия $ \ log_ {10} (\ Delta t) $ против $ \ log_ {10} (Ошибка) $ покажет наклон 4, как и ожидалось из теории (см. График ниже). Более подробную информацию можно найти here, here или here.

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