Aiuto con integratori simplettici

Question

Sto provando a sviluppare una simulazione fisica e voglio implementare un metodo symplectic integration di quarto ordine. Il problema è che devo essere in errore, dal momento che la mia simulazione non funziona affatto quando si utilizza l'integratore simplettico (rispetto a un integratore Runge-Kutta del quarto ordine che funziona abbastanza bene per la simulazione). Ho cercato su Google per sempre e tutto quello che posso trovare sono articoli scientifici sull'argomento. Ho provato ad adattare il metodo usato negli articoli, ma non ho fortuna. Voglio sapere se qualcuno ha il codice sorgente per una simulazione che fa uso di integratori simplettici, preferibilmente per simulare un campo gravitazionale, ma qualsiasi integratore simplettico farebbe. La lingua in cui si trova la sorgente non ha importanza, ma apprezzerei un linguaggio usando la sintassi in stile C. Grazie!

Answer 1

Come richiesto per il codice sorgente: da HERE è possibile scaricare codice MATLAB e FORTRAN per metodi simplettici per sistemi Hamiltoniani e metodi simmetrici per problemi reversibili. E molti altri metodi per affrontare anche le equazioni di diff.

E nella carta THIS è possibile trovare la descrizione degli algoritmi.

Modifica

Se si usa Mathematica this può aiutare anche.

Answer 2

Sono nel campo della fisica degli acceleratori (sorgenti di luce di sincrotrone) e, nel modellare gli elettroni che si muovono attraverso i campi magnetici, usiamo integratori simplettici su base regolare. Il nostro cavallo di battaglia di base è un integratore simplettico del 4 ° ordine. Come notato nei commenti sopra, sfortunatamente questi codici non sono così standardizzati o disponibili facilmente.

Un codice di tracciamento basato su Matlab open source è denominato Accelerator Toolbox. C'è un progetto Sourceforge chiamato atcollab. Vedere una wiki disordinato qui https://sourceforge.net/apps/mediawiki/atcollab/index.php?title=Main_Page

Per gli integratori, si può guardare qui: https://sourceforge.net/p/atcollab/code-0/235/tree/trunk/atintegrators/ Gli integratori sono scritti in C con MEX link per Matlab. Poiché gli elettroni sono relativistici, i termini cinetici e potenziali sembrano leggermente diversi rispetto al caso non relativistico, ma si può ancora scrivere l'Hamiltoniano come H = H1 + H2 dove H1 è una deriva e H2 è un calcio (diciamo dal quadrupolo magneti o altri campi magnetici).

Answer 3

Ecco il codice sorgente per un metodo di composizione del quarto ordine basato sullo schema di Verlet. Una regressione lineare di $ \ log_ {10} (\ Delta t) $ contro $ \ log_ {10} (Errore) $ mostrerà una pendenza di 4, come previsto dalla teoria (vedere il grafico sottostante). Ulteriori informazioni possono essere trovate here, here o 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); 
}