diff --git a/4.1 Methodes directes.ipynb b/4.1 Methodes directes.ipynb
new file mode 100644
index 0000000..b515ed9
--- /dev/null
+++ b/4.1 Methodes directes.ipynb
@@ -0,0 +1,538 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Résolution de systèmes linéaires\n",
+ "\n",
+ "## Méthodes Directes\n",
+ "\n",
+ "D'abord quelques exemples d'utilisations de la factorisation LU et de Cholesky, ensuite on passe aux methodes itératives"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 4.1.1 Exemple: réseau de capillaires\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 4.1.2 Résolution de systèmes linéaires\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 4.1.3 Factorisation LU\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 4.1.4 Critères pour la factorisation LU\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Exercice 1\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 notebook\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "# scipy.linalg.eigh : eigenvalues for symmetric matrix\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",
+ "# scipy.linalg.hilbert : Hilbert matrix\n",
+ "from scipy.linalg import hilbert\n",
+ "\n",
+ "import pprint\n",
+ "\n",
+ "np.set_printoptions(precision=4, suppress=True, linewidth=120)"
+ ]
+ },
+ {
+ "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": [
+ "### 4.1.5 Factorisation LU avec Pivoting\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": [
+ "### 4.1.6 Factorisation de Choleski\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Exercice 2\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 \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 3.\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",
+ "\n",
+ "# What are the eigenvlues of A ? (usually use eig, but here A is symmetric so we can use eigh )\n",
+ "lk, v = eigh(A)\n",
+ "print(f'The eigenvalues of A are {lk}')\n",
+ "\n",
+ "# Cholesky factorisation: A = L L^T \n",
+ "L = cholesky(A, lower=True) # lower=True -> returns the lower-triangular matrix L\n",
+ "print(\"\\nA = L^T L = \")\n",
+ "pprint.pprint(L.dot(L.T))\n",
+ "\n",
+ "print (f\"\\nL = \")\n",
+ "pprint.pprint(L)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Partie 4"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# to complete\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 4.1.7 Problèmes de précision\n",
+ ""
+ ]
+ },
+ {
+ "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 $A$\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 \\mathbin{/} \\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 une échelle logarithmique pour l'axe $y$.\n",
+ " \n",
+ "4. Sur le même graphique, réportez le conditionnement de la matrice $A$ (commande: `np.linalg.cond(A)`)\n",
+ "\n",
+ "Répétez 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"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(Nrange, err, 'b:.', label=\"Erreur\")\n",
+ "plt.plot(Nrange, cond, 'g:.', label=\"Conditionemment\")\n",
+ "\n",
+ "plt.xlabel('n'); plt.ylabel('err');\n",
+ "# plt.xscale('log') \n",
+ "plt.yscale('log')\n",
+ "plt.grid(True)\n",
+ "plt.legend()\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": [
+ "plt.plot(Nrange, err, 'b:.', label=\"Erreur\")\n",
+ "plt.plot(Nrange, cond, 'g:.', label=\"Conditionemment\")\n",
+ "\n",
+ "plt.xlabel('n'); plt.ylabel('err');\n",
+ "# plt.xscale('log') \n",
+ "plt.yscale('log')\n",
+ "plt.grid(True)\n",
+ "plt.legend()\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 de son 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": []
+ }
+ ],
+ "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.7.16"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}