Linear Equation Solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# System of Linear Equations Solver\n",
    "_By Dhruv Jain_\n",
    "\n",
    "### **Objective: Solve a fully determined system of linear equations**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Key libraries: Numpy(for mathematical procedures) and matplotlib(to create plots)\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "import copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Elimination with back subsitution without pivoting and with scaled partial pivoting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ge_lin_equ(a,b,opt):\n",
    "    \"\"\"Dhruv Jain, 19 Sept, 2021\n",
    "    Obj: To solve system of linear equations using Gaussian Elimination with No Pivoting(opt = 1) and Scaled Partial Pivoting (opt = 2)\n",
    "    Args:\n",
    "        a: (n x n) matrix, Coefficient Matrix\n",
    "        b: (n x 1) vector, Right-hand side of the systems of equations\n",
    "        opt: integer, 1 => GE with No Pivoting, 2 => GE with Scaled Partial Pivoting\n",
    "    Returns:\n",
    "        x: (n x 1) vector, computed solution of the system  \n",
    "    \"\"\"\n",
    "    \n",
    "    n = len(b) # Number of columns = Number of rows = length of b vector\n",
    "    r = np.arange(n) # r = [1, 2, .... n-1, n] to track the location of equation \n",
    "    a_orig = a\n",
    "    a = np.hstack((a, np.transpose(np.array([b])))) # Augmented matric [A | b]\n",
    "    \n",
    "    if opt == 1: # GE with No Pivoting\n",
    "        for i in range(n):\n",
    "            for j in range(i+1,n):\n",
    "                if a[i,i] != 0: # Check if leading column element is not 0, if 0 then swap with row of nonzero leading element\n",
    "                    lam = a[j,i]/a[i,i] # Common factor to make subsequent column elements 0\n",
    "                    a[j,:] = a[j,:] - lam*a[i,:] # ERO3, so that all subsequent columen elements = 0\n",
    "                else: # If leading column element is 0 then we use ERO1, r_i <-> r_j\n",
    "                    k = j\n",
    "                    while a[k,i] == 0 and k < n-1: # Check till non-zero element in the leading column found\n",
    "                        k = k + 1\n",
    "\n",
    "                    if i < n-1 and k == n-1 and a[k,i] == 0: # If all subsequent columns are zero then move to next column\n",
    "                        j = n\n",
    "                    else: # If a non-zero element found then interchange rows and track interchange using r\n",
    "                        temp_a = copy.copy(a[i,:])\n",
    "                        a[i,:] = copy.copy(a[k,:])\n",
    "                        a[k,:] = copy.copy(temp_a)\n",
    "\n",
    "                        temp = copy.copy(r[i])\n",
    "                        r[i] = copy.copy(k)\n",
    "                        r[k] = copy.copy(temp)\n",
    "                        j = j - 1\n",
    "\n",
    "    elif opt == 2: # GE with Scaled Partial Pivoting\n",
    "        s = np.zeros(n)\n",
    "        M = np.zeros(n)\n",
    "        for i in range(n-1):\n",
    "            for k in range(i,n):\n",
    "                    origpos = r[k]\n",
    "                    s[k] = max(np.abs(a_orig[origpos])) # Calculation of scale vector using elements of the original A\n",
    "                    M[k] = np.abs(a[k,i])/s[k] # Calculation of scaled coefficients, to find pivot element\n",
    "            \n",
    "            # Interchange rows with pivot element calculated using scaled coefficients \n",
    "            posmax = np.argmax(M[i:])\n",
    "            temp_a = copy.copy(a[i,:])\n",
    "            a[i,:] = copy.copy(a[i+posmax,:])\n",
    "            a[i+posmax,:] = copy.copy(temp_a)\n",
    "\n",
    "            # Track the row change \n",
    "            temp = copy.copy(r[i])\n",
    "            r[i] = copy.copy(i+posmax)\n",
    "            r[i+posmax] = copy.copy(temp)\n",
    "\n",
    "        # ERO3, to elimiate subsequent leading terms \n",
    "            for j in range(i+1,n):\n",
    "                if a[i,i] != 0: # Check if leading column element is not 0, if 0 then swap with row of nonzero leading element\n",
    "                    lam = a[j,i]/a[i,i] # Common factor to make subsequent column elements 0\n",
    "                    a[j,:] = a[j,:] - lam*a[i,:]\n",
    "\n",
    "    # Back subsitution to obtain x, solution vector from the above simplified upper traingle augmented matrix\n",
    "    x = np.zeros(n)\n",
    "    for i in range(n-1,-1,-1):\n",
    "        if a[i,i] != 0:\n",
    "            x[i]=(a[i,-1]-a[i,i+1:n].dot(x[i+1:n]))/a[i,i] \n",
    "        else:\n",
    "            x[i] = 0\n",
    "    \n",
    "    return x    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complex Gradient Method (Minimization method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Conjugate Gradient Method:\n",
    "def cgm(A,b,x0, tol, tol_div, Nmax):\n",
    "    \"\"\"Dhruv Jain, 4 Oct, 2021, MSAAE Purdue\n",
    "    Obj: Use Complex Gradient Method (Minimization method) to iteratively solve for a solution of linear system of equations of the form Ax = b, x is the solution \n",
    "    Args:\n",
    "        A: (n x n) matrix, Coefficeient matrix\n",
    "        b: (n x 1) vector, right-hand side vector\n",
    "        x0: (n x 1) vecotr, inital guess for x*\n",
    "        tol: ||x_n+1 - x_n||infi, step convergence tolerance \n",
    "        tol_div: ||x_n+1 - x_n||infi, step Divergence tolerance\n",
    "        Nmax: max allowable iterations for the iterative solver\n",
    "    Output:\n",
    "        xfinal: (n_iter x n), contains solved x for each iteration, last row is the converged solution\n",
    "        iterations: # iterations to converge  \n",
    "    \"\"\"\n",
    "    r = np.matmul(A,x0) - b # Initial calculation of the residul vector\n",
    "    d = -r # Setup search direction vector\n",
    "    x = np.zeros((Nmax,len(b))) # Create x to save solution calculated after each iteration\n",
    "    x[0,:] = x0 # Assign 1st row as the Initial Guess\n",
    "    \n",
    "    # Goal is to find minima of f = 1/2x.T*A*x - b.T*x + c\n",
    "    for i in range(Nmax-1): # Run loop for set maximum iterations\n",
    "        lam = np.matmul(-d.T, r)/np.matmul(np.matmul(d.T,A),d) # Calculate step size\n",
    "        x[i+1,:] = x[i,:] + lam*d # Update solution using calculated step size and search direction\n",
    "        r = np.matmul(A,x[i+1,:]) - b # Update residual vector using new calculated solution\n",
    "        \n",
    "        if max(abs(x[i+1,:] - x[i,:])) < tol: # Check if solution converges below a set tolerance using the l-infinity norm\n",
    "            xfinal = x[0:i+2,:] # Save solutions for all the iterations \n",
    "            break # Exit out of the loop\n",
    "        elif max(abs(x[i+1,:] - x[i,:])) > tol_div: # Check if solution diverges beyond a set tolerance, if does then stop the loop and rerun with updated arguments\n",
    "            print('Solution is diverging \\n')\n",
    "            xfinal = 0\n",
    "            break # Exit out of the loop\n",
    "        # If solution has not converged/diverged beyond set tolerance\n",
    "        alpha = np.matmul(np.matmul(r.T,A),d)/np.matmul(np.matmul(d.T,A),d) # Calculate an appropriate scaling scalar\n",
    "        d = -r + alpha*d # Update search direction using the updated scalar alpha\n",
    "    \n",
    "    iterations = i+1\n",
    "    return xfinal, iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauss-Seidel Method (Splitting Method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gauss-Seidel Method:\n",
    "def gsm(A,b,x0, tol, tol_div, Nmax):\n",
    "    \"\"\"Dhruv Jain, 4 Oct, 2021, MSAAE Purdue\n",
    "    Obj: Use Gauss-Seidel Method (Splitting Method) to iteratively solve for a solution of linear system of equations of the form Ax = b, x is the solution \n",
    "    Args:\n",
    "        A: (n x n) matrix, Coefficeient matrix\n",
    "        b: (n x 1) vector, right-hand side vector\n",
    "        x0: (n x 1) vecotr, inital guess for x*\n",
    "        tol: ||x_n+1 - x_n||infi, step convergence tolerance \n",
    "        tol_div: ||x_n+1 - x_n||infi, step Divergence tolerance \n",
    "        Nmax: max allowable iterations for the iterative solver\n",
    "    Output:\n",
    "        xfinal: (niter x n), contains solved x for each iteration, last row is the converged solution\n",
    "        iterations: # iterations to converge  \n",
    "    \"\"\"\n",
    "    x = np.zeros((Nmax,len(b))) # Create x to save solution calculated after each iteration\n",
    "    x[0,:] = x0 # Assign 1st row as the Initial Guess\n",
    "    \n",
    "    for k in range(Nmax-1): # Run loop for set maximum iterations\n",
    "        for i in range(len(b)): # To individually compute each element of x, can be parallelized\n",
    "            x[k+1,i] = b[i]/A[i,i] # c term in x = Tx + x \n",
    "            for j in range(len(b)): # Incorporate the aij(k)*xj(k) for all j except j==i\n",
    "                if j < i: # For j<i, use the updated xi \n",
    "                    x[k+1,i] =  x[k+1,i]-(A[i,j]*x[k+1,j])/A[i,i] \n",
    "                if j > i: \n",
    "                    x[k+1,i] =  x[k+1,i]-(A[i,j]*x[k,j])/A[i,i]\n",
    "                    \n",
    "        if max(abs(x[k+1,:] - x[k,:])) < tol: # Check if solution converges below a set tolerance using the l-infinity norm\n",
    "            xfinal = x[0:k+2,:] # Save solutions for all the iterations \n",
    "            break # Exit out of the loop\n",
    "        elif max(abs(x[k+1,:] - x[k,:])) > tol_div: # Check if solution diverges beyond a set tolerance, if does then stop the loop and rerun with updated arguments\n",
    "            print('Solution is diverging \\n')\n",
    "            break # Exit out of the loop\n",
    "    \n",
    "    iterations = k+1\n",
    "    return xfinal, iterations    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jacobi Method (Splitting Method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Jacobi Method:\n",
    "def jacm(A,b,x0, tol, tol_div, Nmax):\n",
    "    \"\"\"Dhruv Jain, 4 Oct, 2021, MSAAE Purdue\n",
    "    Obj: Use Jacobi Method (Splitting Method) to iteratively solve for a solution of linear system of equations of the form Ax = b, x is the solution \n",
    "    Args:\n",
    "        A: (n x n) matrix, Coefficeient matrix\n",
    "        b: (n x 1) vector, right-hand side vector\n",
    "        x0: (n x 1) vecotr, inital guess for x*\n",
    "        tol: ||x_n+1 - x_n||infi, step convergence tolerance \n",
    "        tol_div: ||x_n+1 - x_n||infi, step Divergence tolerance\n",
    "        Nmax: max allowable iterations for the iterative solver\n",
    "    Output:\n",
    "        xfinal: (niter x n), contains solved x for each iteration, last row is the converged solution\n",
    "        iterations: # iterations to converge  \n",
    "    \"\"\"\n",
    "    x = np.zeros((Nmax,len(b))) # Create x to save solution calculated after each iteration\n",
    "    x[0,:] = x0 # Assign 1st row as the Initial Guess\n",
    "    \n",
    "    for k in range(Nmax-1): # Run loop for set maximum iterations\n",
    "        for i in range(len(b)): # To individually compute each element of x, can be parallelized\n",
    "            x[k+1,i] = b[i]/A[i,i] # c term in x = Tx + x \n",
    "            for j in range(len(b)): \n",
    "                if (i != j): # Incorporate the aij(k)*xj(k) for all j except j==i\n",
    "                    x[k+1,i] =  x[k+1,i]-(A[i,j]*x[k,j])/A[i,i]\n",
    "        \n",
    "        if max(abs(x[k+1,:] - x[k,:])) < tol: # Check if solution converges below a set tolerance using the l-infinity norm\n",
    "            xfinal = x[0:k+2,:] # Save solutions for all the iterations \n",
    "            break # Exit out of the loop\n",
    "        elif max(abs(x[k+1,:] - x[k,:])) > tol_div: # Check if solution diverges beyond a set tolerance, if does then stop the loop and rerun with updated arguments\n",
    "            print('Solution is diverging \\n')\n",
    "            break # Exit out of the loop\n",
    "    \n",
    "    iterations = k+1\n",
    "    return xfinal, iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SOR Method (Splitting Method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SOR Method:\n",
    "def sorm(A,b,x0, tol, tol_div, Nmax, W):\n",
    "    \"\"\"Dhruv Jain, 4 Oct, 2021, MSAAE Purdue\n",
    "    Obj: Use SOR Method (Splitting Method) to iteratively solve for a solution of linear system of equations of the form Ax = b, x is the solution \n",
    "    Args:\n",
    "        A: (n x n) matrix, Coefficeient matrix\n",
    "        b: (n x 1) vector, right-hand side vector\n",
    "        x0: (n x 1) vecotr, inital guess for x*\n",
    "        tol: ||x_n+1 - x_n||infi, step convergence tolerance \n",
    "        tol_div: ||x_n+1 - x_n||infi, step Divergence tolerance\n",
    "        Nmax: max allowable iterations for the iterative solver\n",
    "        W: (n_steps x 1) vector, relaxation parameter\n",
    "    Output:\n",
    "        xfinal: (niter x n), contains solved x for each iteration, last row is the converged solution\n",
    "        wf: Relaxation paramter list for which the solution converged\n",
    "        iterations: # iterations to converge  \n",
    "    \"\"\"\n",
    "    \n",
    "    xfinal = []\n",
    "    wf = []\n",
    "    iterations = []\n",
    "    for w in W:\n",
    "        x = np.zeros((Nmax,len(b))) # Create x to save solution calculated after each iteration\n",
    "        x[0,:] = x0 # Assign 1st row as the Initial Guess\n",
    "        check = 0\n",
    "        for k in range(Nmax-1): # Run loop for set maximum iterations\n",
    "            for i in range(len(b)): # To individually compute each element of x, can be parallelized\n",
    "                x[k+1,i] = w*b[i]/A[i,i] + (1-w)*x[k,i] # c term in x = Tx + x \n",
    "                for j in range(len(b)): \n",
    "                    if j < i:\n",
    "                        x[k+1,i] =  x[k+1,i]-w*(A[i,j]*x[k+1,j])/A[i,i] \n",
    "                    if j > i: # Incorporate the aij(k)*xj(k) for all j except j==i\n",
    "                        x[k+1,i] =  x[k+1,i]-w*(A[i,j]*x[k,j])/A[i,i]\n",
    "\n",
    "            if max(abs(x[k+1,:] - x[k,:])) < tol: # Check if solution converges below a set tolerance using the l-infinity norm\n",
    "                xfinal.append(x[0:k+1,:]) # Save solutions for all the iterations \n",
    "                wf.append(w)\n",
    "                iterations.append(k+1)\n",
    "                check = 1\n",
    "                break # Exit out of the loop\n",
    "            elif max(abs(x[k+1,:] - x[k,:])) > tol_div: # Check if solution diverges beyond a set tolerance, if does then stop the loop and rerun with updated arguments\n",
    "                print('Solution is diverging for w = ',round(w,2))\n",
    "                check = 1\n",
    "                break # Exit out of the loop\n",
    "    \n",
    "        if check == 0:\n",
    "            print('No solution for given convergence tolerance and w = ',round(w,2),' for set Nmax found.')\n",
    "            temp = np.zeros(len(b))\n",
    "            temp[:] = np.nan\n",
    "            x[1,:] = temp\n",
    "            xfinal.append(x[0:2,:])\n",
    "            wf.append(w)\n",
    "            iterations.append(k+1)\n",
    "    return xfinal, wf, iterations    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution of the system using Gaussian Elimination with No Pivoting: [0.13104672 0.19990201 0.25839976 0.30809799 0.35032038 0.38619147]\n",
      "Solution of the system using Gaussian Elimination with Scaled Partial Pivoting: [0.13104672 0.19990201 0.25839976 0.30809799 0.35032038 0.38619147]\n",
      "Solution of the system using Complex Gradient Method: [0.13107739 0.19976316 0.25835912 0.30799941 0.35048153 0.38637662]\n",
      "Solution of the system using Gauss-Sidel Method: [0.12867346 0.19636241 0.25474647 0.30507419 0.34831331 0.38526955]\n",
      "Solution of the system using Jacobi Method: [0.12862484 0.19624009 0.25377642 0.30421854 0.34717047 0.38472467]\n",
      "Solution of the system using SOR Method: [0.12810361 0.19551249 0.25386926 0.30434812 0.34783137 0.38504819]\n"
     ]
    }
   ],
   "source": [
    "# A and b matrix for the linear systems of equations\n",
    "A = np.zeros((6,6))\n",
    "for i in range(6):\n",
    "    A[i,i] = -88.824\n",
    "    if i < 5:\n",
    "        A[i,i+1] = 48.024\n",
    "    if i > 0:\n",
    "        A[i,i-1] = 40.8\n",
    "b = np.array([-2.04,0,0,0,0,-20.01])    \n",
    "x0 = np.zeros(6)\n",
    "\n",
    "#x0 = np.array([1,1,1,1])\n",
    "tol = 10e-4\n",
    "tol_div = 10e4\n",
    "Nmax = 50\n",
    "\n",
    "# Gaussian-Elimination with no Pivoting\n",
    "x_a = ge_lin_equ(A,b,opt=1)\n",
    "print('Solution of the system using Gaussian Elimination with No Pivoting:',x_a)\n",
    "\n",
    "# Gaussian-Elimination with Sacled Partial Pivoting\n",
    "x_b = ge_lin_equ(A,b,opt=2)\n",
    "print('Solution of the system using Gaussian Elimination with Scaled Partial Pivoting:',x_b)\n",
    "\n",
    "# Complex Gradient Method\n",
    "xfinal_1, iter_1 = cgm(A,b,x0,tol,tol_div, Nmax)\n",
    "print('Solution of the system using Complex Gradient Method:',xfinal_1[-1,:])\n",
    "\n",
    "# Gauss-Sidel\n",
    "xfinal_a, iterations = gsm(A,b,x0, tol, tol_div, Nmax)\n",
    "print('Solution of the system using Gauss-Sidel Method:',xfinal_a[-1,:])\n",
    "\n",
    "# Jacobi Method\n",
    "xfinal_b, iterations = jacm(A,b,x0, tol, tol_div, Nmax)\n",
    "print('Solution of the system using Jacobi Method:',xfinal_b[-1,:])\n",
    "\n",
    "# SOR Method\n",
    "W = [1]\n",
    "xfinal_c, wf, iterations = sorm(A,b,x0, tol, tol_div, Nmax, W)\n",
    "print('Solution of the system using SOR Method:',xfinal_c[0][-1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}