Hello, I''m trying to implement Runge Kutta integration but am having some trouble. I found the formulas; k1 = f''( t , y ) k2 = f''( t + dt/2, y + k1 * dt/2 ) k3 = f''( t + dt/2, y + k2 * dt/2 ) k4 = f''( t + dt, y + k3 * dt ) k = ( k1 + 2*k2 + 2*k3 + k4 ) / 6 I''m just wondering what to sub into the formulas given above. Currently, I add up all my forces and end up with netforce and do euler integration from there. What values do I put in for k1 .. k4?? Thanks a lot.

The exact format of the algorithm depends a bit on what you''re integrating. In particular, what order is your differential equation? Do you have one equation, or several coupled equations? If you''re doing second order N-body mutually interacting motion, I can list my (seemingly functional) C++ code.

Alternately/preferably, check out:

http://www.ap.univie.ac.at/users/ves/cp0102/dx/node57.html

http://www.ap.univie.ac.at/users/ves/cp0102/dx/node64.html

and

http://www.ap.univie.ac.at/users/ves/cp0102/dx/node69.html


If that doesn''t answer your question, can you be more specific about what you need to know? Eg. do you know how to do Runge-Kutta 4th order integration by pen and paper and just not how to convert to your preferred programming language, or do you not have any idea how the algorithm works? Or is it a matter of applying the algorithm to more complicated sorts of differential equations than you''ve used it for previously?

Thanks for the reply. I didn''t know how to convert the formula''s to my desired program. I think I got it now. I didn''t know how to get k2, k3, k4 using the modified velocity. Now, I just recalculate my forces again and whenever I see velocity used in the calculation, I just use ( velocity + kX, where kX = k1, k2, etc...). Is this correct?

http://www.eskimo.com/~jbm/ballistics/rk.html has more explicit and/or useful versions of the formulas. I can''t tell if what you''re doing is correct, as what you''ve written is rather vague.

If you''re using forces however, you probably have 2nd order differential equations. They also probably coupled equations (as in you have x and y positions and dx/dt and dy/dt velocities, the values of which affect eachothers'' differential equations). If this is the case, I think you need to use the 2nd order coupled equations found on the page linked above (in this post). This is equivalent to 4 coupled first order differential equations. You end up having 4 velocity estimates and 4 acceleration estimates for each dimention of your problem (which are the k values). If you''re doing 2D x and y, that''ll be 16 ''k'' values. (24 for 3D, etc) (not always represented by the letter k).

To get all the k values, add appropriate adjustments to the initial conditions for the time step, storing the results and using them to calculate additional accelerations and velocities as required. I found it useful to store the initial conditions and intermediate values in temporary variables. When done that, use the formulas for y(n+1) in terms of y(n) and various k values to figure out what the new positions and velocities are at the end of the time step.

I hope that was coherant and correct. If someone knows better, please say so.

quote:

---

Original post by Anonymous Poster
Thanks for the reply. I didn''t know how to convert the formula''s to my desired program. I think I got it now. I didn''t know how to get k2, k3, k4 using the modified velocity. Now, I just recalculate my forces again and whenever I see velocity used in the calculation, I just use ( velocity + kX, where kX = k1, k2, etc...). Is this correct?

---

Effectively, yes.

The RK methods are usually clouded by mathematical jargon. Assume that f(t,y)==V(t,y) is your velocity at time t, position y. f'' is dV/dt = A(t,y) (acceleration). Now, k1 finds the acceleration at the start position. Next you take k1 and use it to find a position halfway along the time step, and place the acceleration at that point into k2. Next you take the acceleration k2, and uses it to perform another half-step from the start position,also to the middle, placing the acceleration at that point into k3. [Note that k2 and k3 are going to be very similar, which is why they''re each multiplied by two at the end.] Finally we take k3, then perform a whole step from the start point to the end point, and we take the acceleration at that end point and put it in to k4. k is simply a weighted average of the values, and is a much better value of the acceleration you should use.

As Geoff pointed out, the RK methods are only really effective when you need to consider 2nd- (or higher-) order differential equations.

For example, linear gravity (that is, always pointing downwards) is a 1st order ordinary differential equation (ODE), and so it''s unnecessary to use the RK method, since it just breaks down to Euler anyway.

Velocity in a fluid, air resistance, friction, proper 3D gravity (e.g. from a star), and possibly (depending on how you''ve built your physics model) external forces may all benefit from a better model than Euler''s.


