Table of Contents¶

1 Runge-Kutta methods for ODE integration in Julia

1.1 Preliminary

1.2 Runge-Kutta method of order 1, or the Euler method

1.3 Runge-Kutta method of order 2

1.4 Runge-Kutta method of order 4, "RK4"

1.5 Comparisons

1.6 Comparisons on another integration problem

1.7 Small benchmark

1.8 Conclusion

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:
#Pkg.add("DifferentialEquations")
#Pkg.add("PyPlot")
#Pkg.add("Plots")

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]:

In [23]:

t3 = linspace(0, 10, 101);
sol3 = rungekutta2(f_pend, y0, t3);

In [24]:

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

Out[24]:

Runge-Kutta method of order 4, "RK4"¶

The order 4 Runge-Method uses this update: $$ y_{n+1} = y_n + \frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4),$$ if $h = t_{n+1} - t_n$, and $$\begin{cases} k_1 &= f(y_n, t_n), \\ k_2 &= f(y_n + \frac{h}{2} k_1, t_n + \frac{h}{2}), \\ k_3 &= f(y_n + \frac{h}{2} k_2, t_n + \frac{h}{2}), \\ k_4 &= f(y_n + h k_3, t_n + h). \end{cases}$$

In [25]:

function rungekutta4(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]
        k1 = f(y[i,:], t[i])
        k2 = f(y[i,:] + k1 * h/2, t[i] + h/2)
        k3 = f(y[i,:] + k2 * h/2, t[i] + h/2)
        k4 = f(y[i,:] + k3 * h, t[i] + h)
        y[i+1,:] = y[i,:] + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
    end
    return y
end

Out[25]:

rungekutta4 (generic function with 1 method)

For our simple ODE example, this method is even more efficient.

In [26]:

t = linspace(0, 10, 31);
sol = rungekutta4(f_pend, y0, t);

In [27]:

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

Out[27]:

In [28]:

t = linspace(0, 10, 101);
sol = rungekutta4(f_pend, y0, t);

In [29]:

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

Out[29]:

Comparisons¶

In [30]:

methods = [rungekutta1, rungekutta2, rungekutta4]
markers = [:o, :s, :>]

Out[30]:

3-element Array{Symbol,1}:
 :o
 :s
 :>

In [31]:

function test_1(n=101)
    t = linspace(0, 10, n)
    tspan = (0.0, 10.0)
    t1, sol = odeint_1(pend, y0, tspan)
    plt = plot(t1, sol[:, 1], marker=:d, xaxis="Time t", title="Solution to the pendulum ODE with ($n points)", label="odeint")
    for (method, m) in zip(methods, markers)
        sol = method(f_pend, y0, t)
        plot!(t, sol[:, 1], marker=m, label=string(method))
    end
    display(plt)
end

Out[31]:

test_1 (generic function with 2 methods)

In [32]:

test_1(10)

In [33]:

test_1(20)

In [34]:

test_1(100)

In [35]:

test_1(200)

Comparisons on another integration problem¶

