Using Runge–Kutta method for position calc

In summary, you are trying to calculate a position and velocity of a body over time using an integrator, but you have problems with high precision with few iterations. You use an RK4 algorithm to get better results with fewer iterations.
  • #1
Manuel Goucha
9
0
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, that's why, using RK4, I can get better results with fewer iterations of an amount 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
 
Technology news on Phys.org
  • #2
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
Manuel Goucha said:
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
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
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:

I like Serena said:
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
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

  • Fall Test.zip
    205.1 KB · Views: 222
Last edited:
  • #9
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 happened 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
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
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
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 can't find it anywhere explaining for this type of system
 
  • #17
Manuel Goucha said:
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
Manuel Goucha said:
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 to Using Runge–Kutta method for position calc

1. What is the Runge-Kutta method?

The Runge-Kutta method is a numerical method used to solve ordinary differential equations. It is a popular method for approximating the solution of a differential equation by breaking it down into smaller steps.

2. How does the Runge-Kutta method work?

The Runge-Kutta method works by taking an initial value and using a series of calculations to approximate the solution at the next point. This process is repeated until the desired number of steps is reached, resulting in an approximation of the solution to the differential equation.

3. What are the advantages of using the Runge-Kutta method?

The Runge-Kutta method is known for its accuracy and stability. It can handle a wide range of differential equations and can provide a more accurate solution compared to other numerical methods. It is also relatively easy to implement and can handle both stiff and non-stiff equations.

4. Are there any limitations to using the Runge-Kutta method?

While the Runge-Kutta method is a powerful numerical method, it does have some limitations. It can be computationally expensive for equations with a large number of steps, and it may not be suitable for systems with multiple dimensions or complex boundary conditions.

5. How is the Runge-Kutta method used for position calculation?

The Runge-Kutta method can be used for position calculation by using the differential equations that describe the motion of the object and plugging them into the method. The resulting approximation can then be used to calculate the position of the object at different points in time.

Similar threads

  • Programming and Computer Science
Replies
15
Views
2K
  • Programming and Computer Science
Replies
1
Views
2K
  • Programming and Computer Science
Replies
10
Views
12K
  • Calculus and Beyond Homework Help
Replies
14
Views
2K
  • Programming and Computer Science
Replies
19
Views
6K
  • Programming and Computer Science
Replies
2
Views
4K
  • Astronomy and Astrophysics
Replies
9
Views
9K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
9K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
3K
Back
Top