Aiuto con integratori simplettici

Question

Sto cercando di sviluppare una simulazione fisica e voglio implementare un quarto ordine integrazione simplettico metodo. Il problema è che devo essere sempre dalla parte del torto matematica, dal momento che la mia simulazione non funziona affatto quando si utilizza l'integratore simplettica (rispetto a un quarto ordine Runge-Kutta integratore che funziona ragionevolmente bene per la simulazione). Sono stato googling questo per sempre e tutto quello che posso trovare sono articoli scientifici sull'argomento. Ho cercato di adattare il metodo usato negli articoli, ma sto avendo senza fortuna. Vorrei 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. Che lingua la fonte è in non importa troppo, ma gradirei un linguaggio utilizzando la sintassi in stile C. Grazie!

Answer 1

Per quanto hai chiesto per il codice sorgente: Da qui che è possibile scaricare MATLAB e codice FORTRAN per metodi simplettici per sistemi Hamiltoniani e metodi simmetriche per problemi reversibili. E un sacco di altri metodi per trattare con diff equazioni troppo.

E in questo documento è possibile trovare la descrizione di gli algoritmi.

Modifica

Se si usa Mathematica questo può aiutare anche.

Answer 2

Ecco il codice sorgente per un quarto metodo di composizione ordine basato sullo schema Verlet. Una regressione lineare di $ \ log_ {10} (\ Delta t) $ contro $ \ log_ {10} (errore) $ mostrerà una pendenza di 4, come previsto dalla teoria (si veda il grafico sottostante). Maggiori informazioni possono essere trovate qui , qui o qui .

#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

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

Un open source Matlab codice di monitoraggio basato si chiama 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 / albero / trunk / atintegrators / Gli integratori sono scritti in C con MEX link per Matlab. Poiché gli elettroni sono relativistica i termini cinetica e potenziale aspetto un po 'diverso rispetto al caso non relativistico, ma si può ancora scrivere l'Hamiltoniana come H = H1 + H2 dove H1 è una deriva e H2 è un calcio (dicono da magneti quadrupolo o di altri campi magnetici).