Trợ giúp với các nhà tích hợp symplectic

Question

Tôi đang cố gắng để phát triển một mô phỏng vật lý và tôi muốn thực hiện một phương pháp thứ tư thứ tự symplectic integration. Vấn đề là tôi phải nhận được toán sai, vì mô phỏng của tôi không hoạt động chút nào khi sử dụng trình tích hợp đối xứng (so với trình tích hợp Runge-Kutta thứ tư hoạt động tốt cho mô phỏng). Tôi đã googling này mãi mãi và tất cả những gì tôi có thể tìm thấy là những bài báo khoa học về chủ đề này. Tôi đã cố gắng để thích ứng với phương pháp được sử dụng trong các bài báo, nhưng không có may mắn. Tôi muốn biết nếu có ai có mã nguồn cho một mô phỏng mà làm cho việc sử dụng tích hợp symplectic, tốt nhất để mô phỏng một trường hấp dẫn nhưng bất kỳ tích hợp symplectic sẽ làm. Ngôn ngữ nguồn là gì không quan trọng quá nhiều, nhưng tôi sẽ đánh giá cao một ngôn ngữ sử dụng cú pháp kiểu C. Cảm ơn!

Answer 1

Khi bạn yêu cầu mã nguồn: Từ HERE bạn có thể tải xuống mã MATLAB và FORTRAN cho các phương pháp đối xứng cho hệ thống Hamilton và các phương pháp đối xứng cho các vấn đề có thể đảo ngược. Và rất nhiều phương pháp khác để xử lý các phương trình khác.

Và trong giấy THIS bạn có thể tìm thấy mô tả về thuật toán.

Sửa

Nếu bạn sử dụng Mathematica this có thể giúp quá.

Answer 2

Tôi đang trong lĩnh vực vật lý gia tốc (nguồn ánh sáng synchrotron), và trong việc mô hình các electron di chuyển qua các từ trường, chúng ta sử dụng các nhà tích hợp symplectic một cách thường xuyên. Workhorse cơ bản của chúng tôi là một bộ tích hợp đối xứng bậc 4. Như đã lưu ý trong các ý kiến ​​ở trên, thật không may các mã này không được chuẩn hóa tốt hay dễ dàng.

Mã theo dõi dựa trên mã nguồn mở Matlab được gọi là Hộp công cụ gia tốc. Có một dự án Sourceforge được gọi là atcollab. Xem wiki lộn xộn ở đây https://sourceforge.net/apps/mediawiki/atcollab/index.php?title=Main_Page

Đối với các nhà tích hợp, bạn có thể xem ở đây: https://sourceforge.net/p/atcollab/code-0/235/tree/trunk/atintegrators/ Các nhà tích hợp được viết bằng C với MEX liên kết đến Matlab. Bởi vì các electron tương đối tương đối với các thuật ngữ động học và tiềm năng trông hơi khác so với trường hợp không tương đối tính, nhưng ta vẫn có thể viết Hamilton là H = H1 + H2 trong đó H1 là trôi và H2 là một cú đá (nói từ tứ cực nam châm hoặc từ trường khác).

Answer 3

Đây là mã nguồn cho phương thức bố cục thứ tự thứ tư dựa trên lược đồ Verlet. Hồi quy tuyến tính của $ \ log_ {10} (\ Delta t) $ so với $ \ log_ {10} (Lỗi) $ sẽ hiển thị độ dốc 4, như mong đợi từ lý thuyết (xem biểu đồ bên dưới). Thông tin thêm có thể được tìm thấy here, here hoặc 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); 
}