Ayuda con integradores simplécticos

Question

Estoy tratando de desarrollar una simulación de física y quiero implementar un método de cuarto orden symplectic integration. El problema es que debo sacar las matemáticas mal, ya que mi simulación no funciona en absoluto cuando uso el integrador simpléctico (en comparación con un integrador Runge-Kutta de cuarto orden que funciona razonablemente bien para la simulación). He estado buscando en Google esto para siempre y todo lo que puedo encontrar son artículos científicos sobre el tema. Traté de adaptar el método utilizado en los artículos, pero no tengo suerte. Quiero saber si alguien tiene un código fuente para una simulación que hace uso de integradores simplécticos, preferiblemente para simular un campo gravitacional, pero cualquier integrador simpléctico lo haría. En qué idioma está la fuente no importa demasiado, pero agradecería un lenguaje que use la sintaxis estilo C. ¡Gracias!

Answer 1

Como solicitó el código fuente: A partir del HERE puede descargar el código MATLAB y FORTRAN para métodos simplécticos para sistemas hamiltonianos y métodos simétricos para problemas reversibles. Y muchos otros métodos para lidiar con ecuaciones de diferencias también.

Y en el documento THIS puede encontrar la descripción de los algoritmos.

Editar

Si utiliza Mathematica this puede ayudar también.

Answer 2

Estoy en el campo de la física del acelerador (fuentes de luz de sincrotrón), y en el modelado de electrones que se mueven a través de campos magnéticos, utilizamos integradores simplécticos de manera regular. Nuestro caballo de batalla básico es un integrador simpléctico de 4º orden. Como se señaló en los comentarios anteriores, desafortunadamente estos códigos no están tan bien estandarizados ni están fácilmente disponibles.

Un código de seguimiento basado en Matlab basado en código abierto se llama Accelerator Toolbox. Hay un proyecto de Sourceforge llamado atcollab. Ver un wiki desordenado aquí https://sourceforge.net/apps/mediawiki/atcollab/index.php?title=Main_Page

Para los integradores, se puede ver aquí: https://sourceforge.net/p/atcollab/code-0/235/tree/trunk/atintegrators/ Los integradores están escritos en C con enlace MEX a Matlab. Debido a que los electrones son relativistas, los términos cinéticos y potenciales se ven un poco diferentes que en el caso no relativista, pero aún se puede escribir el Hamiltoniano como H = H1 + H2 donde H1 es una deriva y H2 es una patada (por ejemplo, de cuadrupolo imanes u otros campos magnéticos).

Answer 3

Aquí está el código fuente para un método de composición de cuarto orden basado en el esquema de Verlet. Una regresión lineal de $ \ log_ {10} (\ Delta t) $ contra $ \ log_ {10} (Error) $ mostrará una pendiente de 4, como se espera de la teoría (ver el gráfico a continuación). Más información se puede encontrar 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); 
}