System_of_linear_equations_solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Goal of this project"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this project is to decompose one matrix into two L and U matrices using Doolittle and Crout's decomposition algorithms, in order to solve a system of linear equations.\n",
    "\n",
    "The system of linear equations we are trying to solve is:\n",
    "\n",
    "2x + 3y - 1 = 2 <br/>\n",
    "4x + 4y - 1 = -1<br/>\n",
    "-2x -3y + 4z = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing the libraries needed for the project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this project, the only library we need is numpy for manipulating arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next step: defining the functions for the project"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After this, we are going to define the main functions needed for the LU decomposition using Doolittle and Crout's methods.\n",
    "\n",
    "We will have four main functions: one for forward substitution, one for backward substitution, one for Doolittle's algorithm and one for Crout's algorithm. On top of this, we will also have an auxiliary function to compute the final results of each algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function for the forward substitution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_substitution(L, b):\n",
    "\n",
    "    y = np.full_like(b,0) # Creating the y vector the same size as the b vector\n",
    "    \n",
    "    for k in range(len(b)):\n",
    "        \n",
    "        y[k] = b[k]\n",
    "        \n",
    "        for i in range(k):\n",
    "            \n",
    "            y[k] = y[k] - (L[k, i]*y[i])\n",
    "            \n",
    "        y[k] = y[k] / L[k, k] # Using forward substitution to find the value of y\n",
    "    \n",
    "    # Returning the y vector\n",
    "    \n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This functions takes as input a lower triangular matrix L and a right-side vector b. It returns the solution vector y using the equation Ly = b.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function for the backwards substitution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_substitution(U, y):\n",
    "    \n",
    "    x = np.full_like(y,0) # Creating the x vector the same size as the y vector\n",
    "    \n",
    "    for k in range(len(x), 0, -1):\n",
    "        \n",
    "      x[k-1] = (y[k-1] - np.dot(U[k-1, k:], x[k:])) / U[k-1, k-1] # Using backward substitution to find the value of x\n",
    "     \n",
    "    # Returning the solution vector x\n",
    "    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This functions takes as input a lower triangular matrix U and a right-side vector y. It returns the solution vector x using the equation Ux = y."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function for Doolittle's algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "def doolittle(A):\n",
    "    \n",
    "  # Creating two L and U matrices filled with 0s and the same size as A\n",
    "\n",
    "  L = np.matrix([[0,0,0],[0,0,0],[0,0,0]])\n",
    "  U = np.matrix([[0,0,0],[0,0,0],[0,0,0]])\n",
    "  n = len(A)\n",
    "  \n",
    "  # for-loop in order to set the j,j-th entry of U to 1\n",
    "    \n",
    "  for z in range(n):\n",
    "        \n",
    "    L[z, z] = 1\n",
    "    \n",
    "    U[z, z] = (A[z, z] - np.dot(L[z, :z], U[:z, z]))\n",
    "    \n",
    "    for i in range(z+1, n):\n",
    "        \n",
    "      U[z, i] = (A[z, i] - np.dot(L[z, :z], U[:z, i]))\n",
    "    \n",
    "    for k in range(z+1, n):\n",
    "        \n",
    "      L[k, z] = (A[k, z] - np.dot(L[k, :z], U[:z, z])) / U[z, z]\n",
    "\n",
    "  # Returning the matrices L and U i.e. the A matrix decomposed using Doolittle's algorithm\n",
    "\n",
    "  return (L, U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function decomposes the matrix A into the two L and U matrices using Doolittle's algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function for Crout's algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crout(A):\n",
    "        \n",
    "    # Creating two L and U matrices filled with 0s and the same size as A\n",
    "\n",
    "    L = np.matrix([[0.0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]])\n",
    "    U = np.matrix([[0.0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]])\n",
    "    n = len(A)\n",
    "    \n",
    "    # for-loop in order to set the j,j-th entry of U to 1\n",
    "    \n",
    "    for z in range(n):\n",
    "        \n",
    "        U[z,z] = 1             \n",
    "        \n",
    "        # for-loop starting at L[j][j] in order to solve the j-th column of L\n",
    "        \n",
    "        for j in range(z,n):\n",
    "\n",
    "            # Declaring a temporary L to store values and insert them later in the L matrix\n",
    "\n",
    "            temporary_L = float(A[j,z])\n",
    "            \n",
    "            for k in range(z):\n",
    "                \n",
    "                temporary_L -= L[j,k]*U[k,z]\n",
    "                \n",
    "            L[j,z] = temporary_L\n",
    "            \n",
    "        # for-loop starting at U[j][j+1] in order to solve the j-th row of U\n",
    "        \n",
    "        for j in range(z+1, n):\n",
    "            \n",
    "            # Declaring a temporary U to store values and insert them later in the U matrix\n",
    "            \n",
    "            temporary_U = float(A[z,j])\n",
    "            \n",
    "            for k in range(z):\n",
    "                \n",
    "                temporary_U -= L[z,k]*U[k,j]\n",
    "                \n",
    "            U[z,j] = temporary_U / L[z,z]\n",
    "    \n",
    "    # Returning the matrices L and U i.e. the A matrix decomposed using Crout's algorithm\n",
    "    \n",
    "    return (L, U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function decomposes the matrix A into the two L and U matrices using Crout's algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function for computing the final solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "def computing_final_solution(A, b, algorithm_used):\n",
    "    \n",
    "    # Creating the L and U matrices using the specified algorithm\n",
    "    \n",
    "    L, U = algorithm_used(A)\n",
    "    \n",
    "    print(\"L = \" + str(L) + \"\\n\")\n",
    "    print(\"U = \" + str(U) + \"\\n\")\n",
    "    \n",
    "    # Calling forward then backward substitution\n",
    "    \n",
    "    y = forward_substitution(L,b)\n",
    "    x = backward_substitution(U,y)\n",
    "    \n",
    "    # Returning the solution vector x\n",
    "    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function decomposes the A matrix into L and U matrices by calling the specified algorithm (i.e. Doolittle or Crout). After this, it uses forward then backward substitution to find the x vector (i.e. the solution for the linear system of equations)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Declaring the matrices to solve the system of linear equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Declaring the A and b matrices for the final solution\n",
    "\n",
    "A = np.matrix([[2.0,3.0,-1.0],[4.0,4.0,-1.0],[-2.0,-3.0,4.0]])\n",
    "\n",
    "b = np.array([2.0,-1.0,1.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Printing out the final solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The solution using Doolittle's algorithm:\n",
      "\n",
      "L = [[ 1  0  0]\n",
      " [ 2  1  0]\n",
      " [-1  0  1]]\n",
      "\n",
      "U = [[ 2  3 -1]\n",
      " [ 0 -2  1]\n",
      " [ 0  0  3]]\n",
      "\n",
      "x = [-3.  3.  1.]\n",
      "\n",
      "\n",
      "The solution using Crout's algorithm:\n",
      "\n",
      "L = [[ 2.  0.  0.]\n",
      " [ 4. -2.  0.]\n",
      " [-2.  0.  3.]]\n",
      "\n",
      "U = [[ 1.   1.5 -0.5]\n",
      " [ 0.   1.  -0.5]\n",
      " [ 0.   0.   1. ]]\n",
      "\n",
      "x = [-3.  3.  1.]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Printing out the results\n",
    "\n",
    "# For Doolittle \n",
    "\n",
    "print(\"The solution using Doolittle's algorithm:\" + \"\\n\")\n",
    "print( \"x = \" + str(computing_final_solution(A,b, doolittle)) + \"\\n\" )\n",
    "\n",
    "# For Crout \n",
    "\n",
    "print(\"\\n\" + \"The solution using Crout's algorithm:\" + \"\\n\")\n",
    "print( \"x = \" + str(computing_final_solution(A,b, crout)) + \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, the final solution for both Doolittle and Crout algorithms is x = -3, y = 3 and z = 1. This means that those values of x, y and z can solve the system of linear equations."
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}