{
 "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": "en"
   },
   "source": [
    "# Exercise 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": "en"
   },
   "source": [
    "We will increase the radius of convergence the improved Jacobi method made in TD, namely the Gauss-Seidel method.\n",
    "\n",
    "We will study its convergence in different cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Gauss-Seidel\n",
    "\n",
    "When we calculate the following **x** with Jacobi we do not take advantage of the fact that N is triangular\n",
    "and therefore we know the new value of $x_n$ when we calculate $x_{n-1}$. With Gauss-Seidel\n",
    "the last computed value is always used, which accelerates convergence.\n",
    "\n",
    "To summarize Gauss-Seidel from a matrix point of view, we have:\n",
    "\n",
    "* D = the diagonal matrix extracted from A: `D = np.diag(np.diag(A))`* L = the lower strictly triangular matrix of A: `L = np.tril(A, -1)`\n",
    "* U = the upper strictly triangular matrix of A: `U = np.triu(A, 1)`\n",
    "\n",
    "and an iteration is given by the following formula:\n",
    "\n",
    "$$\n",
    "(D + L)\\, {\\bf x}^{k+1} = -U\\; {\\bf x}^k + {\\bf b}\n",
    "$$\n",
    "Where\n",
    "$$\n",
    "{\\bf x}^{k+1} = D^{-1} \\, ( -L\\, {\\bf x}^{k+1} - U\\; {\\bf x}^k + {\\bf b})\n",
    "$$\n",
    "i.e.\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",
    "Note that I can put $L\\, {\\bf x}^{k+1}$ to the right of the equal sign\n",
    "if I solve my system line by line starting from the top since in\n",
    "in this case the ${\\bf x}^{k+1}$ used are known. This is what we did during the last lab.\n",
    "\n",
    "### Gauss-Seidel overrelaxation\n",
    "\n",
    "As we did with Jacobi, we introduce inertia with $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",
    "Check that you arrive at the following matrix entry:\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",
    "Written thus we see that this method consists in splitting elements of the diagonal on both sides of the equality.\n",
    "This can be interpreted as an advantage linked to a better repartition of the information included in our matrix (it needs to be tested)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "### Let's program overrelaxed Gauss-Seidel\n",
    "\n",
    "We will write two functions:\n",
    "\n",
    "* `solve_triangular(L,b)` which returns the solution of L**x** = **b** when L is lower triangular\n",
    "* `gauss_seidel_r(A, b, x0, w, n)` which does `n` Gauss-Seidel iteration starting at `x0` with `w` the given relaxation coefficient in argument.\n",
    "   This function returns an array of calculated **x** values ​​(thus an array in 2D).\n",
    "   \n",
    "As always, be careful to limit `for` and do as many vector and matrix operations as possible."
   ]
  },
  {
   "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": "en"
   },
   "source": [
    "Does the unrelaxed Gauss-Seidel method converge in this case?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "### The good case\n",
    "\n",
    "Find a `seed` which allows to generate a case which does not converge with the basic Gauss-Seidel but which\n",
    "converges with relaxation ($w=0.2$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Plot the convergence curves for the selected case with and without 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": "en"
   },
   "source": [
    "### Study by $w$\n",
    "\n",
    "Still in our chosen case,\n",
    "indicate what is the interval of\n",
    "values ​​of $w$ which guarantees convergence for our matrix system A **x** = **b** with always the same `x0`\n",
    "and a number of iterations to be determined.\n",
    "\n",
    "Find the optimal value of $w$ to converge fastest for this case.\n",
    "\n",
    "The requested precision for the interval and the optimal value is $10^{-2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}