🚀 Initial Commit for my studies work

Analyse Matricielle/TP1_Methode_de_Gauss.ipynb

@@ -0,0 +1,368 @@
+{
+ "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",
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\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,2,0,0,0],[0,0,3,0,0],[0,0,0,4,0],[0,0,0,0,5]])\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",
+    "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",
+    "        raise ValueError(\"A est ni triangulaire supérieure ni triangulaire inférieure\")\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",
+    "        raise ValueError(\"Erreur de dimension : A doit etre carré\")\n",
+    "    if n != b.size:\n",
+    "        raise valueError(\"Erreur de dimension : le nombre de lignes de A doit être égal au nombr ede colonnes de b\")\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
+}