{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f6fd6083-f549-46d6-b1ce-484ab7dbd226",
   "metadata": {},
   "source": [
    "## MT09 - TP2 - Automne 2024\n",
    "### Elimination de Gauss et factorisation LU\n",
    "\n",
    "1. Ecrire une fonction python \n",
    "```\n",
    "x = solsup(U,b)\n",
    "```\n",
    "qui, étant donné $U$ une matrice triangulaire supérieure de taille $n\\times n$ et $b$ vecteur colonne de taille $n$, calcule $x$ solution de $Ux=b$.\n",
    "\n",
    "Application : utiliser la fonction `solsup()` pour déterminer la solution de \n",
    "\n",
    "$$\n",
    "\\begin{pmatrix} 1 & 2 & 3 \\\\ 0 & 4 & 8 \\\\ 0 & 0 & 5 \\end{pmatrix} x = \\begin{pmatrix} 6 \\\\ 16 \\\\ 15 \\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b770ac94-6d48-4363-a9d6-8648dd15ca61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x =  [[ 1.]\n",
      " [-2.]\n",
      " [ 3.]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([[ 6.],\n",
       "        [16.],\n",
       "        [15.]]),\n",
       " array([[ 6.],\n",
       "        [16.],\n",
       "        [15.]]),\n",
       " array([[0.],\n",
       "        [0.],\n",
       "        [0.]]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "# import matplotlib.pyplot as plt\n",
    "\n",
    "def solsup(U, b):\n",
    "    \"\"\" Résolution du système triangulaire supérieur Ux=b \"\"\"\n",
    "# A COMPLETER\n",
    "    n = np.size(b);\n",
    "    x = np.zeros((n,1));\n",
    "    for i in np.arange(0,n)[::-1]:\n",
    "        x[i] = (b[i] - U[i, i+1:] @ x[i+1:]) / U[i,i];\n",
    "    return x\n",
    "\n",
    "U = np.array([[1, 2, 3], [0, 4, 8], [0, 0, 5]], dtype=np.float64);\n",
    "b = np.array([[6, 16, 15]], dtype=np.float64).T;\n",
    "x = solsup(U, b);\n",
    "print('x = ', x)\n",
    "# Verif\n",
    "U@x, b, U@x - b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6031d8a-afbd-4f15-9d29-dd4f8de2a167",
   "metadata": {},
   "source": [
    "2. Ecrire une fonction python\n",
    "```\n",
    "U, e = trisup(A,b)\n",
    "```\n",
    "qui, étant donné $A$ matrice carrée de taille $n\\times n$ et $b$ vecteur colonne de taille $n$, déterminer par la méthode d'élimination de Gauss, quand c'est possible, la matrice $U$ triangulaire supérieure et le vecteur colonne $e$ tel que $Ax=b$\n",
    "$\\Leftrightarrow$ $Ux = e$.\n",
    "\n",
    "On rappelle l'algorithme d'élimination de Gauss:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04b59759-9488-4db8-a379-c4d6dee06077",
   "metadata": {},
   "source": [
    "### Algorithme de Gauss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d8d803a1-5b37-4915-a875-60ba4a438015",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "def trisup(A, b):\n",
    "# Méthode d'elimination de Gauss\n",
    "    Am = np.copy(A);\n",
    "    bm = np.copy(b) \n",
    "    n = np.size(b);\n",
    "    for k in np.arange(0,n):\n",
    "        if np.abs(Am[k,k])<1e-15:\n",
    "            print(\"ERROR : Pivot nul, k=\", k, '\\n');\n",
    "            sys.exit(1)\n",
    "        for i in np.arange(k+1, n):\n",
    "            c = Am[i,k] / Am[k,k];\n",
    "            bm[i] -= c * bm[k];\n",
    "            Am[i,k] = 0;\n",
    "            Am[i,k+1:n] -= c * Am[k,k+1:n];\n",
    "    if np.abs(Am[n-1,n-1])<1e-15:\n",
    "        print(\"ERROR : Pivot A_nn nul\\n\");\n",
    "        sys.exit(1);\n",
    "    return Am, bm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5309d9f0-d2ca-4d76-962a-69105d902acd",
   "metadata": {},
   "source": [
    "Application : utiliser la fonction `trisup()` pour déterminer un système d'équations équivalent au système\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix} 3 & 1 & 2 \\\\ 3 & 2 & 6 \\\\ 6 & 1 & -1 \\end{pmatrix} x = \\begin{pmatrix} 2 \\\\ 1 \\\\ 4 \\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "768fa758-c431-4e69-80e5-231c7297edf9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[ 3.  1.  2.]\n",
      " [ 0.  1.  4.]\n",
      " [ 0.  0. -1.]] \n",
      "e= [[ 2.]\n",
      " [-1.]\n",
      " [-1.]]\n",
      "x= [[ 1.66666667]\n",
      " [-5.        ]\n",
      " [ 1.        ]]\n",
      "Ax-b= [[0.]\n",
      " [0.]\n",
      " [0.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[3, 1, 2], [3, 2, 6], [6, 1, -1]], dtype=np.float64);\n",
    "b = np.array([[2, 1, 4]], dtype=np.float64).T;\n",
    "U, e = trisup(A,b);\n",
    "print(\"U=\", U, \"\\ne=\", e);\n",
    "x = solsup(U, e);\n",
    "print(\"x=\",x);\n",
    "# Verification\n",
    "print( \"Ax-b=\", A@x-b);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08e6b6f0-ce46-447e-8d73-8506f9651ec7",
   "metadata": {},
   "source": [
    "3. Ecrire une fonction \n",
    "```\n",
    "x = resolG(A, b)\n",
    "```\n",
    "qui, étant donné $A$ matrice carrée de taille $n\\times n$ et $b$ vecteur colonne de taille $n$, calcule et retourne \n",
    "(quand c'est possible) la solution de $Ax=b$ en utilisant la méthode d'élimination de Gauss. Cette fonction doit\n",
    "utiliser les 2 fonctions précédentes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "148ad5ba-57bd-40b2-8a42-d213f5aca3cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def resolG(A,b):\n",
    "    U, e = trisup(A,b);\n",
    "    x = solsup(U, e);\n",
    "    return x;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4aec09e1-8df2-4b78-80fd-526dbe8d4f83",
   "metadata": {},
   "source": [
    "Application : utiliser la fonction `resolG()`pour déterminer la solution de\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix} 1 & 2 & 3 \\\\ 5 & 2 & 1 \\\\ 3 & -1 & 1 \\end{pmatrix} x = \\begin{pmatrix} 5 \\\\ 5 \\\\ 6 \\end{pmatrix},\n",
    "\\quad \\text{puis}\\quad\n",
    "\\begin{pmatrix} 2 & 1 & 5 \\\\ 1 & 2 & 4 \\\\ 3 & 4 & 10 \\end{pmatrix} x = \\begin{pmatrix} 5 \\\\ 5 \\\\ 6 \\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cb3d18d2-1c4a-422d-b9cd-2cf90df62e2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 1., -1.,  2.]), array([5., 5., 6.]), array([5., 5., 6.]))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1, 2, 3], [5, 2, 1], [3, -1, 1]], dtype=np.float64);\n",
    "b = np.array([5, 5, 6], dtype=np.float64);\n",
    "x = resolG(A, b);\n",
    "x, A@x, b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "940b8d0b-af75-4178-9d0b-bb7c6cda392e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ERROR : Pivot nul, k= 2 \n",
      "\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "1",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[1;31mSystemExit\u001b[0m\u001b[1;31m:\u001b[0m 1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fdevu\\anaconda3\\Lib\\site-packages\\IPython\\core\\interactiveshell.py:3534: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.\n",
      "  warn(\"To exit: use 'exit', 'quit', or Ctrl-D.\", stacklevel=1)\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[2, 1, 5], [1, 2, 4], [3, 4, 10]], dtype=np.float64);\n",
    "b = np.array([5, 5, 6], dtype=np.float64);\n",
    "x = resolG(A, b);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c321205a-b078-42b3-92c7-106de867daae",
   "metadata": {},
   "source": [
    "4. Ecrire une fonction\n",
    "```\n",
    "L, U = LU(A)\n",
    "```\n",
    "qui, étant donné $A$ de taille $n\\times n$, détermine une matrice triangulaire supérieure $U$ et une matrice \n",
    "triangulaire inférieure $L$ vérifiant $A=LU$.\n",
    "\n",
    "NB : on poura en particulier reprendre la fonction `trisup()` que l''on modifiera pour obtenir la fonction `LU()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2ddce40e-bffa-4e96-842d-5ece0a835e37",
   "metadata": {},
   "outputs": [],
   "source": [
    "def LU(A):\n",
    "    U = np.copy(A);\n",
    "    n = np.size(A[0,:]);\n",
    "    L = np.eye(n, dtype=np.float64);\n",
    "    for k in np.arange(0,n):\n",
    "        if np.abs(U[k,k])<1e-15:\n",
    "            print(\"Pivot nul\\n\");\n",
    "            break\n",
    "        for i in np.arange(k+1, n):\n",
    "            c = U[i,k] / U[k,k];\n",
    "            L[i,k] = c;\n",
    "            U[i,k] = 0;\n",
    "            U[i,k+1:n] -= c * U[k,k+1:n];\n",
    "    if np.abs(U[n-1,n-1])<1e-15:\n",
    "        print(\"Pivot A_nn nul\\n\");\n",
    "    return L, U"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f57d63c-7cd8-4492-9b3c-7dd3d0c279ff",
   "metadata": {},
   "source": [
    "Application : utiliser la fonction `LU()` pour déterminer $L$ et $U$ dans le cas où\n",
    "\n",
    "$$\n",
    "A = \\begin{pmatrix} 3 & 1 & 2 \\\\ 3 & 2 & 6 \\\\ 6 & 1 & -1 \\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "43887836-0299-4450-a74f-b771fd6f84a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 1.,  0.,  0.],\n",
       "        [ 1.,  1.,  0.],\n",
       "        [ 2., -1.,  1.]]),\n",
       " array([[ 3.,  1.,  2.],\n",
       "        [ 0.,  1.,  4.],\n",
       "        [ 0.,  0., -1.]]),\n",
       " array([[ 3.,  1.,  2.],\n",
       "        [ 3.,  2.,  6.],\n",
       "        [ 6.,  1., -1.]]))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[3, 1, 2], [3, 2, 6], [6, 1, -1]], dtype=np.float64);\n",
    "L, U = LU(A);\n",
    "L, U, L@U"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "907d99c0-a738-4a39-a3e5-c1f4d5d7f750",
   "metadata": {},
   "source": [
    "5. Ecrire une fonction \n",
    "```\n",
    "x = solinf(L,b)\n",
    "```\n",
    "qui, étant donné $L$ matrice triangulaire inférieure de taille $n\\times n$ et $b$ vecteur colonne de taille $n$, calcule $x$ solution de $Lx=b$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f87ff95f-935c-4bce-b206-83801e5be620",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solinf(L, b):\n",
    "    n = np.size(b);\n",
    "    x = np.zeros(n);\n",
    "    for i in np.arange(0, n):\n",
    "        x[i] = (b[i] - L[i, 0:i] @ x[0:i]) / L[i,i];\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72cdec62-7be3-4613-829d-123e180e7317",
   "metadata": {},
   "source": [
    "Application : utiliser la fonction `solinf()` pour déterminer la solution de\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix} 1 & 0 & 0 \\\\ 2 & 3 & 0 \\\\ 1 & 4 & -1 \\end{pmatrix} x = \\begin{pmatrix} 1 \\\\ 8 \\\\ 10 \\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "694129ed-a678-4020-9a99-283387bc01b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 1.,  2., -1.]), array([ 1.,  8., 10.]), array([ 1.,  8., 10.]))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = np.array([[1, 0, 0], [2, 3, 0], [1, 4, -1]], dtype=np.float64);\n",
    "b = np.array([1, 8, 10], dtype=np.float64);\n",
    "x = solinf(L, b);\n",
    "x, L@x, b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19ed15df-2741-4a02-88ea-2536d99cbb6a",
   "metadata": {},
   "source": [
    "6. Ecrire une fonction \n",
    "```\n",
    "x = resolLU(A,b)\n",
    "```\n",
    "qui calcule, quand c'est possible, la solution de $Ax=b$ en utilisant la factorisation $A=LU$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "748cf5f5-7b30-4403-8df0-0683326b8d56",
   "metadata": {},
   "outputs": [],
   "source": [
    "def resolLU(A, b):\n",
    "    L, U = LU(A);\n",
    "    y = solinf(L, b);\n",
    "    x = solsup(U, y);\n",
    "    return x;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5acf9564-0db0-4b7c-8456-82ef7cdb8616",
   "metadata": {},
   "source": [
    "Application : utiliser la fonction `resolLU()` pour déterminer la solution de\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix} 1 & 2 & 3 \\\\ 5 & 2 & 1 \\\\ 3 & -1 & 1 \\end{pmatrix} x = \\begin{pmatrix} 5 \\\\ 5 \\\\ 6 \\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f55ede69-48d1-440e-ac59-4523eb5ee46c",
   "metadata": {},
   "source": [
    "7. Ecrire une fonction\n",
    "```\n",
    "B = inverse(A)\n",
    "```\n",
    "qui calcule la matrice $B$ inverse de $A$ en utilisant `LU()`, `solinf()` et `solsup()` et en minimisant le \n",
    "nombre d'opérations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "79c06b48-1ac5-4b27-b7b6-fdcfba4153b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def inverse(A):\n",
    "    n = np.size(A[0,:]);\n",
    "    Id = np.eye(n, dtype=np.float64);\n",
    "    B = np.zeros((n,n), dtype=np.float64);\n",
    "    L, U = LU(A);\n",
    "    for i in np.arange(0, n):\n",
    "        y = solinf(L, Id[:, i]);\n",
    "        B[:,i] = solsup(U, y);\n",
    "    return B;    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8abace67-6cb6-403e-b1ce-77f034f7537f",
   "metadata": {},
   "source": [
    "Application : utiliser la fonction `inverse()` pour déterminer la matrice $B$ inverse de\n",
    "\n",
    "$$\n",
    "A = \\begin{pmatrix} 1 & 2 & 3 \\\\ 5 & 2 & 1 \\\\ 3 & -1 & 1 \\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "6288f5a8-de7e-4c0c-9aa0-c02f24e6a281",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-0.08823529,  0.14705882,  0.11764706],\n",
       "        [ 0.05882353,  0.23529412, -0.41176471],\n",
       "        [ 0.32352941, -0.20588235,  0.23529412]]),\n",
       " array([[ 1.00000000e+00, -2.22044605e-16, -3.33066907e-16],\n",
       "        [ 1.11022302e-16,  1.00000000e+00,  5.55111512e-17],\n",
       "        [ 5.55111512e-17,  5.55111512e-17,  1.00000000e+00]]),\n",
       " 4.80740671595891e-16)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1, 2, 3], [5, 2, 1], [3, -1, 1]], dtype=np.float64);\n",
    "B = inverse(A);\n",
    "C = B@A;\n",
    "B, C, la.norm(C - np.eye(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "5bf5f3ab-c527-4768-a35c-f40ec3b2dd2b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "C[np.abs(C)<1e-15] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b46e8c89-f9be-4f9b-aa37-6bfae33a568c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 0., 1.]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b59cd4e5-300e-4ff2-ac81-e2cd43af8160",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  },
  "toc-autonumbering": false,
  "toc-showcode": true,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}