2014-01-21 1 views
1

저는 4 차 Runge-Kutta 솔버에서 작업 중이며 다소 어려움이 있습니다. 필자는 on gafferongames 문서를 기반으로 솔버를 작성했지만, 작은 예제를 실행하면 간단한 중력이라도 단순한 오일러 통합으로 얻는 것보다 훨씬 오류가 있습니다. 필자는 자체 포함 된 예제 (인쇄를 포함하여 ~ 60 행의 코드)로 정리했지만 GLM을 실행해야합니다.
그것은 전체적으로 내 문제점을 보여줍니다. 55 행은 분석 솔루션과 RK4 솔루션의 차이점을 인쇄합니다. 이것은 상대적으로 작아야하지만 10 단계의 단계를 거친 후에도 폭발합니다.Runge-Kutta 4 차 적분기 오류

#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; 
} 

누구든지 나를 도울 수 있다면, 나는 이것으로 나의 밧줄 끝에 있습니다.

답변

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;