# Runge Kutta

This topic is 4429 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hey, how you' all doing? Well i followed your advice and started on the Runge kutta method for modeling, so i got these equation of the web site that i got from here: k_1 = hf(x_n,y_n) (4) k_2 = hf(x_n+1/2h,y_n+1/2k_1) (5) k_3 = hf(x_n+1/2h,y_n+1/2k_2) (6) k_4 = hf(x_n+h,y_n+k_3) (7) y_(n+1) = y_n+1/6k_1+1/3k_2+1/3k_3+1/6k_4+O(h^5) and i would like you to check out my logic on what they mean: k_1 = hf(x_n,y_n) (4) this is the first equation that we do, and what we got here is: h=step size in second. f(xn,yn) is the function of position (x) from the initial value of x_n ,from the initial time yn.So my question is, the function itself, is it a basic Euler's calcualtion that is done on there or not? and why would i need the the inital time as opposed to the time step in the function itself? Thx

##### Share on other sites
I only know very little about this integration method and its inner workings. (plz correct me if I'm mistaken somewhere)

What I'm almost sure of, is that it calculates the value of the function at infitesimal distances away from a known value f(x0) at x=x0, using as many orders of derivatives available at x0.
This is merely a one dimensional Taylor series. It allows you to calculate the value of a function at one point, given the values of the function and ALL its derivatives at another -known- point x0. A Taylor series is of the form:

f(x) = f(x0) + f '(x0)*(x-x0) + f ''(x0)*(x-x0)2/2! + f '''(x0)*(x-x0)3/3! + ... + f(n)(x0)*(x-x0)n/n! + ...

Suppose that you know all the values of the function and its derivatives at x=x0. Replacing x with x0+dx...
f(x0+dx) = f(x0) + f '(x0)*dx + f ''(x0)*dx2/2! + f '''(x0)*dx3/3! + ...

Therefore, you can calculate the "behaviour" of the function around a known point x0. The more derivatives you know at x0, the higher the accuracy of the result. Also, keeping dx as small as possible is a good thing to keep in mind.

I hope it's a little more clear now, how this method actually works.
Why it works, is another subject!

edit:
there's a bug with the apostrophes in superscripts, be careful there... The order of derivative in each term must be equal to the number in the factorial in the denominator (and the exponent of dx)

[Edited by - someusername on January 4, 2006 9:52:54 PM]

##### Share on other sites
so all the Runge Kutta method is euler's method wrapped up in Taylors series?
is that correct?
And what is O(h^5) exectly in runge kutta? waht is this o function?
Thx

##### Share on other sites
I don't think I'm qualified to answer that. It should be far more accurate than the euler method, though...

##### Share on other sites
Quote:
 And what is O(h^5) exectly in runge kutta? waht is this o function?

I suppose it's the usual O from algorithm theory. If you state that f(n)=O( g(n)), it means that f(n)<=g(n) for all 'n' in the domain. So, this O(h^5) should be a warning that you should also add a quantity that is smaller in order of magnitude than h^5. That quantity is the rest of the (infinite) derivatives

##### Share on other sites
well, the calculations for eulers and runge kuttas seems to be almost the same, the function being the same, then u just do more steps in there and that's it, therefor i can adampt my Euler's method to fit Runge Kutta method?
But i am not sure what the O(h^5) in teh formula is. I thinks thats an error factor, can somone confirmt the above statements about the same function, and the O function, or deny them?
Thx
And Thx for ur help u post above was great.

##### Share on other sites
The O(h5) is Big-oh notation, a constituent of asymptotic notation (which you'll see used often in algorithmic analysis, for example).

It formally denotes the set of functions

O(f(x)) = {g(x) : f(x) <= c*g(x), x > x0}

for some constant c and x0. Intuitively, this says that for sufficiently large x, c*g(x) is always greater than or equal to f(x), i.e., g(x) is an upper-bound on f(x).

Often, however, the notation is slightly abused and interpreted thusly:

f(x) = O(g(x)) if there exists c and x0 such that f(x) <= c*g(x > x0).

Hence, in the case of O(h5), you may interpret this to be shorthand for c*h5.

The notation is often used to get rid of the details of a function and keep only the important stuff, i.e. the stuff which really dictates it's growth rate at the asymptotic scale (i.e., as the argument goes to infinity). In this sense, Big-oh notation specifies an anonymous function with only it's asymptotic growth rate (the function in the Big-oh notation parantheses) given.

##### Share on other sites
hey, i tried to implemet Runge Kutta method in 1 dimension, but i think i made a mistake, not sure where since i followed the formulas exectly.Here is the relevant part of the source
General data that will be used. What i am modeling here is a 1 kg object moving up at 10 m/s from Earth;s surface.
Set Up:
point[0].dPosX=0;         point[0].dPosY=0;         point[0].dPosZ=0;         point[0].dVeloX=0;         point[0].dVeloY=0;         point[0].dVeloZ=0;         point[0].dMass=5.98e24;         point[1].dPosX=6.38e6;         point[1].dPosY=0;         point[1].dPosZ=0;         point[1].dVeloX=10;         point[1].dVeloY=0;         point[1].dVeloZ=0;         point[1].dMass=1;         DistanceFactor=6.38e4;          TimeFactor=0.1;

Main function for calculations:
void  CDraw::Runge(){int x=0;for(;;){      double k1=0;      double k2=0;      double k3=0;      double k4=0;      k1=TimeFactor*Step(x,iTime);        k2=TimeFactor*Step(x+k1/2,iTime +TimeFactor/2);      k3=TimeFactor*Step(x+k2/2,iTime +TimeFactor/2);         k4=TimeFactor*Step(x+k3,iTime+TimeFactor);       x=x+k1/6+k2/3+k3/3+k4/6;      iTime+=TimeFactor;      cout<<"Time:"<<setfill('0')<<setw(2)<<iTime<<"Velocity x:"<<setfill('0')<<setw(2)<<point[1].dVeloX<<endl;      if(iTime>=1)system("pause");}}double CDraw::Step(double xp,double dt){     point[1].dForceX=0;     point[1].dForceX+=ForceX(point[1],point[0]);     point[1].dVeloX=Velocity(point[1].dVeloX,Acceleration(point[1].dForceX,point[1].dMass),dt);      xp+=point[1].dVeloX*dt;      return xp;}

Supplemetrary stuff that worked b4 on Euler;s simulations, therefor most likely correct, but just in case:
//Do the time scaling using Time factordouble CDraw::Time(long double dTimeMs){     return (dTimeMs*TimeFactor);}//calculate the new velocity. newVelocity=oldVelocioty +AccelerationTimedouble CDraw::Velocity(double dLastVelocity,double dAcceleration,double dTime){     if(bFirst){          return dLastVelocity+dAcceleration*0.5*dTime;     }     return dLastVelocity+dAcceleration*dTime;}//Do acceleration from f=madouble CDraw::Acceleration(double dForce,double dMass){     return dForce/dMass;}//Calculate force in the X direction//Fx=forceGravity*changeinX/radiusdouble CDraw::ForceX(gPOINT point,gPOINT point2){     return (Gravity(point,point2)*(point2.dPosX-point.dPosX)/Pytheg(point,point2));} //Calculate force in the X direction//Fy=forceGravity*changeiny/radiusdouble CDraw::ForceY(gPOINT point,gPOINT point2){     return (Gravity(point,point2)*(point2.dPosY-point.dPosY)/Pytheg(point,point2));}double CDraw::ForceZ(gPOINT point,gPOINT point2){     return (Gravity(point,point2)*(point2.dPosZ-point.dPosZ)/Pytheg(point,point2));}//Gravity, Fg=Gmm/r^2double CDraw::Gravity(gPOINT point,gPOINT point2){     return ((6.67e-11*point.dMass*point2.dMass)/pow(Pytheg(point,point2),2));}//r^2=a^2+b^2//r=sqrt(a^2+b^2)double CDraw::Pytheg(gPOINT point,gPOINT point2){     return sqrt(pow(point2.dPosX-point.dPosX,2)+pow(point2.dPosY-point.dPosY,2)+pow(point2.dPosZ-point.dPosZ,2));}

Thx alot for your time, the responses that i got are well appriciated and are very informative, sry for many question , just trying to get a hang of this stuff.

##### Share on other sites
Quote:
 Original post by someusernameI only know very little about this integration method and its inner workings. (plz correct me if I'm mistaken somewhere)What I'm almost sure of, is that it calculates the value of the function at infitesimal distances away from a known value f(x0) at x=x0, using as many orders of derivatives available at x0. This is merely a one dimensional Taylor series. It allows you to calculate the value of a function at one point, given the values of the function and ALL its derivatives at another -known- point x0. A Taylor series is of the form:f(x) = f(x0) + f '(x0)*(x-x0) + f ''(x0)*(x-x0)2/2! + f '''(x0)*(x-x0)3/3! + ... + f(n)(x0)*(x-x0)n/n! + ...Suppose that you know all the values of the function and its derivatives at x=x0. Replacing x with x0+dx...f(x0+dx) = f(x0) + f '(x0)*dx + f ''(x0)*dx2/2! + f '''(x0)*dx3/3! + ... Therefore, you can calculate the "behaviour" of the function around a known point x0. The more derivatives you know at x0, the higher the accuracy of the result. Also, keeping dx as small as possible is a good thing to keep in mind.I hope it's a little more clear now, how this method actually works.Why it works, is another subject!edit:there's a bug with the apostrophes in superscripts, be careful there... The order of derivative in each term must be equal to the number in the factorial in the denominator (and the exponent of dx)

Forgive me, I'm going to be a little pedantic here... too much teaching [wink]

(1) Runge-Kutta methods are not "merely" a Taylor series. The formulae are consistent with a Taylor series, however, the formulae used to derive the runge-kutta method involve a set of coefficients that are not uniquely determined by Taylor's series and this provides extra degrees of freedom that can improve the general performance of the method in terms of stability and accuracy.

(2) You don't use any integration method to calculate the solution an inifinitesimal point from the initial conditions. That is how integration works analytically, but you won't get very far taking infinitesimal steps [smile] Instead an integration method takes a (generally) small but finite step.

(3) It is important to remember that higher order does not imply greater accuracy. This is often the case, but certainly not always. It depends upon the smoothness of the function you are integrating.

(4) As nilkn explained, big-O notation is asymptotic; It doesn't tell you the size of the error that you get due to truncating the Taylor series.

Have fun! [smile]

##### Share on other sites
thanks for the clarifications. As I said, I only know a little about this method, I just thought I could point out the basic idea and why you need the known point.

Quote:
 (3) It is important to remember that higher order does not imply greater accuracy. This is often the case, but certainly not always. It depends upon the smoothness of the function you are integrating.

Indeed, in many areas, increasing the order introduces error. (approximating with polynomials i.e.)

##### Share on other sites
Quote:
 Original post by reverhey, i tried to implemet Runge Kutta method in 1 dimension, but i think i made a mistake, not sure where since i followed the formulas exectly.Here is the relevant part of the sourceGeneral data that will be used. What i am modeling here is a 1 kg object moving up at 10 m/s from Earth;s surface.Set Up:*** Source Snippet Removed ***Main function for calculations:*** Source Snippet Removed ***Supplemetrary stuff that worked b4 on Euler;s simulations, therefor most likely correct, but just in case:*** Source Snippet Removed ***Thx alot for your time, the responses that i got are well appriciated and are very informative, sry for many question , just trying to get a hang of this stuff.

You are having trouble with this simulation, but it may have nothing to do with the integrator. It could be that the problem is poorly scaled. When that is the case you run the risk of the integrator becoming unstable and producing garbage... but the integrator could still be coded correctly. The difference in the scales of this problem may be the problem so I would suggest you try something simpler first.

Given the size and speed of the mass you are modeling you can treat gravity as a constant force acting downwards without much loss of accuracy. That will make the model much easier to work with for debugging. Also it has an analytic solution. The height of the mass will be

x(t) = x(0) + v(0)*t - g*t2/2

(g is gravitational acceleration = 9.81)

Yay for math!

Hmmm... I took another look at your code and you are doing something really weird there. I think the problem is that you should be modeling this problem as a system but you're not. I can explain that if you want but I'm trying to keep my rambling to a minimum [wink]

In short, try to solve

x'(t) = v(t)
v'(t) = -g

for the initial conditions x(0) = x0 and v(0) = v0.

##### Share on other sites
thx alot for that, i tried to model a simpler system but i am missing something her is simplified code:
void  CDraw::Runge(){int x=0;for(;;){      double k1=0;      double k2=0;      double k3=0;      double k4=0;      k1=TimeFactor*Step(x,iTime);        k2=TimeFactor*Step(x+k1/2,iTime +TimeFactor/2);      k3=TimeFactor*Step(x+k2/2,iTime +TimeFactor/2);         k4=TimeFactor*Step(x+k3,iTime+TimeFactor);       x=x+k1/6+k2/3+k3/3+k4/6;      iTime+=TimeFactor;      cout<<"Velocity x:"<<setfill('0')<<setw(2)<<point[1].dVeloX<<endl;      if(iTime>=1)system("pause");}}double CDraw::Step(double xp,double dt){     point[1].dForceX=-9.81;     point[1].dVeloX=point[1].dVeloX-9.81*dt;     xp+=point[1].dVeloX*dt;      return xp;}

I got something wrong in the method itself since analytical solution don;t work, thats what i have been going from all the time.

BTW, sry i am taking such a long time to get this, i really appriciate ur responses.

I think the Error is in the ITime stuff that results in usage of the new time as opposed of time step upon the calculations, therefor i think this is how it should looke:
void  CDraw::Runge(){int x=0;for(;;){      double k1=0;      double k2=0;      double k3=0;      double k4=0;      k1=TimeFactor*Step(x,0);        k2=TimeFactor*Step(x+k1/2,TimeFactor/2);      k3=TimeFactor*Step(x+k2/2,TimeFactor/2);         k4=TimeFactor*Step(x+k3,TimeFactor);       x=x+k1/6+k2/3+k3/3+k4/6;      iTime+=TimeFactor;      cout<<"Velocity x:"<<setfill('0')<<setw(2)<<point[1].dVeloX<<endl;      if(iTime>=1)system("pause");}}double CDraw::Step(double xp,double dt){     point[1].dForceX=-9.81;     point[1].dVeloX=point[1].dVeloX-9.81*dt;     xp+=point[1].dVeloX*dt;      return xp;}

The answer i get is still wrong (-8 m), but it is not as wrong as the above one (-100 m).
Any ideas?

[Edited by - rever on January 5, 2006 10:11:25 AM]

##### Share on other sites
Quote:
 Original post by reverthx alot for that, i tried to model a simpler system but i am missing something her is simplified code:*** Source Snippet Removed ***I got something wrong in the method itself since analytical solution don;t work, thats what i have been going from all the time.BTW, sry i am taking such a long time to get this, i really appriciate ur responses.I think the Error is in the ITime stuff that results in usage of the new time as opposed of time step upon the calculations, therefor i think this is how it should looke:*** Source Snippet Removed ***The answer i get is still wrong (-8 m), but it is not as wrong as the above one (-100 m).Any ideas?

I think I see the problem. You were close, the problem is with the time but the problem is in the Step() routine. The RK code (before the change) looks correct to me. You should be passing the total time to Step(), not the timestep. You were passing the total time in Runge() (originally) but then treating it like a timestep in Step().

The problem is

point[1].dVeloX=point[1].dVeloX-9.81*dt;

You are not actually integrating the velocity in your code. You are treating it like a known solution. Replace this with something like

point[1].dVeloX=point[1].initialVelocityX-9.81*dt;

Note that 'point[1].initialVelocityX' is a constant and I'm using 'dt' because it is in your code, it would be clearer to rename that variable 't' in Step(). Let me know how that works.

##### Share on other sites
should the velocisty be constant only for a single step and then change in the next step to the value calculated in the previous step?

##### Share on other sites
This website has a good picture.

##### Share on other sites
Quote:
 Original post by revershould the velocisty be constant only for a single step and then change in the next step to the value calculated in the previous step?

If you mean 'should point[1].initialVelocityX change' then no. I would not recommend this way of approaching the problem in general, but it will help you to make progress towards debugging the code.

Sorry, that wasn't clear. point[1].initialVelocityX should never change.

##### Share on other sites
if it should never change then all we are doing is calculating velocity at a time by formula and there is no recursion here at all, that measn that i will run Runge Kutta method only once?the 4 steps i will run only once?

##### Share on other sites
Quote:
 Original post by reverif it should never change then all we are doing is calculating velocity at a time by formula and there is no recursion here at all, that measn that i will run Runge Kutta method only once?the 4 steps i will run only once?

No. The problem is with how you are trying to solve this problem. There are two differential equations that you should be trying to solve, each using RK. However, you are trying to solve for the position using RK but not solving for the velocity using RK. The advise I'm giving you is designed to work within the way you have formulated the problem. I would not approach this problem the way you are doing it.

Those four steps in the RK method are perform for each timestep. If you want to take 1 step, you will perform those four calculations only once. If you want to take 10 steps, you will perform them 10 times. In other words, if you want to use a timstep of 0.1 second and integrate for a total of 1.0 second, you will have to call those four calculation a total of 10 times.

Have you made the change I suggested?

##### Share on other sites
well if i am aproaching the problem wrong then i should change that. So what will do now is i will use RK method to get the veolocities and then use velocities to use the position.Brb in a few mins, i will rewrite the method to accomodate the above requriment.
Thx alot, i think i am getting this now.I will also re read some of the posts to check for any more info that i might ahve misunderstood b4.

##### Share on other sites
Quote:
 Original post by reverwell if i am aproaching the problem wrong then i should change that. So what will do now is i will use RK method to get the veolocities and then use velocities to use the position.Brb in a few mins, i will rewrite the method to accomodate the above requriment.Thx alot, i think i am getting this now.I will also re read some of the posts to check for any more info that i might ahve misunderstood b4.

I hate trying to explain this stuff in text - I'm not very good at it :-( I wish I could show you on a white board. I'm coding up RK4 now to provide a brief demo. I'll post it when it's done. Hopefully, that'll be useful.

##### Share on other sites
that will be very usefull, and if u could can it be something simple like a motion in 1 dimension or something that is easyly understoddd by idiots like me.
Thx alot, this is more help then i hoped i would get:D

##### Share on other sites
Quote:
 Original post by reverthat will be very usefull, and if u could can it be something simple like a motion in 1 dimension or something that is easyly understoddd by idiots like me.Thx alot, this is more help then i hoped i would get:D

No problem [smile]

Here is the demo.

##### Share on other sites
Quote:
 Original post by jjdI hate trying to explain this stuff in text - I'm not very good at it :-( I wish I could show you on a white board. I'm coding up RK4 now to provide a brief demo. I'll post it when it's done. Hopefully, that'll be useful.

Haha and modest too. It was actually a better explanation than what I've been reading in some "good" books and it was helpful for me as well.

Thanks!

##### Share on other sites
Quote:
Quote:
 Original post by jjdI hate trying to explain this stuff in text - I'm not very good at it :-( I wish I could show you on a white board. I'm coding up RK4 now to provide a brief demo. I'll post it when it's done. Hopefully, that'll be useful.

Haha and modest too. It was actually a better explanation than what I've been reading in some "good" books and it was helpful for me as well.

Thanks!

Thanks! Glad you found it useful [smile][smile][smile]

##### Share on other sites
Thx alot man, like the above person said this is prolly one of the best demos of RK4 out there, u should definitly write an article and post it in the archive on gamedev.