init: initial commit

1 files changed, 417 insertions, 0 deletions

diff --git a/PVCM/cama/fr/ma40 Méthodes itératives.ipynb b/PVCM/cama/fr/ma40 Méthodes itératives.ipynb
new file mode 100644
index 0000000..ae66062
--- /dev/null
+++ b/PVCM/cama/fr/ma40 Méthodes itératives.ipynb
@@ -0,0 +1,417 @@
+{
+ "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.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4,
+ "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)\n",
+    "np.random.seed(123)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "## La simulation numérique\n",
+    "\n",
+    "Si le prix des bananes est important, le secteur qui a les plus gros besoins en résolution de systèmes\n",
+    "matriciels est la simulation numérique.\n",
+    "Cela comprend <a href=\"https://www.nas.nasa.gov/SC18/demos/demo10.html\">ceci</a>\n",
+    "\n",
+    "<center><img alt=\"simulation du X-57 par la Nasa\" src=\"images/nasa-x57-pression.jpg\"/></center>\n",
+    "\n",
+    "et tout ce qu'on simule numériquement à savoir tout ce qui touche le transport, l'énergie, la construction, tout\n",
+    "ce qu'on fabrique, qui coûte cher et qui doit avoir un comportement physique bien précis. Mais ce n'est pas tout,\n",
+    "comprendre notre environnement (météo, réchauffement de la planète, chimie...) revient aussi à résoudre des\n",
+    "systèmes matriciels ainsi que d'autres domaines. Cela étant il existe des méthodes de simulation numérique\n",
+    "qui ne passent pas par des systèmes matriciels.\n",
+    "\n",
+    "Pour faire l'image ci-dessus, on transforme des équations physique comme [Navier-Stokes](https://fr.wikipedia.org/wiki/%C3%89quations_de_Navier-Stokes) en un système matriciel où les inconnues sont définies en chaque point\n",
+    "d'un maillage à définir. Dans notre cas l'inconnue est la pression et le maillage est l'intérieur d'une\n",
+    "boite imaginaire qui comprend l'avion et l'air qui circule autour.\n",
+    "\n",
+    "Si la boite est un cube avec 1000 points dans chaque direction, on a 1 milliard de points dans la boite (moins\n",
+    "ce qui est à l'intérieur de l'avion mais restons sur 1 milliards). Alors la matrice a 1 trillion d'éléments (un milliard au carré).\n",
+    "\n",
+    "$$\n",
+    "\\begin{bmatrix}\n",
+    "a_{11} \\; a_{12} \\ldots a_{1,10^9} \\\\\n",
+    "a_{21} \\; a_{22} \\ldots a_{2,10^9} \\\\\n",
+    " \\vdots \\\\\n",
+    "a_{10^9,1} a_{n2} \\ldots a_{10^9,10^9} \\\\\n",
+    "\\end{bmatrix}\n",
+    "\\;\n",
+    "\\begin{bmatrix}\n",
+    "p_{1} \\\\\n",
+    "p_{2} \\\\\n",
+    "\\vdots \\\\\n",
+    "p_{10^9} \\\\\n",
+    "\\end{bmatrix}\n",
+    "=\n",
+    "\\begin{bmatrix}\n",
+    "f_{1} \\\\\n",
+    "f_{2} \\\\\n",
+    "\\vdots \\\\\n",
+    "f_{10^9} \\\\\n",
+    "\\end{bmatrix}\n",
+    "$$\n",
+    "\n",
+    "Inverser une matrice se fait en $O(n³)$ opérations soit $O(10^{27})$ dans notre cas. \n",
+    "\n",
+    "Si vous avez [l'ordinateur le plus puissant du monde](https://www.top500.org/)\n",
+    "vous pouvez faire 1 Eflops soit $10^{18}$ opérations par seconde en pic (en 2024). Aussi il vous faudra $10^{9}$ secondes soit un peu plus de 30 ans. C'est trop.\n",
+    "\n",
+    "**Inverser une matrice ou résoudre par une méthode directe n'est pas la bonne solution pour résoudre un grand système matriciel.**\n",
+    "\n",
+    "#### Exercice 1\n",
+    "Pour une telle simulation il faut aussi calculer la vitesse de l'air en chacun des 1 milliard de points de notre maillage. Une vitesse c'est 3 variables qu'il faut ajouter à la pression. Combien de temps faut-il pour inverser\n",
+    "la matrice qui intègre aussi la vitesse de l'air ?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ceci est une calculatrice, vous pouvez donc écrire la réponse ici\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "Estimer le temps d'un gros calcul avant de le lancer est une bonne idée. \n",
+    "\n",
+    "D'ailleur avec des temps si longs, il ne\n",
+    "faut pas rester aux ordres de grandeurs mais préciser avec la formule exacte. Ainsi avec Choleski c'est $n^3/3$\n",
+    "opérations, donc 3 fois moins, mais bon..."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "## Méthodes itératives\n",
+    "\n",
+    "Les méthodes itératives sont des méthodes qui s'approchent pas à pas de la solution recherchée. Elles permettent\n",
+    "de trouver une approximation de ${\\bf x}$, la solution de $A\\, {\\bf x} = b$.\n",
+    "\n",
+    "En général on\n",
+    "arrête le calcul lorsqu'on estime qu'on est à une distance choisie de la solution (qu'on appelle l'erreur) plutôt\n",
+    "que d'attendre d'être à la solution exacte.\n",
+    "De toute facon la solution exacte sur un ordinateur qui est limité en nombre de chiffres après la virgule peut\n",
+    "être inatteignable. Donc **on cherchera jamais à avoir une erreur plus petite que notre précision maximale**.\n",
+    "\n",
+    "L'idée fondamentale de l'algorithme itératif est d'avoir une formule du type $\\; {\\bf x}^{t+1} = B \\, {\\bf x}^t + {\\bf c}\\;$ ou en Python:\n",
+    "\n",
+    "```\n",
+    "x = np.random(size = c.size)\n",
+    "while np.square(x - old_x) > seuil:\n",
+    "    old_x = x\n",
+    "    x = B @ x + c\n",
+    "```\n",
+    "\n",
+    "La magie c'est si **x** converge. Dans ce cas on a atteint un point fixe c.a.d. que ${\\bf x}^{t+1} = {\\bf x}^t$\n",
+    "et donc \n",
+    "\n",
+    "$${\\bf x}^t = B \\, {\\bf x}^t + {\\bf c} \\quad \\textrm{c.a.d.} \\quad (Id -B) \\, {\\bf x}^t = {\\bf c}$$\n",
+    "\n",
+    "On retrouve notre $A \\; {\\bf x} = {\\bf b}$ cela étant en pratique on ne découpe pas A en Id et B car\n",
+    "ca ne converge pas (sauf cas particuliers). Regardons comment fait la méthode de Jacobi."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "## Méthode de Jacobi\n",
+    "\n",
+    "La méthode de Jacobi découpe la matrice A en M et N avec\n",
+    "\n",
+    "* M la matrice diagonale qui comprend les éléments de la diagonale de A\n",
+    "* N = M - A  (donc 0 sur la diagonale et l'opposé des éléments de A ailleurs)\n",
+    "\n",
+    "Ainsi le système à résoudre est $\\; (M - N) \\, {\\bf x} = {\\bf b}$.\n",
+    "\n",
+    "La formule itérative est donc pour l'itération k+1 exprimée en fonction de l'itération k :\n",
+    "\n",
+    "$$\n",
+    "M \\; {\\bf x}^{k+1} =  N\\; {\\bf x}^k + {\\bf b}\n",
+    "$$\n",
+    "\n",
+    "et comme M est une matrice diagonale on a :\n",
+    "\n",
+    "$$\n",
+    "\\begin{array}{l}\n",
+    "a_{11} x_1^{k+1} = \\qquad -a_{12} \\, x_2^k - a_{13} \\, x_3^k  \\ldots - a_{1n} \\, x_n^k + b_1 \\\\\n",
+    "a_{22} x_2^{k+1} = -a_{21} \\, x_1^k \\qquad - a_{23} \\, x_3^k  \\ldots - a_{2n} \\, x_n^k + b_2 \\\\\n",
+    "a_{33} x_3^{k+1} = -a_{31} \\, x_1^k - a_{32} \\, x_2^k  \\qquad \\ldots - a_{3n} \\, x_n^k + b_3 \\\\\n",
+    " \\vdots \\\\\n",
+    "a_{nn} x_n^{k+1} = -a_{n1} \\, x_1^k - a_{n2} \\, x_2^k  \\ldots - a_{n-1,n-1} \\, x_{n-1}^k \\qquad + b_n \\\\\n",
+    "\\end{array}\n",
+    "$$\n",
+    "\n",
+    "On voit en filigramme $A \\; {\\bf x} = {\\bf b}$.\n",
+    "\n",
+    "Pour calculer $x_1^{k+1}$ il faut diviser le terme de droite de la première ligne par $a_{11}$ il faut donc que  $a_{11} \\ne 0$.\n",
+    "\n",
+    "Pour calculer les termes suivants $x_i^{k+1}$ il faut aussi que $a_{ii}$ soient non nul donc **il\n",
+    "faut que A n'ait pas pas de zéro sur sa diagonale**.\n",
+    "\n",
+    "Cela se retrouve dans l'écriture suivante qui reprend la formule initiale :\n",
+    "$ {\\bf x}^{k+1} =  M^{-1} \\;(N\\; {\\bf x}^k + {\\bf b})$ \n",
+    "\n",
+    ">   On voit que pour être efficace, il faut que M soit facilement inversible, sinon autant inverser A directement.\n",
+    "Ici c'est bien le cas puisque M est diagonale.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "#### Programmons Jacobi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A:\n",
+      " [[2 2 6 1]\n",
+      " [3 9 6 1]\n",
+      " [0 1 9 0]\n",
+      " [0 9 3 4]] \n",
+      "b:\n",
+      " [11 19 10 16] \n",
+      "\n",
+      "M:\n",
+      " [[2 0 0 0]\n",
+      " [0 9 0 0]\n",
+      " [0 0 9 0]\n",
+      " [0 0 0 4]]\n",
+      "N:\n",
+      " [[ 0 -2 -6 -1]\n",
+      " [-3  0 -6 -1]\n",
+      " [ 0 -1  0  0]\n",
+      " [ 0 -9 -3  0]]\n",
+      "\n",
+      "x_0 = [0.398 0.738 0.182 0.175]\n",
+      "x_1 = [4.127 1.837 1.029 2.203]\n",
+      "x_2 = [-0.526 -0.195  0.907 -0.906]\n",
+      "x_3 = [3.427 1.782 1.133 3.759]\n",
+      "x_4 = [-1.56  -0.204  0.913 -0.86 ]\n",
+      "x_5 = [3.395 2.118 1.134 3.775]\n",
+      "x_6 = [-1.907 -0.196  0.876 -1.616]\n",
+      "x_7 = [3.877 2.342 1.133 3.784]\n",
+      "x_8 = [-2.133 -0.357  0.851 -2.12 ]\n",
+      "x_9 = [4.364 2.49  1.151 4.165]\n",
+      "x_10 = [-2.525 -0.574  0.834 -2.467]\n",
+      "x_11 = [4.804 2.671 1.175 4.665]\n",
+      "x_12 = [-3.027 -0.792  0.814 -2.89 ]\n",
+      "x_13 = [5.293 2.898 1.199 5.17 ]\n",
+      "x_14 = [-3.581 -1.027  0.789 -3.421]\n",
+      "x_15 = [5.87  3.159 1.225 5.719]\n",
+      "x_16 = [-4.194 -1.298  0.76  -4.026]\n",
+      "x_17 = [6.531 3.45  1.255 6.35 ]\n",
+      "x_18 = [-4.891 -1.608  0.728 -4.704]\n",
+      "x_19 = [7.277 3.779 1.29  7.073]\n"
+     ]
+    }
+   ],
+   "source": [
+    "A = np.random.randint(10, size=(4,4))\n",
+    "b = A.sum(axis=1)                     # ainsi la solution est [1,1,1,1]\n",
+    "print('A:\\n', A, \"\\nb:\\n\", b, \"\\n\")\n",
+    "\n",
+    "M = np.diag(A)        # attention, c'est un vecteur\n",
+    "N = np.diag(M) - A    # np.diag d'une matrice donne un vecteur, np.diag d'un vecteur donne une matrice\n",
+    "print(f\"M:\\n {np.diag(M)}\\nN:\\n {N}\\n\")\n",
+    "\n",
+    "x = np.random.random(4)\n",
+    "for i in range(20):\n",
+    "    print(f\"x_{i} = {x}\")\n",
+    "    x = (N @ x + b) / M"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "Ca ne converge pas vraiment... (si vous relancez et voyez des `NaN` c'est qu'il y a un zéro sur la diagonale)\n",
+    "\n",
+    "2e essai :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A:\n",
+      " [[24  2  4  8]\n",
+      " [ 0 24  9  3]\n",
+      " [ 4  6 16  5]\n",
+      " [ 6  2  1 32]] \n",
+      "b:\n",
+      " [38 36 31 41] \n",
+      "\n",
+      "M:\n",
+      " [[24  0  0  0]\n",
+      " [ 0 24  0  0]\n",
+      " [ 0  0 16  0]\n",
+      " [ 0  0  0 32]]\n",
+      "N:\n",
+      " [[ 0 -2 -4 -8]\n",
+      " [ 0  0 -9 -3]\n",
+      " [-4 -6  0 -5]\n",
+      " [-6 -2 -1  0]]\n",
+      "\n",
+      "x_0 = [0.766 0.777 0.028 0.174]\n",
+      "x_1 = [1.456 1.468 1.4   1.088]\n",
+      "x_2 = [0.865 0.839 0.683 0.873]\n",
+      "x_3 = [1.109 1.135 1.134 1.045]\n",
+      "x_4 = [0.951 0.944 0.908 0.967]\n",
+      "x_5 = [1.031 1.039 1.043 1.015]\n",
+      "x_6 = [0.984 0.982 0.973 0.99 ]\n",
+      "x_7 = [1.009 1.011 1.014 1.005]\n",
+      "x_8 = [0.995 0.994 0.992 0.997]\n",
+      "x_9 = [1.003 1.003 1.004 1.002]\n",
+      "x_10 = [0.998 0.998 0.998 0.999]\n",
+      "x_11 = [1.001 1.001 1.001 1.   ]\n",
+      "x_12 = [1.    0.999 0.999 1.   ]\n",
+      "x_13 = [1. 1. 1. 1.]\n",
+      "x_14 = [1. 1. 1. 1.]\n",
+      "x_15 = [1. 1. 1. 1.]\n",
+      "x_16 = [1. 1. 1. 1.]\n",
+      "x_17 = [1. 1. 1. 1.]\n",
+      "x_18 = [1. 1. 1. 1.]\n",
+      "x_19 = [1. 1. 1. 1.]\n"
+     ]
+    }
+   ],
+   "source": [
+    "A = np.random.randint(10, size=(4,4))\n",
+    "A = A + np.diag(A.sum(axis=0))\n",
+    "b = A.sum(axis=1)                     # ainsi la solution est [1,1,1,1]\n",
+    "print('A:\\n', A, \"\\nb:\\n\", b, \"\\n\")\n",
+    "\n",
+    "\n",
+    "M = np.diag(A)        # attention, c'est un vecteur\n",
+    "N = np.diag(M) - A    # np.diag d'une matrice donne un vecteur, np.diag d'un vecteur donne une matrice\n",
+    "print(f\"M:\\n {np.diag(M)}\\nN:\\n {N}\\n\")\n",
+    "\n",
+    "x = np.random.random(4)\n",
+    "for i in range(20):\n",
+    "    print(f\"x_{i} = {x}\")\n",
+    "    x = (N @ x + b) / M"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "Magie !"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "### Pourquoi le 2e cas marche ?\n",
+    "\n",
+    "Pour qu'une méthode itérative du type  ${\\bf x} = B\\; {\\bf x} + {\\bf c}$  converge il faut au choix\n",
+    "\n",
+    "* $\\rho(B) < 1\\quad$ avec $\\rho$ le rayon spectral (qui est la plus grande valeur propre en valeur absolue)\n",
+    "* $||B|| < 1\\quad$ avec une norme matricielle est subordonnée à une norme vectorielle.\n",
+    "             \n",
+    "\n",
+    "#### Cas de la méthode de Jacobi\n",
+    "\n",
+    "Pour la matrice B de Jacobi, on respecte ces conditions si A est à **diagonale dominante** à savoir que chaque\n",
+    "élément de la diagonale est plus grand en module que la sommes tous les autres éléments en module de sa ligne ou colonne ($|a_{i,i}| \\ge \\sum_{j \\ne i} |a_{i,j}|$).\n",
+    "\n",
+    "Jacobi converge aussi si A est symétrique, réelle et définie positive \n",
+    "(c.a.d. $\\forall {\\bf x}, \\; {\\bf x}^T \\, A \\, {\\bf x} > 0$)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "### Temps calcul\n",
+    "\n",
+    "Si on converge en 10 itérations alors\n",
+    "on utilise 10 multiplications matricielles, 10 additions vectorielles et 10 divisions vectorielles soit 180 opérations\n",
+    "à comparer aux $4^3 / 3 = 22$ opérations d'une méthode directe. Zut !\n",
+    "\n",
+    "Quelques raisons qui font qu'on y perd sont\n",
+    "\n",
+    "* La matrice A est trop petite, les méthodes itératives sont intéressantes pour les grandes matrice\n",
+    "* La méthode de Jacobi n'est pas la meilleure (mais elle est très facilement parallélisable)\n",
+    "* Le rayon spectral de la matrice est trop grand. Plus le rayon spectral est\n",
+    "  petit et plus on converge vite."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ]
+}
\ No newline at end of file