{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f6fd6083-f549-46d6-b1ce-484ab7dbd226",
   "metadata": {},
   "source": [
    "## MT09 - TP2 - Automne 2025\n",
    "### Élimination de Gauss et factorisation LU\n",
    "\n",
    "1. Écrire 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": 25,
   "id": "b770ac94-6d48-4363-a9d6-8648dd15ca61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x =  [ 1. -2.  3.]\n",
      "Ux-b =  [0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "# import matplotlib.pyplot as plt\n",
    "\n",
    "def solsup_nonvectorisee(U,b):\n",
    "    \"\"\" Résolution du système triangulaire supérieur Ux=b \"\"\"\n",
    "    # Version non vectorisée\n",
    "    n = np.size(b);\n",
    "    x = np.zeros(n);\n",
    "    for i in np.arange(0,n)[::-1]:\n",
    "        sum = 0.\n",
    "        for j in range(i+1, n):\n",
    "            sum += U[i,j]*x[j]\n",
    "        x[i] = (b[i] - sum) / U[i,i];\n",
    "    return x\n",
    "\n",
    "def solsup(U, b):\n",
    "    \"\"\" Résolution du système triangulaire supérieur Ux=b \"\"\"\n",
    "    # Version vectorisee\n",
    "    n = np.size(b);\n",
    "    x = np.zeros(n);\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);\n",
    "x = solsup_nonvectorisee(U, b);\n",
    "print('x = ', x)\n",
    "# Verif\n",
    "print(\"Ux-b = \", U@x - b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6031d8a-afbd-4f15-9d29-dd4f8de2a167",
   "metadata": {},
   "source": [
    "2. Écrire 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": 31,
   "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",
    "            raise SystemExit(\"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",
    "            #for j in range(k+1,n):\n",
    "            #    Am[i,j] -= c * Am[k,j];\n",
    "    if np.abs(Am[n-1,n-1])<1e-15:\n",
    "        raise SystemExit(\"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": 34,
   "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",
      "\n",
      "e= [ 2. -1. -1.] \n",
      "\n",
      "x= [ 1.66666667 -5.          1.        ] \n",
      "\n",
      "Ax-b=\n",
      " [0. 0. 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);\n",
    "U, e = trisup(A,b);\n",
    "print(\"U=\", U, \"\\n\\ne=\", e, '\\n');\n",
    "x = solsup(U, e);\n",
    "print(\"x=\",x, '\\n');\n",
    "# Vérification\n",
    "print( \"Ax-b=\\n\", A@x-b);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08e6b6f0-ce46-447e-8d73-8506f9651ec7",
   "metadata": {},
   "source": [
    "3. Écrire 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": 72,
   "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": 75,
   "id": "cb3d18d2-1c4a-422d-b9cd-2cf90df62e2b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U = \n",
      " [[  1.     2.     3.  ]\n",
      " [  0.    -8.   -14.  ]\n",
      " [  0.     0.     4.25]] \n",
      "e = \n",
      " [  5.  -20.    8.5] \n",
      "\n",
      "b =  [5. 5. 6.] \n",
      "x =  [ 1. -1.  2.] \n",
      "Ax-b =  [0. 0. 0.]\n"
     ]
    }
   ],
   "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",
    "U, e = trisup(A, b)\n",
    "print(\"U = \\n\", U, \"\\ne = \\n\", e, \"\\n\")\n",
    "\n",
    "x = resolG(A, b);\n",
    "print('b = ', b, '\\nx = ', x, \"\\nAx-b = \", A@x-b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "940b8d0b-af75-4178-9d0b-bb7c6cda392e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "det(A) =  -1.0177044392397329e-15\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "('ERROR : Pivot nul, k=', 2, '\\n')",
     "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 ('ERROR : Pivot nul, k=', 2, '\\n')\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",
    "import numpy.linalg as la\n",
    "print('det(A) = ', la.det(A))\n",
    "x = resolG(A, b);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "2b51522e-9f29-404f-9708-dbf7ee7d68ab",
   "metadata": {},
   "outputs": [
    {
     "ename": "SystemExit",
     "evalue": "('ERROR : Pivot nul, k=', 2, '\\n')",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mSystemExit\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[77], line 5\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01mnumpy\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mlinalg\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01mla\u001b[39;00m\n\u001b[0;32m      4\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdet(A) = \u001b[39m\u001b[38;5;124m'\u001b[39m, la\u001b[38;5;241m.\u001b[39mdet(A))\n\u001b[1;32m----> 5\u001b[0m x \u001b[38;5;241m=\u001b[39m resolG(A, b)\n",
      "Cell \u001b[1;32mIn[72], line 2\u001b[0m, in \u001b[0;36mresolG\u001b[1;34m(A, b)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21mresolG\u001b[39m(A,b):\n\u001b[1;32m----> 2\u001b[0m     U, e \u001b[38;5;241m=\u001b[39m trisup(A,b);\n\u001b[0;32m      3\u001b[0m     x \u001b[38;5;241m=\u001b[39m solsup(U, e);\n\u001b[0;32m      4\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m x\n",
      "Cell \u001b[1;32mIn[31], line 9\u001b[0m, in \u001b[0;36mtrisup\u001b[1;34m(A, b)\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m np\u001b[38;5;241m.\u001b[39marange(\u001b[38;5;241m0\u001b[39m,n):\n\u001b[0;32m      8\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m np\u001b[38;5;241m.\u001b[39mabs(Am[k,k])\u001b[38;5;241m<\u001b[39m\u001b[38;5;241m1e-15\u001b[39m:\n\u001b[1;32m----> 9\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mSystemExit\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mERROR : Pivot nul, k=\u001b[39m\u001b[38;5;124m\"\u001b[39m, k, \u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m'\u001b[39m);\n\u001b[0;32m     10\u001b[0m         sys\u001b[38;5;241m.\u001b[39mexit(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m     11\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m np\u001b[38;5;241m.\u001b[39marange(k\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m, n):\n",
      "\u001b[1;31mSystemExit\u001b[0m: ('ERROR : Pivot nul, k=', 2, '\\n')"
     ]
    }
   ],
   "source": [
    "%tb"
   ]
  },
  {
   "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": 82,
   "id": "2ddce40e-bffa-4e96-842d-5ece0a835e37",
   "metadata": {},
   "outputs": [],
   "source": [
    "def LU(A):\n",
    "    \"\"\"Factorisation LU : A=LU\"\"\"\n",
    "    U = np.copy(A);\n",
    "    m, n = A.shape;\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",
    "            raise SystemExit(\"LU: Oops, 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",
    "        raise SystemExit(\"LU: Oops, 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": 85,
   "id": "43887836-0299-4450-a74f-b771fd6f84a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L=\n",
      " [[ 1.  0.  0.]\n",
      " [ 1.  1.  0.]\n",
      " [ 2. -1.  1.]] \n",
      "U=\n",
      " [[ 3.  1.  2.]\n",
      " [ 0.  1.  4.]\n",
      " [ 0.  0. -1.]]\n",
      "LU-A = \n",
      " [[0. 0. 0.]\n",
      " [0. 0. 0.]\n",
      " [0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[3, 1, 2], [3, 2, 6], [6, 1, -1]], dtype=np.float64);\n",
    "L, U = LU(A);\n",
    "print('L=\\n', L, '\\nU=\\n', U)\n",
    "print('LU-A = \\n', L@U-A)"
   ]
  },
  {
   "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": 88,
   "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": 91,
   "id": "694129ed-a678-4020-9a99-283387bc01b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x =  [ 1.  2. -1.] , Lx-b =  [0. 0. 0.]\n"
     ]
    }
   ],
   "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",
    "print('x = ', x, ', Lx-b = ', 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": 106,
   "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 L, U, y, 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": "code",
   "execution_count": 110,
   "id": "77a5c595-7969-4420-80f2-d75a3190857b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L =  [[1.    0.    0.   ]\n",
      " [5.    1.    0.   ]\n",
      " [3.    0.875 1.   ]] \n",
      "\n",
      "U =  [[  1.     2.     3.  ]\n",
      " [  0.    -8.   -14.  ]\n",
      " [  0.     0.     4.25]] \n",
      "\n",
      "y =  [  5.  -20.    8.5] \n",
      "\n",
      "x =  [ 1. -1.  2.]\n",
      "x =  [ 1. -1.  2.] , Ax-b =  [0. 0. 0.]\n"
     ]
    }
   ],
   "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",
    "L,U,y,x = resolLU(A, b)\n",
    "print('L = ',L,'\\n\\nU = ', U,'\\n\\ny = ',y,'\\n\\nx = ',x)\n",
    "print('x = ', x, ', Ax-b = ', A@x-b )"
   ]
  },
  {
   "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": 126,
   "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": 129,
   "id": "6288f5a8-de7e-4c0c-9aa0-c02f24e6a281",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = \n",
      " [[ 1.  2.  3.]\n",
      " [ 5.  2.  1.]\n",
      " [ 3. -1.  1.]] \n",
      "B = \n",
      " [[-0.08823529  0.14705882  0.11764706]\n",
      " [ 0.05882353  0.23529412 -0.41176471]\n",
      " [ 0.32352941 -0.20588235  0.23529412]] \n",
      "C = AB = \n",
      " [[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n"
     ]
    }
   ],
   "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",
    "eps = 2 * np.finfo(np.float64).eps # 2 = facteur de securité\n",
    "C[np.abs(C) < eps] = 0 # Pour enlever le bruit numerique \n",
    "print(\"A = \\n\", A, \"\\nB = \\n\", B, \"\\nC = AB = \\n\", C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "4b6dffec-a4f7-4356-900a-e533be0ce9a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.220446049250313e-16"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.finfo(np.float64).eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "fb08f39b-750d-4848-9000-6b761ed185c0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[1;31mInit signature:\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfinfo\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
       "\u001b[1;31mDocstring:\u001b[0m     \n",
       "finfo(dtype)\n",
       "\n",
       "Machine limits for floating point types.\n",
       "\n",
       "Attributes\n",
       "----------\n",
       "bits : int\n",
       "    The number of bits occupied by the type.\n",
       "dtype : dtype\n",
       "    Returns the dtype for which `finfo` returns information. For complex\n",
       "    input, the returned dtype is the associated ``float*`` dtype for its\n",
       "    real and complex components.\n",
       "eps : float\n",
       "    The difference between 1.0 and the next smallest representable float\n",
       "    larger than 1.0. For example, for 64-bit binary floats in the IEEE-754\n",
       "    standard, ``eps = 2**-52``, approximately 2.22e-16.\n",
       "epsneg : float\n",
       "    The difference between 1.0 and the next smallest representable float\n",
       "    less than 1.0. For example, for 64-bit binary floats in the IEEE-754\n",
       "    standard, ``epsneg = 2**-53``, approximately 1.11e-16.\n",
       "iexp : int\n",
       "    The number of bits in the exponent portion of the floating point\n",
       "    representation.\n",
       "machep : int\n",
       "    The exponent that yields `eps`.\n",
       "max : floating point number of the appropriate type\n",
       "    The largest representable number.\n",
       "maxexp : int\n",
       "    The smallest positive power of the base (2) that causes overflow.\n",
       "min : floating point number of the appropriate type\n",
       "    The smallest representable number, typically ``-max``.\n",
       "minexp : int\n",
       "    The most negative power of the base (2) consistent with there\n",
       "    being no leading 0's in the mantissa.\n",
       "negep : int\n",
       "    The exponent that yields `epsneg`.\n",
       "nexp : int\n",
       "    The number of bits in the exponent including its sign and bias.\n",
       "nmant : int\n",
       "    The number of bits in the mantissa.\n",
       "precision : int\n",
       "    The approximate number of decimal digits to which this kind of\n",
       "    float is precise.\n",
       "resolution : floating point number of the appropriate type\n",
       "    The approximate decimal resolution of this type, i.e.,\n",
       "    ``10**-precision``.\n",
       "tiny : float\n",
       "    An alias for `smallest_normal`, kept for backwards compatibility.\n",
       "smallest_normal : float\n",
       "    The smallest positive floating point number with 1 as leading bit in\n",
       "    the mantissa following IEEE-754 (see Notes).\n",
       "smallest_subnormal : float\n",
       "    The smallest positive floating point number with 0 as leading bit in\n",
       "    the mantissa following IEEE-754.\n",
       "\n",
       "Parameters\n",
       "----------\n",
       "dtype : float, dtype, or instance\n",
       "    Kind of floating point or complex floating point\n",
       "    data-type about which to get information.\n",
       "\n",
       "See Also\n",
       "--------\n",
       "iinfo : The equivalent for integer data types.\n",
       "spacing : The distance between a value and the nearest adjacent number\n",
       "nextafter : The next floating point value after x1 towards x2\n",
       "\n",
       "Notes\n",
       "-----\n",
       "For developers of NumPy: do not instantiate this at the module level.\n",
       "The initial calculation of these parameters is expensive and negatively\n",
       "impacts import times.  These objects are cached, so calling ``finfo()``\n",
       "repeatedly inside your functions is not a problem.\n",
       "\n",
       "Note that ``smallest_normal`` is not actually the smallest positive\n",
       "representable value in a NumPy floating point type. As in the IEEE-754\n",
       "standard [1]_, NumPy floating point types make use of subnormal numbers to\n",
       "fill the gap between 0 and ``smallest_normal``. However, subnormal numbers\n",
       "may have significantly reduced precision [2]_.\n",
       "\n",
       "This function can also be used for complex data types as well. If used,\n",
       "the output will be the same as the corresponding real float type\n",
       "(e.g. numpy.finfo(numpy.csingle) is the same as numpy.finfo(numpy.single)).\n",
       "However, the output is true for the real and imaginary components.\n",
       "\n",
       "References\n",
       "----------\n",
       ".. [1] IEEE Standard for Floating-Point Arithmetic, IEEE Std 754-2008,\n",
       "       pp.1-70, 2008, http://www.doi.org/10.1109/IEEESTD.2008.4610935\n",
       ".. [2] Wikipedia, \"Denormal Numbers\",\n",
       "       https://en.wikipedia.org/wiki/Denormal_number\n",
       "\n",
       "Examples\n",
       "--------\n",
       ">>> np.finfo(np.float64).dtype\n",
       "dtype('float64')\n",
       ">>> np.finfo(np.complex64).dtype\n",
       "dtype('float32')\n",
       "\u001b[1;31mFile:\u001b[0m           c:\\users\\fdevu\\anaconda3\\lib\\site-packages\\numpy\\__init__.py\n",
       "\u001b[1;31mType:\u001b[0m           type\n",
       "\u001b[1;31mSubclasses:\u001b[0m     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "?np.finfo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "id": "ceb57aaf-8241-4018-808b-551de25f9fd3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False,  True,  True],\n",
       "       [ True, False,  True],\n",
       "       [ True,  True, False]])"
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.abs(C) < eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93d0ec2e-ef17-46d1-8359-17522bac53f9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "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
}