{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Runge-Kutta methods for ODE integration in OCaml\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 [OCaml programming language](https://www.ocaml.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": [ "The OCaml toplevel, version 4.04.2\n" ] }, { "data": { "text/plain": [ "- : int = 0\n" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Sys.command \"ocaml -version\";;" ] }, { "cell_type": "code", "execution_count": 147, "metadata": {}, "outputs": [], "source": [ "#thread ;;\n", "#require \"jupyter.notebook\" ;;\n", "#require \"jupyter.archimedes\" ;;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I don't want to try to find a reference implementation of an Ordinary Differential Equations solver in OCaml. I'm sure there is, just don't care." ] }, { "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": 11, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "val b : float = 0.25\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "val c : float = 5.\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "let b = 0.25 ;;\n", "let c = 5.0 ;;" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "val y0 : float array = [|3.04156; 0.|]\n" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "let y0 = [| 3.14156 -. 0.1; 0.0 |];;" ] }, { "cell_type": "code", "execution_count": 148, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "val f_pend : float array -> float -> float array =