Hilfe mit symplektischen Integratoren

Question

Ich versuche, eine Physik-Simulation zu entwickeln, und ich möchte eine Methode der vierten Ordnung symplectic integration implementieren. Das Problem ist, dass ich die Mathematik falsch interpretieren muss, da meine Simulation bei Verwendung des symplektischen Integrators überhaupt nicht funktioniert (im Vergleich zu einem Runge-Kutta-Integrator vierter Ordnung, der für die Simulation einigermaßen gut funktioniert). Ich habe das für immer gegooglet und alles, was ich finden kann, sind wissenschaftliche Artikel zu diesem Thema. Ich habe versucht, die in den Artikeln verwendete Methode anzupassen, habe aber kein Glück. Ich möchte wissen, ob irgendjemand einen Quellcode für eine Simulation hat, die symplektische Integratoren verwendet, vorzugsweise um ein Gravitationsfeld zu simulieren, aber jeder symplektische Integrator würde dies tun. In welcher Sprache die Quelle ist, ist nicht so wichtig, aber ich würde eine Sprache schätzen, die C-artige Syntax verwendet. Vielen Dank!

Answer 1

Wie Sie nach Quellcode gefragt: Von HERE können Sie MATLAB und FORTRAN-Code für symplektische Methoden für Hamilton-Systeme und symmetrische Methoden für reversible Probleme herunterladen. Und viele andere Methoden, um auch mit Differentialgleichungen umzugehen.

Und in THIS Papier können Sie die Beschreibung der Algorithmen finden.

bearbeiten

Wenn Sie Mathematica verwenden this helfen.

Answer 2

Ich bin auf dem Gebiet der Beschleunigerphysik (Synchrotronlichtquellen), und bei der Modellierung von Elektronen, die sich durch Magnetfelder bewegen, verwenden wir regelmäßig symplektische Integratoren. Unser grundlegendes Arbeitspferd ist ein symplektischer Integrator 4. Ordnung. Wie in den obigen Bemerkungen erwähnt, sind diese Codes leider nicht so gut standardisiert oder einfach verfügbar.

Ein Open-Source-Matlab-basierter Tracking-Code heißt Accelerator Toolbox. Es gibt ein Sourceforge-Projekt namens atcollab. Hier finden Sie eine unordentliche Wiki hier https://sourceforge.net/apps/mediawiki/atcollab/index.php?title=Main_Page

Für die Integratoren, können Sie sich hier: https://sourceforge.net/p/atcollab/code-0/235/tree/trunk/atintegrators/ Die Integratoren sind geschrieben in C mit MEX-Link zu Matlab. Da die Elektronen relativistisch sind, sehen die kinetischen und potentiellen Terme ein wenig anders aus als im nicht-relativistischen Fall, aber man kann den Hamiltonian immer noch als H = H1 + H2 schreiben, wobei H1 eine Drift und H2 ein Kick ist (etwa vom Quadrupol) Magnete oder andere magnetische Felder).

Answer 3

Hier ist der Quellcode für eine Kompositionsmethode vierter Ordnung basierend auf dem Verlet-Schema. Eine lineare Regression von $ \ log_ {10} (\ Delta t) $ vs. $ \ log_ {10} (Fehler) $ zeigt eine Steigung von 4, wie von der Theorie erwartet (siehe Grafik unten). Weitere Informationen finden Sie here, here oder 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); 
}