c++ - Runge-Kutta 四阶粒子平流代码示例

Question

我一直在尝试将 RK4 集成到我正在做的模拟中。下面的函数是我基于第 12 页 this 上的方程式使用 RK4 对 3 维力场进行积分的最佳尝试。网站。

在我的代码中，粒子类主要存储速度和位置列表，并且可以计算给定位置的力(力与速度无关)。另外，我知道我的函数又长又长，可以使用 for 循环减少，但我(目前)想要匹配我链接的论文中使用的结构。

当我尝试使用这种方法模拟粒子时，错误比我使用跳蛙积分方法时严重得多。因此我认为我的 RK4 实现有问题。在使用 RK4 求解耦合微分方程时，如果我误解了 RK4 的工作原理，请告诉我。

// 4th Order Runge-Kutta
void Update(Particle * p, double dt) {

    double * v   = p->getVel();
    double * pos = p->getPos();

    double initPos[3] = {pos[0], pos[1], pos[2]};
    double initVel[3] = {v[0], v[1], v[2]};
    double mass = 0.01;

    double k[4][3]; // related to dv
    double l[4][3]; // related to dr

    p->findForce();

    k[0][0] = dt*p->force[0]/mass;
    k[0][1] = dt*p->force[1]/mass;
    k[0][2] = dt*p->force[2]/mass;

    l[0][0] = dt*v[0];
    l[0][1] = dt*v[1];
    l[0][2] = dt*v[2];

    // Set position to midpoint, using l[0]
    pos[0] = initPos[0] + l[0][0]/2;
    pos[1] = initPos[1] + l[0][1]/2;
    pos[2] = initPos[2] + l[0][2]/2;

     p->findForce();

    k[1][0] = dt*p->force[0]/mass;
    k[1][1] = dt*p->force[1]/mass;
    k[1][2] = dt*p->force[2]/mass;

    l[1][0] = dt*(v[0]+k[0][0]/2);
    l[1][1] = dt*(v[1]+k[0][1]/2);
    l[1][2] = dt*(v[2]+k[0][2]/2);

    // Set position to midpoint, using l[1]
    pos[0] = initPos[0] + l[1][0]/2;
    pos[1] = initPos[1] + l[1][1]/2;
    pos[2] = initPos[2] + l[1][2]/2;

    p->findForce();

    k[2][0] = dt*p->force[0]/mass;
    k[2][1] = dt*p->force[1]/mass;
    k[2][2] = dt*p->force[2]/mass;

    l[2][0] = dt*(v[0]+k[1][0]/2);
    l[2][1] = dt*(v[1]+k[1][1]/2);
    l[2][2] = dt*(v[2]+k[1][2]/2);

    // Set position to endpoint, using l[2]
    pos[0] = initPos[0] + l[2][0];
    pos[1] = initPos[1] + l[2][1];
    pos[2] = initPos[2] + l[2][2];

    p->findForce();

    k[3][0] = dt*p->force[0]/mass;
    k[3][1] = dt*p->force[1]/mass;
    k[3][2] = dt*p->force[2]/mass;

    l[3][0] = dt*(v[0]+k[2][0]);
    l[3][1] = dt*(v[1]+k[2][1]);
    l[3][2] = dt*(v[2]+k[2][2]);

    // Finalize pos and v
    pos[0] = initPos[0] + (l[0][0] + 2*l[1][0] + 2*l[2][0] + l[3][0])/6;
    pos[1] = initPos[1] + (l[0][1] + 2*l[1][1] + 2*l[2][1] + l[3][1])/6;
    pos[2] = initPos[2] + (l[0][2] + 2*l[1][2] + 2*l[2][2] + l[3][2])/6;

    v[0]   = initVel[0] + (k[0][0] + 2*k[1][0] + 2*k[2][0] + k[3][0])/6;
    v[1]   = initVel[1] + (k[0][1] + 2*k[1][1] + 2*k[2][1] + k[3][1])/6;
    v[2]   = initVel[2] + (k[0][2] + 2*k[1][2] + 2*k[2][2] + k[3][2])/6;
}

Answer 1

您不能一次积分一个粒子，这将导致粒子集合采用 1 阶方法，从而在适用于 4 阶方法的步长上出现较大的漂移。

您必须一次整合所有粒子，即计算所有粒子的第 0 阶段，为所有粒子的第 1 阶段设置状态 1，计算力和速度，即 k量，一次从状态 1 计算所有粒子的阶段 1。然后计算阶段 2 的状态 2，一次计算所有粒子的那些 k vector 等。