Using Runge–Kutta method for position calc

  • #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:
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:
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
 

Answers and Replies

  • #2
I like Serena
Homework Helper
6,577
176
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.
 
  • #3
jtbell
Mentor
15,735
3,892
I don't know how to apply Runge–Kutta metthod to a 3 equation physics system
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.
 
  • #4
AlephZero
Science Advisor
Homework Helper
6,994
292
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.
 
  • #5
jtbell
Mentor
15,735
3,892
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.
 
  • #6
Thanks for all the exaplanations but as you said to use:

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.
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:
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
 
  • #7
I like Serena
Homework Helper
6,577
176
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.
 
  • #8
When I run a simulation on a simple particle falling at V0 = 0, I use this conditions:

Code:
Particle = record
  P: TVect2;
  V: TVect2;
  A: TVect2;
  F: TVect2;
end;
Code:
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
 

Attachments

Last edited:
  • #9
I like Serena
Homework Helper
6,577
176
You're not saying what you think is going wrong.
 
  • #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
 
  • #11
I like Serena
Homework Helper
6,577
176
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?
 
  • #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
 
  • #13
I like Serena
Homework Helper
6,577
176
Looks as if you are working with uninitialized values.

Looking at your program you have:
Code:
NewDerivative: Derivative;
which is apparently never initialized.
 
  • #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
 
  • #15
I like Serena
Homework Helper
6,577
176
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.
 
  • #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
 
  • #17
rcgldr
Homework Helper
8,728
546
I don't know much about 4 order derivative expressions.
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 + 1/2 Δt a2
a3 = f(v3, p3)

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

p[i+1] = p + 1/6 Δt (v1 + 2 v2 + 2 v3 + v4)
v[i+1] = v + 1/6 Δt (a1 + 2 a2 + 2 a3 + a4)
a[i+1] = f(v[i+1], p[i+1])
 
Last edited:
  • #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
 
  • #19
rcgldr
Homework Helper
8,728
546
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.
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])
 
  • #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
 

Related Threads on Using Runge–Kutta method for position calc

Replies
7
Views
1K
Replies
5
Views
966
Replies
15
Views
31K
  • Last Post
Replies
10
Views
17K
  • Last Post
Replies
2
Views
8K
Replies
8
Views
845
Replies
3
Views
2K
  • Last Post
Replies
8
Views
3K
  • Last Post
Replies
2
Views
3K
  • Last Post
Replies
5
Views
40K
Top