path: root/PVCM/cama/en/ma40 Méthodes itératives.ipynb

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