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.
In [1]:
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:
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.
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]:
In [5]:
y0 = [pi - 0.1; 0.0]
Out[5]:
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]:
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]:
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]:
In [9]:
function odeint(f, y0, tspan)
t, sol = odeint_1(f, y0, tspan)
return sol
end
Out[9]:
In [10]:
t, sol = odeint_1(pend, y0, tspan)
Out[10]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
In [30]:
methods = [rungekutta1, rungekutta2, rungekutta4]
markers = [:o, :s, :>]
Out[30]:
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]:
In [32]:
test_1(10)
In [33]:
test_1(20)
In [34]:
test_1(100)
In [35]:
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:
In [36]:
y0 = [0; 1; 0.1]
Out[36]:
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]:
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]:
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]:
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]:
In [48]:
function f(y, t)
return [y[2], y[3], y[4], y[1]^(-5/3)]
end
Out[48]:
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]:
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]:
In [55]:
test_3(10)
In [56]:
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).
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
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
DifferentialEquations implementation is much faster than our manual implentations!That's it for today, folks! See my other notebooks, available on GitHub.