BTW, the best bang-for-buck is the RK2 method:
k1 = f''(t,y)
k2 = f''(t+dt/2, y+k1*dt/2)
k = (k1+k2)/2

It''s vastly more accurate than the Euler method (though less accurate than the RK4 method outlined by revearz), and only twice as slow

[teamonkey]

quote:

---

Original post by teamonkey
[...]
[Note that k2 and k3 are going to be very similar, which is why they''re each multiplied by two at the end.]

---

http://mathworld.wolfram.com/Runge-KuttaMethod.html says that Runge Kutta integration is "A method of numerically integrating ordinary differential equations by using a trial step at the midpoint of an interval to cancel out lower-order error terms." In other words, the k2 and k3 values are multiplied by 2 because that''s the correct weighting (it''s actually 1/3 since you divide the whole sum of k terms by 6) such that error terms from the various estimates cancel eachother out, which is why the integration is a better approximation of the true solution.

quote:

---

Original post by teamonkey
As Geoff pointed out, the RK methods are only really effective when you need to consider 2nd- (or higher-) order differential equations.

---

I actually said (meant to say?) that the particular version of the 4th order equations which were designed for 2nd order coupled differential equations on the page I liked to should be used. I disagree with your further statements:

quote:

---

Original post by teamonkey
For example, linear gravity (that is, always pointing downwards) is a 1st order ordinary differential equation (ODE), and so it''s unnecessary to use the RK method, since it just breaks down to Euler anyway.

---

Linear (constant, independent of height, I assume you meant) gravity is a 2nd order differential equation. Just because the value of the acceleration (the 2nd derivative of the position) is constant, doesn''t mean the differential equation is somehow lower order. The 2nd derivative of the position appears in the equation, thus it is a 2nd order equation. First order, or linear differential equations, have only first deriviates in their equations, which would be the velocity in this case.

What makes you say RK2 or RK4 break down to Euler for first order differential equations (be they your version of what that means, or mine)? Higher order integration methods almost always produce better approximations to the true solution than lower order methods (see below for qualifications of that statement). I don''t know of any reason to assume Euler and an RK method are equivalent for certain type of differential equations...

quote:

---

Original post by teamonkey
BTW, the best bang-for-buck is the RK2 method:

---

The RK2 algorithm is 2nd order, meaning that if you reduce the step size by half, the error in the integration decreases roughly proportional to (0.5)^2. RK4 is 4th order, meaning a half reduction in step size decreases the error roughly proportional to (0.5)^4. Single integration steps done with RK2 occur about twice as fast as those done with RK4, assuming that the force/acceleration cacluation takes significantly longer than the actual integration calculations. If this is the case, then significantly better integration speed may be obtained with RK4 than with RK2 for some accuracies / step sizes. (this assumes speed is determined by ''distance'' or ''elapsed time'' (or whatever the appropriate variable names are) in terms of the problem, not number of steps taken) If the acceleration calculation is trivial, this may not be the case. The actual accuracy may be limited my the accuracy of the computer''s representation of the numbers being integrated, especially for very small step sizes.

quote:

---

Original post by Geoff the Medio
In other words, the k2 and k3 values are multiplied by 2 because that's the correct weighting

---

A weighting. There are a number of different weightings you can use, particularly those discovered by Merson and Fehlberg. You can use any variation of constants (even negative ones), although they should almost always add up to 1. I found a good number of different RK variations in "Solving Ordinary Differential Equations: Nonstiff Problems" by E. Hairer.

quote:

---

What makes you say RK2 or RK4 break down to Euler for first order differential equations (be they your version of what that means, or mine)?

---

It effectively breaks down (not necessarily mathematically). Think of a ball falling in a linear gravitational field. If you treat the acceleration as a constant A that is independent of position, velocity and time it makes the differential equations a lot simpler. Performing the RK4 method in this case makes no sense: k1 = k2 = k3 = k4 = k = A, which is simply the result Euler's method would have given you anyway.

I realise this is simply a special case, however.

quote:

---

The RK2 algorithm is 2nd order, meaning that if you reduce the step size by half, the error in the integration decreases roughly proportional to (0.5)^2.

---

IIRC, the average error of a nth-order algorithm scales with h^(n+1), where h is your time step [Numerical Recipes in C]. So although the RK2 method is 2nd order, the average error is 3rd order. The error for RK4 is 5th-order.

