init: initial commit

1 files changed, 799 insertions, 0 deletions

diff --git a/PVCM/cama/en/ma31 Système d'équations.ipynb b/PVCM/cama/en/ma31 Système d'équations.ipynb
new file mode 100644
index 0000000..2d5453a
--- /dev/null
+++ b/PVCM/cama/en/ma31 Système d'équations.ipynb
@@ -0,0 +1,799 @@
+{
+ "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.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2,
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import numpy.linalg as lin\n",
+    "\n",
+    "np.set_printoptions(precision=3, linewidth=150, suppress=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Systèmes matriciels\n",
+    "\n",
+    "Un système de plusieurs équations à autant d'inconnues peut s'écrire sous forme d'un système matriciel de la forme $A {\\bf x} = {\\bf b}$ c.a.d.\n",
+    "\n",
+    "$$\n",
+    "\\begin{bmatrix}\n",
+    "a_{11} a_{12} \\ldots a_{1n} \\\\\n",
+    "a_{21} a_{22} \\ldots a_{2n} \\\\\n",
+    " \\vdots \\\\\n",
+    "a_{n1} a_{n2} \\ldots a_{nn} \\\\\n",
+    "\\end{bmatrix}\n",
+    "\\;\n",
+    "\\begin{bmatrix}\n",
+    "x_{1} \\\\\n",
+    "x_{2} \\\\\n",
+    "\\vdots \\\\\n",
+    "x_{n} \\\\\n",
+    "\\end{bmatrix}\n",
+    "=\n",
+    "\\begin{bmatrix}\n",
+    "b_{1} \\\\\n",
+    "b_{2} \\\\\n",
+    "\\vdots \\\\\n",
+    "b_{n} \\\\\n",
+    "\\end{bmatrix}\n",
+    "$$\n",
+    "\n",
+    "Ainsi si votre épicier ne vous donne pas le prix des pommes, poires et bananes mais que vous en achetez 3 fois des quantités différentes, vous pouvez retrouver le prix de chaque fruit. Il suffit de résoudre les 3 équations\n",
+    "qui correspondent aux 3 achats ou le système matriciel 3x3 correpondant. Pour résourdre ce dernier la facon intuitive est d'inverser A  et de calculer $A^{-1}\\, {\\bf b}$ :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([0.8, 0.9, 0.6])"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# A est la quantité de chaque fruit achetée\n",
+    "# x est le prix de chaque fruit\n",
+    "# b est la somme qu'on a payé pour chaque course\n",
+    "\n",
+    "A = np.array([[6,5,4], [5,3,2], [7,3,2]])\n",
+    "b = np.array([11.7, 7.9, 9.5])\n",
+    "\n",
+    "x = lin.inv(A) @ b\n",
+    "x"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Les systèmes matriciels sont probablement l'application la plus importante des matrices car quoi de plus important que de connaitre le prix des bananes ?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Résolution d'un système matriciel\n",
+    "\n",
+    "On regardera les méthodes suivantes :\n",
+    "\n",
+    "* Méthode du pivot de Gauss (rendre la matrice triangulaire)\n",
+    "* Décomposition LU (effet de bord de la résolution par le pivot de Gauss)\n",
+    "* Méthode de Gauss-Jordan (rendre la matrice diagonale)\n",
+    "\n",
+    "**Note :** Si vous trouvez la suite difficile, faites une pause et regardez cette [explication sur une matrice 3x3](https://www.youtube.com/watch?v=VoMsTiJOUE4). Elle commence par triangulariser la matrice comme le fait\n",
+    "la méthode du pivot de Gauss, puis elle continue pour diagonaliser la matrice comme le fait Gauss-Jordan et\n",
+    "finalement elle arrive à l'inverse de la matrice. C'est très bien pour comprendre ou faire les calculs à la main\n",
+    "sur une petite matrice, c'est horrible à programmer. Ce que je vous montre c'est la transformation de l'intuition\n",
+    "qu'offre cette vidéo en un calcul matriciel rapide et simple à programmer."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Méthode du pivot de Gauss\n",
+    "\n",
+    "On transforme A en une matrice triangulaire qui permet ensuite de résoudre le système en O(n²) opérations.\n",
+    "\n",
+    "Pour cela on commence à mettre des 0 sur la première colonne en dessous de la diagonale. Pour cela il\n",
+    "suffit de multiplier A par la matrice $E_1$ suivante :\n",
+    "\n",
+    "\n",
+    "$$\n",
+    "E_1 = \n",
+    "\\begin{bmatrix}\n",
+    "\\;1 \\quad 0\\; 0 \\ldots 0 \\\\\n",
+    "\\frac{-a_{21}}{a_{11}} \\, 1\\; 0  \\ldots 0 \\\\\n",
+    "\\frac{-a_{31}}{a_{11}} \\, 0\\; 1  \\ldots 0 \\\\\n",
+    "\\vdots \\\\\n",
+    "\\frac{-a_{n1}}{a_{11}}\\; 0\\; 0  \\ldots 1 \\\\\n",
+    "\\end{bmatrix}\n",
+    "$$\n",
+    "\n",
+    "$E_2$ sera la matrice équivalente avec des termes $\\frac{-a_{k2}}{a_{22}}$ sous la diagonales afin des 0 dans la 2e colonne de A sous la diagonale, etc.\n",
+    "\n",
+    "Si vous ne connaissez pas cette méthode, je vous conseille de vérifier à la main sur une matrice 4x4 que\n",
+    "$E_1 \\, A$ met bien des zéros là où il faut. Ensuite faites le travail sur la 2e colonne.\n",
+    "\n",
+    "Bien sur si on multiplie A par $E_1$, il faut aussi multiplier ${\\bf b}$ par $E_1$ pour que le système matriciel soit\n",
+    "équivalent. Cette multiplication $E_1 \\, \\bf b$ peut se calculer nettement plus rapidement qu'avec un \n",
+    "produit matrice vecteur. Donnez la bonne formule et vérifiez dans le code ci-dessous.\n",
+    "\n",
+    "#### Système matriciel avec matrice triangulaire\n",
+    "\n",
+    "Pour finir il faut résoudre U x = c avec U une matrice triangulaire supérieur. Cela se fait en partant\n",
+    "de la dernière ligne qui donne directement `x[-1] = c[-1] / U[-1,-1]`. Une fois `x[-1]` connu, on en\n",
+    "déduit la valeur de `x[-2]` avec l'avant dernière ligne, etc.\n",
+    "\n",
+    "On faisant les calculs on constate qu'il faut environ n²/2 additions et multiplications.\n",
+    "\n",
+    "Là encore, si vous ne vous rappelez pas comment on fait, faites le à la main avec une matrice 3x3. Vous pouvez aussi ajouter des `print` dans le code qui suit pour voir les matrices intermédiaires."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Cf analyse du code ci-dessous\n",
+    "\n",
+    "def solve_gauss(A, b):\n",
+    "    A = A.astype(np.float64)   # si A a des entiers, on va avoir des erreurs de calculs\n",
+    "    for i in range(len(A)-1):\n",
+    "        coefs = - A[i+1:,i] / A[i,i]  # Normalement il faut tester que A[i,i] != 0\n",
+    "        A[i+1:, i:] += np.outer(coefs, A[i, i:]) # ou np.einsum('i,j -> ij', coefs, A[i, i:])\n",
+    "        b[i+1:] += coefs * b[i]\n",
+    "    # A est maintenant triangulaire surpérieur\n",
+    "    res = np.zeros(len(b), dtype=b.dtype)\n",
+    "    res[-1] = b[-1] / A[-1,-1]\n",
+    "    for i in range(len(A)-1)[::-1]:\n",
+    "        res[i] = (b[i] - A[i,i+1:] @ res[i+1:]) / A[i,i]\n",
+    "    return res\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Analyse du code\n",
+    "\n",
+    "Comme **E** = **Id** + **Coefs**, on a **E A** = **A** + **Coefs A** et donc **A** += **Coef A**. Ici **Coefs** est la matrice **E** ci-dessus sans les 1 sur la diagonale.\n",
+    "\n",
+    "Le problème est que `Coef @ A` fait n³ opérations ce qui est beaucoup trop et surtout inutile sachant que **Coef** est un vecteur dans une matrice de zéros. Si ce vecteur est dans la colonne $i$ alors seuls les termes de la ligne $i$ de A serviront pour générer ce produit matriciel (à faire sur un papier). De plus on a `res[j,k] = Coef[j,i] * A[i,k]` ce qui s'écrit en Numpy `res = np.outer(coefs, A[i,:])` si `coefs` est la colonne $i$ de `Coefs`.\n",
+    "\n",
+    "Comme la matrice **Coefs** n'a des valeurs non nulles qu'en dessous de la diagonale, cela veut dire que le résultat `res` de la multiplication n'a que des zéros au dessus de la ligne $i$. Donc autant ne faire que le calcul pour ce qui est en\n",
+    "dessus de la ligne $i$ : `A[i+1:, i:] += np.outer(coefs, A[i, i:])` avec `coefs` un vecteur de taille $n-i$. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A\n",
+      " [[2.878 0.086 7.04  4.388 7.731]\n",
+      " [1.997 0.405 9.546 0.681 3.954]\n",
+      " [9.664 5.849 3.799 5.096 8.049]\n",
+      " [7.136 3.182 6.452 1.741 3.594]\n",
+      " [1.064 3.951 8.733 7.093 9.787]] \n",
+      "b\n",
+      " [22.122 16.584 32.457 22.105 30.627]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array([1., 1., 1., 1., 1.])"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A = 10 * np.random.random(size=(5,5))\n",
+    "b = A.sum(axis=1)\n",
+    "print(f\"A\\n {A} \\nb\\n {b}\")\n",
+    "solve_gauss(A, b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Complexité du pivot de Gauss\n",
+    "\n",
+    "#### Exercice 10.1\n",
+    "\n",
+    "Calculer la complexité de `solve_gauss`. Vous devez trouver n³/3 multiplications et autant d'additions\n",
+    "plus des termes de degrés inférieurs. Pour compter les opérations je vous suggère de faire un dessin \n",
+    " et d'imaginer votre dessin en 3 dimensions."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Décomposition LU (Lower, Upper)\n",
+    "\n",
+    "Si on a besoin de résoudre plusieurs systèmes matriciels avec A, alors il est préférable de décomposer\n",
+    "A en un produit d'une matriciel triangulaire inférieur et d'une matrice triangulaire supérieure.\n",
+    "\n",
+    "$A = L\\, U$\n",
+    "\n",
+    "Pour cela il suffit de relancer le pivot de Gauss mais au lieu de modifier **b** à chaque itération, on\n",
+    "fabrique à la volée l'inverse de la matrices $E_n \\ldots E_2\\, E_1$ : \n",
+    "\n",
+    "$$\n",
+    "E_n \\ldots E_2\\, E_1 \\, A\\, x = E_n \\ldots E_2\\, E_1 \\, b \\quad \\textrm{donc} \\\\\n",
+    "(E_n \\ldots E_2\\, E_1)^{-1} \\, E_n \\ldots E_2\\, E_1 \\, A\\, x =  b \\\\\n",
+    "E_1^{-1} \\, E_2^{-1} \\ldots E_n^{-1} \\; E_n \\ldots E_2\\, E_1 \\, A\\, x =  b \\\\\n",
+    "$$\n",
+    "\n",
+    "Ce calcul est très simple car ces matrices inverses ont la\n",
+    "même forme que les $E_k$ avec simplement le valeurs opposées sous la diagonale :\n",
+    "\n",
+    "$$\n",
+    "E_1^{-1} = \n",
+    "\\begin{bmatrix}\n",
+    "\\;1 \\quad 0\\; 0 \\ldots 0 \\\\\n",
+    "\\frac{a_{21}}{a_{11}} \\, 1\\; 0  \\ldots 0 \\\\\n",
+    "\\frac{a_{31}}{a_{11}} \\, 0\\; 1  \\ldots 0 \\\\\n",
+    "\\vdots \\\\\n",
+    "\\frac{a_{n1}}{a_{11}}\\; 0\\; 0  \\ldots 1 \\\\\n",
+    "\\end{bmatrix}\n",
+    "$$\n",
+    "\n",
+    "et le produit $E^{-1} = E_1^{-1} \\,E_1^{-1} \\,E_2^{-1} \\,E_3^{-1} \\, \\ldots \\,E_n^{-1}$ est simplement la concaténation\n",
+    "des colonnes (je vous laisse vérifier) :\n",
+    "\n",
+    "$$\n",
+    "E^{-1} = \n",
+    "\\begin{bmatrix}\n",
+    "\\;1 \\quad 0\\; \\; 0 \\ldots 0 \\\\\n",
+    "\\frac{a_{21}}{a_{11}} \\; \\; 1\\; \\; 0  \\ldots 0 \\\\\n",
+    "\\frac{a_{31}}{a_{11}} \\, \\frac{a_{32}}{a_{22}}  \\; 1  \\ldots 0 \\\\\n",
+    "\\vdots \\\\\n",
+    "\\frac{a_{n1}}{a_{11}}\\; \\frac{a_{n2}}{a_{22}}\\; \\frac{a_{n3}}{a_{33}}  \\ldots 1 \\\\\n",
+    "\\end{bmatrix}\n",
+    "$$\n",
+    "\n",
+    "et ainsi on a $L = E^{-1}$ et $U = E_n \\ldots E_2\\, E_1 \\, A$ puisque tout le travail fait sur A est\n",
+    "justement pour mettre des 0 sous la diagonale et donc en faire une matrice triangulaire supérieure.\n",
+    "\n",
+    "Une fois la décomposition LU faites, résoudre $L\\, U \\, {\\bf x} = {\\bf b}$ se fait en 2 étapes :\n",
+    "\n",
+    "* on résoud $L\\, {\\bf y}  = {\\bf b}$ ce qui donne ${\\bf y}$ \n",
+    "* puis on résoud $U\\, {\\bf x} = {\\bf y}$ et on a notre solution ${\\bf x}$\n",
+    "\n",
+    "en O(n²) opérations à chaque fois puisqu'il s'agit de matrices triangulaires."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def LU(A):\n",
+    "    L = np.diag([1.,] * len(A))\n",
+    "    for i in range(len(A)-1):\n",
+    "        E = np.diag([1.,] * len(A))\n",
+    "        E[i+1:,i] = - A[i+1:,i] / A[i,i]\n",
+    "        L[i+1:,i] = -E[i+1:,i]\n",
+    "        A[i:, i:] = E[i:,i:] @ A[i:,i:]\n",
+    "    return L, A"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A\n",
+      " [[7.512 7.513 5.909 4.879 0.641]\n",
+      " [8.468 6.44  4.16  1.435 4.431]\n",
+      " [6.661 4.823 2.936 3.073 4.768]\n",
+      " [5.927 6.654 1.168 9.273 0.706]\n",
+      " [0.631 1.075 2.517 6.229 7.547]]\n",
+      "L\n",
+      "[[  1.      0.      0.      0.      0.   ]\n",
+      " [  1.127   1.      0.      0.      0.   ]\n",
+      " [  0.887   0.906   1.      0.      0.   ]\n",
+      " [  0.789  -0.358 122.023   1.      0.   ]\n",
+      " [  0.084  -0.219 -40.943  -0.357   1.   ]]\n",
+      "U\n",
+      "[[   7.512    7.513    5.909    4.879    0.641]\n",
+      " [   0.      -2.029   -2.502   -4.066    3.708]\n",
+      " [   0.       0.      -0.036    2.432    0.838]\n",
+      " [   0.       0.       0.    -292.814 -100.761]\n",
+      " [   0.       0.       0.       0.       6.662]]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 0.,  0.,  0.,  0.,  0.],\n",
+       "       [ 0.,  0.,  0., -0.,  0.],\n",
+       "       [ 0.,  0.,  0.,  0.,  0.],\n",
+       "       [ 0.,  0.,  0., -0., -0.],\n",
+       "       [ 0.,  0.,  0., -0., -0.]])"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A = 10 * np.random.random(size=(5,5))\n",
+    "print(\"A\\n\",A)\n",
+    "L,U = LU(A.copy())  # Attention, notre fonction modifie A donc si on veut le réutiliser il faut une copie\n",
+    "print(f\"L\\n{L}\\nU\\n{U}\")\n",
+    "A - (L @ U)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Gauss Jordan\n",
+    "\n",
+    "L'idée est de diagonaliser A par des multiplications matricielles similaires à celles de Gauss mais qui\n",
+    "annulent aussi les termes au dessus de la diagonale :\n",
+    "\n",
+    "$$\n",
+    "E_3 = \n",
+    "\\begin{bmatrix}\n",
+    "1 \\; 0\\; \\frac{-a_{13}}{a_{33}} \\;  0 \\ldots 0 \\\\\n",
+    "0 \\; 1\\; \\frac{-a_{23}}{a_{33}} \\;  0 \\ldots 0 \\\\\n",
+    "0 \\; 0\\quad 1 \\quad  0 \\ldots 0 \\\\\n",
+    "0 \\; 0\\; \\frac{-a_{43}}{a_{33}}  1 \\ldots 0 \\\\\n",
+    "\\vdots \\\\\n",
+    "0 \\; 0\\; \\frac{-a_{n3}}{a_{33}}  0 \\ldots 1 \\\\\n",
+    "\\end{bmatrix}\n",
+    "$$\n",
+    "\n",
+    "\n",
+    "Cette méthode est plus compacte à écrire que Gauss car il n'y a pas remonté finale, mais elle est plus lente car on modifie tout A et b à\n",
+    "chaque itération."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def solve_gauss_jordan(A,b):\n",
+    "    for i in range(len(A)):\n",
+    "        d1 = np.diag([1.,] * len(A))\n",
+    "        d1[:,i] = - A[:,i] / A[i,i]\n",
+    "        A = d1 @ A\n",
+    "        b = d1 @ b\n",
+    "    return b / np.diag(A)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([1., 1., 1., 1., 1.])"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A = 10 * np.random.random(size=(5,5))\n",
+    "b = A.sum(axis=1)\n",
+    "solve_gauss_jordan(A, b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Comparaison de la vitesse de méthodes"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n = 500\n",
+    "A = 10 * np.random.random(size=(n,n))\n",
+    "b = A.sum(axis=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "2.63 s ± 218 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit solve_gauss_jordan(A, b)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "101 ms ± 3.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit solve_gauss(A, b)   # 2.5 fois plus rapide que Gauss-Jordan"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "16.4 ms ± 227 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit B = lin.inv(A)  # 6 fois plus rapide que mon pivot de Gauss"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "14.2 ms ± 277 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit lin.solve(A, b) # un peu plus rapide"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Pour n=5000, `lin.solve` est 3 fois plus rapide que `lin.inv` sur ma machine. Et plus on augmente $n$ et plus\n",
+    "l'écart se creuse \n",
+    "\n",
+    "**Note**: Les calculs utilisent tous les processeurs de ma machine, Numpy sait paralléliser et utilise\n",
+    "les bonnes bibliothèques (BLAS et Lapack) s'il est bien installé.\n",
+    "\n",
+    "#### Morale\n",
+    "\n",
+    "* Une bonne méthode est mieux qu'une mauvaise (et oui !).\n",
+    "* Une bonne bibliothèque fait toute la différence (vous ne ferez pas mieux que BLAS et Lapack)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Erreurs d'arrondi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "finfo(resolution=1e-06, min=-3.4028235e+38, max=3.4028235e+38, dtype=float32)"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.set_printoptions(precision=20, linewidth=150, suppress=True)\n",
+    "\n",
+    "np.finfo(np.float32)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A\n",
+      " [[0.000001 1.      ]\n",
+      " [1.       2.      ]] \n",
+      "b\n",
+      " [1. 3.]\n"
+     ]
+    }
+   ],
+   "source": [
+    "e = 1E-6\n",
+    "A = np.array([[e, 1], [1, 2]], dtype='float32')\n",
+    "b = np.array([1., 3.], dtype='float32')\n",
+    "print(f\"A\\n {A} \\nb\\n {b}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "x = [1.013279 0.999999]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array([1.      , 3.013277], dtype=float32)"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "x = solve_gauss(A.copy(),b)\n",
+    "print(\"x =\", x)\n",
+    "A @ x                # on ne retrouve pas exactement b"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Il y a un problème d'arrondi évident sachant que la solution est [1 + 2*e, 1 - e] :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "La solution est : [1.000002 0.999999]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array([1., 3.], dtype=float32)"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "x = np.array([1 + 2*e, 1 - e], dtype='float32')\n",
+    "print(\"La solution est :\", x)\n",
+    "A @ x             # on retrouve b"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Le problème est qu'on divise les valeurs de la première colonne de A par le premier pivot : 1E-6.\n",
+    "\n",
+    "Aussi on multiplie par 1E6 ce qui va faire des grandes valeurs or la précision est au niveau du 6e chiffre,\n",
+    "on pert les suivants comme le montre la simple addition qui suit :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([1000000.], dtype=float32)"
+      ]
+     },
+     "execution_count": 18,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "x = np.array([1E6], dtype='float32')   # 1 million\n",
+    "y = np.array([1E-2], dtype='float32')  # 0.01\n",
+    "x+y                                    # résultat : 1 million sans rien après la virgule"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**À mélanger des grands nombres avec des petits, on fait des erreurs d'arrondi.**\n",
+    "\n",
+    "C'est vrai pour les calculs matriciels qu'on fait, c'est toujours vrai et c'est pour cela qu'on aime normaliser\n",
+    "les données sur lesquelles on travaille lorsqu'on le peut."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Solution au problème d'arrondi dans le cas du pivot de Gauss\n",
+    "\n",
+    "Il suffit d'échanger les lignes (ou colonnes mais c'est plus compliqué) pour avoir la plus grande valeur,\n",
+    "en valeur absolue, comme pivot. C'est la méthode du pivot partiel.\n",
+    "\n",
+    "La méthode du pivot total consiste à prendre comme pivot la plus grande valeur parmis la ligne et la colonne concernées.\n",
+    "\n",
+    "#### Exercice pour le TP\n",
+    "\n",
+    "Programmer la méthode du pivot partiel et vérifier ce que cela donne sur l'exemple ci-dessus qui pose\n",
+    "des problèmes d'arrondi."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "A = np.array([1,2,3])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.isin?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array(False)"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.isin(5,A)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "B = np.array([1,5,6])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([ True, False, False])"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.isin(B,A)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ]
+}
\ No newline at end of file