diff --git a/4.1 Methodes directes.ipynb b/4.1 Methodes directes.ipynb new file mode 100644 index 0000000..2720ce7 --- /dev/null +++ b/4.1 Methodes directes.ipynb @@ -0,0 +1,475 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Méthodes Directes\n", + "\n", + "D'abord quelques exemples d'utilisations de la factorisation LU et de Choleski, ensuite on passe aux methodes itéreatives" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercice \n", + "\n", + "On considère le système linéaire $A\\mathbf{x} = \\mathbf{b}$ où :\n", + "\n", + "$$A = \\left(\n", + " \\begin{array}{ccc}\n", + " 3 & 6 & 7 \\\\[0.5em]\n", + " 1 & 1 & 4 \\\\[0.5em]\n", + " 2 & 4 & 8\n", + " \\end{array} \\right), \\quad\n", + " \\mathbf{b} = \\left( \\begin{array}{c} 4 \\\\[0.5em] 5 \\\\[0.5em] 6 \\end{array} \n", + "\\right) \\, .$$\n", + " \n", + "1. Calculer la factorisation $LU$ de la matrice $A$ avec Python:\n", + "\n", + "2. Résoudre le système linéaire $A \\mathbf{x} = \\mathbf{b}$ en utilisant la factorisation trouvée au point\n", + "précédent (Ne plus utiliser Python. Mais ici on va le faire)\n", + "\n", + "3. Calculer le déterminant de la matrice $A$ en utilisant sa factorisation $LU$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# importing libraries used in this book\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# scipy.linalg.eigh : eigenvalues for symmetric matric\n", + "from scipy.linalg import eigh\n", + "# scipy.linalg.eig : eigenvalues of a matrix\n", + "from scipy.linalg import eig\n", + "# scipy.linalg.det : determinant of a matrix\n", + "from scipy.linalg import det\n", + "# scipy.linalg.lu : LU decomposition\n", + "from scipy.linalg import lu\n", + "# scipy.linalg.cholesky : Cholesky decomposition\n", + "from scipy.linalg import cholesky\n", + "from scipy.linalg import hilbert\n", + "\n", + "import pprint" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = np.array([[3, 6, 7],\n", + " [1, 1, 4],\n", + " [2, 4, 8] ])\n", + "\n", + "# LU factorisation with pivoting\n", + "P, L, U = lu(A)\n", + "\n", + "print(\"A = P L U\")\n", + "pprint.pprint(P.dot(L.dot(U)) )\n", + "\n", + "print(\"P:\")\n", + "pprint.pprint(P)\n", + "\n", + "print (\"L:\")\n", + "pprint.pprint(L)\n", + "\n", + "print (\"U:\")\n", + "pprint.pprint(U)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 2\n", + "\n", + "Si $P$ est l'identité, $A = LU$.\n", + "\n", + "1. résoudre pour $\\mathbf y$ tel que $L \\mathbf y = \\mathbf{b}$ par substitution progressive\n", + "\n", + "2. résoudre pour $\\mathbf x$ tel que $U \\mathbf x = \\mathbf{y}$ par substitution retrograde\n", + "\n", + "Si $P$ n'est pas l'identité, $A = PLU$. \n", + "\n", + "Attention, une matrice de permutation est une matrice orthogonale car les colonnes sont orthonormée.\n", + "Donc $P^{-1} = P^T$. Du coup : $P^TA = LU$ et \n", + "$$A\\mathbf x = \\mathbf b \\Leftrightarrow P^TA\\mathbf x = P^T \\mathbf b\\Leftrightarrow LU\\mathbf x = P^T \\mathbf b.$$\n", + "Alors il faut modifier les calculs précedants comme suit:\n", + "\n", + "1. résoudre pour $\\mathbf y$ tel que $L \\mathbf y = P^T\\mathbf{b}$ par substitution progressive\n", + "\n", + "2. résoudre pour $\\mathbf x$ tel que $U \\mathbf x = \\mathbf{y}$ par substitution retrograde\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# to complete\n", + "\n", + "# check the residual of the equation\n", + "print(b - A.dot(x))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 3\n", + "\n", + "$$\\det A = \\det P \\, \\det L \\, \\det U$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# to complete\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*Les mineurs principaux** d'une matrice $A \\in \\mathbb{R}^{n \\times n}$ \n", + "sont les déterminants des matrices $A_p = (a_{i,j})_{1 \\leq i, j \\leq p}$, $p = 1, ..., n$.\n", + "\n", + "**Critère de Sylvester:** une matrice symétrique $A \\in \\mathbb{R}^{n\\times n}$ \n", + "est définie positive si et seulement si les mineurs\n", + "principaux de $A$ sont tous positifs." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercice \n", + "\n", + "On considère le système linéaire $A\\mathbf{x} = \\mathbf{b}$\n", + "où\n", + "\n", + "$$A = \\left(\n", + " \\begin{array}{ccc}\n", + " \\varepsilon & 1 & 2\\\\\n", + " 1 & 3 & 1 \\\\\n", + " 2 & 1 & 3 \\\\\n", + " \\end{array}\n", + " \\right).$$\n", + " \n", + "1. Déterminer pour quelles valeurs du paramètre réel\n", + "$\\varepsilon \\in \\mathbb{R}$, la matrice $A$ est symétrique\n", + "définie positive \n", + "\n", + "2. Soit maintenant $\\varepsilon=0$. On veut résoudre le système\n", + "$A\\mathbf{x}=\\mathbf{b}$ par une méthode directe; quelle\n", + "factorisation de la matrice $A$ envisageriez-vous? Justifier\n", + "votre réponse.\n", + "\n", + "3. En considérant $\\varepsilon = 2$, vérifier que dans ce cas la\n", + "matrice $A$ est définie positive et calculer sa\n", + "factorisation de Cholesky $A = L L^T$.\n", + "\n", + "4. En supposant que $\\mathbf{b} = (1,1,1)^T$, résoudre le\n", + "système linéaire $A \\mathbf{x} = \\mathbf{b}$ en utilisant la\n", + "factorisation de Cholesky calculée au point c).\n", + "\n", + "Reference Python pour la factorisation de Cholesky :\n", + "[scipy.linalg.cholesky](https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.linalg.cholesky.html)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 1 et 2\n", + "\n", + "cf corrigé de la série" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "epsilon = 2\n", + "A = np.array([[epsilon, 1, 2],\n", + " [1, 3, 1],\n", + " [2, 1, 3] ])\n", + "\n", + "# Is A symmetric ?\n", + "print(f'max(abs(A-AT)) = {np.max(np.abs(A-A.T))} ')\n", + "\n", + "# What are the eigenvlues of $A$ ? (usually, use eig, but here A is symmetric, we can use eigh )\n", + "lk, v = eigh(A)\n", + "print(f'The eigenvalues of A are {lk}')\n", + "\n", + "# Cholesky factorisation: lower : return lower-triangular matrix, A = L L^T\n", + "L = cholesky(A, lower=True) \n", + "print(f\"\\n A = L^T L = {L.dot(L.T)}\\n\")\n", + "\n", + "print (f\"L = {L}\")\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 4" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# to complete\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Problèmes de précision\n", + "\n", + "Les erreurs d'arrondis peuvent causer des différences importantes\n", + "entre la solution calculée par la méthode d'élimination de Gauss (MEG)\n", + "et la solution exacte. Cela arrive si le\n", + "*conditionnement* de la matrice du système est très grand.\n", + "\n", + "La matrice de Hilbert de taille $n\\times n$ est une matrice\n", + "symétrique définie par\n", + "\n", + "$$a_{ij} = \\frac{1}{i+j-1}, \\qquad i,j=1,\\ldots,n$$\n", + "\n", + "On peut construire une matrice de Hilbert de taille $n$\n", + "quelconque en utilisant la commande `A = scipy.linalg.hilbert(n)`. Par\n", + "example, pour $n = 4$, on a:\n", + "$$A = \\begin{bmatrix}\n", + " 1 & \\frac{1}{2} & \\frac{1}{3} & \\frac{1}{4} \\\\\n", + " \\frac{1}{2} & \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} \\\\\n", + " \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} \\\\\n", + " \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7}\n", + "\\end{bmatrix}$$\n", + "\n", + "On considère les systèmes linéaires $A_n \\mathbf {x}_n=\\mathbf {b}_n$ ou $A_n$ est la matrice de Hilbert de taille\n", + "$n$ avec $n=4,6,8,10,12,14,\\ldots, 20$ tandis que $\\mathbf {b}_n$ est choisi de sorte que la solution exacte soit $\\mathbf{x}_n= (1,1,\\cdots, 1)^T$.\n", + "\n", + "1. Pour chaque $n$, calculez le conditionnement de la matrice\n", + "\n", + "2. Résolvez le système linéaire par la factorisation $LU$ et notez $\\vec{x}^{LU}_n$ la solution calculée.\n", + "\n", + "3. Dessinez le graphique avec le conditionnement obtenu ainsi que l'erreur rélatif\n", + " $\\Vert \\mathbf x_n -\\mathbf x^{LU}_n\\Vert/\\Vert \\mathbf x_n\\Vert$ (où $\\Vert\\cdot\\Vert$\n", + " est la norme euclidienne d'un vecteur, \n", + " $\\Vert\\mathbf {x}\\Vert = \\sqrt{\\mathbf {x}^T\\cdot\\mathbf {x}}$). Utilisez ue échelle logarithmique pour l'axe $y$.\n", + " \n", + "4. Sur le même graphique, reportez le conditionnement de la matrice $A$, `np.linalg.cond(A)`\n", + "\n", + "Repetez la même chose avec la factorisation de Cholesky `L = cholesky(A, lower=True)` pour $n=4,6,8,10,12$. Que se passe-t-il si $n=14$\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Nrange = range(2,20,2)\n", + "err = []\n", + "cond = []\n", + "\n", + "for n in Nrange :\n", + " A = hilbert(n)\n", + " P, L, U = lu(A)\n", + "\n", + " x = np.ones([n,1])\n", + "\n", + " b = A.dot(x)\n", + "\n", + " # to complete (2-3 lines)\n", + " \n", + "\n", + " err.append( np.linalg.norm(x-xLU) / np.linalg.norm(x) )\n", + " cond.append( np.linalg.cond(A) )\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "plt.plot(Nrange, err, 'b:.',Nrange, cond, 'g:.')\n", + "\n", + "plt.xlabel('n'); plt.ylabel('err');\n", + "# plt.xscale('log') \n", + "plt.yscale('log')\n", + "plt.grid(True)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Nrange = range(2,13,2)\n", + "err = []\n", + "cond = []\n", + "\n", + "for n in Nrange :\n", + " A = hilbert(n)\n", + " L = cholesky(A, lower=True)\n", + "\n", + " x = np.ones([n,1])\n", + "\n", + " b = A.dot(x)\n", + "\n", + " # to complete (2 lines)\n", + "\n", + " err.append( np.linalg.norm(x-xCho) / np.linalg.norm(x) )\n", + " cond.append( np.linalg.cond(A) )\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "plt.plot(Nrange, err, 'b:.',Nrange, cond, 'g:.')\n", + "\n", + "plt.xlabel('n'); plt.ylabel('err');\n", + "# plt.xscale('log') \n", + "plt.yscale('log')\n", + "plt.grid(True)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n = 14\n", + "A = hilbert(n)\n", + "# L = cholesky(A, lower=True)\n", + "lk, v = eigh(A)\n", + "\n", + "print(f'lk[0] = {lk[0]:0.5e} , the matrix is not positive definite. Let''s verify with the first eigenvector')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "La première valeur propre est négative, cela signifie que $(\\mathbf v, A\\mathbf v) = (\\mathbf v, \\lambda\\mathbf v) = \\lambda (\\mathbf v, \\mathbf v) = \\lambda \\|\\mathbf v\\|^2<0$, où $\\mathbf v$ est le vecteur propre associé à cette valeur propre $\\lambda$ négative.\n", + "\n", + "Vérifions:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print ( A.dot(v[0]).T.dot(v[0]) )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Pourquoi ? En vérité, le calcul des valeurs et vecteurs propres de $A$ a aussi un problème d'arrondi à cause du conditionnement. On peut en effet vérifier que le rapport entre les composantes de $\\mathbf v$ et $A\\mathbf v$ n'est ni égale à $\\lambda$, ni une autre constante :" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(v[0]/A.dot(v[0]) )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/4.2 Methodes iteratives.ipynb b/4.2 Methodes iteratives.ipynb new file mode 100644 index 0000000..240af56 --- /dev/null +++ b/4.2 Methodes iteratives.ipynb @@ -0,0 +1,646 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Méthodes itératives" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# importing libraries used in this book\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "\n", + "# scipy.linalg.eig : eigenvalues of a matric\n", + "from scipy.linalg import eig\n", + "# scipy.linalg.lu : LU decomposition\n", + "from scipy.linalg import lu\n", + "# scipy.linalg.cholesky : Cholesky decomposition\n", + "from scipy.linalg import cholesky\n", + "from scipy.linalg import hilbert" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exemple\n", + "\n", + "$$A=\\begin{pmatrix}\n", + "1 & 2 & 3 & 4 \\\\\n", + "5 & 6 & 7 & 8\\\\\n", + "9 & 10 & 11 & 12\\\\\n", + "13 & 14 & 15 & 16\n", + "\\end{pmatrix}.$$\n", + "\n", + "On a alors\n", + "$$D=\\begin{pmatrix}\n", + "1 & 0 & 0 & 0 \\\\\n", + "0 & 6 & 0 & 0\\\\\n", + "0 & 0 & 11 & 0\\\\\n", + "0 & 0 & 0 & 16\n", + "\\end{pmatrix};$$\n", + "la matrice d'itération de la méthode de Jacobi est\n", + "$$B_J = D ^{-1}(D - A) = I - D^{-1}A =\n", + "\\begin{pmatrix}\n", + "0 & -2 & -3 & -4 \\\\\n", + "-5/6 & 0 & -7/6 & -4/3\\\\\n", + "-9/11 & -10/11 & 0 & -12/11\\\\\n", + "-13/16 & -14/16 & -15/16 & 0\n", + "\\end{pmatrix}.$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Jacobi\n", + "# first diag extract the diagonal of A, the second builds a diagonal matrix starting from a vector\n", + "D = np.diag(np.diag(A)) \n", + "\n", + "Bj = np.linalg.solve(D,(D-A))\n", + "# What are the eigenvlues of $A$ ? \n", + "lk, v = eig(Bj)\n", + "print(f'The eigenvalues of Bj are {lk}\\n')\n", + "print(f'The spectral radius of Bj is {np.max(np.abs(lk)):.4f} > 1, therefore Jacobi does not converge')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Gauss-Seidel\n", + "# Return a copy of an array with elements above the k-th diagonal zeroed. \n", + "Pgs = np.tril(A,k=0)\n", + "\n", + "Bgs = np.linalg.solve(Pgs,(Pgs-A))\n", + "# What are the eigenvlues of $A$ ? \n", + "lk, v = eig(Bgs)\n", + "print(f'The eigenvalues of Bgs are {lk}\\n')\n", + "print(f'The spectral radius of Bgs is {np.max(np.abs(lk)):.4f} > 1, therefore Gauss-Seidel does not converge')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Exercice \n", + "reprenez cet exemple avec la matrice \n", + "```python\n", + "A = np.array([[3,2,1,0],[3,4,3,2],[3,2,5,2],[1,2,3,4]])\n", + "```\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exemple\n", + "\n", + "$$A = \\left(\n", + " \\begin{array}{ccc}\n", + " 3 & 2 & 1 \\\\[0.5em]\n", + " 1 & 4 & 1 \\\\[0.5em]\n", + " 2 & 4 & 8\n", + " \\end{array} \\right), \\quad\n", + " \\mathbf{b} = \\left( \\begin{array}{c} 4 \\\\ 5 \\\\ 6 \\end{array} \n", + "\\right) \\, .$$\n", + "\n", + "Exprimez (sans Calculer) $P$, $B$, et $\\rho(B)$ dans le cas de Jacobi et Gauss-Seidel.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = np.array([[3, 2, 1],\n", + " [1, 4, 1],\n", + " [2, 4, 8] ])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Jacobi\n", + "\n", + "# To complete\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Gauss-Seidel\n", + "\n", + "# To complete\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Méthode de Richardson\n", + "\n", + "$A$ et $\\mathbf b$ donnés, on cherche à approximer la solution $\\mathbf x$ de $A\\mathbf x = \\mathbf b$.\n", + "\n", + "Soit $\\mathbf x^{(0)}$ donné, $\\mathbf r^{(0)}= \\mathbf b - A\\mathbf x^{(0)}$ pour $k=0,1,2,...$ :\n", + "\n", + "\\begin{eqnarray*}\n", + "\\text{trouver }\\mathbf{z}^{(k)} \\text{ tel que } && P \\mathbf{z}^{(k)} = \\mathbf{r}^{(k)}\\\\\n", + "\\text{choisir }\\alpha_k && \\\\\n", + "&& \\mathbf{x}^{(k+1)} = \\mathbf{x}^{(k)} + \\alpha_k \\mathbf{z}^{(k)} \\\\\n", + "&& \\mathbf{r}^{(k+1)} = \\mathbf{r}^{(k)} - \\alpha_k A \\mathbf{z}^{(k)}.\n", + "\\end{eqnarray*}\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Exemples: Jacobi et Gauss-Seidel avec relaxation\n", + "\n", + "Nous avons vu que les méthodes de Jacobi appliquée à la matrice\n", + "$$A=\\begin{pmatrix}\n", + "3 &2&1&0\\\\3&4&3&2\\\\3&2&5&2\\\\1&2&3&4\n", + "\\end{pmatrix}$$\n", + "ne converge pas, alors que Gauss-Seidel converge.\n", + "\n", + "Nous allons utiliser la méthode de Richardson stationnaire avec un paramètre de relaxation $\\alpha$ (donc constant) constant pour voir si on peut améliorer la convergence.\n", + "\n", + "#### Partie 1 matrices d'itérations\n", + "\n", + "La matrice d'itérations pour la méthod de Richardson est $B(\\alpha_k) = I - P^{-1}A$. En particulier, pour $\\alpha$ constant, nous pouvons faire un graphique du rayon spectral $\\rho(B(\\alpha))$ dans les cas de Jacobi et Gauss-Seidel:\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = np.array([[3,2,1,0],[3,4,3,2],[3,2,5,2],[1,2,3,4]])\n", + "D = np.diag(np.diag(A)) \n", + "Pgs = np.tril(A,k=0)\n", + "\n", + "Alphas = np.linspace(0, 2, 3001)\n", + "\n", + "rhoBj = []\n", + "rhoBgs = []\n", + "\n", + "\n", + "for alpha in Alphas :\n", + " Bgs = np.linalg.solve(Pgs,(Pgs-alpha*A))\n", + " Bj = np.linalg.solve(D,(D-alpha*A))\n", + "\n", + " lk, v = eig(Bj)\n", + " rhoBj.append( np.max(np.abs(lk)) )\n", + "\n", + " lk, v = eig(Bgs)\n", + " rhoBgs.append( np.max(np.abs(lk)) )\n", + "\n", + " \n", + "plt.plot(Alphas, rhoBj, 'b:', Alphas, rhoBgs, 'g:')\n", + "\n", + "plt.xlabel('alpha'); plt.ylabel('rho');\n", + "plt.title('Spectral radius')\n", + "plt.grid(True)\n", + "plt.legend(['rho(B_J(alpha))','rho(B_GS(alpha))'])\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "On constate que:\n", + "\n", + "1. En utilisant $P=D$, Richardson converge pour $\\alpha\\lessapprox 0.8$ and the optimal one is for $\\alpha \\approx 0.75$\n", + "2. En utilisant $P$ égale au triangle inférieur de $A$, Richardson converge pour $\\alpha\\lessapprox 2$ and the optimal one is for $\\alpha \\approx 1.5$\n", + "\n", + "#### Partie 2 Implémentation de la Méthode de Richardson\n", + "\n", + "Définissez une fonction Python qui implémente la méthode de Richardson stationnaire et qui a la structure suivante:\n", + "```python\n", + "def Richardson(A, b, x0, P, alpha, maxIterations, tolerance) :\n", + " # Richardson method to approximate the solution of Ax=b\n", + " # x0 : initial guess\n", + " # P : préconditioner\n", + " # alpha : relaxiation parameter\n", + " # maxIterations : maximum number of iterations\n", + " # tolerance : tolerance for relative residual\n", + "```\n", + "\n", + "La méthode de Richardson, comme toute méthode iterative a besoin d'un critere d'arrêt. Dans ce cas, le plus simple est de poser les critères suivants:\n", + "\n", + "1. On fixe le nombre maximal d'itérations\n", + "2. Le résidu rélatif est plus petit qu'une tolérance $\\epsilon$ : \n", + " $$ \\frac{\\| \\mathbf b - A \\mathbf x^{(k)} \\|}{\\|\\mathbf b \\|} < \\epsilon$$\n", + "On s'arrête dès que l'un des ces critières est rempli\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def Richardson(A, b, x0, P, alpha, maxIterations, tolerance) :\n", + " # Richardson method to approximate the solution of Ax=b\n", + " # x0 : initial guess\n", + " # P : préconditioner\n", + " # alpha : relaxiation parameter\n", + " # maxIterations : maximum number of iterations\n", + " # tolerance : tolerance for relative residual\n", + " # Outputs :\n", + " # xk : approximate solution\n", + " # vector of relative norm of the residuals\n", + " \n", + "\n", + " # To complete\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 3 Utilisation de la Méthode de Richardson\n", + "\n", + "Uliliser la fonction `Richardson` pour approcher la solution de $A\\mathbf x = \\mathbf b$ pour $\\mathbf b = (0, 1, -1, 1)^T$ avec un tolerance sur le residu rélatif de $10^{-6}$ et\n", + "le vecteur null comme point initial, et un nombre maximum d'itération de $200$\n", + "\n", + "1. En utilisant $P=D$, avec : $\\alpha = 0.7, 0.75, 0.79,0.81, 0.9$\n", + "2. En utilisant $P$ égale au triangle inférieur de $A$, avec : $\\alpha = 1, 1.5, 1.95,2.05$ \n", + "\n", + "Comment intérprétez=vous ces résultats ?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "b = np.array([0,1,-1,1])\n", + "# Define the initial guess\n", + "x0 = 0*b # trick to have te right size for x0\n", + "\n", + "tolerance = 1e-6\n", + "maxIter = 200" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Jacobi-like precondtioner\n", + "P = np.diag(np.diag(A))\n", + "legend = []\n", + "\n", + "for alpha in [0.7, 0.75, 0.79,0.81, 0.9] :\n", + " x,relRes = Richardson(A,b,x0,P,alpha,maxIter,tolerance)\n", + " print(relRes[-1])\n", + " plt.plot( relRes, ':')\n", + " legend.append(str(alpha))\n", + "\n", + "plt.legend(legend)\n", + "\n", + "plt.xlabel('n'); plt.ylabel('res');\n", + "plt.title('Relative residuals')\n", + "plt.yscale('log')\n", + "plt.show()\n", + "\n", + " \n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Gauss-Seidel-like precondtioner\n", + "P = np.tril(A,k=0)\n", + "legend = []\n", + "\n", + "for alpha in [1, 1.5, 1.95,2.05] :\n", + " x,relRes = Richardson(A,b,x0,P,alpha,maxIter,tolerance)\n", + " print(relRes[-1])\n", + " plt.plot( relRes, ':')\n", + " legend.append(str(alpha))\n", + "\n", + "plt.legend(legend)\n", + "\n", + "plt.xlabel('n'); plt.ylabel('res');\n", + "plt.title('Relative residuals')\n", + "plt.yscale('log')\n", + "plt.show()\n", + "\n", + " \n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 4 Implémentation de la Méthode du Gradient Préconditionné\n", + "\n", + "Définissez une fonction Python qui implémente la méthode de Gradient Préconditionné et qui a la structure suivante:\n", + "```python\n", + "def PrecGradient(A, b, x0, P, maxIterations, tolerance) :\n", + " # Preconditioned Gradient method to approximate the solution of Ax=b\n", + " # x0 : initial guess\n", + " # P : préconditioner\n", + " # maxIterations : maximum number of iterations\n", + " # tolerance : tolerance for relative residual\n", + "```\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def PrecGradient(A, b, x0, P, maxIterations, tolerance) :\n", + " # Preconditioned Gradient method to approximate the solution of Ax=b\n", + " # x0 : initial guess\n", + " # P : préconditioner\n", + " # maxIterations : maximum number of iterations\n", + " # tolerance : tolerance for relative residual\n", + " # Outputs :\n", + " # xk : approximate solution\n", + " # vector of relative norm of the residuals\n", + " \n", + " # To complete. Start from Richardson, it is easier\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 5 Implémentation de la Méthode du Gradient Conjugué Préconditionné\n", + "\n", + "Définissez une fonction Python qui implémente la méthode de Gradient Préconditionné et qui a la structure suivante:\n", + "```python\n", + "def PrecConjugateGradient(A, b, x0, P, maxIterations, tolerance) :\n", + " # Preconditionate Conjugate Gradient method to approximate the solution of Ax=b\n", + " # x0 : initial guess\n", + " # P : préconditioner\n", + " # maxIterations : maximum number of iterations\n", + " # tolerance : tolerance for relative residual\n", + "```\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def PrecConjugateGradient(A, b, x0, P, maxIterations, tolerance) :\n", + " # Preconditionate Conjugate Gradient method to approximate the solution of Ax=b\n", + " # x0 : initial guess\n", + " # P : préconditioner\n", + " # maxIterations : maximum number of iterations\n", + " # tolerance : tolerance for relative residual\n", + " # Outputs :\n", + " # xk : approximate solution\n", + " # vector of relative norm of the residuals\n", + " \n", + " # we do not keep track of all the sequence, just the last two entries\n", + " xk = x0\n", + " rk = b - A.dot(xk)\n", + " zk = np.linalg.solve(P,rk)\n", + " pk = zk\n", + " \n", + " # To complete. Start from PrecGradient, it is easier" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Partie 6 Utilisation de la Méthode de Gradient Preconditioné et du Gradient Conjugué Préconditionné\n", + "\n", + "Uliliser les fonctions `PrecGradient` et `PrecConjugateGradient` pour approcher la solution de $A\\mathbf x = \\mathbf b$ pour $\\mathbf b = (0, 1, -1, 1)^T$ avec un tolerance sur le residu rélatif de $10^{-6}$ et\n", + "le vecteur null comme point initial, et un nombre maximum d'itération de $200$\n", + "\n", + "1. En utilisant $P=D$\n", + "2. En utilisant $P$ égale au triangle inférieur de $A$\n", + "\n", + "Comment intérprétez=vous ces résultats ?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "legend = []\n", + "\n", + "# Jacobi-like precondtioner\n", + "P = np.diag(np.diag(A))\n", + "\n", + "x,relRes = PrecGradient(A,b,x0,P,maxIter,tolerance)\n", + "print(relRes[-1])\n", + "plt.plot( relRes, ':')\n", + "legend.append('Gradient, P=D')\n", + "\n", + "repeat the same for\n", + "# PrecConjugateGradient, then\n", + "# Gauss-Seidel-like precondtioner (both Gradient and Conjugate Gradient)\n", + "\n", + "\n", + "\n", + "plt.legend(legend)\n", + "\n", + "plt.xlabel('n'); plt.ylabel('res');\n", + "plt.title('Relative residuals')\n", + "plt.yscale('log')\n", + "plt.show()\n", + "\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Autres Exemples\n", + "\n", + "Repetez les expériences numériques ci dessous pour approcher les solution de\n", + "\n", + "$$\\left(\n", + " \\begin{array}{ccc}\n", + " 2 & 1 \\\\[0.5em]\n", + " 1 & 3 \n", + " \\end{array} \\right)\n", + " \\mathbf x = \n", + " \\left( \\begin{array}{c} 1 \\\\ 0 \\end{array} \\right) \\quad,\\quad \n", + " \\left(\n", + " \\begin{array}{ccc}\n", + " 2 & 2 \\\\[0.5em]\n", + " -1 & 3 \n", + " \\end{array} \\right)\n", + " \\mathbf x = \n", + " \\left( \\begin{array}{c} 1 \\\\ 0 \\end{array} \\right) \\quad\\text{et}\\quad \n", + " \\left(\n", + " \\begin{array}{ccc}\n", + " 5 & 7 \\\\[0.5em]\n", + " 7 & 10 \n", + " \\end{array} \\right)\n", + " \\mathbf x = \n", + " \\left( \\begin{array}{c} 1 \\\\ 0 \\end{array} \\right) .$$\n", + " \n", + "Pour quelles matrices et méthodes vous vous attendez convergence ? Pour la dernière matrice, le résultat n'est pas celui attendu pourquoi ?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "b = np.array([1,0])\n", + "# Define the initial guess\n", + "x0 = np.array([1,1])\n", + "\n", + "tolerance = 1e-6\n", + "maxIter = 10\n", + "A = np.array([[2,1],[1,3]])\n", + "# A = np.array([[2,2],[-1,3]])\n", + "# A = np.array([[5,7],[7,10]])\n", + "\n", + "legend = []\n", + "\n", + "# To complete: similarly to previous example:\n", + "# Jacobi, Gauss-Seidel (simple)\n", + "# PrecGrandient and ConjugateGradient, with\n", + "# both Jacobi-like and Gauss-seidel-like\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "plt.legend(legend)\n", + "\n", + "plt.xlabel('n'); plt.ylabel('res');\n", + "plt.title('Relative residuals')\n", + "plt.yscale('log')\n", + "plt.show()\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Matrice de Hilbert\n", + "\n", + "Peut-on utiliser une la méthode de Richardson pour resoudre le problème du mauvais conditionnement de la matrice de Hilbert ?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n = 10\n", + "A = hilbert(n)\n", + "xexact = np.ones(n)\n", + "b = A.dot(xexact)\n", + "\n", + "# Define the initial guess\n", + "x0 = b\n", + "\n", + "tolerance = 1e-6\n", + "maxIter = 10\n", + "\n", + "legend = []\n", + "\n", + "# To complete: similarly to previous example:\n", + "# Jacobi, Gauss-Seidel (simple)\n", + "# PrecGrandient and ConjugateGradient, with\n", + "# both Jacobi-like and Gauss-seidel-like\n", + "\n", + "\n", + "\n", + "\n", + "plt.legend(legend)\n", + "\n", + "plt.xlabel('n'); plt.ylabel('res');\n", + "plt.title('Relative residuals')\n", + "plt.yscale('log')\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}