Consider the following ODE on $t\in[0, 1]$: $$ \begin{cases} y'''(t) = 12 y(t)^{4/5} + \cos(y'(t))^3 - \sin(y''(t)) \\ y(0) = 0, y'(0) = 1, y''(0) = 0.1 \end{cases} $$

It can be written in a vectorial form like the first one:

In [36]:

y0 = [0; 1; 0.1]

Out[36]:

3-element Array{Float64,1}:
 0.0
 1.0
 0.1

In [37]:

function f(y, t)
    return [y[2], y[3], 12 * y[1]^(4/5) + cos(y[2])^3 - sin(y[3])]
end

Out[37]:

f (generic function with 1 method)

In [38]:

function f_2(t, y, dy)
    dy[1] = y[2]
    dy[2] = y[3]
    dy[3] = 12 * y[1]^(4/5) + cos(y[2])^3 - sin(y[3])
end

Out[38]:

f_2 (generic function with 1 method)

In [39]:

function test_2(n=101)
    t = linspace(0, 1, n)
    tspan = (0.0, 1.0)
    t1, sol = odeint_1(f_2, y0, tspan)
    plt = plot(t1, sol[:, 1], marker=:d, xaxis="Time t", title="Solution to an ODE with ($n points)", label="odeint")
    for (method, m) in zip(methods, markers)
        sol = method(f, y0, t)
        plot!(t, sol[:, 1], marker=m, label=string(method))
    end
    display(plt)
end

Out[39]:

test_2 (generic function with 2 methods)

In [40]:

test_2(10)

In [41]:

test_2(50)

Consider the following ODE on $t\in[0, 3]$: $$ \begin{cases} y''''(t) = y(t)^{-5/3} \\ y(0) = 10, y'(0) = -3, y''(0) = 1, y'''(0) = 1 \end{cases} $$

It can be written in a vectorial form like the first one:

In [47]:

y0 = [10.0, -3.0, 1.0, 1.0]

Out[47]:

4-element Array{Float64,1}:
 10.0
 -3.0
  1.0
  1.0

In [48]:

function f(y, t)
    return [y[2], y[3], y[4], y[1]^(-5/3)]
end

Out[48]:

f (generic function with 1 method)

In [53]:

function f_2(t, y, dy)
    dy[1] = y[2]
    dy[2] = y[3]
    dy[3] = y[4]
    dy[4] = y[1]^(-5/3)
end

Out[53]:

f_2 (generic function with 1 method)

In [54]:

function test_3(n=101)
    t = linspace(0, 1, n)
    tspan = (0.0, 1.0)
    t1, sol = odeint_1(f_2, y0, tspan)
    plt = plot(t1, sol[:, 1], marker=:d, xaxis="Time t", title="Solution to an ODE with ($n points)", label="odeint")
    for (method, m) in zip(methods, markers)
        sol = method(f, y0, t)
        plot!(t, sol[:, 1], marker=m, label=string(method))
    end
    display(plt)
end

Out[54]:

test_3 (generic function with 2 methods)

In [55]:

test_3(10)

In [56]:

test_3(50)

WARNING: both OrdinaryDiffEq and StochasticDiffEq export "NLSOLVEJL_SETUP"; uses of it in module DifferentialEquations must be qualified

search: end endof endswith ENDIAN_BOM send append! backend backends backend_name

Our hand-written Runge-Kutta method of order 4 seems to be as efficient as the odeint method from scipy... and that's because odeint basically uses a Runge-Kutta method of order 4 (with smart variants).

Small benchmark¶

We can also compare their speed:

In [59]:

function benchmark(n=101)
    t = linspace(0, 1, n)
    tspan = (0.0, 1.0)

    print("Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\n")
    @time t1, sol = odeint_1(f_2, y0, tspan)
    for method in methods
        print("Time of solving an ODE with the $method method for $n points ...\n")
        @time sol = method(f, y0, t)
    end
end

for n in [20, 100, 1000]
    print("\nFor $n points...\n")
    benchmark(n)
end

For 20 points...
Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...
  0.000516 seconds (634 allocations: 44.906 KiB)
Time of solving an ODE with the rungekutta1 method for 20 points ...
  0.000022 seconds (231 allocations: 13.266 KiB)
Time of solving an ODE with the rungekutta2 method for 20 points ...
  0.000032 seconds (421 allocations: 25.141 KiB)
Time of solving an ODE with the rungekutta4 method for 20 points ...
  0.000066 seconds (877 allocations: 57.203 KiB)

For 100 points...
Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...
  0.000435 seconds (634 allocations: 44.906 KiB)
Time of solving an ODE with the rungekutta1 method for 100 points ...
  0.000090 seconds (1.19 k allocations: 68.297 KiB)
Time of solving an ODE with the rungekutta2 method for 100 points ...
  0.000139 seconds (2.18 k allocations: 130.172 KiB)
Time of solving an ODE with the rungekutta4 method for 100 points ...
  0.000317 seconds (4.56 k allocations: 297.234 KiB)

For 1000 points...
Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...
  0.000468 seconds (634 allocations: 44.906 KiB)
Time of solving an ODE with the rungekutta1 method for 1000 points ...
  0.000748 seconds (13.46 k allocations: 709.891 KiB)
Time of solving an ODE with the rungekutta2 method for 1000 points ...
  0.001345 seconds (23.93 k allocations: 1.310 MiB)
Time of solving an ODE with the rungekutta4 method for 1000 points ...
  0.012061 seconds (48.89 k allocations: 2.972 MiB, 79.27% gc time)

Using BenchmarkTools.jl is also interesting as it is more precise than the builtin @time benchmark macro.

In [61]:

using BenchmarkTools, Compat

In [64]:

function benchmark(n=101)
    t = linspace(0, 1, n)
    tspan = (0.0, 1.0)

    print("Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\n")
    @btime t1, sol = $odeint_1($f_2, $y0, $tspan)
    for method in methods
        print("Time of solving an ODE with the $method method for $n points ...\n")
        @btime sol = $method($f, $y0, $t)
    end
end

for n in [20, 100, 1000]
    print("\nFor $n points...\n")
    benchmark(n)
end

For 20 points...
Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...
  164.017 μs (635 allocations: 44.94 KiB)
Time of solving an ODE with the rungekutta1 method for 20 points ...
  5.765 μs (232 allocations: 13.33 KiB)
Time of solving an ODE with the rungekutta2 method for 20 points ...
  10.680 μs (422 allocations: 25.20 KiB)
Time of solving an ODE with the rungekutta4 method for 20 points ...
  23.242 μs (878 allocations: 57.27 KiB)

For 100 points...
Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...
  164.906 μs (635 allocations: 44.94 KiB)
Time of solving an ODE with the rungekutta1 method for 100 points ...
  29.808 μs (1192 allocations: 68.36 KiB)
Time of solving an ODE with the rungekutta2 method for 100 points ...
  57.445 μs (2182 allocations: 130.23 KiB)
Time of solving an ODE with the rungekutta4 method for 100 points ...
  117.637 μs (4558 allocations: 297.30 KiB)

For 1000 points...
Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...
  165.346 μs (635 allocations: 44.94 KiB)
Time of solving an ODE with the rungekutta1 method for 1000 points ...
  296.730 μs (13458 allocations: 709.95 KiB)
Time of solving an ODE with the rungekutta2 method for 1000 points ...
  580.526 μs (23936 allocations: 1.31 MiB)
Time of solving an ODE with the rungekutta4 method for 1000 points ...
  1.233 ms (48888 allocations: 2.97 MiB)

The order 1 method is simpler and so faster than the order 2, which itself is simpler and faster than the order 4 method.
And we can check that the DifferentialEquations implementation is much faster than our manual implentations!

Conclusion¶

That's it for today, folks! See my other notebooks, available on GitHub.