As you say, the limitations of the number format come in to play. If you have an integration time step in your game of 1/30second, then the average error in the RK2 method is about 3E-5, which is of the same order of accuracy expected by 32-bit floats (I think - I could be wrong). The RK4 method would have an error in the order of 4E-8, which means that you're wasting processor time using it if you're using floats, since any advantage of using RK4 is drowned out by the inaccuracy of the floats. Conversely, using the RK2 method negates the accuracy you think you have when using doubles. For comparison, the Euler method would have an average error as large as 0.001! So I would say that the RK2 method is pretty good.

Actually, the best method to use is any of the methods combined with adaptive step-size control. That is, if the error's too great, repeat it with a smaller time steps (you can work out the size of the largest step you need to use to get your desired accuracy). I haven't quite worked out a way that you can do this well in real-time games though.

[teamonkey]

[edited by - teamonkey on July 3, 2003 8:33:49 AM]

quote:

---

Original post by teamonkey
A weighting. There are a number of different weightings you can use [...]

---

Ok, I agree with that. I should have better qualified my statement. And I forgot.

quote:

---

Original post by teamonkey
It effectively breaks down (not necessarily mathematically). Think of a ball falling in a linear gravitational field. If you treat the acceleration as a constant A that is independent of position, velocity and time it makes the differential equations a lot simpler. Performing the RK4 method in this case makes no sense: k1 = k2 = k3 = k4 = k = A, which is simply the result Euler''s method would have given you anyway.

I realise this is simply a special case, however.

---

You''re forgetting the fact that this is a second order differential equation, which is the same as two coupled first order equations. Check out the "Simultaneous Equations of Second Order" on http://www.eskimo.com/~jbm/ballistics/rk.html , which talks about coupled second order quations, which break down to 4 coupled first order equations, but is equilvalent for what I''m talking about. You''re right that all the m and p values (estimates of change in velocity during the time step, what you called k) will be the same. For constant acceleration, this means you could (significantly) simplify the algorithm to improve performance (as you said).

It''s also important to consider the l and k values (estimates of change in position during the time step, which you seem to have forgotten?) which are determined by adding different multiples (0, 1/2 or 1) of the m and p values to the positions at the start of the time step. An average of these (l or k values) is then taken in order to get the change in the position, but it doesn''t simplify to just a simple Euler step. (If I added right, it comes to x(n+1) = x(n) + u(n)*h + m/2 , where x is the position, u is the velocity, and m is the (constant) acceleration).

Obviously a much more efficient version that the full explicit version linked to could be found... but it''s not quite the same as Euler.

You know, unless I screwed up...

quote:

---

Original post by teamonkey
IIRC, the average error of a nth-order algorithm scales with h^(n+1), where h is your time step [Numerical Recipes in C]. So although the RK2 method is 2nd order, the average error is 3rd order. The error for RK4 is 5th-order.

---

I think you''re confusing global order and local order (might have the names wrong...). RK4 is locally 5th order, but globally 4th order. Similarly, RK2 is 3rd local and 2nd global.

I''m a bit fuzzy on this, but I think that essentially means it scales like I said when you use it repeatedly (more than one integration step). I think this arisis due to errors from previous steps increasing the error for the current step, in addition to the current step adding it''s own error that still would have been there even if all steps prior to the current one had been perfectly accurate (true solution).

quote:

---

Original post by teamonkey
Actually, the best method to use is any of the methods combined with adaptive step-size control. That is, if the error''s too great, repeat it with a smaller time steps (you can work out the size of the largest step you need to use to get your desired accuracy). I haven''t quite worked out a way that you can do this well in real-time games though.

---

The "main" problem is that if you start changing the step sizes, the perceived rate of integration will change. It doesn''t matter how often you make step calculations (frame rate), if the steps are really tiny sometimes, and really big other times, then actual real time rate the ''object'' will update its position will vary wildly. The good thing about constant step size algorithms is that you can change the step size due to external factors... like how long it''s been since the previous frame was updated.

I think we agreed on that.

I did mess up, sorta... The formula should have been:

x(n+1) = x(n) + u(n)*h + h*m/2

where m is h*F where F is the acceleration, and h is the step size.

This, oddly enough, simplifies to

x(n+1) = x(n) + u(n)*h + h*h*F/2

which I believe you might recognize as the formula for constant acceleration motion.

So, RK4 actually integrates perfectly, up to machine percision.

Which I''m fairly certain is better than Euler.

Thanks a lot for all your help. I think I get it now.