Runge-Kutta 4阶粒子平流代码示例

Question

我一直在试图将RK4集成到我正在做的模拟中。下面的功能是我最好的尝试，使用RK4根据第12页上的方程在this网站上的三维力场进行积分。

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

当我尝试使用此方法模拟一个粒子时，该错误显着比我使用跳跃积分方法时更糟糕。因此我认为我的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的状态1，包含所有粒子，计算力和速度，即阶段1的k数量所有的粒子一次从状态1.然后计算阶段2的状态2，一次计算所有粒子的那些k向量等。