{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 10 - Initial value problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "nitial value problems of ordinary differential equations, explicit and implicit Euler method, Runge-Kutta methods, Leap-frog, Adams-Bashford, Adams-Moulton, predictor-corrector, Bulirsch-Stoer algorithm, stiff equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate, optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Runge-Kutta methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [Runge–Kutta methods](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) are a family of [explicit and implicit](https://en.wikipedia.org/wiki/Explicit_and_implicit_methods) [iterative](https://en.wikipedia.org/wiki/Iterative_method) methods for the numerical solution of the [initial value problem](https://en.wikipedia.org/wiki/Initial_value_problem)\n",
    "\n",
    "$$\n",
    "\\frac{dy}{dx} = f(x, y), \\quad y(x_0) = y_0.\n",
    "$$\n",
    "\n",
    "Here $ y $ is an unknown function (scalar or vector) of $ x $. The function $ f $ and the [initial conditions](https://en.wikipedia.org/wiki/Initial_condition) $ x_0 $, $ y_0 $ are given. The Runge–Kutta methods take the form\n",
    "\n",
    "$$\n",
    "y_{i + 1} = y_{i} + h \\sum_{j = 1}^{r} b_{j} k_{j}(h), \\quad i = 0, \\dots, n - 1.\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "k_j = f(x_i + c_j h, y_i + h \\sum_{l = 1}^{j - 1} a_{jl} k_l)\n",
    "$$\n",
    "\n",
    "in the case of explicit methods and\n",
    "\n",
    "$$\n",
    "k_j = f(x_i + c_j h, y_i + h \\sum_{l = 1}^{r}  a_{jl} k_l)\n",
    "$$\n",
    "\n",
    "in the case of implicit methods. Here, $ y_i $ and $ y_{i + 1} $ are the approximations of $ y(x_i)$ and $ y(x_{i + 1}) $, respectively, and $ h > 0 $ is a step size. To specify a particular method, one needs to provide the integer $ r $ (the number of stages), and the coefficients $ a_{ij} $, $ b_i $, and $ c_i $ (see the [list of Runge-Kutta methods](https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods)). The matrix $ (a_{ij}) $ is called the Runge–Kutta matrix, while the $ b_i $ and $ c_i $ are known as the weights and the nodes. These data can be arranged in a Butcher tableau:\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "{c|cccc}\n",
    "c_1 & a_{11} & a_{12} & \\cdots & a_{1r} \\\\\n",
    "c_2 & a_{21} & a_{22} & \\cdots & a_{2r} \\\\\n",
    "\\vdots & \\vdots & \\vdots  & \\ddots & \\vdots \\\\\n",
    "c_r & a_{r1} & a_{r2} & \\cdots & a_{rr}\\\\\n",
    "\\hline\n",
    "& b_1 & b_2 & \\cdots & b_{r}\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "**Note:** In the case of explicit methods the Runge-Kutta matrix is lower triangular. \n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forward (explicit) Euler method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [forward Euler method](https://en.wikipedia.org/wiki/Euler_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.1:** Implement the forward Euler method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_euler(f, x_0, x_n, y_0, n): # Runge-Kutta 1st order method\n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(forward_euler(f, 0, 10, 2, 10000)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.2:** The fundamental law of radioactive decay can be mathematically expressed as\n",
    "       \n",
    "$$\n",
    "\\frac{dN(t)}{dt} = - \\lambda N(t), \\qquad N(t_0) = N_0,\n",
    "$$\n",
    "\n",
    "where $ N $ is the number of radioactive nuclei, $ t $ is time, $ \\lambda = \\ln(2) \\, / \\, t_{1/2} $ is the decay constant and $ t_{1/2} $ is the half-life of the decaying quantity. The analytical solution of the decay equation is\n",
    "    \n",
    "$$\n",
    "N(t) = N_0 e^{-\\lambda t}.\n",
    "$$\n",
    "    \n",
    "Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the forward Euler method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Runge-Kutta $2^{\\mathrm{nd}}$ order method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.3:** Implement the Runge-Kutta $ 2^{\\mathrm{nd}} $ order method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def runge_kutta_2(f, x_0, x_n, y_0, n):\n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(runge_kutta_2(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.4:** Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the Runge-Kutta $ 2^{\\mathrm{nd}} $ order method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Runge-Kutta $3^{rd}$ order method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.5:** Implement the Runge-Kutta $ 3^{\\mathrm{rd}} $ order method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def runge_kutta_3(f, x_0, x_n, y_0, n):\n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(runge_kutta_3(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.6:** Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the Runge-Kutta $ 3^{\\mathrm{rd}} $ order method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Runge-Kutta $4^{th}$ order method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.7:** Implement the Runge-Kutta $ 4^{\\mathrm{th}} $ order method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def runge_kutta_4(f, x_0, x_n, y_0, n):\n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(runge_kutta_4(f, 0, 10, 2, 10)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.8:** Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the Runge-Kutta $ 4^{\\mathrm{th}} $ order method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Backward (implicit) Euler method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.9:** Implement the backward Euler method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_euler(f, x_0, x_n, y_0, n): # implicit method\n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(backward_euler(f, 0, 10, 2, 1000)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.10:** Solve the equation\n",
    "    \n",
    "$$\n",
    "\\frac{dx}{dy} = - 100 y + 100, \\quad y(0) = 5\n",
    "$$\n",
    "    \n",
    "using the backward Euler method. Compare the numerical solution with the analytical solution \n",
    "\n",
    "$$\n",
    "y(x) = (y_0 - 1) e^{-100x} + 1\n",
    "$$\n",
    "    \n",
    "and with the solution obtained by the forward Euler method.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multistep methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Leap-frog method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.11:** Implement the Leap-frog method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def leap_frog(f, x_0, x_n, y_0, n): \n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(leap_frog(f, 0, 10, 2, 10000)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.12:** Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the Leap-frog method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adams-Bashforth method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [Adams-Bashforth method](https://en.wikipedia.org/wiki/Linear_multistep_method#Adams%E2%80%93Bashforth_methods)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.13:** Implement the Adams-Bashforth method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def adams_bashforth(f, x_0, x_n, y_0, n): \n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(adams_bashforth(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.14:** Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the Adams-Bashforth method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adams-Moulton method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [Adams-Moulton method](https://en.wikipedia.org/wiki/Linear_multistep_method#Adams%E2%80%93Moulton_methods)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.15:** Implement the Adams-Moulton method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def adams_moulton(f, x_0, x_n, y_0, n): \n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(adams_moulton(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.16:** Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the Adams-Moulton method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predictor-corrector methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [predictor-corrector method](https://en.wikipedia.org/wiki/Predictor%E2%80%93corrector_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.17:** Implement the predictor-corrector method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predictor_corrector(f, x_0, x_n, y_0, n):\n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(predictor_corrector(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.18:** Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the predictor-corrector method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bulirsch-Stoer method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [Bulirsch-Stoer algorithm](https://en.wikipedia.org/wiki/Bulirsch%E2%80%93Stoer_algorithm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.19:** Implement the Bulirsch-Stoer method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bulirsch_stoer(f, x_0, x_n, y_0, n, m=5): \n",
    "\n",
    "    # add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y): \n",
    "    return -0.5 * y\n",
    "try:\n",
    "    np.testing.assert_almost_equal(bulirsch_stoer(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \\\n",
    "        decimal=4)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 10.20:** Solve the decay equation for [Ra-226](https://en.wikipedia.org/wiki/Isotopes_of_radium#Radium-226) using the Bulirsch-Stoer method and compare the numerical and analytical solutions.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}