{
 "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": null,
   "id": "b770ac94-6d48-4363-a9d6-8648dd15ca61",
   "metadata": {},
   "outputs": [],
   "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",
    "    n = np.size(b);\n",
    "    x = np.zeros((n,1));\n",
    "    # VOIR POLY\n",
    "    # A COMPLETER ...\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",
    "# Verification\n",
    "print(U@x, b)\n",
    "print(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": null,
   "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",
    "        # VOIR COURS\n",
    "        # A COMPLETER\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": null,
   "id": "768fa758-c431-4e69-80e5-231c7297edf9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "id": "148ad5ba-57bd-40b2-8a42-d213f5aca3cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def resolG(A,b):\n",
    "    # A COMPLETER\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": null,
   "id": "cb3d18d2-1c4a-422d-b9cd-2cf90df62e2b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "940b8d0b-af75-4178-9d0b-bb7c6cda392e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "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",
    "        # A COMPLETER\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": null,
   "id": "43887836-0299-4450-a74f-b771fd6f84a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[3, 1, 2], [3, 2, 6], [6, 1, -1]], dtype=np.float64);\n",
    "L, U = LU(A);\n",
    "print(L, U, A-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": null,
   "id": "f87ff95f-935c-4bce-b206-83801e5be620",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solinf(L, b):\n",
    "    # A COMPLETER\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": null,
   "id": "694129ed-a678-4020-9a99-283387bc01b3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "id": "748cf5f5-7b30-4403-8df0-0683326b8d56",
   "metadata": {},
   "outputs": [],
   "source": [
    "def resolLU(A, b):\n",
    "    # A COMPLETER\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": null,
   "id": "79c06b48-1ac5-4b27-b7b6-fdcfba4153b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def inverse(A):\n",
    "    n = np.size(A[0,:]);\n",
    "    # A COMPLETER\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": null,
   "id": "6288f5a8-de7e-4c0c-9aa0-c02f24e6a281",
   "metadata": {},
   "outputs": [],
   "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",
    "print( B, C, la.norm(C - np.eye(3)) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bf5f3ab-c527-4768-a35c-f40ec3b2dd2b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b46e8c89-f9be-4f9b-aa37-6bfae33a568c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "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
}