Search code examples
c++vectorrunge-kutta

Right Runge Kutta 4th method approach?


I have this runge kutta code. However, one mentioned my approach is wrong. And I couldn't really understand why from him, so anyone here, who could give a hint on why this way is wrong?

    Vector3d r = P.GetAcceleration();

    Vector3d s = P.GetAcceleration() + 0.5*m_dDeltaT*r;

    Vector3d t = P.GetAcceleration() + 0.5*m_dDeltaT*s;

    Vector3d u = P.GetAcceleration() + m_dDeltaT*t;

    P.Velocity += m_dDeltaT * (r + 2.0 * (s + t) + u) / 6.0);

====EDIT====

Vector3d are storing the coordinates, x, y, z.

The GetAcceleration returns the acceleration for each x, y, and z.


Solution

  • You have some acceleration function

    a(p,q) where p=(x,y,z) and q=(vx,vy,vz)
    

    Your order 1 system that can be solved via RK4 is

    dotp = q
    dotq = a(p,q)
    

    The stages of the RK method involve an offset of the state vector(s)

    k1p = q
    k1q = a(p,q)
    
    p1 = p + 0.5*dt*k1p
    q1 = q + 0.5*dt*k1q
    k2p = q1
    k2q = a(p1,q1)
    
    p2 = p + 0.5*dt*k2p
    q2 = p + 0.5*dt*k2q
    k3p = q2
    k3q = a(p2,q2)
    

    etc. You can either adjust the state vectors of the point P for each step, saving the original coordinates, or use a temporary copy of P to compute k2, k3, k4.