init: initial commit

1 files changed, 416 insertions, 0 deletions

diff --git a/PVCM/cama/en/ma40 Méthodes itératives.ipynb b/PVCM/cama/en/ma40 Méthodes itératives.ipynb
new file mode 100644
index 0000000..c1d93e0
--- /dev/null
+++ b/PVCM/cama/en/ma40 Méthodes itératives.ipynb
@@ -0,0 +1,416 @@
+{
+ "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": "en"
+   },
+   "source": [
+    "## Numerical simulation\n",
+    "\n",
+    "If the price of bananas is important, the sector with the greatest system resolution needs\n",
+    "matrices is numerical simulation.\n",
+    "This includes <a href=\"https://www.nas.nasa.gov/SC18/demos/demo10.html\">this</a>\n",
+    "\n",
+    "<center><img alt=\"simulation du X-57 par la Nasa\" src=\"images/nasa-x57-pression.jpg\"/></center>\n",
+    "\n",
+    "and everything that is digitally simulated, namely everything related to transport, energy, construction, everything\n",
+    "what we manufacture, which is expensive and which must have a very specific physical behavior. But that's not all,\n",
+    "understanding our environment (weather, global warming, chemistry, etc.) also means solving\n",
+    "matrix systems as well as other areas. However, there are numerical simulation methods\n",
+    "that do not go through matrix systems.\n",
+    "\n",
+    "To make the image above, we transform physical equations like [Navier-Stokes](https://fr.wikipedia.org/wiki/%C3%89quations_de_Navier-Stokes) into a matrix system where the unknowns are defined at each point\n",
+    "of a mesh to be defined. In our case the unknown is the pressure and the mesh is the interior of a\n",
+    "imaginary box that includes the plane and the air that circulates around it.\n",
+    "\n",
+    "If the box is a cube with 1000 points in each direction, we have 1 billion points in the box (minus\n",
+    "what is inside the plane but let's stay on 1 billion). Then the matrix has 1 quadrillion elements (one billion squared).\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",
+    "Inverting a matrix is ​​done in $O(n³)$ operations or $O(10^{27})$ in our case.\n",
+    "\n",
+    "If you have [the most powerful computer in the world](https://www.top500.org/) in 2024\n",
+    "you can do 1 Eflops or $10^{18}$ operations per second. Also you will need $10^{9}$ seconds or a little over 30 years. It's too much.\n",
+    "\n",
+    "**Inverting a matrix or solving by a direct method is not the right solution to solve a large matrix system.**\n",
+    "\n",
+    "#### Exercise 12.1\n",
+    "For such a simulation it is also necessary to calculate the speed of the air in each of the 1 billion points of our grid. A speed is 3 variables that must be added to the pressure. How long does it take to reverse\n",
+    "the matrix which also integrates the speed of the 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": "en"
+   },
+   "source": [
+    "It is a good idea to estimate the time of a large calculation before launching it.\n",
+    "\n",
+    "Besides, with such long times, he\n",
+    "should not stay with orders of magnitude but specify with the exact formula. So with Choleski it's $n^3/3$\n",
+    "operations therefore it runs 3 times faster, but it is still too much."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "en"
+   },
+   "source": [
+    "## Iterative Methods\n",
+    "\n",
+    "Iterative methods are methods that approach the desired solution step by step. They allow\n",
+    "to find an approximation of ${\\bf x}$, the solution of $A\\, {\\bf x} = b$.\n",
+    "\n",
+    "In general we\n",
+    "stops the calculation when we estimate that we are at a chosen distance from the solution (which we call the error) rather\n",
+    "than waiting to be at the exact solution.\n",
+    "Anyway the exact solution on a computer which is limited in number of decimal places can\n",
+    "be unattainable. So **we will never seek to have an error smaller than our maximum precision**.\n",
+    "\n",
+    "The fundamental idea of ​​the iterative algorithm is to have a formula like $\\; {\\bf x}^{t+1} = B \\, {\\bf x}^t + {\\bf c}\\;$ or in Python:\n",
+    "\n",
+    "```\n",
+    "x = np.random(size = c.size)\n",
+    "while np.square(x - old_x) > threshold:\n",
+    "    old_x = x\n",
+    "    x = B @ x + c\n",
+    "```\n",
+    "\n",
+    "The magic is if **x** converges. In this case we have reached a fixed point i.e. that ${\\bf x}^{t+1} = {\\bf x}^t$\n",
+    "and so\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",
+    "We find our $A \\; {\\bf x} = {\\bf b}$ that being in practice we do not cut A into Id and B because\n",
+    "ca does not converge (except particular cases). Let's see how the Jacobi method does it."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "en"
+   },
+   "source": [
+    "## Jacobi method\n",
+    "\n",
+    "The Jacobi method cuts the matrix A into M and N with\n",
+    "\n",
+    "*M the diagonal matrix that includes the elements of the diagonal of A* N = M - A (so 0 on the diagonal and the opposite of the elements of A elsewhere)\n",
+    "\n",
+    "So the system to solve is $\\; (M - N) \\, {\\bf x} = {\\bf b}$.\n",
+    "\n",
+    "The iterative formula is therefore for iteration k+1 expressed as a function of iteration k:\n",
+    "\n",
+    "$$\n",
+    "M \\; {\\bf x}^{k+1} =  N\\; {\\bf x}^k + {\\bf b}\n",
+    "$$\n",
+    "\n",
+    "and since M is a diagonal matrix, we have:\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_3^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_3^k  \\ldots - a_{n-1,n-1} \\, x_{n-1}^k \\qquad + b_n \\\\\n",
+    "\\end{array}\n",
+    "$$\n",
+    "\n",
+    "We see in watermark $A \\; {\\bf x} = {\\bf b}$.\n",
+    "\n",
+    "To calculate $x_1^{k+1}$ it is necessary to divide the right term of the first line by $a_{11}$ it is thus necessary that $a_{11} \\ne 0$.\n",
+    "\n",
+    "To calculate the following terms $x_i^{k+1}$ it is also necessary that $a_{ii}$ be non-zero so **it\n",
+    "A must not have zero on its diagonal**.\n",
+    "\n",
+    "This can be found in the following writing, which takes up the initial formula:\n",
+    "$ {\\bf x}^{k+1} =  M^{-1} \\;(N\\; {\\bf x}^k + {\\bf b})$\n",
+    "\n",
+    "> We see that in order to be effective, M must be easily invertible, otherwise we might as well invert A directly.\n",
+    "Here it is indeed the case since M is diagonal."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "en"
+   },
+   "source": [
+    "#### Let's program 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": "en"
+   },
+   "source": [
+    "It doesn't really converge... (if you rerun and see `NaN` it means there is a zero on the diagonal)\n",
+    "\n",
+    "2nd try:"
+   ]
+  },
+  {
+   "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": "en"
+   },
+   "source": [
+    "Magic!"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "en"
+   },
+   "source": [
+    "### Why does the 2nd case work?\n",
+    "\n",
+    "For an iterative method of the ${\\bf x} = B\\; {\\bf x} + {\\bf c}$ type to converge, it is necessary to choose\n",
+    "\n",
+    "* $\\rho(B) < 1\\quad$ with $\\rho$ the spectral radius (which is the largest eigenvalue in absolute value)\n",
+    "* $||B|| < 1\\quad$ with a matrix norm is subordinate to a vector norm.\n",
+    "             \n",
+    "\n",
+    "#### Case of the Jacobi method\n",
+    "\n",
+    "For the Jacobi matrix B, these conditions are met if A is at **dominant diagonal**, namely that each\n",
+    "diagonal element is greater in modulus than the sum of all other modulus elements in its row or column ($|a_{i,i}| \\ge \\sum_{j \\ne i} |a_{i,j}|$).\n",
+    "\n",
+    "Jacobi also converges if A is symmetric, real and positive definite\n",
+    "(i.e. $\\forall {\\bf x}, \\; {\\bf x}^T \\, A \\, {\\bf x} > 0$)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "en"
+   },
+   "source": [
+    "### Calculation time\n",
+    "\n",
+    "If we converge in 10 iterations then\n",
+    "10 matrix multiplications, 10 vector additions and 10 vector divisions are used, i.e. 180 operations\n",
+    "to be compared to the $4^3 / 3 = 22$ operations of a direct method. Damn !\n",
+    "\n",
+    "Some reasons why we lose are\n",
+    "\n",
+    "* Matrix A is too small, iterative methods are interesting for large matrices\n",
+    "* Jacobi's method is not the best (but it is very easily parallelizable)\n",
+    "* The spectral radius of the matrix is too large. The greater the spectral radius\n",
+    "  smaller and the faster we converge."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ]
+}
\ No newline at end of file