{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": "true"
},
"source": [
"# Table of Contents\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Runge-Kutta methods for ODE integration in Julia\n",
"\n",
"- I want to implement and illustrate the [Runge-Kutta method](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (actually, different variants), in the [Julia programming language](https://www.julialang.org/).\n",
"\n",
"- The Runge-Kutta methods are a family of numerical iterative algorithms to approximate solutions of [Ordinary Differential Equations](https://en.wikipedia.org/wiki/Ordinary_differential_equation). I will simply implement them, for the mathematical descriptions, I let the interested reader refer to the Wikipedia page, or [any](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#References) [good](https://www.directtextbook.com/isbn/9780521007948) [book](https://www.decitre.fr/livres/analyse-numerique-et-equations-differentielles-9782868838919.html) or [course](https://courses.maths.ox.ac.uk/node/4294) on numerical integration of ODE.\n",
"- I will start with the order 1 method, then the order 2 and the most famous order 4.\n",
"- They will be compared on different ODE."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Preliminary"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia Version 0.6.0\n",
"Commit 9036443 (2017-06-19 13:05 UTC)\n",
"Platform Info:\n",
" OS: Linux (x86_64-pc-linux-gnu)\n",
" CPU: Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz\n",
" WORD_SIZE: 64\n",
" BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)\n",
" LAPACK: libopenblas64_\n",
" LIBM: libopenlibm\n",
" LLVM: libLLVM-3.9.1 (ORCJIT, skylake)\n"
]
}
],
"source": [
"versioninfo()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For comparison, let's use this [mature and fully featured package `DifferentialEquations`](http://docs.juliadiffeq.org/latest/) that provides a `solve` function to numerically integrate ordinary different equations, and the `Plots` package with `PyPlot` backend for plotting:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# If needed:\n",
"#Pkg.add(\"DifferentialEquations\")\n",
"#Pkg.add(\"PyPlot\")\n",
"#Pkg.add(\"Plots\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"gr()\n",
"using DifferentialEquations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I will use as a first example the one included in [the scipy (Python) documentation for this `odeint` function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html).\n",
"\n",
"$$\\theta''(t) + b \\theta'(t) + c \\sin(\\theta(t)) = 0.$$\n",
"\n",
"If $\\omega(t) := \\theta'(t)$, this gives\n",
"$$ \\begin{cases}\n",
"\\theta'(t) = \\omega(t) \\\\\n",
"\\omega'(t) = -b \\omega(t) - c \\sin(\\theta(t))\n",
"\\end{cases} $$\n",
"\n",
"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))]$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We assume the values of $b$ and $c$ to be known, and the starting point to be also fixed:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
"5.0"
],
"text/plain": [
"5.0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b = 0.25\n",
"c = 5.0"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 3.04159\n",
" 0.0 "
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y0 = [pi - 0.1; 0.0]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"f_pend (generic function with 1 method)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function pend(t, y, dy)\n",
" dy[1] = y[2]\n",
" dy[2] = (-b * y[2]) - (c * sin(y[1]))\n",
"end\n",
"\n",
"function f_pend(y, t)\n",
" return [y[2], (-b * y[2]) - (c * sin(y[1]))]\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `solve` function from `DifferentialEquations` will be used to solve this ODE on the interval $t \\in [0, 10]$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0, 10.0)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tspan = (0.0, 10.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is used like this, and our implementations will follow this signature."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"odeint_1 (generic function with 1 method)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function odeint_1(f, y0, tspan)\n",
" prob = ODEProblem(f, y0, tspan)\n",
" sol = solve(prob)\n",
" return sol.t, hcat(sol.u...)'\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"odeint (generic function with 1 method)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function odeint(f, y0, tspan)\n",
" t, sol = odeint_1(f, y0, tspan)\n",
" return sol\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"([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])"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t, sol = odeint_1(pend, y0, tspan)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t, sol[:, 1], xaxis=\"Time t\", title=\"Solution to the pendulum ODE\", label=\"\\\\theta (t)\")\n",
"plot!(t, sol[:, 2], label=\"\\\\omega (t)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----\n",
"## Runge-Kutta method of order 1, or the Euler method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The approximation is computed using this update:\n",
"$$y_{n+1} = y_n + (t_{n+1} - t_n) f(y_n, t_n).$$\n",
"\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)$."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"rungekutta1 (generic function with 1 method)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function rungekutta1(f, y0, t)\n",
" n = length(t)\n",
" y = zeros((n, length(y0)))\n",
" y[1,:] = y0\n",
" for i in 1:n-1\n",
" h = t[i+1] - t[i]\n",
" y[i+1,:] = y[i,:] + h * f(y[i,:], t[i])\n",
" end\n",
" return y\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"t = linspace(0, 10, 101);\n",
"sol = rungekutta1(f_pend, y0, t);"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t, sol[:, 1], xaxis=\"Time t\", title=\"Solution to the pendulum ODE with Runge-Kutta 1\", label=\"\\\\theta (t)\")\n",
"plot!(t, sol[:, 2], label=\"\\\\omega (t)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"t2 = linspace(0, 10, 1001);\n",
"sol2 = rungekutta1(f_pend, y0, t2);"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t2, sol2[:, 1], xaxis=\"Time t\", title=\"Solution to the pendulum ODE with Runge-Kutta 1\", label=\"\\\\theta (t)\")\n",
"plot!(t2, sol2[:, 2], label=\"\\\\omega (t)\")"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"t3 = linspace(0, 10, 2001);\n",
"sol3 = rungekutta1(f_pend, y0, t3);"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t3, sol3[:, 1], xaxis=\"Time t\", title=\"Solution to the pendulum ODE with Runge-Kutta 1\", label=\"\\\\theta (t)\")\n",
"plot!(t3, sol3[:, 2], label=\"\\\\omega (t)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----\n",
"## Runge-Kutta method of order 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The order 2 Runge-Method uses this update:\n",
"$$ y_{n+1} = y_n + h f(t + \\frac{h}{2}, y_n + \\frac{h}{2} f(t, y_n)),$$\n",
"if $h = t_{n+1} - t_n$."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"rungekutta2 (generic function with 1 method)"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function rungekutta2(f, y0, t)\n",
" n = length(t)\n",
" y = zeros((n, length(y0)))\n",
" y[1,:] = y0\n",
" for i in 1:n-1\n",
" h = t[i+1] - t[i]\n",
" y[i+1,:] = y[i,:] + h * f(y[i,:] + f(y[i,:], t[i]) * h/2, t[i] + h/2)\n",
" end\n",
" return y\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For our simple ODE example, this method is already quite efficient."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"t3 = linspace(0, 10, 21);\n",
"sol3 = rungekutta2(f_pend, y0, t3);"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t3, sol3[:, 1], xaxis=\"Time t\", title=\"Solution to the pendulum ODE with Runge-Kutta 2 (21 points)\", label=\"\\\\theta (t)\")\n",
"plot!(t3, sol3[:, 2], label=\"\\\\omega (t)\")"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"t3 = linspace(0, 10, 101);\n",
"sol3 = rungekutta2(f_pend, y0, t3);"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t3, sol3[:, 1], xaxis=\"Time t\", title=\"Solution to the pendulum ODE with Runge-Kutta 2 (101 points)\", label=\"\\\\theta (t)\")\n",
"plot!(t3, sol3[:, 2], label=\"\\\\omega (t)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----\n",
"## Runge-Kutta method of order 4, *\"RK4\"*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The order 4 Runge-Method uses this update:\n",
"$$ y_{n+1} = y_n + \\frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4),$$\n",
"if $h = t_{n+1} - t_n$, and\n",
"$$\\begin{cases}\n",
"k_1 &= f(y_n, t_n), \\\\\n",
"k_2 &= f(y_n + \\frac{h}{2} k_1, t_n + \\frac{h}{2}), \\\\\n",
"k_3 &= f(y_n + \\frac{h}{2} k_2, t_n + \\frac{h}{2}), \\\\\n",
"k_4 &= f(y_n + h k_3, t_n + h).\n",
"\\end{cases}$$"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"rungekutta4 (generic function with 1 method)"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function rungekutta4(f, y0, t)\n",
" n = length(t)\n",
" y = zeros((n, length(y0)))\n",
" y[1,:] = y0\n",
" for i in 1:n-1\n",
" h = t[i+1] - t[i]\n",
" k1 = f(y[i,:], t[i])\n",
" k2 = f(y[i,:] + k1 * h/2, t[i] + h/2)\n",
" k3 = f(y[i,:] + k2 * h/2, t[i] + h/2)\n",
" k4 = f(y[i,:] + k3 * h, t[i] + h)\n",
" y[i+1,:] = y[i,:] + (h/6) * (k1 + 2*k2 + 2*k3 + k4)\n",
" end\n",
" return y\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For our simple ODE example, this method is even more efficient."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"t = linspace(0, 10, 31);\n",
"sol = rungekutta4(f_pend, y0, t);"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t, sol[:, 1], xaxis=\"Time t\", title=\"Solution to the pendulum ODE with Runge-Kutta 4 (31 points)\", label=\"\\\\theta (t)\")\n",
"plot!(t, sol[:, 2], label=\"\\\\omega (t)\")"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"t = linspace(0, 10, 101);\n",
"sol = rungekutta4(f_pend, y0, t);"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t, sol[:, 1], xaxis=\"Time t\", title=\"Solution to the pendulum ODE with Runge-Kutta 4 (101 points)\", label=\"\\\\theta (t)\")\n",
"plot!(t, sol[:, 2], label=\"\\\\omega (t)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----\n",
"## Comparisons"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Symbol,1}:\n",
" :o\n",
" :s\n",
" :>"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"methods = [rungekutta1, rungekutta2, rungekutta4]\n",
"markers = [:o, :s, :>]"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"test_1 (generic function with 2 methods)"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function test_1(n=101)\n",
" t = linspace(0, 10, n)\n",
" tspan = (0.0, 10.0)\n",
" t1, sol = odeint_1(pend, y0, tspan)\n",
" plt = plot(t1, sol[:, 1], marker=:d, xaxis=\"Time t\", title=\"Solution to the pendulum ODE with ($n points)\", label=\"odeint\")\n",
" for (method, m) in zip(methods, markers)\n",
" sol = method(f_pend, y0, t)\n",
" plot!(t, sol[:, 1], marker=m, label=string(method))\n",
" end\n",
" display(plt)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"test_1(10)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"test_1(20)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"test_1(100)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"test_1(200)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparisons on another integration problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the following ODE on $t\\in[0, 1]$:\n",
"$$\n",
"\\begin{cases}\n",
" y'''(t) = 12 y(t)^{4/5} + \\cos(y'(t))^3 - \\sin(y''(t)) \\\\\n",
" y(0) = 0, y'(0) = 1, y''(0) = 0.1\n",
"\\end{cases}\n",
"$$\n",
"\n",
"It can be written in a vectorial form like the first one:"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 0.0\n",
" 1.0\n",
" 0.1"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y0 = [0; 1; 0.1]"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"f (generic function with 1 method)"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function f(y, t)\n",
" return [y[2], y[3], 12 * y[1]^(4/5) + cos(y[2])^3 - sin(y[3])]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"f_2 (generic function with 1 method)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function f_2(t, y, dy)\n",
" dy[1] = y[2]\n",
" dy[2] = y[3]\n",
" dy[3] = 12 * y[1]^(4/5) + cos(y[2])^3 - sin(y[3])\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"test_2 (generic function with 2 methods)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function test_2(n=101)\n",
" t = linspace(0, 1, n)\n",
" tspan = (0.0, 1.0)\n",
" t1, sol = odeint_1(f_2, y0, tspan)\n",
" plt = plot(t1, sol[:, 1], marker=:d, xaxis=\"Time t\", title=\"Solution to an ODE with ($n points)\", label=\"odeint\")\n",
" for (method, m) in zip(methods, markers)\n",
" sol = method(f, y0, t)\n",
" plot!(t, sol[:, 1], marker=m, label=string(method))\n",
" end\n",
" display(plt)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"test_2(10)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"test_2(50)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the following ODE on $t\\in[0, 3]$:\n",
"$$\n",
"\\begin{cases}\n",
" y''''(t) = y(t)^{-5/3} \\\\\n",
" y(0) = 10, y'(0) = -3, y''(0) = 1, y'''(0) = 1\n",
"\\end{cases}\n",
"$$\n",
"\n",
"It can be written in a vectorial form like the first one:"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Float64,1}:\n",
" 10.0\n",
" -3.0\n",
" 1.0\n",
" 1.0"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y0 = [10.0, -3.0, 1.0, 1.0]"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"code_folding": []
},
"outputs": [
{
"data": {
"text/plain": [
"f (generic function with 1 method)"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function f(y, t)\n",
" return [y[2], y[3], y[4], y[1]^(-5/3)]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"code_folding": []
},
"outputs": [
{
"data": {
"text/plain": [
"f_2 (generic function with 1 method)"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function f_2(t, y, dy)\n",
" dy[1] = y[2]\n",
" dy[2] = y[3]\n",
" dy[3] = y[4]\n",
" dy[4] = y[1]^(-5/3)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"test_3 (generic function with 2 methods)"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function test_3(n=101)\n",
" t = linspace(0, 1, n)\n",
" tspan = (0.0, 1.0)\n",
" t1, sol = odeint_1(f_2, y0, tspan)\n",
" plt = plot(t1, sol[:, 1], marker=:d, xaxis=\"Time t\", title=\"Solution to an ODE with ($n points)\", label=\"odeint\")\n",
" for (method, m) in zip(methods, markers)\n",
" sol = method(f, y0, t)\n",
" plot!(t, sol[:, 1], marker=m, label=string(method))\n",
" end\n",
" display(plt)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"test_3(10)"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: both OrdinaryDiffEq and StochasticDiffEq export \"NLSOLVEJL_SETUP\"; uses of it in module DifferentialEquations must be qualified\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"search: \u001b[1me\u001b[22m\u001b[1mn\u001b[22m\u001b[1md\u001b[22m \u001b[1me\u001b[22m\u001b[1mn\u001b[22m\u001b[1md\u001b[22mof \u001b[1me\u001b[22m\u001b[1mn\u001b[22m\u001b[1md\u001b[22mswith \u001b[1mE\u001b[22m\u001b[1mN\u001b[22m\u001b[1mD\u001b[22mIAN_BOM s\u001b[1me\u001b[22m\u001b[1mn\u001b[22m\u001b[1md\u001b[22m app\u001b[1me\u001b[22m\u001b[1mn\u001b[22m\u001b[1md\u001b[22m! back\u001b[1me\u001b[22m\u001b[1mn\u001b[22m\u001b[1md\u001b[22m back\u001b[1me\u001b[22m\u001b[1mn\u001b[22m\u001b[1md\u001b[22ms back\u001b[1me\u001b[22m\u001b[1mn\u001b[22m\u001b[1md\u001b[22m_name\n",
"\n"
]
}
],
"source": [
"test_3(50)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Small benchmark\n",
"We can also compare their speed:"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"For 20 points...\n",
"Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\n",
" 0.000516 seconds (634 allocations: 44.906 KiB)\n",
"Time of solving an ODE with the rungekutta1 method for 20 points ...\n",
" 0.000022 seconds (231 allocations: 13.266 KiB)\n",
"Time of solving an ODE with the rungekutta2 method for 20 points ...\n",
" 0.000032 seconds (421 allocations: 25.141 KiB)\n",
"Time of solving an ODE with the rungekutta4 method for 20 points ...\n",
" 0.000066 seconds (877 allocations: 57.203 KiB)\n",
"\n",
"For 100 points...\n",
"Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\n",
" 0.000435 seconds (634 allocations: 44.906 KiB)\n",
"Time of solving an ODE with the rungekutta1 method for 100 points ...\n",
" 0.000090 seconds (1.19 k allocations: 68.297 KiB)\n",
"Time of solving an ODE with the rungekutta2 method for 100 points ...\n",
" 0.000139 seconds (2.18 k allocations: 130.172 KiB)\n",
"Time of solving an ODE with the rungekutta4 method for 100 points ...\n",
" 0.000317 seconds (4.56 k allocations: 297.234 KiB)\n",
"\n",
"For 1000 points...\n",
"Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\n",
" 0.000468 seconds (634 allocations: 44.906 KiB)\n",
"Time of solving an ODE with the rungekutta1 method for 1000 points ...\n",
" 0.000748 seconds (13.46 k allocations: 709.891 KiB)\n",
"Time of solving an ODE with the rungekutta2 method for 1000 points ...\n",
" 0.001345 seconds (23.93 k allocations: 1.310 MiB)\n",
"Time of solving an ODE with the rungekutta4 method for 1000 points ...\n",
" 0.012061 seconds (48.89 k allocations: 2.972 MiB, 79.27% gc time)\n"
]
}
],
"source": [
"function benchmark(n=101)\n",
" t = linspace(0, 1, n)\n",
" tspan = (0.0, 1.0)\n",
"\n",
" print(\"Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\\n\")\n",
" @time t1, sol = odeint_1(f_2, y0, tspan)\n",
" for method in methods\n",
" print(\"Time of solving an ODE with the $method method for $n points ...\\n\")\n",
" @time sol = method(f, y0, t)\n",
" end\n",
"end\n",
"\n",
"for n in [20, 100, 1000]\n",
" print(\"\\nFor $n points...\\n\")\n",
" benchmark(n)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using [`BenchmarkTools.jl`](https://github.com/JuliaCI/BenchmarkTools.jl) is also interesting as it is more precise than the builtin `@time` benchmark macro."
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [],
"source": [
"using BenchmarkTools, Compat"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"For 20 points...\n",
"Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\n",
" 164.017 μs (635 allocations: 44.94 KiB)\n",
"Time of solving an ODE with the rungekutta1 method for 20 points ...\n",
" 5.765 μs (232 allocations: 13.33 KiB)\n",
"Time of solving an ODE with the rungekutta2 method for 20 points ...\n",
" 10.680 μs (422 allocations: 25.20 KiB)\n",
"Time of solving an ODE with the rungekutta4 method for 20 points ...\n",
" 23.242 μs (878 allocations: 57.27 KiB)\n",
"\n",
"For 100 points...\n",
"Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\n",
" 164.906 μs (635 allocations: 44.94 KiB)\n",
"Time of solving an ODE with the rungekutta1 method for 100 points ...\n",
" 29.808 μs (1192 allocations: 68.36 KiB)\n",
"Time of solving an ODE with the rungekutta2 method for 100 points ...\n",
" 57.445 μs (2182 allocations: 130.23 KiB)\n",
"Time of solving an ODE with the rungekutta4 method for 100 points ...\n",
" 117.637 μs (4558 allocations: 297.30 KiB)\n",
"\n",
"For 1000 points...\n",
"Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\n",
" 165.346 μs (635 allocations: 44.94 KiB)\n",
"Time of solving an ODE with the rungekutta1 method for 1000 points ...\n",
" 296.730 μs (13458 allocations: 709.95 KiB)\n",
"Time of solving an ODE with the rungekutta2 method for 1000 points ...\n",
" 580.526 μs (23936 allocations: 1.31 MiB)\n",
"Time of solving an ODE with the rungekutta4 method for 1000 points ...\n",
" 1.233 ms (48888 allocations: 2.97 MiB)\n"
]
}
],
"source": [
"function benchmark(n=101)\n",
" t = linspace(0, 1, n)\n",
" tspan = (0.0, 1.0)\n",
"\n",
" print(\"Time of solving an ODE with the 'solve' method from 'DifferentialEquations' ...\\n\")\n",
" @btime t1, sol = $odeint_1($f_2, $y0, $tspan)\n",
" for method in methods\n",
" print(\"Time of solving an ODE with the $method method for $n points ...\\n\")\n",
" @btime sol = $method($f, $y0, $t)\n",
" end\n",
"end\n",
"\n",
"for n in [20, 100, 1000]\n",
" print(\"\\nFor $n points...\\n\")\n",
" benchmark(n)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- The order 1 method is simpler and so faster than the order 2, which itself is simpler and faster than the order 4 method.\n",
"- And we can check that the `DifferentialEquations` implementation is much faster than our manual implentations!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusion\n",
"\n",
"> *That's it for today, folks!* See my other notebooks, [available on GitHub](https://github.com/Naereen/notebooks/)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.0",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.0"
},
"toc": {
"colors": {
"hover_highlight": "#DAA520",
"running_highlight": "#FF0000",
"selected_highlight": "#FFD700"
},
"moveMenuLeft": true,
"nav_menu": {
"height": "193px",
"width": "251px"
},
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 4,
"toc_cell": true,
"toc_section_display": "block",
"toc_window_display": true
}
},
"nbformat": 4,
"nbformat_minor": 2
}