init: initial commit

1 files changed, 336 insertions, 0 deletions

diff --git a/PVCM/cama/fr/ma42 Surrelaxation pour Gauss-Seidel -- Exercice.ipynb b/PVCM/cama/fr/ma42 Surrelaxation pour Gauss-Seidel -- Exercice.ipynb
new file mode 100644
index 0000000..bb2a484
--- /dev/null
+++ b/PVCM/cama/fr/ma42 Surrelaxation pour Gauss-Seidel -- Exercice.ipynb
@@ -0,0 +1,336 @@
+{
+ "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": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "# Exercice ma21"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import numpy.linalg as lin\n",
+    "import matplotlib.pylab as plt\n",
+    "\n",
+    "%matplotlib inline\n",
+    "\n",
+    "np.set_printoptions(precision=3, linewidth=150, suppress=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "On va augmenter le rayon de convergence la méthode Jacobi amméliorée faite en TD à savoir la méthode de Gauss-Seidel.\n",
+    "\n",
+    "On étudiera sa convergence dans différents cas."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "## Gauss-Seidel\n",
+    "\n",
+    "Lorsqu'on calcul le **x** suivant avec Jacobi on ne profite pas du fait que N est triangulaire\n",
+    "et donc qu'on connait la nouvelle valeur de $x_n$ lorsqu'on calcule $x_{n-1}$. Avec Gauss-Seidel\n",
+    "on utilise toujours la dernière valeur calculée ce qui accélère la convergence.\n",
+    "\n",
+    "Pour résumer Gauss-Seidel d'un point de vu matriciel on a :\n",
+    "\n",
+    "* D = la matrice diagonale extraite de A : `D = np.diag(np.diag(A))`\n",
+    "* L = la matrice stritecement triangulaire inférieure de A : `L = np.tril(A, -1)`\n",
+    "* U = la matrice stritecement triangulaire supérieure de A : `U = np.triu(A, 1)`\n",
+    "\n",
+    "et une itération est donnée par la formule suivante :\n",
+    "\n",
+    "$$\n",
+    "(D + L)\\, {\\bf x}^{k+1} = -U\\; {\\bf x}^k + {\\bf b}\n",
+    "$$\n",
+    "ou\n",
+    "$$\n",
+    "{\\bf x}^{k+1} = D^{-1} \\, ( -L\\, {\\bf x}^{k+1} - U\\; {\\bf x}^k + {\\bf b})\n",
+    "$$\n",
+    "c.a.d.\n",
+    "$$\n",
+    "\\begin{bmatrix}\n",
+    "x_{1}^{k+1} \\\\\n",
+    "x_{2}^{k+1} \\\\\n",
+    "\\vdots \\\\\n",
+    "x_{n}^{k+1} \\\\\n",
+    "\\end{bmatrix}\n",
+    "=\n",
+    "\\begin{bmatrix}\n",
+    "1/a_{11} \\quad 0 \\quad \\ldots \\quad 0 \\\\\n",
+    "0 \\quad 1/a_{22} \\quad \\ldots \\quad 0 \\\\\n",
+    " \\vdots \\\\\n",
+    "0 \\quad 0  \\quad \\ldots \\quad 1/a_{nn} \\\\\n",
+    "\\end{bmatrix}\n",
+    "\\;\n",
+    "\\left(\n",
+    "\\;\n",
+    "-\n",
+    "\\begin{bmatrix}\n",
+    "0 \\quad 0 \\quad \\ldots \\quad 0 \\\\\n",
+    "a_{21} \\; 0 \\quad \\ldots \\quad 0 \\\\\n",
+    " \\vdots \\\\\n",
+    "a_{n1} \\, a_{n2}  \\; \\ldots \\quad 0 \\\\\n",
+    "\\end{bmatrix}\n",
+    "\\;\n",
+    "\\begin{bmatrix}\n",
+    "x_{1}^{k+1} \\\\\n",
+    "x_{2}^{k+1} \\\\\n",
+    "\\vdots \\\\\n",
+    "x_{n}^{k+1} \\\\\n",
+    "\\end{bmatrix}\n",
+    "-\n",
+    "\\begin{bmatrix}\n",
+    "0 \\; a_{12} \\; \\ldots \\; a_{1n} \\\\\n",
+    "0 \\quad 0 \\; \\ldots \\; a_{2n}  \\\\\n",
+    " \\vdots \\\\\n",
+    "0 \\quad 0  \\; \\ldots \\; 0 \\\\\n",
+    "\\end{bmatrix}\n",
+    "\\;\n",
+    "\\begin{bmatrix}\n",
+    "x_{1}^k \\\\\n",
+    "x_{2}^k \\\\\n",
+    "\\vdots \\\\\n",
+    "x_{n}^k \\\\\n",
+    "\\end{bmatrix}\n",
+    "+\n",
+    "\\begin{bmatrix}\n",
+    "b_{1} \\\\\n",
+    "b_{2} \\\\\n",
+    "\\vdots \\\\\n",
+    "b_{n} \\\\\n",
+    "\\end{bmatrix}\n",
+    "\\; \\right)\n",
+    "$$\n",
+    "\n",
+    "Notons que je peux mettre $L\\, {\\bf x}^{k+1}$ à droite du signe égal\n",
+    "si je résoud mon système ligne par ligne en commencant par le haut puisque dans\n",
+    "ce cas les ${\\bf x}^{k+1}$ utilisés sont connus. C'est ce qu'on a fait lors du dernier TP.\n",
+    "\n",
+    "### Surrelaxation de Gauss-Seidel\n",
+    "\n",
+    "Comme on a fait avec Jacobi, on introduit de l'inertie avec $w$ :\n",
+    "\n",
+    "$$\n",
+    "{\\bf x}^{k+1} = w \\, D^{-1} \\, ( -L\\, {\\bf x}^{k+1} - U\\; {\\bf x}^k + {\\bf b}) + (1-w) \\; {\\bf x}^k\n",
+    "$$\n",
+    "\n",
+    "Vérifiez que l'on arrive à l'écriture matricielle suivante :\n",
+    "\n",
+    "$$\n",
+    "\\left(\\frac{D}{w} + L\\right)\\, {\\bf x}^{k+1} = \\left(\\frac{1-w}{w} \\, D - U\\right)\\; {\\bf x}^k + {\\bf b}\n",
+    "$$\n",
+    "\n",
+    "Écrit ainsi on voit que cette méthode consiste à avoir les éléments de la diagonale des 2 cotés de l'égalité.\n",
+    "On peut interpréter cela comme un avantage lié à une meilleure répartition de l'information contenue dans la matrice (à tester pour savoir)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "### Programmons Gauss-Seidel surrelaxé\n",
+    "\n",
+    "On écrira deux fonctions :\n",
+    "\n",
+    "* `solve_triangular(L,b)` qui retourne la solution de L **x** = **b** lorsque L est triangulaire inférieure\n",
+    "* `gauss_seidel_r(A, b, x0, w, n)` qui fait `n` iteration de Gauss-Seidel en démarrant à `x0` avec `w` le coefficient de relaxation donné en argument. \n",
+    "   Cette fonction retourne un tableau des valeurs de **x** calculées (donc tableau en 2D).\n",
+    "   \n",
+    "Comme toujours, attention à limiter les `for` et à faire le plus possible d'opérations vectorielles et matricielles."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.random.seed(123)\n",
+    "A = np.random.randint(10, size=(4,4))\n",
+    "b = A.sum(axis=1)\n",
+    "x0 = np.random.random(4)\n",
+    "\n",
+    "res = gauss_seidel_r(A, b, x0, w=0.2, n=100)\n",
+    "print(res[-1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_convergences(values, result):\n",
+    "    error = np.square(values - result).sum(axis = -1) / np.square(result).sum(axis=-1)\n",
+    "    error2 = np.square(np.diff(values)).sum(axis = -1) / np.square(values).sum(axis=-1)\n",
+    "    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14,4))\n",
+    "    ax1.plot(range(len(error)), error)\n",
+    "    ax1.set_title('Erreur absolue normalisée')\n",
+    "    ax1.semilogy();\n",
+    "    ax2.plot(range(len(error2)), error2)\n",
+    "    ax2.set_title('Erreur relative normalisée')\n",
+    "    ax2.semilogy()\n",
+    "    print(\"Itération du minimum :\",np.argmin(error), np.argmin(error2))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_convergences(res, np.ones(4))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "Est-ce que la méthode de Gauss-Seidel non relaxée converge dans ce cas ?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "### Le bon cas\n",
+    "\n",
+    "Trouver un `seed` qui permet de générer un cas qui ne converge pas avec Gauss-Seidel de base mais qui \n",
+    "converge avec la relaxation ($w=0.2$)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "Tracer les courbes de convergence pour le cas retenu avec et sans relaxation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "### Étude de $w$\n",
+    "\n",
+    "Toujours dans notre cas retenu,\n",
+    "indiquer quel est l'intervale de\n",
+    "valeurs de $w$ qui garantit la convergence pour notre système matriciel A **x** = **b** avec toujours le même `x0` \n",
+    "et un nombre d'itérations à déterminer.\n",
+    "\n",
+    "Trouver la valeur optimiale de $w$ pour converger le plus rapidement pour ce cas. \n",
+    "\n",
+    "La précision demandée pour l'intervale et la valeur optimale est de $10^{-2}$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ]
+}
\ No newline at end of file