part1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: Gaussian Elimination\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "2 & -1 & 3 & 7\\\\\n",
    "4 & 4 & 0 & 7 \\\\\n",
    "2 & 1 & 1 & 3 \\\\\n",
    "6 & 5 & 4 & 17\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$b = \\begin{bmatrix}\n",
    "15 \\\\\n",
    "11 \\\\\n",
    "7 \\\\\n",
    "31\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "$$Ax = \\begin{bmatrix}\n",
    "2 & -1 & 3 & 7\\\\\n",
    "4 & 4 & 0 & 7 \\\\\n",
    "2 & 1 & 1 & 3 \\\\\n",
    "6 & 5 & 4 & 17\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x_1 \\\\\n",
    "x_2 \\\\\n",
    "x_3 \\\\\n",
    "x_4\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "15 \\\\\n",
    "11 \\\\\n",
    "7 \\\\\n",
    "31\n",
    "\\end{bmatrix} = b\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to solve the system above with Gaussian elimination s.t. $Ax = b$. First we import our dependencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from typing import Any"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create the function gauss, which rewrites the matrix into a upper diagonal matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss(matrix: np.ndarray, answer: np.ndarray):\n",
    "    number_of_variables = matrix.shape[0]\n",
    "    index_vector: np.ndarray = np.arange(number_of_variables)\n",
    "    maximum_vector: np.ndarray[int, Any] = np.zeros(len(index_vector))\n",
    "    matrix = np.column_stack((matrix, answer))\n",
    "    number_of_variables = matrix.shape[0]\n",
    "    j = 0\n",
    "    \n",
    "    for i in range(number_of_variables):\n",
    "        index_vector[i] = i\n",
    "        smax = 0\n",
    "        \n",
    "        for j in range(number_of_variables):\n",
    "            smax = max(smax, abs(matrix[i, j]))\n",
    "\n",
    "        maximum_vector[i] = smax\n",
    "\n",
    "    # scaled partial pivoting\n",
    "    for k in range(number_of_variables - 1):\n",
    "        rmax = 0\n",
    "        for i in range(k, number_of_variables):\n",
    "            r = abs(matrix[index_vector[i], k] / maximum_vector[index_vector[i]])\n",
    "\n",
    "            if(r > rmax):\n",
    "                rmax = r\n",
    "                j = i\n",
    "\n",
    "        switch = index_vector[j]\n",
    "        index_vector[j] = index_vector[k]\n",
    "        index_vector[k] = switch\n",
    "\n",
    "        # elimination\n",
    "        for i in range(k + 1, number_of_variables):\n",
    "            xmult = matrix[index_vector[i], k]/matrix[index_vector[k], k]\n",
    "            matrix[index_vector[i]] = matrix[index_vector[i]] - xmult * matrix[index_vector[k]] # very cool matrix function\n",
    "\n",
    "        print(f\"Step: {k + 1}\\n\", matrix, \"\\n\")\n",
    "\n",
    "    # sort afterwards\n",
    "    sorted_matrix = np.zeros(matrix.shape)\n",
    "    for i in range(number_of_variables):\n",
    "        sorted_matrix[i] = matrix[index_vector[i]]\n",
    "\n",
    "    return sorted_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we create the function solve, which solves this matrix using backwards substitution. Below you can clearly see how the Gaussian substition happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step: 1\n",
      " [[ 0 -2  2  4  8]\n",
      " [ 0  2 -2  1 -3]\n",
      " [ 2  1  1  3  7]\n",
      " [ 0  2  1  8 10]] \n",
      "\n",
      "Step: 2\n",
      " [[ 0  0  0  5  5]\n",
      " [ 0  2 -2  1 -3]\n",
      " [ 2  1  1  3  7]\n",
      " [ 0  0  3  7 13]] \n",
      "\n",
      "Step: 3\n",
      " [[ 0  0  0  5  5]\n",
      " [ 0  2 -2  1 -3]\n",
      " [ 2  1  1  3  7]\n",
      " [ 0  0  3  7 13]] \n",
      "\n",
      "Upper diagonal matrix:\n",
      " [[ 2.  1.  1.  3.  7.]\n",
      " [ 0.  2. -2.  1. -3.]\n",
      " [ 0.  0.  3.  7. 13.]\n",
      " [ 0.  0.  0.  5.  5.]]\n",
      "\n",
      "Solution:  [1. 0. 2. 1.]\n"
     ]
    }
   ],
   "source": [
    "def solve(matrix: np.ndarray):\n",
    "    # backwards substitution\n",
    "    number_of_variables = matrix.shape[0]\n",
    "    answer = np.zeros(number_of_variables)\n",
    "    for i in range(number_of_variables - 1, -1, -1):\n",
    "        for j in range(matrix.shape[1] - 2, i, -1):\n",
    "            matrix[i, -1] = matrix[i, -1] - matrix[i, j] * answer[j]\n",
    "        answer[i] = answer[i] + matrix[i, -1]/matrix[i,i]\n",
    "\n",
    "    return answer\n",
    "\n",
    "def main():\n",
    "    matrix = np.array([[2, -1, 3, 7], \n",
    "                  [4, 4, 0, 7],\n",
    "                  [2, 1, 1, 3],\n",
    "                  [6, 5, 4, 17]])\n",
    "\n",
    "    answer_vector = np.array([15, 11, 7, 31])\n",
    "    \n",
    "    matrix = gauss(matrix, answer_vector)\n",
    "    print(\"Upper diagonal matrix:\\n\", matrix)\n",
    "    print(\"\\nSolution: \", solve(matrix))\n",
    "\n",
    "main()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "1c7434731747d88815ea98a4470429bf1ebbf9287f8433fda1e0ecccc300fe0c"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}