init: initial commit

1 files changed, 335 insertions, 0 deletions

diff --git a/PVCM/cama/en/ma42 Surrelaxation pour Gauss-Seidel -- Exercice.ipynb b/PVCM/cama/en/ma42 Surrelaxation pour Gauss-Seidel -- Exercice.ipynb
new file mode 100644
index 0000000..3798257
--- /dev/null
+++ b/PVCM/cama/en/ma42 Surrelaxation pour Gauss-Seidel -- Exercice.ipynb
@@ -0,0 +1,335 @@
+{
+ "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": []
+  }
+ ]
+}
\ No newline at end of file