Lectures/Lecture04_linear_systems_Part1.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyOfqZV0gyl0EK134uVmQWL0",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/hamidrezanorouzi/numericalMethods/blob/main/Lectures/Lecture04_linear_systems_Part1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# **System of Linear Equations - Part 1**\n",
        "\n",
        "&nbsp;\n",
        "\n",
        "Lecturer: **Hamidreza Norouzi**\n",
        "\n",
        "&nbsp;"
      ],
      "metadata": {
        "id": "-ueWCmps_jP8"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### **Notes**❕\n",
        "This content forms a part of the instructional presentations for the **`numerical methods in chemical engineering`** course designed for undergraduate chemical engineering students at Amirkabir University of Technology.\n",
        "\n",
        "Feel free to utilize the information and source codes provided in this material, ensuring appropriate acknowledgment of the original document.\n",
        "\n",
        "The visual elements featured in this document are either original or have been obtained from the following sources, unless specified otherwise:\n",
        "\n",
        "* Steven C. Chapra, Applied Numerical Methods with Matlab for Engineers and Scientists, 3rd edition, McGraw-Hill (2012).\n",
        "* Amos Gilat and Vish Subramanian, Numerical Methods for Engineers and Scientists, 3rd edition, Wiley (2014).\n",
        "\n",
        "<div align=\"center\">\n",
        "🟧 🟧 🟧\n",
        "</dive>\n",
        "\n",
        "---"
      ],
      "metadata": {
        "id": "J5VWgZjnGDwn"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 🔴 1) Where do we face solving a set of linear equations?"
      ],
      "metadata": {
        "id": "_eI4X3NHGUra"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* **Chemical Systems (Tens to Thousands of Equations):** M&E balances are employed in diverse chemical systems such as distillation columns, adsorption columns, and chemical processes. These systems often entail solving systems of equations ranging from tens to thousands, reflecting the intricate interplay of chemical reactions and mass/energy transfers.\n",
        "\n",
        "* **Reactor Engineering (Tens of Equations):** Reactor engineering deals with systems characterized by numerous elementary reactions. Typically, tens of equations are formulated to describe reaction kinetics and species concentrations within the reactor.\n",
        "\n",
        "* **Curve Fitting (Usually Less than 10 Equations):** Curve fitting scenarios involve the determination of linear systems, typically containing fewer than ten equations. These are essential for fitting experimental data to mathematical models or functions.\n",
        "\n",
        "* **Solving ODEs and PDEs (Thousands to Millions of Equations):** Solving ordinary differential equations (ODEs) and partial differential equations (PDEs), often encountered in boundary value problems, necessitates techniques like linearization and discretization. These processes transform complex continuous problems into sets of linear equations that may range from thousands to millions in scale.\n"
      ],
      "metadata": {
        "id": "zyudnRtXG972"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 🔵 2) Cramer's rule\n",
        "* It is best suited for small systems, typically a set with  2 to 4 equations.\n",
        "\n",
        "\n",
        "\\begin{cases}\n",
        "a_{1,1}x_1 + a_{1,2}x_2 + ... + a_{1,n}x_n = b_1 \\\\\n",
        "a_{2,1}x_1 + a_{2,2}x_2 + ... + a_{2,n}x_n = b_2 \\\\\n",
        ". \\\\\n",
        ". \\\\\n",
        ". \\\\\n",
        "a_{n,1}x_1 + a_{n,2}x_2 + ... + a_{n,n}x_n = b_n \\\\\n",
        "\\end{cases}\n",
        "\n",
        "* In matrix form it is:\n",
        "$$\n",
        "[A]\\{x\\}=\\{B\\}\n",
        "$$\n"
      ],
      "metadata": {
        "id": "jiCFlMsSDwAc"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* The solution to the above set is given by:\n",
        "$$\n",
        "x_1 = \\frac{D_1}{D}, x_2=\\frac{D_2}{D}, ..., x_n = \\frac{D_n}{D}  \\tag{2-1}\n",
        "$$\n",
        "\n",
        "&nbsp;\n",
        "\n",
        "\\begin{align}\n",
        "D_1 &=\n",
        "\\begin{vmatrix}\n",
        "b_{1} & a_{1,2} & ... & a_{1,n} \\\\\n",
        "b_{2} & a_{2,2} & ... & a_{2,n} \\\\\n",
        ". \\\\ . \\\\ . \\\\\n",
        "b_{n} & a_{n,2} & ... & a_{n,n}\n",
        "\\end{vmatrix} \\notag\n",
        "&,\n",
        "D_n =\n",
        "\\begin{vmatrix}\n",
        "a_{1,1} & a_{1,2} & ... & b_{1} \\\\\n",
        "a_{2,1} & a_{2,2} & ... & b_{2} \\\\\n",
        ". \\\\ . \\\\ . \\\\\n",
        "a_{n,1} & a_{n,2} & ... & b_{n}\n",
        "\\end{vmatrix} \\tag{2-2}\n",
        "&&\n",
        "\\end{align}\n",
        "\n",
        "\n",
        "\\begin{align}\n",
        "D &=\n",
        "\\begin{vmatrix}\n",
        "a_{1,1} & a_{1,2} & ... & a_{1,n} \\\\\n",
        "a_{2,1} & a_{2,2} & ... & a_{2,n} \\\\\n",
        ". \\\\ . \\\\ . \\\\\n",
        "a_{n,1} & a_{n,2} & ... & a_{n,n}\n",
        "\\end{vmatrix} \\tag{2-3}\n",
        "&&\n",
        "\\end{align}"
      ],
      "metadata": {
        "id": "mQRv-XuUkgo3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Arithmetic operations in Cramer's rule\n",
        "* When the number of equations increases, the required operations grow exponentially and the Cramer's rule becomes infeasible.\n",
        "\n",
        "| No. equations | No. determinantes | No. multiplication | Total operations|\n",
        "|:---|:---:|:---:|:---:|\n",
        "|2   | 3 | 2 | 6|\n",
        "|3 | 4 | 12 | 48 |\n",
        "|n | n+1 | (n-1)n! | (n+1)(n-1)n!|"
      ],
      "metadata": {
        "id": "5cCyy7CYZZ4y"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 🟢 3) Gauss elimination"
      ],
      "metadata": {
        "id": "9MusQFbYNmW9"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3-1) Gauss elimination without pivoting\n",
        "\n",
        "* Consider the following set of equations:\n",
        "$$\n",
        "[A]\\{x\\} = \\{B\\}\n",
        "$$\n",
        " where $[A]$ is the matrix of coefficient and $\\{B\\}$ is the vector of known variables and $\\{x\\}$ is the vector of unknowns.\n",
        "* Gauss elimination is based on the operations (on rows of the system) to convert the coefficient matrix into a upper-triangle matrix:\n",
        "$$\n",
        "[U]\\{x\\} = \\{B^{'}\\}\n",
        "$$\n",
        "\n",
        "* Then using back-substitution from the last equation to obtain the solution of the system."
      ],
      "metadata": {
        "id": "Wbao-_GRP7xw"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<div align=\"center\">\n",
        "<img src=\"https://drive.google.com/uc?id=1f2csAYPmDogI2vuas9KbmJ7JOyo6GcHq\" width = \"400\">\n",
        "</div>"
      ],
      "metadata": {
        "id": "q3qYCtY36k1c"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Rules of for triangulation:\n",
        "\n",
        "* Any equation can be multiplied by a nonzero scalar without affecting the solution.\n",
        "* Any equation can be added to (or subtracted from) another equation without affecting the solution.\n",
        "* Any two equations can interchange positions within the set without affecting the solution.\n"
      ],
      "metadata": {
        "id": "FjKHe2p4UzvD"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### ❓ **Example 1:**\n",
        "Solve the following system of equations using gauss elimination.\n",
        "\n",
        "\\begin{cases}\n",
        "3x_1 +18 x_2+ 9x_3= 18\\\\\n",
        "2x_1 +3 x_2 +3 x_3=117\\\\\n",
        "4x_1 +1 x_2 +2 x_3=283\n",
        "\\end{cases}\n",
        "\n",
        "&nbsp;\n"
      ],
      "metadata": {
        "id": "fjcfhraJDFj1"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "💡 *Solution*\n",
        "\n",
        "\n",
        "The above set in matrix form becomes:\n",
        "$$ \\begin{bmatrix}\n",
        "3 & 18 & 9 \\\\\n",
        "2 & 3 & 3 \\\\\n",
        "4 & 1 & 2\n",
        "\\end{bmatrix}\n",
        "\\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3\\end{bmatrix} =\n",
        "\\begin{bmatrix} 18 \\\\ 117\\\\ 283 \\end{bmatrix}\n",
        "$$"
      ],
      "metadata": {
        "id": "xHWu_lQTPIem"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Step 1:** Write the augmented matrix:\n",
        "\n",
        "\n",
        "\\begin{array}{ccc|c}\n",
        "3 & 18 & 9 & 18 \\\\\n",
        "2 & 3 & 3 & 117 \\\\\n",
        "4 & 1 & 2 & 283 \\\\\n",
        "\\end{array}"
      ],
      "metadata": {
        "id": "5cokIKJoOZL7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Step 2:** $a_{1,1}$ is chosen as pivot element to make elements under the pivot element zero:\n",
        "\n",
        "1) Row2 - (2/3) * Row1 → Row2\n",
        "\n",
        "\\begin{array}{ccc|c}\n",
        "3 & 18 & 9 & 18 \\\\\n",
        "0 & -9 & -3 & 105 \\\\\n",
        "4 & 1 & 2 & 283 \\\\\n",
        "\\end{array}\n",
        "\n",
        "2) Row3 - (4/3) * Row1 → Row3\n",
        "\n",
        "\\begin{array}{ccc|c}\n",
        "3 & 18 & 9 & 18 \\\\\n",
        "0 & -9 & -3 & 105 \\\\\n",
        "0 & -23 & -10 & 259 \\\\\n",
        "\\end{array}\n",
        "\n"
      ],
      "metadata": {
        "id": "Jrj0MQPrFDZt"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Step 3:** $a_{2,2}$ is chosen as pivot element to make elements under the pivot element zero:\n",
        "\n",
        "1) Row3 - (23/7) * Row2 → Row3\n",
        "\n",
        "\\begin{array}{ccc|c}\n",
        "3 & 18 & 9 & 18 \\\\\n",
        "0 & -9 & -3 & 105 \\\\\n",
        "0 &  0 & \\frac{-7}{3} & \\frac{-28}{3} \\\\\n",
        "\\end{array}\n",
        "\n"
      ],
      "metadata": {
        "id": "2O4X3daSOeMG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Step 4:** Back substitution:\n",
        "\n",
        "&nbsp;\n",
        "\n",
        "\\begin{cases}\n",
        "3x_1+18x_2+9x_3 = 18 \\\\\n",
        "-9x_2 - 3x_3 = 105 \\\\\n",
        "-\\frac{7}{3}x_3 = \\frac{-28}{3}\n",
        "\\end{cases}\n",
        "\n",
        "&nbsp;\n",
        "\n",
        "\n",
        "$$\n",
        "x_3 = 4,  x_2 = -13, x_1 = 117\n",
        "$$"
      ],
      "metadata": {
        "id": "Sp1d0__9OgB_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3-2) Code for Gauss elimination without pivoting"
      ],
      "metadata": {
        "id": "8wwEzeaM1QXD"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "def gaussEliminationNoPivoting(A, B):\n",
        "\n",
        "  n = len(A)\n",
        "  augMat = np.concatenate([A,B],1)\n",
        "\n",
        "  # forward pass\n",
        "  for k in range(n-1):\n",
        "\n",
        "    # making elements below the pivot element zero\n",
        "    for m in range(k+1,n):\n",
        "      augMat[m,:] -= (augMat[m,k]/augMat[k,k]) * augMat[k,:]\n",
        "\n",
        "  # back substitution\n",
        "  x = np.zeros((n,1))\n",
        "  x[n-1] = augMat[n-1,n]/augMat[n-1,n-1]\n",
        "\n",
        "  for k in range(n-2, -1, -1):\n",
        "    x[k] = (augMat[k,n] - np.dot(augMat[k,k+1:n],x[k+1:n]) )/augMat[k,k]\n",
        "\n",
        "  return x"
      ],
      "metadata": {
        "id": "TYAwEWAH10GK"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "A = np.array([[3.0, 18 , 9],[2,3,3],[4,1,2]]);\n",
        "B = np.array([[18, 117, 283]]).T\n",
        "\n",
        "x = gaussEliminationNoPivoting(A,B)\n",
        "print(x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "sDWzCL5n28Kp",
        "outputId": "1f7ede50-7d4c-4659-d0fe-515eace2f8a8"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[ 72.]\n",
            " [-13.]\n",
            " [  4.]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3-3) Gauss elimination with partial pivoting\n",
        "* Consider the following system of equations:\n",
        "  \n",
        "\\begin{cases}\n",
        "2x_2 + 5x_3 = 10\\\\\n",
        "x_1+ 3x_2 - x_3 = 4\\\\\n",
        "3x_1+x_2+x_3 = -2\n",
        "\\end{cases}\n",
        "\n",
        "* This will give the following form:\n",
        "\n",
        "$$ \\begin{bmatrix}\n",
        "0 & 2 & 5 \\\\\n",
        "1 & 3 & -1 \\\\\n",
        "3 & 1 & 1\n",
        "\\end{bmatrix}\n",
        "\\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3\\end{bmatrix} =\n",
        "\\begin{bmatrix} 10 \\\\ 4\\\\ -2 \\end{bmatrix}\n",
        "$$\n"
      ],
      "metadata": {
        "id": "i9PMpiJ_u-rO"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* In each step, before eliminating the elements below the pivot element, it is advantageous to determine the coefficient with the **largest absolute value** in the column **below the pivot element**. The rows can then be **switched** so that the largest element is the pivot element.\n",
        "\n",
        "* By rearranging rows of the augmented matrix, we can obtain a diagonal dominant matrix (partial pivoting).\n",
        "\n",
        " * This reduces the possibility of division by zero.\n",
        "\n",
        " * Increases the accuracy of calculation by dividing the elements by larger values (decreasing round-off error).\n"
      ],
      "metadata": {
        "id": "q0jugQ07yKbp"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* Lets consider the following system of two equations. The exact solution is $x_1 = 1/3$ and $x_2 = 2/3$\n",
        "\\begin{cases}\n",
        "3.0\\times10^{-12}x_1 + 3x_2 = 1 \\\\\n",
        "2.000000000001x_1 + x_2 = 1\n",
        "\\end{cases}\n",
        "* Now, lets test the solution of Gauss elemination without pivoting:"
      ],
      "metadata": {
        "id": "XVCe-AAc5h6r"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "5ZqlZZH4GDKM",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "1a58ddc7-83f6-4701-b78b-7e0565ac97b3"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Answer is \n",
            " [[0.33336297]\n",
            " [0.66666667]]\n"
          ]
        }
      ],
      "source": [
        "A = np.array([[3.0e-12, 3.0],[1.0,1.0]]);\n",
        "B = np.array([[2.000000000001, 1]]).T\n",
        "\n",
        "x = gaussEliminationNoPivoting(A,B)\n",
        "print(\"Answer is \\n\",x)"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3-4) Code for Gauss elimination with pivoting\n"
      ],
      "metadata": {
        "id": "M5Rep1To6ANA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "def gaussElimination(A, B):\n",
        "  tol = 0.000001\n",
        "  #number of rows and columns\n",
        "  nr,nc = A.shape;\n",
        "  if nr!=nc:\n",
        "    print('input matrix A is not square!')\n",
        "    return None\n",
        "\n",
        "  #number of equations in the set\n",
        "  n,m = B.shape\n",
        "  if n!= nr:\n",
        "    print('A and B size mismatch')\n",
        "    return None\n",
        "\n",
        "  Aug = np.concatenate([A,B],1);\n",
        "  detA = 1;\n",
        "\n",
        "  #main loop\n",
        "  for k in range(n-1):\n",
        "\n",
        "    #first: partial pivoting\n",
        "    pRow = np.argmax(abs(Aug[k:, k])) + k\n",
        "\n",
        "    #interchanges the rows, if necessary\n",
        "    if pRow != k:\n",
        "      temp = np.array(Aug[k,:])\n",
        "      Aug[k,:] = Aug[pRow,:]\n",
        "      Aug[pRow,:] = temp\n",
        "      detA = -detA  #change of sign\n",
        "\n",
        "\n",
        "    if abs(Aug[k,k]) < tol:\n",
        "      print('Singular matrix!')\n",
        "      return None\n",
        "\n",
        "\n",
        "    # making elements below the pivot element zero\n",
        "    for m in range(k+1,n):\n",
        "      Aug[m,:] -= Aug[m,k]/Aug[k,k] * Aug[k,:]\n",
        "\n",
        "  #the last equation never checked\n",
        "  if abs(Aug[n-1,n-1]) < tol:\n",
        "    print('Singular matrix!')\n",
        "    return None\n",
        "\n",
        "  # back substitution\n",
        "  x = np.zeros((n,1))\n",
        "  x[n-1] = Aug[n-1,n]/Aug[n-1,n-1];\n",
        "  detA = detA*Aug[n-1,n-1];\n",
        "\n",
        "  for k in range(n-2, -1, -1):\n",
        "    x[k] = (Aug[k,n] - np.dot(Aug[k,k+1:n],x[k+1:n]) )/Aug[k,k]\n",
        "    detA = detA * Aug[k,k]\n",
        "\n",
        "  return x, detA;"
      ],
      "metadata": {
        "id": "9ZiNq0Xb6t22"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "A = np.array([[3.0e-12, 3.0],[1.0,1.0]]);\n",
        "B = np.array([[2.00000000001, 1]]).T\n",
        "\n",
        "x, det = gaussElimination(A,B)\n",
        "print(\"Answer is \\n\",x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "6jzZPAfa7EfR",
        "outputId": "0a55e72b-4671-46f4-c7f5-86cae95a3337"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Answer is \n",
            " [[0.33333333]\n",
            " [0.66666667]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3-5) Determinant evaluation with Gauss elimination\n",
        "* The determinant of a **triangular matrix** can be simply computed as the product of its diagonal elements:\n",
        "\n",
        "$$\n",
        "D = a_{1,1}.a_{2,2}.a_{3,3} . . . a_{n,n} \\tag{3-1}\n",
        "$$\n",
        "\n",
        "* If the programs uses **partial pivoting** for obtaining the triangular  matrix, the sign of the determinant switches with every row interchange. So, Eq. (3-1) should be changed to:\n",
        "\n",
        "$$\n",
        "D = a_{1,1}.a_{2,2}.a_{3,3} . . . a_{n,n}(-1)^p \\tag{3-2}\n",
        "$$\n",
        " where $p$ is the number of times that rows are interchanged/pivoted."
      ],
      "metadata": {
        "id": "SMlHdt3JmuJk"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 4) Tridiagonal systems\n",
        "* Some chemical engineering problems lead to a set of equations whose most of its elements in the coefficient matrix is zero: sparse matrix.\n",
        "* If nonzero elements are clustered about the diagonal elements, the matrix is called banded. As an example, the following matrix is tri-diagonal matrix.\n",
        "\n",
        "<div align=\"center\">\n",
        "<img src=\"https://drive.google.com/uc?id=1plVWLLmIAOiFkV7iuTBiqQ_wQ52ctQOt\" width=\"500\">\n",
        "</div>"
      ],
      "metadata": {
        "id": "78WTsppWoQK5"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* Memory requirement:\n",
        " * Consider $n×n$ matrix, Gauss elimination requires $n^2$ elements to store the values of the coefficients.\n",
        "\n",
        " * But with the tri-diagonal notation it only requires to save $e_i$,  $f_i$ and $g_i$, each has $(n-1) + n + (n-1)$ elements $→ 3n-2$ elements in total.\n",
        "* Arithmetic operations:\n",
        " * Applying Gauss elimination method requires $O(n^3)$ operations [explained in section 9.2.2] while tri-diagonal method only requires $O(n)$ operations."
      ],
      "metadata": {
        "id": "Id19FQdGp3Se"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## ❓ **Example 2**\n",
        "Solve the following tri-diagonal system using Gauss elimination.\n",
        "\n",
        "$$ \\begin{bmatrix}\n",
        "2.04 & -1 & 0 & 0 \\\\\n",
        "-1 & 2.04 & -1& 0 \\\\\n",
        "0 & -1 & 2.04& -1 \\\\\n",
        "0 & 0 & -1 & 2.04 \\\\\n",
        "\\end{bmatrix}\n",
        "\\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4\\end{bmatrix} =\n",
        "\\begin{bmatrix} 40.8 \\\\ 0.8\\\\ 0.8 \\\\ 200.8 \\end{bmatrix}\n",
        "$$"
      ],
      "metadata": {
        "id": "GXSMFWboxdef"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "💡 *solution*\n",
        "\n",
        "1) Similar to Gauss elimination, we start zeroing the elements under each pivot element.\n",
        "$$\n",
        "f_2 = f_2 - \\frac{e_2}{f_1}g_1 = 2.04 - \\frac{-1}{2.04}(-1)=1.55 \\\\\n",
        "r_2 = r_2 - \\frac{e_2}{f_1}r_1 = 0.8 - \\frac{-1}{2.04}(40.8) = 20.8\n",
        "$$\n",
        "2) Notice that $g_2$ is unmodified, since the elemenet above it is zero.\n",
        "\n",
        "3) Perforoming a similar calculation will result the following system:\n",
        "\n",
        "$$ \\begin{bmatrix}\n",
        "2.04 & -1 & 0 & 0 \\\\\n",
        "0 & 1.550 & -1& 0 \\\\\n",
        "0 & 0 & 1.395& -1 \\\\\n",
        "0 & 0 & 0 & 1.323 \\\\\n",
        "\\end{bmatrix}\n",
        "\\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4\\end{bmatrix} =\n",
        "\\begin{bmatrix} 40.8 \\\\ 20.8\\\\ 14.221 \\\\ 210.996 \\end{bmatrix}\n",
        "$$"
      ],
      "metadata": {
        "id": "TrBsApBcyqzy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "4) Back substitution:\n",
        "\n",
        "\n",
        "\\begin{aligned}\n",
        "&x_4 = \\frac{r_4}{f_4} = \\frac{210.996}{1.323} =159.480 \\\\\n",
        "&x_3= \\frac{r_3-g_3x_4}{f_3} = \\frac{14.221-(-1)159.480}{1.395} = 124.538 \\\\\n",
        "&x_2= \\frac{r_2-g_2x_3}{f_2} = \\frac{20.800-(-1)124.538}{1.550} = 93.778 \\\\\n",
        "&x_1= \\frac{r_1-g_1x_2}{f_1} = \\frac{80.800-(-1)93.778}{2.04} = 65.970\\\\\n",
        "\\end{aligned}\n"
      ],
      "metadata": {
        "id": "3D6fb3mu0X0k"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### General formula for Tri-diagonal system\n",
        "* It is called **Thomas** algorithm and consists of two steps:\n",
        " * Forward swip:\n",
        " $$\n",
        " \\begin{aligned}\n",
        " g^{'}_i &=\n",
        " \\begin{cases}  \n",
        " \\frac{g_i}{f_i} & i=1 \\\\ \\frac{g_i}{f_i-e_i g^{'}_{i-1}} & i=2,3,...,n-1\n",
        " \\end{cases} \\notag\n",
        " \\\\\n",
        " r^{'}_i &=\n",
        " \\begin{cases}\n",
        " \\frac{r_i}{f_i} & i=1 \\\\\n",
        " \\frac{r_i-e_i r^{'}_{i-1}}{f_i-e_i g^{'}_{i-1}} & i=2,3,...,n\n",
        " \\end{cases}\n",
        " \\end{aligned}\n",
        " $$\n",
        "\n",
        " * Back substitution:\n",
        "  \n",
        "  $$\n",
        "  \\begin{aligned}\n",
        "  x_n &= r^{'}_n \\\\\n",
        "  x_i &= r^{'}_i - g^{'}_i x_{i+1} \\ \\ \\ \\  \\text{for} \\ \\ \\ i= n-1, n-2, ..., 1\n",
        "  \\end{aligned}\n",
        "  $$"
      ],
      "metadata": {
        "id": "BxssmqgO48D0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "def triDiagonal(e, f, g, r):\n",
        "\n",
        "    #number of equations in the set\n",
        "    n = len(f);\n",
        "\n",
        "    gg = np.zeros((n,));\n",
        "    rr = np.zeros((n,));\n",
        "\n",
        "    #forward sweeping\n",
        "    gg[0] = g[0]/f[0];\n",
        "\n",
        "    for i in range(1,n-1):\n",
        "        gg[i] = g[i]/(f[i]-e[i]*gg[i-1]);\n",
        "\n",
        "    rr[0] = r[0]/f[0];\n",
        "    for i in range(1,n):\n",
        "        rr[i] = (r[i]-e[i]*rr[i-1])/(f[i]-e[i]*gg[i-1]);\n",
        "\n",
        "    # back substitution\n",
        "    X = np.zeros((n,));\n",
        "    X[n-1] = rr[n-1];\n",
        "    for i in range(n-2,-1,-1):\n",
        "        X[i] = rr[i] - gg[i]*X[i+1];\n",
        "\n",
        "    return X;"
      ],
      "metadata": {
        "id": "PcLdvFgdCfup"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "e = np.array([0,-1,-1,-1])\n",
        "f = np.array([2.04, 2.04, 2.04, 2.04])\n",
        "g = np.array([-1,-1,-1,0])\n",
        "r = np.array([40.8, 0.8, 0.8, 200.8])\n",
        "\n",
        "x= triDiagonal(e,f,g,r)\n",
        "\n",
        "print(x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "6jyPpmehBspb",
        "outputId": "c5e8c176-9a99-4013-fbc7-ee572eeee406"
      },
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[ 65.96983437  93.77846211 124.53822833 159.47952369]\n"
          ]
        }
      ]
    }
  ]
}