# Runge-Kutta methods for ODE integration in Julia¶

• I want to implement and illustrate the Runge-Kutta method (actually, different variants), in the Julia programming language.

• The Runge-Kutta methods are a family of numerical iterative algorithms to approximate solutions of Ordinary Differential Equations. I will simply implement them, for the mathematical descriptions, I let the interested reader refer to the Wikipedia page, or any good book or course on numerical integration of ODE.

• I will start with the order 1 method, then the order 2 and the most famous order 4.
• They will be compared on different ODE.

## Preliminary¶

In [1]:
versioninfo()

Julia Version 0.6.0
Commit 9036443 (2017-06-19 13:05 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz
WORD_SIZE: 64
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
LAPACK: libopenblas64_
LIBM: libopenlibm
LLVM: libLLVM-3.9.1 (ORCJIT, skylake)


For comparison, let's use this mature and fully featured package DifferentialEquations that provides a solve function to numerically integrate ordinary different equations, and the Plots package with PyPlot backend for plotting:

In [2]:
# If needed:

In [3]:
using Plots
gr()
using DifferentialEquations


I will use as a first example the one included in the scipy (Python) documentation for this odeint function.

$$\theta''(t) + b \theta'(t) + c \sin(\theta(t)) = 0.$$

If $\omega(t) := \theta'(t)$, this gives $$\begin{cases} \theta'(t) = \omega(t) \\ \omega'(t) = -b \omega(t) - c \sin(\theta(t)) \end{cases}$$

Vectorially, if $y(t) = [\theta(t), \omega(t)]$, then the equation is $y' = f(t, y)$ where $f(t, y) = [y_2(t), -b y_2(t) - c \sin(y_1(t))]$.

We assume the values of $b$ and $c$ to be known, and the starting point to be also fixed:

In [4]:
b = 0.25
c = 5.0

Out[4]:
5.0
In [5]:
y0 = [pi - 0.1; 0.0]

Out[5]:
2-element Array{Float64,1}:
3.04159
0.0    
In [6]:
function pend(t, y, dy)
dy[1] = y[2]
dy[2] = (-b * y[2]) - (c * sin(y[1]))
end

function f_pend(y, t)
return [y[2], (-b * y[2]) - (c * sin(y[1]))]
end

Out[6]:
f_pend (generic function with 1 method)

The solve function from DifferentialEquations will be used to solve this ODE on the interval $t \in [0, 10]$.

In [7]:
tspan = (0.0, 10.0)

Out[7]:
(0.0, 10.0)

It is used like this, and our implementations will follow this signature.

In [8]:
function odeint_1(f, y0, tspan)
prob = ODEProblem(f, y0, tspan)
sol = solve(prob)
return sol.t, hcat(sol.u...)'
end

Out[8]:
odeint_1 (generic function with 1 method)
In [9]:
function odeint(f, y0, tspan)
t, sol = odeint_1(f, y0, tspan)
return sol
end

Out[9]:
odeint (generic function with 1 method)
In [10]:
t, sol = odeint_1(pend, y0, tspan)

Out[10]:
([0.0, 0.00200268, 0.0220295, 0.0813368, 0.17913, 0.306465, 0.472717, 0.676887, 0.920686, 1.19691  …  6.29546, 6.72503, 7.20537, 7.57654, 8.0385, 8.41262, 8.86583, 9.2365, 9.64917, 10.0], [3.04159 0.0; 3.04159 -0.000999424; … ; -0.500599 1.27096; 0.0236952 1.5641])
In [11]:
plot(t, sol[:, 1], xaxis="Time t", title="Solution to the pendulum ODE", label="\\theta (t)")
plot!(t, sol[:, 2], label="\\omega (t)")

Out[11]:

## Runge-Kutta method of order 1, or the Euler method¶

The approximation is computed using this update: $$y_{n+1} = y_n + (t_{n+1} - t_n) f(y_n, t_n).$$

The math behind this formula are the following: if $g$ is a solution to the ODE, and so far the approximation is correct, $y_n \simeq g(t_n)$, then a small step $h = t_{n+1} - t_n$ satisfy $g(t_n + h) \simeq g(t_n) + h g'(t_n) \simeq y_n + h f(g(t_n), t_n) + \simeq y_n + h f(y_n, t_n)$.

In [12]:
function rungekutta1(f, y0, t)
n = length(t)
y = zeros((n, length(y0)))
y[1,:] = y0
for i in 1:n-1
h = t[i+1] - t[i]
y[i+1,:] = y[i,:] + h * f(y[i,:], t[i])
end
return y
end

Out[12]:
rungekutta1 (generic function with 1 method)
In [14]:
t = linspace(0, 10, 101);
sol = rungekutta1(f_pend, y0, t);

In [15]:
plot(t, sol[:, 1], xaxis="Time t", title="Solution to the pendulum ODE with Runge-Kutta 1", label="\\theta (t)")
plot!(t, sol[:, 2], label="\\omega (t)")

Out[15]:

With the same number of points, the Euler method (i.e. the Runge-Kutta method of order 1) is less precise than the reference solve method. With more points, it can give a satisfactory approximation of the solution:

In [16]:
t2 = linspace(0, 10, 1001);
sol2 = rungekutta1(f_pend, y0, t2);

In [17]:
plot(t2, sol2[:, 1], xaxis="Time t", title="Solution to the pendulum ODE with Runge-Kutta 1", label="\\theta (t)")
plot!(t2, sol2[:, 2], label="\\omega (t)")

Out[17]:
In [18]:
t3 = linspace(0, 10, 2001);
sol3 = rungekutta1(f_pend, y0, t3);

In [19]:
plot(t3, sol3[:, 1], xaxis="Time t", title="Solution to the pendulum ODE with Runge-Kutta 1", label="\\theta (t)")
plot!(t3, sol3[:, 2], label="\\omega (t)")

Out[19]:

## Runge-Kutta method of order 2¶

The order 2 Runge-Method uses this update: $$y_{n+1} = y_n + h f(t + \frac{h}{2}, y_n + \frac{h}{2} f(t, y_n)),$$ if $h = t_{n+1} - t_n$.

In [20]:
function rungekutta2(f, y0, t)
n = length(t)
y = zeros((n, length(y0)))
y[1,:] = y0
for i in 1:n-1
h = t[i+1] - t[i]
y[i+1,:] = y[i,:] + h * f(y[i,:] + f(y[i,:], t[i]) * h/2, t[i] + h/2)
end
return y
end

Out[20]:
rungekutta2 (generic function with 1 method)

For our simple ODE example, this method is already quite efficient.

In [21]:
t3 = linspace(0, 10, 21);
sol3 = rungekutta2(f_pend, y0, t3);

In [22]:
plot(t3, sol3[:, 1], xaxis="Time t", title="Solution to the pendulum ODE with Runge-Kutta 2 (21 points)", label="\\theta (t)")
plot!(t3, sol3[:, 2], label="\\omega (t)")

Out[22]: