Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Using Runge–Kutta method for position calc

  1. Nov 24, 2011 #1
    Hi, I'm trying to calculate a postion and velocity of a body over time using an integrator at each time step, I've only used simple integrators so far but I wanted to a better one that I've seen, RK4 - Runge–Kutta method, to calculate new values of position equation.

    I've been using the simple SUVAT and euler method:

    a = F/m
    v = v0 + a*dt
    x = x0 + v0*dt + 0.5*a*dt*dt

    but this gives me bad results when needing high precision with few iterations, thats why, using RK4, I can get better results with fewer iterations of an amout of calculues leaving more free cpu.

    my problem is that I don't know how to apply Runge–Kutta metthod to a 3 equation physics system, because Runge–Kutta algorithm only uses v and x, so how can I adapt it to my system?

    I have this algorithm:

    Code (Text):
    EulerIntegrate(Pos, Vel, Acc, For, Mass, dt: Double);
    var
      MDiv1: TValue;
      DTDiv: TValue;
      i: Integer;
    begin
      Pos := Pos + (Vel * dt);
      Vel  := Vel + (Acc * dt);
      Acc := Acc + F / M;
    end;
    and the RK4 algotithm I'm trying to use:

    Code (Text):
    RK4Integrate(var Pos, Vel, Acc, dt: double);

    type
      State = record
        Pos: TVect2;          // position
        Vel: TVect2;          // velocity
      end;

      Derivative = record
        DPos: TVect2;          // derivative of position: velocity
        DVel: TVect2;          // derivative of velocity: acceleration
      end;

      function Evaluate(Initial: State; Acc, dt: Double; D: Derivative): Derivative;
      var
        State: State;
        Output: Derivative;
      begin
        State.Pos := Initial.Pos + D.DPos * DT;
        State.Vel  := Initial.Vel + D.DVel * DT;

        Output.DPos := State.Vel;
        Output.DVel := Acc;

        Result := Output;
      end;

      procedure Integrate(var State: State; Acc, dt: Double);
      var
        A, B, C, D: Derivative;
        NewDerivative: Derivative;
        DXDT, DVDT: Double;
      begin
        A := Evaluate(State, Acc, 0,         NewDerivative);
        B := Evaluate(State, Acc, DT*0.5, a);
        C := Evaluate(State, Acc, DT*0.5, b);
        D := Evaluate(State, Acc, DT,       c);

        DXDT := (A.DPos + ((B.DPos + C.DPos)* 2) + D.DPos) * 1/6;
        DVDT := (A.DVel + ((B.DVel + C.DVel)* 2) + D.DVel) * 1/6;

        State.Pos := (State.Pos + DXDT) * DT;
        State.Vel := (State.Vel + DVDT) * DT;
      end;

    var
      Start:    State;
      Derivate: Derivative;

    begin
      Start.Pos := P;
      Start.Vel := V;

      Integrate(Start, A, DT);

      P := Start.Pos;
      V := Start.Vel;
    end;
    but this doens't work, and it doesn't gives me any change at the Acceleration over time, I can't find any better informations on this on the internet and cna't find a way to implement like the first system example, where I give 3 variables - Pos, Vel, Acc, and 2 known values, Force and Mass

    I really apreciate any help solving this
     
  2. jcsd
  3. Nov 24, 2011 #2

    I like Serena

    User Avatar
    Homework Helper

    Welcome to PF, Manuel Goucha! :smile:

    Extend your State and Derivative records with an acceleration, and add formulas equivalent to the speed formulas.
    That should do the trick.
    Oh, and add an DADT in your Integrate procedure.
     
  4. Nov 24, 2011 #3

    jtbell

    User Avatar

    Staff: Mentor

    Do a Google search for something like "runge kutta system of equations".

    Basically, your x, v and t become components of a "vector" and you write your equations in the form of a vector function. The Runge-Kutta method generalizes from a scalar function to a vector function.
     
  5. Nov 24, 2011 #4

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    The usual way to solve a second order equation with RK is to convert it into a system of two first order equations. If you call the variables x1 and x2, you have

    dx1/dt = x2
    dx2/dt = f/m

    where x1 corresponds to your "x" and x2 to your "v".

    As jtbell said, it is probably neater to store x1 and x2 in an array of length 2.
     
  6. Nov 24, 2011 #5

    jtbell

    User Avatar

    Staff: Mentor

    For a small system like this, it's feasible to write out the equations individually, in non-vector form. I've done this for solving them in an Excel spreadsheet. But textbooks and other references generally present the R-K equations in vector form, so you have to start with that, and derive the component (non-vector) equations that apply for your specific situation.
     
  7. Nov 25, 2011 #6
    Thanks for all the exaplanations but as you said to use:

    I've already done that but for some reasson it still doesn't work correctly:

    P = pos vector
    V = vel Vector
    A = Acc vector
    F = force Vector

    Code (Text):
    procedure RK4Integrate(var P, V, A, F: TVect2; DT: Single);

    type
      State = record
        Pos: Vect2;          // position
        Vel: Vect2;          // velocity
        Acc: Vect2;          // acceleration
      end;

      Derivative = record
        DPos: Vect2;          // derivative of position: velocity
        DVel: Vect2;          // derivative of velocity: acceleration
        DAcc: Vect2;          // derivative of acceleration
      end;

      function Evaluate(Initial: State; For: Vect2; DT: Single; D: Derivative): Derivative;
      var
        State: TState;
        Output: TDerivative;
      begin
        State.Pos := VectorAdd(Initial.Pos, VectorScale(D.DPos, DT));
        State.Vel := VectorAdd(Initial.Vel, VectorScale(D.DVel, DT));
        State.Acc := VectorAdd(Initial.Acc, VectorScale(D.DAcc, DT));

        Output.DPos := State.Vel;
        Output.DVel := State.Acc;
        Output.DAcc := For;

        Result := Output;
      end;

      procedure Integrate(var State: State; For: Vect2; DT: Single);
      var
        A, B, C, D: Derivative;
        NewDerivative: Derivative;
        DXDT, DVDT, DADT: Vect2;
      begin
        A := Evaluate(State, For, 0,         NewDerivative);
        B := Evaluate(State, For, DT*0.5, a);
        C := Evaluate(State, For, DT*0.5, b);
        D := Evaluate(State, For, DT,       c);

        DXDT := VectorScale(VectorAdd(A.DPos, VectorAdd(VectorScale(VectorAdd(B.DPos, C.DPos), 2), D.DPos)), 1/6);
        DVDT := VectorScale(VectorAdd(A.DVel, VectorAdd(VectorScale(VectorAdd(B.DVel, C.DVel), 2), D.DVel)), 1/6);
        DADT := VectorScale(VectorAdd(A.DAcc, VectorAdd(VectorScale(VectorAdd(B.DAcc, C.DAcc), 2), D.DAcc)), 1/6);


        State.Pos := VectorScale(VectorAdd(State.Pos, DXDT), DT);
        State.Vel := VectorScale(VectorAdd(State.Vel, DVDT), DT);
        State.Acc := VectorScale(VectorAdd(State.Acc, DADT), DT);
      end;

    var
      Start:    State;
      Derivate: Derivative;

    begin
      Start.Pos := P;
      Start.Vel := V;
      Start.Acc := A;

      Integrate(Start, F, DT);

      P := Start.Pos;
      V := Start.Vel;
      A := Start.Acc;
    end;
    is it something like this?


    thanks again all for helping
     
  8. Nov 25, 2011 #7

    I like Serena

    User Avatar
    Homework Helper

    So what is it that does not work properly?

    I can see that you put in a constant force for all points, which can't be right.
    It defeats the purpose of adding acceleration to the solution.
     
  9. Nov 25, 2011 #8
    When I run a simulation on a simple particle falling at V0 = 0, I use this conditions:

    Code (Text):

    Particle = record
      P: TVect2;
      V: TVect2;
      A: TVect2;
      F: TVect2;
    end;
     
    Code (Text):

    procedure SimStart;
    begin
      Particle.P := Vector(0, 0); <- Start its Position
      Particle.A := Vector(0, -9.8); <- Start its Acceleration as normal Gravity acceleration
    end;

    procedure TimeStep;
    begin
      if Particle.P[1] + Particle.V[1] < -500 <- do a simple collision check
      then Particle.V[1] := -A.V[1];

      RK4Integrate(Particle.P, Particle.V, Particle.A, Particle.F, 0.1, 1); <- integrate with RK4
      DrawParticle(Particle); <- Draw to view particle position;
    end;
     
    if I change the integrator to the euler it works and I see the simulation running, but if I change it to the RK4, it doens't even start, I get numeric erors, The problem is that I ont even know what I using in this expression to check and try to find the errors of sintax

    I'm using delphi to test this, I've sent the code as attachment

    thanks again so far
     

    Attached Files:

    Last edited: Nov 25, 2011
  10. Nov 25, 2011 #9

    I like Serena

    User Avatar
    Homework Helper

    You're not saying what you think is going wrong.
     
  11. Nov 25, 2011 #10
    Thats my problem, I don't know much about 4 order derivative expressions, and I think I'm close but somehow or a numeric error appears or the results get wrong like geting very large numbers or oscilating from positive to negative very quickly, so if I need to have a falling body to get more and more speed at each timestep, it just doesn't happend like that

    I don't know who else to run to, since none of my teachers know about this method and none of the persons I know even likes physics
     
  12. Nov 25, 2011 #11

    I like Serena

    User Avatar
    Homework Helper

    Okay, could you show then which numerical errors you get?

    And perhaps you can print out the values used and calculated in each iteration?
    In particular the first 2 iterations?
     
  13. Nov 25, 2011 #12
    The only results I'm getting now is this:

    Step 0:
    Position.X = 100;
    Position.Y = 0;
    Velocity.X = 0;
    Velocity.Y = 0;
    Acceleration.X = 0;
    Acceleration.Y = 9.8;

    Step1:
    Position.X = 1;
    Position.Y = -0.0081;
    Velocity.X = 0;
    Velocity.Y = 0;
    Acceleration.X = 0;
    Acceleration.Y = -0.097;

    Step2:
    Position.X = 0.00999;
    Position.Y = 0;
    Velocity.X = 0;
    Velocity.Y = 0;
    Acceleration.X = 0;
    Acceleration.Y = -0.000097;

    and after that the simulation can't proces because the numbers are too small
    So I conclude that it cant't be right, when it changes the X value of the position
     
  14. Nov 25, 2011 #13

    I like Serena

    User Avatar
    Homework Helper

    Looks as if you are working with uninitialized values.

    Looking at your program you have:
    Code (Text):
    NewDerivative: Derivative;
    which is apparently never initialized.
     
  15. Nov 25, 2011 #14
    Yes, I've seen that, but this is an adaptation from another code I had here, is that the only problem? if so, what should I use insted? A 0 Value Derivative?

    I've tried to use this:

    NewDerivative.DPos := Vector(0, 0);
    NewDerivative.DVel := Vector(0, 0);
    NewDerivative.DAcc := Vector(0, 0);

    but same error
     
  16. Nov 25, 2011 #15

    I like Serena

    User Avatar
    Homework Helper

    You need to set the derivative based on your input values.
    DPos = Vel
    DVel = Acc
    DAcc = <some formula for the change in acceleration>

    You could set DAcc to zero if you're using a constant force.
     
  17. Nov 25, 2011 #16
    I think I'm trying to do something impossible here, or I might be cursed with this code, because everything I try to change it always gives me bad results or impossible solving systems.

    I've changes the NewDerivative to the starting values but it is still giving me wrong values:

    NewDerivative.DPos := State.Pos;
    NewDerivative.DVel := State.Vel;
    NewDerivative.DAcc := State.Acc;

    I think my code could be all incorrect, could you plase explain me if 4th order Runge–Kutta could de used to calculate a system with 3 variables, Pos, Vel, Acc, and 2 constants Force and Mass? I cant find it anywhere explaining for this type of system
     
  18. Nov 25, 2011 #17

    rcgldr

    User Avatar
    Homework Helper

    The term "4th order" just refers to the fact that RK4 errors are related to Δt4 (where Δt is the time step).

    Here's an example RK4 double integration. I found similar algorithms doing a web search. a1 is acceleration at the start of the time step, a2 and a3 are acceleration at middle of the time step ("average" acceleration), and a4 is acceleration at end of time step. v1 is velocity at start of time step, v2 and v3 are velocity at middle of time step ("average" velocity), and v4 is velocity at end of time step.

    equation for acceleration a = f(v, p)

    p = current position
    v = current velocity
    a = current acceleration = f(v, p)

    rk4 iteration (position based on current estimated velocity):

    v1 = v
    p1 = p

    a1 = f(v1, p1)
    v2 = v + 1/2 Δt a1
    p2 = p + 1/2 Δt v2

    a2 = f(v2, p2)
    v3 = v + 1/2 Δt a2
    p3 = p + 1/2 Δt v3

    a3 = f(v3, p3)
    v4 = v + Δt a3
    p4 = p + Δt v4

    a4 = f(v4, p4)

    v[i+1] = v + 1/6 Δt (a1 + 2 a2 + 2 a3 + a4)
    p[i+1] = p + 1/6 Δt (v1 + 2 v2 + 2 v3 + v4)
    a[i+1] = f(v[i+1], p[i+1])

    rk4 iteration (position based on previous estimated velocity):

    p1 = p
    v1 = v
    a1 = f(v1, p1)

    p2 = p + 1/2 Δt v1
    v2 = v + 1/2 Δt a1
    a2 = f(v2, p2)

    p3 = p + 1/2 Δt v2
    v3 = v[i] + 1/2 Δt a2
    a3 = f(v3, p3)

    p4 = p[i] + Δt v3
    v4 = v[i] + Δt a3
    a4 = f(v4, p4)

    p[i+1] = p[i] + 1/6 Δt (v1 + 2 v2 + 2 v3 + v4)
    v[i+1] = v[i] + 1/6 Δt (a1 + 2 a2 + 2 a3 + a4)
    a[i+1] = f(v[i+1], p[i+1])[/i][/i][/i][/i][/i]
     
    Last edited: Nov 25, 2011
  19. Nov 25, 2011 #18
    Thanks. My problem was to find a method that uses the system of 3 vars, and that doesn't need to have a function to get the acceleration, I needed to give the aceleration instead and it will calculate those 3 vars only by having the force as a given constant, not a function.


    but it seems to work after I adapted it
     
  20. Nov 25, 2011 #19

    rcgldr

    User Avatar
    Homework Helper

    If you can directly integrate acceleration into velocity and position, there's no point in using a method like Runge-Kutta. If acceleration is not a function of velocity or position, there's no point in using the feedback steps used in RK4. Instead of Euler, you could use the trapezoidal method which uses the average acceleration (or velocity) for each time step:

    ...
    a[i+1] = f(i+1)
    v[i+1] = v + 1/2 Δt (a + a[i+1])
    p[i+1] = p + 1/2 Δt (v + v[i+1])
     
  21. Nov 26, 2011 #20
    Yes, I 've figured that Heun's method was better that RK for this type of code but I wanted to use the more precise and exact method there is, so I thoug as RK was a forth order, that I would be better than Suvat, Heun or verlet
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Using Runge–Kutta method for position calc
Loading...