# Using Runge–Kutta method for position calc

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

I like Serena
Homework Helper
Welcome to PF, Manuel Goucha! Extend your State and Derivative records with an acceleration, and add formulas equivalent to the speed formulas.
That should do the trick.

jtbell
Mentor
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.

AlephZero
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

As jtbell said, it is probably neater to store x1 and x2 in an array of length 2.

jtbell
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.

Thanks for all the exaplanations but as you said to use:

Welcome to PF, Manuel Goucha! Extend your State and Derivative records with an acceleration, and add formulas equivalent to the speed formulas.
That should do the trick.

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

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;
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);

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

I like Serena
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.

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 + Particle.V < -500 <- do a simple collision check
then Particle.V := -A.V;

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

• Fall Test.zip
205.1 KB · Views: 176
Last edited:
I like Serena
Homework Helper
You're not saying what you think is going wrong.

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

I like Serena
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?

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

I like Serena
Homework Helper
Looks as if you are working with uninitialized values.

Looking at your program you have:
Code:
NewDerivative: Derivative;
which is apparently never initialized.

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

I like Serena
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.

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

rcgldr
Homework Helper
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:
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

rcgldr
Homework Helper
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])

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