Hilfe für symplektischer Integratoren

Question

Ich versuche, eine Physik-Simulation zu entwickeln, und ich mag eine vierte Ordnung symplektischer Integration implementieren Methode. Das Problem ist, dass ich muß die Mathematik falsch bekommen, da meine Simulation gar nicht funktioniert, wenn der symplektischer Integrator mit (im Vergleich zu einer vierten Ordnung Runge-Kutta Integrator, der einigermaßen gut für die Simulation funktioniert). Ich habe dies für immer war googeln und alles, was ich wissenschaftliche Artikel sind zu diesem Thema finden kann. Ich versuchte, das Verfahren in den Artikeln verwendet anzupassen, aber kein Glück habe. Ich möchte wissen, ob jemand Quellcode für eine Simulation hat das macht symplektischer Integratoren verwenden, vorzugsweise ein Gravitationsfeld zu simulieren, aber jeder symplektischen Integrator tun würde. Welche Sprache ist die Quelle in nicht allzu wichtig, aber ich würde eine Sprache mit C-Stil-Syntax zu schätzen wissen. Dank!

Answer 1

Wie Sie für den Quellcode gestellt: Von hier können Sie herunterladen MATLAB und FORTRAN Code für Methoden zur symplektischer Hamilton-Operator-Systemen und Methoden für symmetrische reversible Probleme. Und viele andere Methoden für zu mit Diff Gleichungen zu tun.

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

Bearbeiten

Wenn Sie Mathematica verwenden diese kann auch helfen.

Answer 2

Hier ist der Quellcode für eine vierte Methode, um Zusammensetzung auf dem Schema Verlet basiert. Eine lineare Regression von log_ $ \ {10} (\ Delta t) $ vs. $ \ log_ {10} (Error) zeigt $ eine Steigung von 4, wie von der Theorie zu erwarten (die Grafik unten). Weitere Informationen finden Sie hier werden, hier oder hier .

#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

Ich bin auf dem Gebiet der Beschleunigerphysik (Synchrotron-Lichtquellen) und Elektronen in der Modellierung durch Magnetfelder bewegen, verwenden wir symplektischer Integratoren auf einer regelmäßigen Basis. Unsere Basis Arbeitspferd ist ein vierte Ordnung symplektischer Integrator. Wie oben in den Kommentaren erwähnt, leider sind diese Codes nicht so gut standardisiert oder easilly zur Verfügung.

Ein Open-Source-Matlab-basierte Tracking-Code wird Accelerator Toolbox genannt. Es gibt ein Projekt namens Source 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 hier sehen: https://sourceforge.net/p/atcollab/code- 0/235 / Baum / trunk / atintegrators / Die Integratoren sind in C mit MEX Link in Matlab geschrieben. Da die Elektronen relativistische die kinetische und potentielle Begriffe sind ein wenig anders aussehen als in der nicht-relativistischen Fall, aber man kann immer noch den Hamilton-Operator als H = H1 + H2 schreiben, wo H1 eine Drift und H2 ist ein Kick (sagen wir von Quadrupol-Magnete oder andere Magnetfelder).