The Runge-Kutta algorithm is the magic formula behind most of the physics simulations shown on this web site. The Runge-Kutta algorithm lets us solve a differential equation numerically (that is, approximately); it is known to be very accurate and well-behaved for a wide range of problems.
Consider the single variable problem
x' = f (t, x)
with initial condition x(0) = x0 . Suppose that xn is the value of the variable at time tn . The Runge-Kutta formula takes xn and tn and calculates an approximation for xn+1 at a brief time later, tn+h . It uses a weighted average of approximated values of f (t, x) at several times within the interval (tn, tn+h) . The formula is given by
xn+1 = xn + h⁄6 (a + 2 b + 2 c + d)
where
a = f (tn, xn)
b = f (tn + h⁄2, xn + h⁄2 a)
c = f (tn + h⁄2, xn + h⁄2 b)
d = f (tn + h, xn + h c)
To run the simulation, we start with x0 and find x1 using the formula above. Then we plug in x1 to find x2 and so on.
With multiple variables, the Runge-Kutta algorithm looks similar to the previous equations, except that the variables become vectors.
The Runge-Kutta Algorithm is fairly simple, but to describe it precisely we need to develop some notation. Suppose there are m variables x1, x2, ..., xm each of which vary over time. For example, in the single spring simulation, x1 is position, x2 is velocity. Suppose further that there are m differential equations for these m variables
x1' = f1(x1, x2, ..., xm)
x2' = f2(x1, x2, ..., xm)
...
xm' = fm(x1, x2, ..., xm)
Notice there are no derivatives on the right hand side of any of those equations, and there are only first derivatives on the left hand side. These equations can be summarized in vector form as
x' = f (x)
where x = (x1, x2, ..., xm) and we allow some loose "vector of functions" concept where f = (f1, f2, ..., fm) . Next we label our time states xn, xn+1 which are separated by time interval of length h . That is, xn is the value of the m variables at time tn . And x1,n is the value of the first variable x1 at time tn .
xn = (x1,n, x2,n, ..., xm,n)
xn+1 = (x1,n+1, x2,n+1, ..., xm,n+1)
Suppose we have the state of the simulation at time tn as xn . To compute the state a short time h later and put the results into xn+1 , the Runge-Kutta algorithm does the following:
an = f(xn)
bn = f(xn + h⁄2 an)
cn = f(xn + h⁄2 bn)
dn = f(xn + h cn)
xn+1 = xn + h⁄6 (an
+ 2 bn + 2 cn + dn)
The new vector xn+1 gives you the state of the simulation after the small time h has elapsed. To spell out the above in more detail, we can drop the vector notation and write the Runge-Kutta algorithm like this:
xj, n = fj(x1,n, x2,n, . . . , xm,n)
xj, n = fj( (x1, n + h⁄2 x1, n), (x2, n
+ h⁄2 x2, n), . . . , (xm, n
+ h⁄2 xm, n) )
xj, n = fj( (x1,n + h⁄2 x1,n), (x2,n
+ h⁄2 x2,n), . . . , (xm,n
+ h⁄2 xm,n) )
xj, n = fj( (x1,n + h x1,n), (x2,n
+ h x2,n), . . . , (xm,n
+ h xm,n) )
xj, n+1 = xj, n + h⁄6 (xj, n + 2 xj, n
+ 2 xj, n + xj, n)
The above equations are applied for each variable j=(1, ..., m) to get the full set of variables in the vector xn+1 .
Most of the simulations shown on this website do not have differential equations that depend explicitly on time. That is, you won't see the variable t on the right-hand side of the differential equations. One simulation that does depend on time is the Chaotic Driven Pendulum because the driving force (which applies the twist to the pendulum) varies over time according to cos(k t) .
When time appears explicitly in the differential equations we can add a time variable t to the state vector x . Suppose we assign this role to the variable x2 . This new variable has the extremely simple differential equation
x2' = 1
That says that the rate of change of the variable x2 is a constant. Since we are taking derivatives with respect to time we can also write the above equation as
$$x_2' = \frac{d}{d t} x_2 = 1$$
This integrates very easily to give x2 = t , which is what we wanted: time as a variable. Suppose that in the driven pendulum simulation we set up x2 in this way. Then the driving force is given by cos(k x2) .
You may ask: Why have time as a variable? We already know the value of t at each time step! The Runge-Kutta algorithm works by averaging the predicted rates at various points in the time interval from t to t+h . Therefore, when the rates (differential equations) depend explictly on t , we also need to know the value of t at those points within the time interval. Putting time in as a variable makes for nicer cleaner computer code.
If you want to, you can avoid keeping time as an additional variable. The following is an equivalent formulation of the Runge-Kutta algorithm where time t is passed in as a variable to each function in f .
an = f(t, xn)
bn = f(t + h⁄2, xn + h⁄2 an)
cn = f(t + h⁄2, xn + h⁄2 bn)
dn = f(t + h, xn + h cn)
xn+1 = xn + h⁄6 (an
+ 2 bn + 2 cn + dn)
This is equivalent to the formulation where time is kept as one of the variables in x . Whether you use this formulation or the earlier (cleaner) one is entirely up to you.
The source code for the Runge-Kutta algorithm is available.
This web page was first published April 2001.
