{ "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": [ "
\n", "\n", "**Note:** In the case of explicit methods the Runge-Kutta matrix is lower triangular. \n", " \n", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.1:** Implement the forward Euler method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.3:** Implement the Runge-Kutta $ 2^{\\mathrm{nd}} $ order method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.5:** Implement the Runge-Kutta $ 3^{\\mathrm{rd}} $ order method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.7:** Implement the Runge-Kutta $ 4^{\\mathrm{th}} $ order method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.9:** Implement the backward Euler method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.11:** Implement the Leap-frog method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.13:** Implement the Adams-Bashforth method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.15:** Implement the Adams-Moulton method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.17:** Implement the predictor-corrector method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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": [ "
\n", "\n", "**Exercise 10.19:** Implement the Bulirsch-Stoer method.\n", " \n", "
" ] }, { "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": [ "
\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", "
" ] }, { "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 }