Doolitles Method.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Doolittle(M):\n",
    "    M = M*1.0\n",
    "    \n",
    "    s = len(M)\n",
    "    \n",
    "    I = np.identity(len(M))\n",
    "\n",
    "    for h in range(s-1):\n",
    "        for i in range(h,s-1):\n",
    "            if M[i+1][h] != 0:\n",
    "                num = float(M[i+1][h]/M[h][h])\n",
    "                I[i+1][h] = num\n",
    "                for j in range(s):\n",
    "                    temp = M[i+1][j]\n",
    "                    M[i+1][j] = temp - (num*M[h][j])\n",
    "    return I,M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def triangular(M):\n",
    "    n = len(M)\n",
    "    if n != len(M[0]):\n",
    "        return 0\n",
    "    upper = True\n",
    "    for i in range(n):\n",
    "        for j in range(i):\n",
    "            if np.round(M[i][j], 3) !=0:\n",
    "                upper = False\n",
    "    lower = True\n",
    "    for i in range(n):\n",
    "        a = n-i-1\n",
    "        for j in range(a):\n",
    "            b = n-j-1\n",
    "            if np.round(M[i][b], 3) !=0:\n",
    "                lower = False\n",
    "    if lower:\n",
    "        return 1\n",
    "    elif upper:\n",
    "        return 2\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solver(L,U,b):\n",
    "    # Forward substitution\n",
    "    if triangular(L) != 1:\n",
    "        return \"L is not a lower triangular matrix\"\n",
    "    n = len(b)\n",
    "    z = np.zeros(n)\n",
    "    for i in range(n):\n",
    "        z[i] = b[i]\n",
    "        for j in range(i):\n",
    "            z[i] = z[i]-L[i][j]*z[j]    \n",
    "        z[i] = z[i]/L[i][i]\n",
    "    print(\"z = \", z)\n",
    "    \n",
    "    if triangular(U) != 2:\n",
    "        return \"U is not an upper triangular matrix\"\n",
    "    \n",
    "    # Backward substitution\n",
    "    y = np.zeros(n)\n",
    "    for i in range(n):\n",
    "        a = n-1-i\n",
    "        y[a] = z[a]\n",
    "        for j in range(i):\n",
    "            b = n-1-j\n",
    "            y[a] = y[a]-U[a][b]*y[b]    \n",
    "        y[a] = y[a]/U[a][a]\n",
    "    print(\"x = \",y)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solvin a system of equations using LU decomposition\n",
    "A = np.array([[10,-1,2,0], \n",
    "              [-1,11,-1,3],\n",
    "              [2,-1,10,-1],\n",
    "              [0,3,-1,8]])\n",
    "A = A*1.0 # to make it a float\n",
    "b = np.array([6,25,-11,15])\n",
    "b= b*1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.          0.          0.          0.        ]\n",
      " [-0.1         1.          0.          0.        ]\n",
      " [ 0.2        -0.0733945   1.          0.        ]\n",
      " [ 0.          0.27522936 -0.08173077  1.        ]]\n",
      "[[ 1.00000000e+01 -1.00000000e+00  2.00000000e+00  0.00000000e+00]\n",
      " [ 0.00000000e+00  1.09000000e+01 -8.00000000e-01  3.00000000e+00]\n",
      " [ 0.00000000e+00  0.00000000e+00  9.54128440e+00 -7.79816514e-01]\n",
      " [ 0.00000000e+00 -4.44089210e-16  0.00000000e+00  7.11057692e+00]]\n"
     ]
    }
   ],
   "source": [
    "L,U = Doolittle(A)\n",
    "print(L)\n",
    "print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z =  [  6.          25.6        -10.32110092   7.11057692]\n",
      "x =  [ 1.  2. -1.  1.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 1.,  2., -1.,  1.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solver(L,U,b)"
   ]
  },
  {
   "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.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}