{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 03 - Linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basic linear algebra operations, direct methods for solving linear equation systems, forward and backward substitution, Gaussian elimination, LU decomposition, Thomas algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scalar product\n",
    "\n",
    "The [scalar product](https://en.wikipedia.org/wiki/Dot_product) of two vectors $ x \\in \\mathbb{R}^{n} $ and $ y \\in \\mathbb{R}^{n} $ is defined as \n",
    "\n",
    "$$\n",
    "x \\cdot y = \\sum_{i=0}^{n-1} x_i y_i.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.1:** Write a function that calculates the scalar product of two vectors (do not use the [`numpy.dot`](https://numpy.org/doc/stable/reference/generated/numpy.dot.html) function nor the `@` operator).\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scalar_product(x, y):\n",
    "    \"\"\"\n",
    "    Calculates scalar product of two vectors.\n",
    "    Args:\n",
    "        x (array_like): Vector of size n\n",
    "        y (array_like): Vector of size n\n",
    "    Returns:\n",
    "        numpy.float: scalar product of x and y\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`numpy.dot`](https://numpy.org/doc/stable/reference/generated/numpy.dot.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.rand(3)\n",
    "y = np.random.rand(3)\n",
    "try:\n",
    "    np.testing.assert_array_almost_equal(scalar_product(x, y), np.dot(x, y), decimal=7)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix-vector product\n",
    "\n",
    "The product of matrix $ \\mathbb{A} \\in \\mathbb{R}^{m \\ \\times \\ n} $ and vector $ x \\in \\mathbb{R}^{n} $ is defined as\n",
    "\n",
    "$$\n",
    "(\\mathbb{A} \\cdot x)_i = \\sum_{j = 0}^{n-1} a_{ij} x_j, \\quad i = 0, 1, \\dots, m-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.2:** Write a function that calculates the product of a matrix and a vector (do not use the [`numpy.dot`](https://numpy.org/doc/stable/reference/generated/numpy.dot.html) function nor the `@` operator).\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_vector_product(A, x):\n",
    "    \"\"\"\n",
    "    Calculates matrix-vector product.\n",
    "    Args:\n",
    "        A (array_like): A m-by-n matrix\n",
    "        x (array_like): Vector of size n\n",
    "    Returns:\n",
    "        numpy.ndarray: Matrix-vector product\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`numpy.dot`](https://numpy.org/doc/stable/reference/generated/numpy.dot.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3)\n",
    "x = np.random.rand(3)\n",
    "try:\n",
    "    np.testing.assert_array_almost_equal(matrix_vector_product(A, x), np.dot(A, x), decimal=7)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix-matrix product\n",
    "\n",
    "The [product of two matrices](https://en.wikipedia.org/wiki/Matrix_multiplication) $ \\mathbb{A} \\in \\mathbb{R}^{m \\ \\times \\ n}  $ and $ \\mathbb{B} \\in \\mathbb{R}^{n \\ \\times \\ p}  $ is defined as\n",
    "\n",
    "$$\n",
    "(\\mathbb{A} \\cdot \\mathbb{B})_{ij} = \\sum_{k = 0}^{n-1} a_{ik} b_{kj}, \\qquad i = 0, 1, \\dots, m-1, \\quad j = 0, 1, \\dots, p-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.3:** Write a function that calculates the product of two matrices (do not use the `numpy.dot` function nor the `@` operator).\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_matrix_product(A, B):\n",
    "    \"\"\"\n",
    "    Calculates matrix-matrix product.\n",
    "    Args:\n",
    "        A (array_like): A m-by-n matrix\n",
    "        B (array_like): A n-by-p matrix\n",
    "    Returns:\n",
    "        numpy.ndarray: Matrix-matrix product\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`numpy.dot`](https://numpy.org/doc/stable/reference/generated/numpy.dot.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3)\n",
    "B = np.random.rand(3, 3)\n",
    "try:\n",
    "    np.testing.assert_array_almost_equal(matrix_matrix_product(A, B), np.dot(A, B), decimal=7)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct methods for linear algebraic equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forward substitution\n",
    "\n",
    "A system of linear algebraic equations $ \\mathbb{A} x = b $ with lower triangular coefficient matrix $ \\mathbb{A} \\in \\mathbb{R}^{n \\ \\times \\ n} $ can be solved by the [forward substitution](https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution) algorithm,\n",
    "\n",
    "$$\n",
    "x_i := \\frac{1}{a_{ii}} \\left( b_i - \\sum_{j = 0}^{i - 1} a_{ij} x_j \\right), \\qquad i = 0, 1, \\dots, n-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.4:** Implement the forward substitution algorithm.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_substitution(A, b):\n",
    "    \"\"\"\n",
    "    Solves a system of linear equations with lower triangular matrix.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n lower triangular matrix\n",
    "        b (array_like): RHS vector of size n\n",
    "    Returns:\n",
    "        numpy.ndarray: Vector of solution\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.tril(np.random.rand(3, 3)) # create a lower triangular matrix from a random matrix\n",
    "b = np.random.rand(3)\n",
    "try:\n",
    "    np.testing.assert_array_almost_equal(forward_substitution(A, b), la.solve(A, b), decimal=7)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Backward substitution\n",
    "\n",
    "A system of linear algebraic equations $ \\mathbb{A} x = b $ with upper triangular coefficient matrix $ \\mathbb{A} \\in \\mathbb{R}^{n \\ \\times \\ n} $ can be solved by the [backward substitution](https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution) algorithm,\n",
    "\n",
    "$$\n",
    "x_i := \\frac{1}{a_{ii}} \\left( b_i - \\sum_{j = i + 1}^{n-1} a_{ij} x_j \\right), \\qquad i = n-1, n-2, \\dots, 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.5:** Implement the backward substitution algorithm.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_substitution(A, b):\n",
    "    \"\"\"\n",
    "    Solves a system of linear equation with upper triangular matrix.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n upper triangular matrix\n",
    "        b (array_like): RHS vector of size n\n",
    "    Returns:\n",
    "        numpy.ndarray: Vector of solution\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.triu(np.random.rand(3, 3)) # create an upper triangular matrix from a random matrix\n",
    "b = np.random.rand(3)\n",
    "try:\n",
    "    np.testing.assert_array_almost_equal(backward_substitution(A, b), la.solve(A, b), decimal=7)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian elimination\n",
    "\n",
    "A matrix $ \\mathbb{A} \\in \\mathbb{R}^{n \\ \\times \\ n} $ can be transformed into an upper triangular form using the [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) algorithm,\n",
    "\n",
    "$$\n",
    "a_{jk} := a_{jk} - \\frac{a_{ji}}{a_{ii}} a_{ik}, \\qquad i = 0, 1, \\dots, n-1, \\quad j = i + 1, i + 2, \\dots, n-1, \\quad k = 0, 1, \\dots, n-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.6:** Implement the Gaussian elimination algorithm.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian_elimination(A):\n",
    "    \"\"\"\n",
    "    Transforms given matrix into an upper triangular form using the Gaussian elimination algorithm.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n matrix\n",
    "    Returns:\n",
    "        numpy.ndarray: Upper triangular matrix\n",
    "    \"\"\"\n",
    "\n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the function above to transform the following matrix $ \\mathbb{A} $ into an upper triangular form,\n",
    "\n",
    "$$\n",
    "\\mathbb{A} = \\begin{pmatrix}\n",
    "0 & 1 \\\\\n",
    "1 & 1\n",
    "\\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "**Warning:** The Gaussian elimination algorithm itself is unable to transform the matrix above into an upper triangular form because it leads to division by zero.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Pivoting:** The [pivot or pivot element](https://en.wikipedia.org/wiki/Pivot_element) is the element of a matrix, or an array, which is selected by an algorithm to do certain calculations. In the case of Gaussian elimination, the algorithm requires that pivot elements are distinct from zero. Therefore, interchanging rows or columns in the case of a zero pivot element is necessary. Furthermore, it is generally desirable to choose a pivot element with large absolute value, which dramatically reduces the round-off error.\n",
    "\n",
    "In partial pivoting, the algorithm selects the entry with largest absolute value from the column of the matrix that is currently being considered as the pivot element. In complete pivoting, the algorithm interchanges both rows and columns in order to use the largest (by absolute value) element in the matrix as the pivot. For Gaussian elination, partial pivoting is usually sufficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.7:** Implement pivoting into the Gaussian elimination algorithm.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian_elimination_with_pivoting(A):\n",
    "    \"\"\"\n",
    "    Transforms given matrix into an upper triangular form using the Gaussian elimination algorithm with pivoting.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n matrix\n",
    "    Returns:\n",
    "        numpy.ndarray: Upper triangular matrix\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the function above to transform the following matrix $ \\mathbb{A} $ into an upper triangular form,\n",
    "\n",
    "$$\n",
    "\\mathbb{A} = \\begin{pmatrix}\n",
    "0 & 1 \\\\\n",
    "1 & 1\n",
    "\\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.8:** Pass the right-hand side vector into the Gaussian elimination algorithm and perform the identical operations on it.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian_elimination_with_pivoting_vector(A, b):\n",
    "    \"\"\"\n",
    "    Transforms given matrix into an upper triangular form using the Gaussian elimination algorithm with pivoting, \n",
    "    performs identical operations on RHS vector.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n regular matrix\n",
    "        b (array_like): RHS vector of size n\n",
    "    Returns:\n",
    "        numpy.ndarray: Upper triangular matrix\n",
    "        numpy.ndarray: RHS vector corresponding to upper triangular matrix\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3)\n",
    "b = np.random.rand(3)\n",
    "A_upper, bb = gaussian_elimination_with_pivoting_vector(A, b)\n",
    "try:\n",
    "    np.testing.assert_array_almost_equal(backward_substitution(A_upper, bb), la.solve(A, b), decimal=7)\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-info\">\n",
    "\n",
    "**Note:** The code above uses Gaussian elimination to transform a given random matrix into an upper triangular form, solves a system of linear algebraic equations with the transformed matrix using the backward substitution, and compare the result with the one obtained by [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LU decomposition\n",
    "\n",
    "[LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition) factors a matrix as the product of a lower triangular matrix and an upper triangular matrix. LU decomposition can be achieved by [Crout algorithm](https://en.wikipedia.org/wiki/Crout_matrix_decomposition),\n",
    "\n",
    "$$\n",
    "u_{ij} = a_{ij} - \\sum_{k = 1}^{i - 1} l_{ik} u_{kj}, \n",
    "$$\n",
    "\n",
    "$$\n",
    "l_{ji} = \\frac{1}{u_{ii}} \\left( a_{ji} - \\sum_{k = 1}^{i - 1} l_{jk} u_{ki} \\right), \n",
    "$$\n",
    "\n",
    "$$\n",
    "i = 1, 2, \\dots, n, \\quad j = i, i + 1, \\dots, n.\n",
    "$$\n",
    "\n",
    "A system of linear algebraic equations $ \\mathbb{A} x = b $ can be solved using LU decomposition as follows,\n",
    "\n",
    "$$\n",
    "\\mathbb{A} x = (\\mathbb{L} \\mathbb{U}) x = \\mathbb{L} ( \\mathbb{U} x) = b, \n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbb{L} y = b, \n",
    "$$\n",
    "$$\n",
    "\\mathbb{U} x = y.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "**Exercise 03.9:** Implement the LU decomposition algorithm.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lu_decomposition(A):\n",
    "    \"\"\"\n",
    "    Factors given matrix as the product of a lower and an upper triangular matrix using LU decomposition.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n matrix\n",
    "    Returns:\n",
    "        numpy.ndarray: Lower triangular matrix\n",
    "        numpy.ndarray: Upper triangular matrix\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3)\n",
    "b = np.random.rand(3)\n",
    "L, U = lu_decomposition(A)\n",
    "try:\n",
    "    np.testing.assert_array_almost_equal(backward_substitution(U, forward_substitution(L, b)), la.solve(A, b), decimal=7)\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-info\">\n",
    "\n",
    "**Note:** The code above decomposes a given random matrix using the LU decomposition, solves the systems of algebraic equations using the forward and backward substitutions, and compare the result with the one obtained by [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Thomas algorithm\n",
    "\n",
    "A system of linear algebraic equations $ \\mathbb{A} x = b $ with tridiagonal coefficient matrix $ \\mathbb{A} $ can be efficiently solved by the [Thomas algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm). Tridiagonal matrix $ \\mathbb{A} $ can be represented by three vectors $ p $, $ q $, and $ r $, as follows,\n",
    "\n",
    "$$\n",
    "\\mathbb{A} = \\begin{pmatrix}\n",
    "q_0 & p_0 & 0 & 0 & \\cdots & 0 & 0 & 0 \\\\\n",
    "r_1 & q_1 & p_1 & 0 & \\cdots & 0 & 0 & 0 \\\\\n",
    "0 & r_2 & q_2 & p_2 & \\cdots & 0 & 0 & 0 \\\\\n",
    "\\vdots  & \\vdots  & \\ddots & \\ddots & \\ddots & \\vdots & \\vdots & \\vdots \\\\\n",
    "\\vdots  & \\vdots  & \\vdots & \\ddots & \\ddots & \\ddots & \\vdots & \\vdots \\\\\n",
    "\\vdots  & \\vdots  & \\vdots & \\vdots & \\ddots & \\ddots & \\ddots & \\vdots \\\\\n",
    "0 & 0 & 0 & 0 & \\cdots & r_{n-2} & q_{n-2} & p_{n-2} \\\\\n",
    "0 & 0 & 0 & 0 & \\cdots & 0 & r_{n-1} & q_{n-1}\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "\n",
    "where $ p_{n-1} = 0 $ and $ r_0 = 0 $. The solution of $ \\mathbb{A} x = b $ is then obtained iteratively as\n",
    "\n",
    "$$\n",
    "x_i = \\mu_i x_{i+1} + \\rho_i, \\qquad i = n-1, n-2, \\dots, 0,\n",
    "$$\n",
    "\n",
    "where the coefficients $ \\mu_i $ and $ \\rho_i $ are\n",
    "\n",
    "$$\n",
    "\\mu_i = -\\frac{p_i}{r_i \\mu_{i-1} + q_i}, \\qquad \\rho_i = \\frac{b_i - r_i \\rho_{i-1}}{r_i \\mu_{i-1} + q_i}, \\qquad i = 1, 2, \\dots, n-1\n",
    "$$\n",
    "\n",
    "with $ \\mu_0 = -p_0 \\, / \\, q_0 $ and $ \\rho_0 = b_0 \\, / \\, q_0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 03.10:** Implement the Thomas algorithm.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def thomas_algorithm(A, b):\n",
    "    \"\"\"\n",
    "    Solves system of linear equations with a tridiagonal matrix using Thomas algorithm.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n regular matrix\n",
    "        b (array_like): RHS vector of size n\n",
    "    Returns:\n",
    "        numpy.ndarray: Vector of solution\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.triu(np.tril(np.random.rand(5, 5), 1), -1) # create a tridiagonal matrix from random matrix\n",
    "b = np.random.rand(5)\n",
    "try:\n",
    "    np.testing.assert_array_almost_equal(thomas_algorithm(A, b), la.solve(A, b), decimal=7)\n",
    "except AssertionError as E:\n",
    "    print(E)\n",
    "else:\n",
    "    print(\"The implementation is correct.\")"
   ]
  }
 ],
 "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
}