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.
versioninfo()
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:
# If needed:
#Pkg.add("DifferentialEquations")
#Pkg.add("PyPlot")
#Pkg.add("Plots")
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:
b = 0.25
c = 5.0
y0 = [pi - 0.1; 0.0]
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
The solve function from DifferentialEquations will be used to solve this ODE on the interval $t \in [0, 10]$.
tspan = (0.0, 10.0)
It is used like this, and our implementations will follow this signature.
function odeint_1(f, y0, tspan)
prob = ODEProblem(f, y0, tspan)
sol = solve(prob)
return sol.t, hcat(sol.u...)'
end
function odeint(f, y0, tspan)
t, sol = odeint_1(f, y0, tspan)
return sol
end
t, sol = odeint_1(pend, y0, tspan)
plot(t, sol[:, 1], xaxis="Time t", title="Solution to the pendulum ODE", label="\\theta (t)")
plot!(t, sol[:, 2], label="\\omega (t)")
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)$.
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
t = linspace(0, 10, 101);
sol = rungekutta1(f_pend, y0, t);
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)")
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:
t2 = linspace(0, 10, 1001);
sol2 = rungekutta1(f_pend, y0, t2);
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)")
t3 = linspace(0, 10, 2001);
sol3 = rungekutta1(f_pend, y0, t3);
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)")
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$.
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
For our simple ODE example, this method is already quite efficient.
t3 = linspace(0, 10, 21);
sol3 = rungekutta2(f_pend, y0, t3);
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)")
t3 = linspace(0, 10, 101);
sol3 = rungekutta2(f_pend, y0, t3);
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)")
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}$$
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
For our simple ODE example, this method is even more efficient.
t = linspace(0, 10, 31);
sol = rungekutta4(f_pend, y0, t);
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)")
t = linspace(0, 10, 101);
sol = rungekutta4(f_pend, y0, t);
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)")
methods = [rungekutta1, rungekutta2, rungekutta4]
markers = [:o, :s, :>]
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
test_1(10)
test_1(20)
test_1(100)
test_1(200)
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:
y0 = [0; 1; 0.1]
function f(y, t)
return [y[2], y[3], 12 * y[1]^(4/5) + cos(y[2])^3 - sin(y[3])]
end
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
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
test_2(10)
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:
y0 = [10.0, -3.0, 1.0, 1.0]
function f(y, t)
return [y[2], y[3], y[4], y[1]^(-5/3)]
end
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
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
test_3(10)
test_3(50)
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).
We can also compare their speed:
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
Using BenchmarkTools.jl is also interesting as it is more precise than the builtin @time benchmark macro.
using BenchmarkTools, Compat
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
DifferentialEquations implementation is much faster than our manual implentations!That's it for today, folks! See my other notebooks, available on GitHub.