المساعدة في التكامل المتماثل

Question

أحاول تطوير محاكاة فيزياء وأريد تنفيذ ترتيب رابع التكامل التوافقي طريقة. المشكلة هي أنني يجب أن أخطأ في الرياضيات ، لأن المحاكاة الخاصة بي لا تعمل على الإطلاق عند استخدام التكامل المتماثل (مقارنة بتكامل Runge-Kutta من الدرجة الرابعة التي تعمل بشكل جيد للمحاكاة). لقد كنت غوغل هذا إلى الأبد وكل ما يمكنني العثور عليه هو مقالات علمية حول هذا الموضوع. حاولت تكييف الطريقة المستخدمة في المقالات ، لكنني لا أحالف. أريد أن أعرف ما إذا كان لدى أي شخص رمز مصدر لمحاكاة تستخدم المتكاملين المتماثلون ، ويفضل أن يحاكي مجال الجاذبية ولكن أي تكامل متكامل سيفعله. ما هي اللغة التي لا يهم المصدر أكثر من اللازم ، لكنني أقدر لغة باستخدام بناء الجملة على غرار C. شكرًا!

Answer 1

كما طلبت رمز المصدر: من هنا يمكنك تنزيل رمز MATLAB و FORTRAN للطرق المتماثلة لأنظمة هاميلتون والأساليب المتماثلة للمشاكل القابلة للانعكاس. والكثير من الطرق الأخرى للتعامل مع معادلات الاختلاف أيضًا.

و في هذه ورقة يمكنك العثور على وصف الخوارزميات.

يحرر

إذا كنت تستخدم Mathematica هذه قد يساعد أيضا.

Answer 2

فيما يلي الرمز المصدري لطريقة تكوين الطلب الرابع بناءً على مخطط Verlet. سيظهر الانحدار الخطي لـ $ log_ {10} ( delta t) $ مقابل $ log_ {10} (خطأ) $ ميلًا قدره 4 ، كما هو متوقع من النظرية (انظر الرسم البياني أدناه). يمكن العثور على مزيد من المعلومات هنا, هنا أو هنا.

#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

أنا في مجال فيزياء التسريع (مصادر ضوء السنكروترون) ، وفي نمذجة الإلكترونات تتحرك عبر الحقول المغناطيسية ، نستخدم المتكاملات المتماثلة بشكل منتظم. لدينا العمود الفقري الأساسي هو تكامل Symplectic من الترتيب الرابع. كما هو مذكور في التعليقات أعلاه ، لسوء الحظ ، فإن هذه الرموز ليست موحدة بشكل جيد أو متاحة Easilly.

يسمى رمز تتبع MATLAB المفتوح المصدر واحد Accelerator Toolbox. هناك مشروع SourceForge يسمى Atcollab. شاهد ويكي فوضوي هناhttps://sourceforge.net/apps/mediawiki/atcollab/index.php؟title=main_page

بالنسبة للتكامل ، يمكنك أن تنظر هنا:https://sourceforge.net/p/atcollab/code-0/235/tree/trunk/atintegrators/يتم كتابة المتكاملين في C مع رابط MEX إلى MATLAB. نظرًا لأن الإلكترونات نسبية ، فإن المصطلحات الحركية والمحتملة تبدو مختلفة قليلاً عما كانت عليه في الحالة غير المعتلية ، ولكن لا يزال بإمكان المرء كتابة Hamiltonian على أنه H = H1+H2 حيث يكون H1 عبارة عن انجراف و H2 عبارة أو غيرها من الحقول المغناطيسية).