{
 "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": []
  }
 ]
}