CODEWITHSUNDEEP
c++simulationrunge-kutta
Simon Ewing
Simon Ewing

Reputation: 125

Runge-Kutta 4th order Particle Advection code sample

I have been trying to implement RK4 integration into a simulation that I am doing. The function below is my best attempt at using RK4 to integrate over a force field in 3 dimensions, based on the equations on page 12 at this site.

In my code, the Particle class essentially stores velocity and position lists, and can calculate the force given a position (the force is independent of velocity). Also, I know that my function is long and drawn out, and can be reduced using a for loop, but I (currently) want to match the structure used in the paper I linked.

When I try to simulate a Particle using this method, the error is significantly worse than when I use a leapfrog integration method. Therefore I think something is wrong with my RK4 implementation. Please let me know if I misunderstand how RK4 works, when using it to solve coupled differential equations.

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

Upvotes: 1

Views: 914

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

You can not integrate one particle at a time, this will result in an order 1 method for the ensemble of particles and thus present a large drift at the step sizes that are appropriate for an order 4 method.

You have to integrate all particles at once, that is, compute stage 0 for all particles, set the state 1 for stage 1 with all particles, compute forces and velocities, that is, the k quantities, for stage 1 for all particles at once from the state 1. Then compute the state 2 for stage 2, compute those k vectors for all particles at once etc.

Upvotes: 1

Related Questions