diff options
Diffstat (limited to 'PVCM/cama/en/ma40 Méthodes itératives.ipynb')
| -rw-r--r-- | PVCM/cama/en/ma40 Méthodes itératives.ipynb | 416 |
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 |