Runge-Kutta 四阶积分器出错答案

【问题标题】：Error with Runge-Kutta 4th order integratorRunge-Kutta 四阶积分器出错
【发布时间】：2014-02-10 09:03:11
【问题描述】：

我一直在研究 4 阶 Runge-Kutta 求解器，但遇到了一些困难。我已经根据文章on gafferongames 编写了求解器，但是当我运行一个包含的小示例时，我得到的错误远比我使用简单的欧拉积分得到的错误要糟糕得多，即使对于简单的重力也是如此。我已将其整理成一个独立的示例（约 60 行代码，包括打印），但它需要 GLM 才能运行。
它完整地显示了我的问题。第 55 行打印出解析解和 RK4 解之间的差异。这应该相对较小，但即使在它采取的 10 个奇数步骤之后它也会爆炸。

#include <iostream>
#include <glm/glm.hpp>
struct State{
    glm::vec3 position, velocity;
};
class Particle{
public:
    glm::vec3 position, velocity, force;
    float mass;

    void solve(float dt);
    glm::vec3 acceleration() const {return force/mass;}
    State evalDerivative(float dt, const State& curr);
    void analytic(float t,  glm::vec3 a);
};
int main(int argc, char* argv[]){
    Particle p;
    p.position = glm::vec3(0.f);
    p.mass = 1.0f;
    for(int i = 1; i <= 10; i++)    {
        p.force = glm::vec3(0.f, -9.81f, 0.f);
        p.solve(.016f);
        p.analytic(i*.016f, glm::vec3(0.f, -9.81f, 0.f));
    }
    getchar();
    return 0;
}
void Particle::solve(float dt){
    State t;t.position = glm::vec3(0.f); t.velocity = glm::vec3(0.f);
    State k1 = evalDerivative(0, t);
    State k2 = evalDerivative(dt*.5f, k1);
    State k3 = evalDerivative(dt*.5f, k2);
    State k4 = evalDerivative(dt, k3);

    position += (k1.position + 2.f*(k2.position + k3.position) + k4.position)/6.f;
    velocity += (k1.velocity + 2.f*(k2.velocity + k3.velocity) + k4.velocity)/6.f;
    force = glm::vec3(0.f);
}
State Particle::evalDerivative(float dt, const State& curr){
    State s;
    s.position = position + curr.position*dt;
    s.velocity = velocity + curr.velocity*dt;

    s.position = s.velocity;
    s.velocity = acceleration();
    return s;
}
void Particle::analytic(float t, glm::vec3 a){
    glm::vec3 tPos = glm::vec3(0.f) + 0.5f*a*t*t;
    glm::vec3 tVel = glm::vec3(0.f) + a*t;

    glm::vec3 posdiff = tPos - position;
    glm::vec3 veldiff = tVel - velocity;
    std::cout << "POSITION: " << posdiff.x << ' ' << posdiff.y << ' ' << posdiff.z << std::endl;
    std::cout << "VELOCITY: " << veldiff.x << ' ' << veldiff.y << ' ' << veldiff.z << std::endl << std::endl;
}

如果有人能帮我解决这个问题，我已经束手无策了。

【问题讨论】：

标签： c++ math numerical-methods runge-kutta

【解决方案1】：

好吧，我觉得我很愚蠢。我已经为此工作了几个小时，但我错过了一小步：

position += (k1.position + 2.f*(k2.position + k3.position) + k4.position)/6.f;
velocity += (k1.velocity + 2.f*(k2.velocity + k3.velocity) + k4.velocity)/6.f;

应该是

position += (k1.position + 2.f*(k2.position + k3.position) + k4.position)*dt/6.f;
velocity += (k1.velocity + 2.f*(k2.velocity + k3.velocity) + k4.velocity)*dt/6.f;

【讨论】：