{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\newcommand{\\nr}[1]{\\|#1\\|}\n",
    "\\newcommand{\\RR}{\\mathbb{R}}\n",
    "\\newcommand{\\N}{\\mathbb{N}}\n",
    "$$\n",
    "### MEU352 2023/2024 - Analyse numérique matricielle et optimisation\n",
    "\n",
    "# TP1 - Résolution de systèmes linéaires triangulaires. Méthode de Gauss.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 0. Manipulation de vecteurs et de matrices.\n",
    "\n",
    "On aura besoin des modules de python ``numpy`` et ``matplotlib.pyplot``. On peut les charger en exécutant les commandes\n",
    "\n",
    "``import numpy as np``\n",
    "\n",
    "``import matplotlib.pyplot as plt``\n",
    "\n",
    "(on désignera alors le module ``numpy`` par ``np`` et ``matplotlib.pyplot`` par ``plt``. \n",
    "\n",
    "**Q1.** Executez les commandes suivantes et affichez le résultat. Essayez de comprendre ce que vous avez obtenu.\n",
    "\n",
    "``\n",
    "u = np.array([1,2,3,4,5])\n",
    "v = np.array([[1,2,3,4,5]])\n",
    "su=u.shape\n",
    "sv=v.shape\n",
    "ut = np.transpose(u)\n",
    "vt = np.transpose(v)\n",
    "vt2 = np.array([[1],[2],[3],[4],[5]])\n",
    "A = np.array([[1,2,0,0,0],[0,0,2,3,1],[0,0,0,2,2],[0,0,0,0,1],[1,1,1,0,0]])\n",
    "B = np.array([[1,2,3,4,5],[2,3,4,5,6],[3,4,5,6,7],[4,5,6,7,8],[5,6,7,8,9]])\n",
    "d=np.diag(A)\n",
    "dd=np.array([np.diag(A)])\n",
    "dt=np.transpose(d)\n",
    "ddt=np.transpose(dd)\n",
    "Ad=np.diag(np.diag(A))``\n",
    "\n",
    "**Q2.** Même question pour les commandes suivantes.\n",
    "\n",
    "``u*v,  u*vt, vt*u, u/v, u/vt, v/v, v/vt, np.vdot(u,v), np.vdot(u,vt)``\n",
    "\n",
    "``A*B, np.dot(A,B)``\n",
    "\n",
    "``np.dot(A,u), np.dot(A,v), np.dot(v,A), np.dot(A,vt), np.linalg.inv(A), np.dot(np.linalg(inv(A)),A))``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "u = np.array([1, 2, 3, 4, 5])\n",
    "v = np.array([[1, 2, 3, 4, 5]])\n",
    "su = u.shape\n",
    "sv = v.shape\n",
    "ut = np.transpose(u)\n",
    "vt = np.transpose(v)\n",
    "vt2 = np.array([[1], [2], [3], [4], [5]])\n",
    "A = np.array(\n",
    "    [\n",
    "        [1, 2, 0, 0, 0],\n",
    "        [0, 2, 0, 0, 0],\n",
    "        [0, 0, 3, 0, 0],\n",
    "        [0, 0, 0, 4, 0],\n",
    "        [0, 0, 0, 0, 5],\n",
    "    ],\n",
    ")\n",
    "B = np.array(\n",
    "    [\n",
    "        [1, 2, 3, 4, 5],\n",
    "        [2, 3, 4, 5, 6],\n",
    "        [3, 4, 5, 6, 7],\n",
    "        [4, 5, 6, 7, 8],\n",
    "        [5, 6, 7, 8, 9],\n",
    "    ],\n",
    ")\n",
    "d = np.diag(A)\n",
    "dd = np.array([np.diag(A)])\n",
    "dt = np.transpose(d)\n",
    "ddt = np.transpose(dd)\n",
    "Ad = np.diag(np.diag(A))\n",
    "\n",
    "print(np.dot(np.linalg.inv(A), A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 1.  Résolution d'un système linéaire triangulaire.\n",
    "\n",
    "Soit $A\\in\\mathcal{M}_n(\\RR)$ une matrice triangulaire inférieure inversible, de taille $n\\in\\N,$ et $b\\in\\RR^n$. Comme $A$ est triangulaire inférieure, on peut résoudre le système $Ax=b$ par une technique dite de *descente* : la solution $x=(x_1,\\dots,x_n)$ est obtenue en calculant successivement ses composantes $x_i$ par les formules\n",
    "$$\n",
    "\\begin{align}\n",
    "x_1&=\\frac{b_1}{A_{11}}\\\\\n",
    "x_2&=\\frac{b_2-A_{21}\\,x_1}{A_{22}}\\\\\n",
    "&\\vdots\\\\\n",
    "x_n&=\\frac{b_n-(A_{n1}\\,x_1+\\cdots +A_{n\\,n-1}\\,x_{n-1})}{A_{nn}}\\\\\n",
    "&\\\\\n",
    "&\\\\\n",
    "\\bigg(\\,\\,x_i&=\\frac{b_i-(A_{i1}\\,x_1+\\cdots +A_{i\\,i-1}\\,x_{i-1})}{A_{ii}}\\,\\, \\bigg)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "**Q1.** Définir une fonction ``descente`` qui prend en argument une matrice $A$ triangulaire inférieure inversible et un vecteur $b$ et qui retourne la solution $x$ du système $Ax=b$. Tester votre fonction sur une matrice $A$ à coefficients aléatoires et un second membre $b$ tel que la solution $x$ de $Ax=b$ soit connue.\n",
    "\n",
    "**Q2.** Écrire la solution $x$ du système $Ax=b$ lorsque $A$ est cette fois-ci triangulaire supérieure, en fonction des coefficients de $A$ et de $b$, en résolvant successivement les équations depuis la dernière jusqu'à la première (on dit qu'on résout le système $Ax=b$ par *remontée*).\n",
    "\n",
    "**Q3.** Modifier votre fonction ``descente`` en une fonction que vous appelerez ``remonte_descente`` qui permet la résolution du système $Ax=b$ lorsque $A$ est triangulaire inférieure ou triangulaire supérieure. Votre fonction devra tester si la matrice $A$ est triangulaire supérieure ou inférieure.\n",
    "\n",
    "*Commandes python : essayez les commandes ``np.tril(A), np.triu(A), np.tril(A,k), np.triu(A,k)``, avec $k=1$ ou $k=-1$, et ``np.random.rand(n,n)``, avec $n\\in\\N$. La somme $(A_{i1}\\,x_1+\\cdots +A_{i\\,i-1}\\,x_{i-1})$ peut être vue comme un produit scalaire entre deux vecteurs, utiliser ``np.vdot`` pour le produit scalaire*. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def remontee_descente(A, b):\n",
    "    x = 0 * b\n",
    "    n = len(b)\n",
    "    if np.allclose(A, np.triu(A)):\n",
    "        for i in range(n - 1, -1, -1):\n",
    "            x[i] = (b[i] - np.dot(A[i, i + 1 :], x[i + 1 :])) / A[i, i]\n",
    "    elif np.allclose(A, np.tril(A)):\n",
    "        for i in range(n):\n",
    "            x[i] = (b[i] - np.dot(A[i, :i], x[:i])) / A[i, i]\n",
    "    else:\n",
    "        msg = \"A est ni triangulaire supérieure ni triangulaire inférieure\"\n",
    "        raise ValueError(msg)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1093356479670479e-31\n"
     ]
    }
   ],
   "source": [
    "n = 5\n",
    "A = np.random.rand(n, n)\n",
    "A = np.tril(A) + np.eye(n) * np.linalg.norm(A)\n",
    "xe = np.array([1] * n)\n",
    "b = np.dot(A, xe)\n",
    "x = remontee_descente(A, b)\n",
    "\n",
    "print(np.dot(x - xe, x - xe))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 2.  La méthode de Gauss.\n",
    "\n",
    "\n",
    "On suppose que $A$ est une matrice carrée inversible et qu'il est\n",
    "possible d'appliquer la méthode de Gauss à la matrice $A$ et donc la transformer en une matrice triangulaire supérieure $U$\n",
    "à coefficients diagonaux non nuls simplement en effectuant\n",
    "successivement des opérations élémentaires sur les lignes du type\n",
    "$L_i$ devient $L_i + \\beta L_j$. On suppose donc que les pivots de\n",
    "la méthode de Gauss sont tous non nuls. \n",
    "\n",
    "**Q1.** Vérifier que l'algorithme suivant permet de transformer une matrice donnée $A$ en une matrice triangulaire supérieure $U$ par la méthode de Gauss :\n",
    "\n",
    "\n",
    "```\n",
    "U = A   # on prend une copie qu'on écrasera\n",
    "pour j = 0 à n-1\n",
    "    pour i = j + 1 à n - 1\n",
    "         beta = U(i,j)/U(j,j)                # U(j,j) est le pivot\n",
    "         pour k = j  à  n -1\n",
    "            U(i, k) = U(i,k) - beta * U(j,k) # ligne i devient ligne i - beta * ligne j\n",
    "         fin k\n",
    "    fin i\n",
    "fin j\n",
    "retourner U\n",
    "```\n",
    "\n",
    "**Q2.** Ecrire une fonction Python de la forme ```met_gauss(A)``` correspondant à cet algorithme.\n",
    "\n",
    "*Remarque : vous pouvez écrire les commandes $U(i, k) = U(i,k) - \\beta U(j,k)$, pour $k=j,\\dots,n-1$, sans utiliser de boucle sur $k$, en écrivant le vecteur $(U(i,j),\\dots,U(i,n-1))$ comme ``U[i,j:]``.*\n",
    "\n",
    "**Q3.** Appliquer cette fonction aux matrices\n",
    "$$\n",
    "A=\\left (\n",
    "\\begin{array}{ccc}\n",
    "9 & 8 & 6 \\\\\n",
    "7 & 6 & 12 \\\\\n",
    "9 & 3 & 9\n",
    "\\end{array}\n",
    "\\right )\n",
    "\\qquad \\mbox{ et } \\qquad \n",
    "B=\\left (\n",
    "\\begin{array}{cccc}\n",
    "11 & 8 & 3 & 13 \\\\\n",
    " 2 & 12 & 7 & 10 \\\\\n",
    " 3 & 3 & 17 & 13 \\\\\n",
    " 11 & 2 & 12 & 7\n",
    "\\end{array}\n",
    "\\right )\n",
    "$$\n",
    "\n",
    "Les réponses attendues sont respectivement\n",
    "\\begin{equation*}\n",
    "A=\\left (\n",
    "\\begin{array}{ccc}\n",
    "9 & 8 & 6 \\\\\n",
    "0 & -0.2222222 & 7.3333333 \\\\\n",
    "0 & 0 & -162\n",
    "\\end{array}\n",
    "\\right )\\qquad\n",
    " \\mbox{ et }\\qquad \n",
    "B=\\left (\n",
    "\\begin{array}{cccc}\n",
    "11 & 8 & 3 & 13 \\\\\n",
    " 0 &   10.5454545 &  6.45454545 &  7.63636364\\\\\n",
    " 0 & 0 & 15.6810345 &  8.86206897 \\\\\n",
    " 0 & 0 & 0 &  -8.81693238\n",
    "\\end{array}\n",
    "\\right )\n",
    "\\end{equation*}\n",
    "\n",
    "**Q4.** Adapter votre fonction ```met_gauss(A)``` en une fonction ```met_gauss_sys(A,b)``` de façon à que l'on puisse l'utiliser pour résoudre un système $Ax=b$, avec $b\\in\\RR^n$ donné. Pour cela, il faut le long de la méthode de Gauss faire les mêmes opérations sur la matrice $A$ et sur le second membre $b$. Cette fonction retournera la solution $x$ du système $Ax=b$ en écrivant le système triangulaire équivalent obtenu par la méthode de Gauss, et en résolvant ce système triangulaire avec la fonction ```remonte_descente```. La tester avec une matrice $A$ et un vecteur $b$ aléatoires par exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def met_gauss(A):\n",
    "    U = A\n",
    "    n = len(A)\n",
    "    for j in range(n):\n",
    "        for i in range(j + 1, n):\n",
    "            beta = U[i, j] / U[j, j]\n",
    "            U[i, j:] = U[i, j:] - beta * U[j, j:]\n",
    "    return U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def met_gauss_sys(A, b):\n",
    "    n, m = A.shape\n",
    "    if n != m:\n",
    "        msg = \"Erreur de dimension : A doit etre carré\"\n",
    "        raise ValueError(msg)\n",
    "    if n != b.size:\n",
    "        msg = \"Erreur de dimension : le nombre de lignes de A doit être égal au nombr ede colonnes de b\"\n",
    "        raise valueError(\n",
    "            msg,\n",
    "        )\n",
    "    U = np.zeros((n, n + 1))\n",
    "    U = A\n",
    "    V = b\n",
    "    for j in range(n):\n",
    "        for i in range(j + 1, n):\n",
    "            beta = U[i, j] / U[j, j]\n",
    "            U[i, j:] = U[i, j:] - beta * U[j, j:]\n",
    "            V[i] = V[i] - beta * V[j]\n",
    "    return remontee_descente(U, V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3096323621833204e-32\n"
     ]
    }
   ],
   "source": [
    "n = 5\n",
    "A = np.random.rand(n, n) + float(n) * np.eye(n)\n",
    "b = np.random.rand(n)\n",
    "x = met_gauss_sys(A, b)\n",
    "\n",
    "print(np.dot(b - A.dot(x), b - A.dot(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## La méthode de Gauss avec stratégie de pivot partiel.\n",
    "\n",
    "En pratique, pour des questions de stabilité numérique, on a intérêt à choisir à l'étape $j$ un pivot $A_{k,j}$, avec $k\\geq j$ tel que $|A_{k,j}|$ est maximal (car cela signifie diviser par la quantité la plus grande possible). À l'étape $j$ de la méthode de Gauss, on commence alors par choisir $p$ tel que $|A_{p,j}|=\\max_{k\\geq j}|A_{k,j}|$ et on échange les lignes $p$ et $j$ de $A$.\n",
    "\n",
    "**Q5.** Créer une fonction ```met_gauss_pivot(A,b)``` qui permet la résolution du système $Ax=b$ en utilisant cette stratégie de choix de pivot. La tester avec le même exemple que dans la question précédente.\n",
    "\n",
    "**Q6. (Comparaison des deux méthodes).** Pour $n=10,20,30,\\dots,200$ :\n",
    "* construire une matrice $A\\in\\mathcal{M}_n(\\RR)$ aléatoire, un vecteur $x_{ex}\\in\\RR^n$ aléatoire et calculer $b=Ax_{ex}$ ;\n",
    "* Résoudre le système $Ax=b$ (dont la solution est $x=x_{ex}$) par la méthode de Gauss avec et sans choix de pivot ;\n",
    "* Calculer la norme $\\|x-x_{ex}\\|$ pour chacune des méthodes.\n",
    "\n",
    "Comparer les résultats obtenus pour les deux méthodes. Vous pouvez représenter $\\|x-x_{ex}\\|$, ou, ce qu'est mieux, $\\mathrm{log}(\\|x-x_{ex}\\|)$ en fonction de la taille $n$ de la matrice.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}