Class 6.1

Integration algorithms: Euler, Runge-Kutta, implicit vs. explicit

In most cases the differential equations that we create for a model are too complex to solve and analyze analytically. There are numerical methods to use to find the solutions. However we should remember that these solutions will be approximate, and the larger the time step we use, the more error we generate.

Consider the single variable problem

x' = f(t,x)

with initial condition x(0)=x₀. Suppose that x_n is the value of the variable at time t_n.

The Euler algorithm is the most straightforward way to solve a difference equation that is used to approximate the differential equation.

x_n+1 = x_n + f(t_n,x_n) dt

The major source of error is that we assume the right-hand-side function f to be constant over the time step dt. Obviously, the smaller the dt, the less the difference between the continuos function f(t,p) and its discrete approximation.

The Runge-Kutta algorithm interpolates the function between t_n and t_n+1 in a tricky way. The Runge-Kutta formula uses a weighted average of approximated values of f(t,x) at several times within the interval (t_n,t_n+dt). The formula is given by

x_n+1 = x_n + (dt/6)(k₁ + 2k₂ + 2k₃ + k₄)

where

k₁ = f(t_n, x_n)
k₂ = f(t_n + dt/2, x_n + (dt/2)k₁)
k₃ = f(t_n + dt/2, x_n + (dt/2)k₂)
k₄ = f(t_n + dt, x_n + dt k₃)

To run the simulation, we simply start with x₀ and find x₁ using the formula above. Then we plug in x₁ to find x₂ and so on.

Note that even though we use 6 intermediate values when calculating x_n+1, we calculate the function f only 4 times. Calculating the right-hand-side function is what usually takes most of the CPU time.

Two-variable Runge-Kutta Algorithm

The Runge-Kutta algorithm can be easily extended to a set of first order differential equations.  You wind up with essentially the same formulas shown above, but all the variables (except for time) are vectors.

To give an idea of how it works, here is an example where we expand the vector notation.  That is, instead of using the highly compact vector notation, we write out all the components of the vector.

Suppose there are only two variables, x, y and two differential equations

x' = f(t,x,y)
y' = g(t,x,y)

Again we let x_n be the value of x at time t_n, and similarly for y_n.  Then the formulas for the Runge-Kutta algorithm are

x_n+1 = x_n + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)
y_n+1 = y_n + (h/6)(j₁ + 2j₂ + 2j₃ + j₄)

where

k₁ = f(t_n, x_n, y_n)
j₁ = g(t_n, x_n, y_n)
k₂ = f(t_n + h/2, x_n + (h/2)k₁, y_n + (h/2)j₁)
j₂ = g(t_n + h/2, x_n + (h/2)k₁, y_n + (h/2)j₁)
k₃ = f(t_n + h/2, x_n + (h/2)k₂, y_n + (h/2)j₂)
j₃ = g(t_n + h/2, x_n + (h/2)k₂, y_n + (h/2)j₂)
k₄ = f(t_n + h, x_n + h k₃, y_n + h j₃)
j₄ = g(t_n + h, x_n + h k₃, y_n + h j₃)

Given starting values x₀,y₀ we can plug them into the formula to find x₁,y₁. Then we can plug in x₁,y₁ to find x₂,y₂ and so on.

Multi-variable Runge-Kutta Algorithm

Suppose we have more than two variables which are x, y, z,.... We then need to convert the differential equations into a system of first order equations, as follows:

x' = f(t, x, y, z, ...)
y' = g(t, x, y, z, ...)
z' = h(t, x, y, z, ...)
...

You will have as many equations as there are variables. The algorithm is similar to that shown for the one-variable version. The details can be confusing and easy to get wrong, so I recommend looking at source code. A simplified version of the source code is available here.