{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8ea8318a",
   "metadata": {},
   "source": [
    "# L3 Analyse Numérique – TP3\n",
    "\n",
    "[Dequay A](mailto:antoine.dequay@ens-rennes.fr) &\n",
    "[Le Barbenchon P](mailto:pierre.le-barbenchon@ens-rennes.fr). TP ENS Rennes\n",
    "\n",
    "[Boutin B](mailto:benjamin.boutin@univ-rennes1.fr). Cours et TP Université de Rennes 1 - UFR Mathématiques  \n",
    "\n",
    "Ce TP a pour objet les méthodes directes et itératives pour résoudre $Ax=b$.\n",
    "\n",
    "Ainsi dans tout ce TP, on n'utilisera jamais la commande $\\verb\"np.linalg.solve\"$ pour résoudre les systèmes linéaires.\n",
    "\n",
    "- Exercice 1. Résolution de systèmes triangulaires\n",
    "- Exercice 2. Algorithme de décomposition LU\n",
    "- Exercice 3. Comparaison des méthodes de Jacobi, Gauss-Seidel et de relaxation.\n",
    "- Exercice 4. Résolution d'un problème de Poisson 2D"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74bb2a34",
   "metadata": {},
   "source": [
    "Nous importons avant tout quelques librairies qui seront utiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "28b598aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "980ee6a8",
   "metadata": {},
   "source": [
    "## Exercice 1. Résolution de systèmes triangulaires\n",
    "\n",
    "<b>Théorème (Décompostion LU)</b> Soit $A\\in\\mathcal{M}_n(\\mathbb{K})$, dont tous les mineurs principaux sont non nuls, c'est-à-dire que pour tout $k\\in\\{1,\\ldots,n\\}$ la sous-matrice extraite $(A_{i,j})_{1\\leq i,j\\leq k}$ est inversible.\n",
    "\n",
    "Alors il existe un unique couple $(L,U)$ de matrices carrées d'ordre $n$ tel que :\n",
    "- $L$ est triangulaire inférieure avec uniquement des 1 sur la diagonale;\n",
    "- $U$ est triangulaire supérieure avec des coefficients diagonaux tous non nuls;\n",
    "- $A=LU$.\n",
    "\n",
    "Une fois cette décomposition connue, la résolution d'un système linéaire $Ax=b$ est équivalente à la résolution de deux systèmes linéaires triangulaires, élémentaires à résoudre :\n",
    "$$\n",
    "(Ax=b) \\Longleftrightarrow (Ly=b\\textrm{ et }Ux=y).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5880772d",
   "metadata": {},
   "source": [
    "### Question 1)\n",
    "\n",
    "Programmer l'algorithme ```Lower(L,b)``` premettant de résoudre des systèmes triangulaires inférieurs ($Ly=b$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4f90b69",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b578ae67",
   "metadata": {},
   "source": [
    "### Question 2)\n",
    "\n",
    "Tester cet algorithme sur une matrice $M$ telle que $M_{i,j} = 1+i+j$ pour tous $0 \\leq i < j \\leq n-1$, et $M_{i,i} = 1$ et pour des $b$ aléatoires."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88d2c720",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5a73912f",
   "metadata": {},
   "source": [
    "### Question 3) \n",
    "\n",
    "Programmer l'algorithme ```Upper(U,b)``` premettant de résoudre des systèmes triangulaires supérieurs ($Uy=b$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1238ac49",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "69d1f5c3",
   "metadata": {},
   "source": [
    "## Exercice 2. Algorithme de décomposition LU\n",
    "\n",
    "L'algorithme de construction des matrices $L$ et $U$ s'appuie sur les opérations d'élimination du pivot de Gauss.\n",
    "Pour réduire de moitié la complexité en espace (<b>forme compacte</b>), l'algorithme agit directement sur la matrice $A$, la matrice $U$ étant construite au fil des itérations dans la partie supérieure de $A$, et la matrice $L$ correspondant aux opérations du pivot étant stockée dans la partie inférieure de $A$. On a ainsi à l'issue de l'algorithme $U_{i,j}=A_{i,j}$ si $i\\leq j$, et $L_{i,j}=A_{i,j}$ si $i>j$, les autres coefficients de $L$ et $U$ étant nuls ou égaux à 1.\n",
    "\n",
    "\n",
    "\\begin{align*}\n",
    "&\\textbf{Algorithme de la factorisation LU.}\\\\\n",
    "&\\textbf{Forme compacte}\\\\\n",
    "&\\text{Pour }k=1,\\ldots,n-1\\\\\n",
    "&\\phantom{.}\\hspace{2em}\\text{Pour } i=k+1,\\ldots, n\\\\\n",
    "&\\phantom{.}\\hspace{4em}A_{i,k}\\leftarrow A_{i,k}/A_{k,k}\\\\\n",
    "&\\phantom{.}\\hspace{4em}\\text{Pour } j=k+1,\\ldots,n\\\\\n",
    "&\\phantom{.}\\hspace{6em}A_{i,j}\\leftarrow A_{i,j}-A_{i,k}A_{k,j}\\\\\n",
    "&\\phantom{.}\\hspace{4em}\\text{Fin pour } j\\\\\n",
    "&\\phantom{.}\\hspace{2em}\\text{Fin pour }i\\\\\n",
    "&\\text{Fin pour }k.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbe5e1f0",
   "metadata": {},
   "source": [
    "Attention ! Quand on définit un ```np.array``` avec des valeurs entières, Python ne fait des opérations que dans les entiers. Ainsi pour la matrice $A = \\begin{pmatrix} 14 & 4 \\\\ 5 & 7 \\end{pmatrix}$, l'algorithme renverra la matrice $\\begin{pmatrix} 14 & 4 \\\\ 0 & 7 \\end{pmatrix}$, alors que la vraie décomposition compacte LU est $\\begin{pmatrix} 14 & 4 \\\\ \\frac{5}{14} & \\frac{39}{7} \\end{pmatrix}$. Pour palier ce problème d'opérations faites dans le monde des entiers, il faut, en définissant la matrice $A$ :\n",
    "- soit mettre des flottants : ```A = np.array([[14.,4.],[5.,7.]])```\n",
    "- soit dire qu'on manipule des flottants : ```A = np.array([[14,4],[5,7]], dtype = 'float64')```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d55c4800",
   "metadata": {},
   "source": [
    "### Question 1)\n",
    "\n",
    "Programmer cet algorithme (on essayera si possible d'obtenir une programmation concise, idéalement avec seulement une boucle sur $k$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c2cbb96",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ac32cd23",
   "metadata": {},
   "source": [
    "### Question 2)\n",
    "\n",
    "Programmer une fonction ```Pascal(n)``` renvoyant la matrice de Pascal. Par exemple la matrice de Pascal de dimension 5 est la suivante \n",
    "$$\n",
    "P = \\begin{pmatrix}\n",
    "1 & 1 & 1 & 1 & 1\\\\\n",
    "1 & 2 & 3 & 4 & 5\\\\\n",
    "1 & 3 & 6 & 10 & 15\\\\\n",
    "1 & 4 & 10 & 20 & 35\\\\\n",
    "1 & 5 & 15 & 35 & 70\n",
    "\\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0278d82",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "682807bf",
   "metadata": {},
   "source": [
    "### Question 3)\n",
    "\n",
    "Tester la décomposition LU sur la matrice de Pascal pour plusieurs $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59ef7d01",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9f5612fa",
   "metadata": {},
   "source": [
    "## Exercice 3. Comparaison des méthodes de Jacobi, Gauss-Seidel et de relaxation\n",
    "\n",
    "Le but de la suite de ce TP est d'introduire et de comparer quelques méthodes itératives pour résoudre le système $Ax=b$.\n",
    "\n",
    "<b>Principe :</b> La matrice $A$ est décomposée sous la forme $A=M-N$ avec $M$ facile et peu coûteuse à inverser et la solution $x$ est approchée par une suite $(x_k)_{k\\in\\mathbb N}$ d'éléments de $\\mathbb R^n$:\n",
    "$$\\begin{align*} \n",
    "x_{0} &\\in \\mathbb R^n,\\\\\n",
    "Mx_{k+1} &= Nx_{k}+b.\n",
    "\\end{align*}$$\n",
    "La matrice $M^{-1}N$ est fondamentale dans l'analyse de convergence. En effet la solution $x$ est l'unique point fixe d'une application affine sur $\\mathbb R^n$ dont la différentielle (constante) a pour matrice $M^{-1}N$. En particulier la convergence itérative de la méthode, requise pour tout choix de l'initialisation $x_0$, est équivalente à la condition suivante sur le rayon spectral:\n",
    "$$\\rho(M^{-1}N)<1.$$\n",
    "\n",
    "Quelques possibilités pour construire $M$ et $N$ : écrivons $A=D-E-F$ avec $D$ la partie diagonale de $A$, $-E$ sa partie triangulaire inférieure stricte et $-F$ sa partie triangulaire supérieure stricte. Autrement dit\n",
    "$$\n",
    "\\begin{aligned}\n",
    "D_{ii}&=A_{ii} && i=1,\\ldots,n,\\\\\n",
    "E_{ij}&=-A_{ij} && 1\\leq j<i\\leq n,\\\\\n",
    "F_{ij}&=-A_{ij} && 1\\leq i<j\\leq n,\n",
    "\\end{aligned}\n",
    "$$\n",
    "tous les autres termes étant nuls. Pour la construction de ces trois parties utilisera la commande ```np.diag```.\n",
    "\n",
    "\n",
    "- Méthode de Jacobi : $M=D$ et $N=E+F$.\n",
    "- Méthode de Gauss-Seidel : $M=D-E$ et $N=F$.\n",
    "- Méthode de relaxation : $M=\\frac 1\\omega D-E$ et $N=\\frac{1-\\omega}\\omega D+F$ avec $\\omega\\in\\mathbb R^+_{*}$.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "722f89de",
   "metadata": {},
   "source": [
    "Soit la matrice carrée $A$ de taille $n$, tridiagonale, intervenant dans le TP précédent, définie par\n",
    "$$\n",
    "  A = \n",
    "  \\begin{bmatrix}\n",
    "    2 & -1 &  0 & \\ldots & 0 \\\\\n",
    "    -1 & 2 & -1 &  \\ddots & \\vdots  \\\\\n",
    "    0  & \\ddots & \\ddots & \\ddots & 0 \\\\\n",
    "    \\vdots & \\ddots & -1 & 2 & -1 \\\\\n",
    "    0 & \\ldots & 0 & -1 & 2\n",
    "  \\end{bmatrix}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dec49ab5",
   "metadata": {},
   "source": [
    "### Question 1)\n",
    "\n",
    "Programmer les trois méthodes précedentes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e550066",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "21b3cce9",
   "metadata": {},
   "source": [
    "### Question 2)\n",
    "\n",
    "\n",
    "Application : utiliser ces trois méthodes sur la matrice $A$ et le second\n",
    "  membre $B = (1,0,\\ldots,0,1)^T$, pour\n",
    "  $n=10$. On effectuera 100 itérations. Dans le cas de la méthode de relaxation,\n",
    "  on choisira $\\omega = \\tfrac{3}{2}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "923304c3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8da413c0",
   "metadata": {},
   "source": [
    "### Question 3)\n",
    "\n",
    "\n",
    "Dans le cas où $n=20$, déterminer le paramètre optimal dans la méthode de \n",
    "  relaxation, en représentant le rayon spectral de la matrice d'itération obtenue en\n",
    "  fonction de $\\omega$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9797b0d9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "96ca742a",
   "metadata": {},
   "source": [
    "\n",
    "### Question 4)\n",
    "\n",
    " Pour différentes valeurs de $n$, comparer le nombre d'itérations\n",
    "  nécessaires pour avoir une précision de $10^{-12}$ entre la solution\n",
    "  approchée et la solution exacte. Quelle relation y a t-il entre le nombre\n",
    "  d'itérations et le rayon spectral ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47dc2e43",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "49d0b17d",
   "metadata": {},
   "source": [
    " ### Question 5) \n",
    "Vérifier que la méthode de Jacobi appliquée à la matrice \n",
    "$$\n",
    "    C =\n",
    "    \\begin{bmatrix}\n",
    "      1 & 2 & -2 \\\\\n",
    "      1 & 1 & 1 \\\\\n",
    "      2 & 2 & 1 \\\\\n",
    "     \\end{bmatrix}\n",
    "$$\n",
    "converge, mais que la méthode de Gauss-Seidel diverge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4588c7a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4009c960",
   "metadata": {},
   "source": [
    "### Question 6)\n",
    "Vérifier que la méthode de Gauss-Seidel appliquée à la matrice \n",
    "$$\n",
    "    D =\n",
    "    \\begin{bmatrix}\n",
    "      2 & -1 & 1 \\\\\n",
    "      2 & 2 & 2 \\\\\n",
    "      -1 & -1 & 2 \\\\\n",
    "     \\end{bmatrix}\n",
    "$$\n",
    "  converge, mais que la méthode de Jacobi diverge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eba2d4fb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f4536108",
   "metadata": {},
   "source": [
    "## Exercice 4. Résolution d'un problème de Poisson 2D\n",
    "\n",
    "On souhaite résoudre le problème de Poisson en dimension 2 d'espace\n",
    "$$\n",
    "\\dfrac{\\partial^2 u}{\\partial x^2} + \\dfrac{\\partial^2 u}{\\partial y^2} = f(x,y), \\quad (x,y)\\in[0,1]^2,\n",
    "$$\n",
    "assorti des conditions de Dirichlet au bord \n",
    "$$ u(x,y)=0,\\quad (x,y)\\in\\partial [0,1]^2.$$\n",
    "\n",
    "Pour ce faire, on approche l'équation par différences finies. Soit $N$ un entier naturel, on pose $h=1/(N+1)$ le pas de discrétisation pour chacune des directions. Pour tout $0\\leq i,j\\leq N+1$, $u_{i,j}$ désigne une approximation de la valeur $u(x_i,y_j)$, les points de la grille étant $x_i= i h$, $y_j = jh$.\n",
    "L'EDP est approchée par des différences centrées :\n",
    "$$\n",
    "\\dfrac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2} + \\dfrac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2} = f(x_i,y_j), \\quad 1\\leq i,j \\leq N,\n",
    "$$\n",
    "et les conditions de bord sont données: $u_{i,j}=0$ lorsque $ij(N+1-i)(N+1-j)=0$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3db500de",
   "metadata": {},
   "source": [
    "### Question 1)\n",
    "\n",
    " Construire la matrice $A$ de taille $N^2\\times N^2$ permettant de définir le vecteur des inconnues $U=(u_{i,j})_{1\\leq i,j\\leq N}$ de $\\mathbb R^{N\\times N}$ comme solution du problème\n",
    "$$AU = F.$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e99ac998",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "46b98ece",
   "metadata": {},
   "source": [
    "### Question 2) \n",
    "En utilisant les méthodes itératives précédemment mises en œuvre, résoudre ce problème pour le second membre $f$ donné dans l'exemple de code ci-dessous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1557fcfa",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4f5d792e",
   "metadata": {},
   "source": [
    "\n",
    "### Question 3)\n",
    "Représenter la solution. On pourra adapter l'exemple suivant pour tracer une surface en 3 dimensions :\n",
    "```python\n",
    "from mpl_toolkits.mplot3d import axes3d\n",
    "from matplotlib import cm\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection = '3d')\n",
    "n = 20\n",
    "X = np.linspace(0,1,n+2)\n",
    "Y = X\n",
    "X, Y = np.meshgrid(X, Y)\n",
    "def f(x,y):\n",
    "    return -(x-0.5)**2 - (y-0.5)**2\n",
    "Z = f(X,Y)\n",
    "surf = ax.plot_surface(X, Y, Z,cmap = cm.coolwarm,linewidth = 0, \n",
    "    \t\tantialiased = False)\n",
    "plt.title(\"La surface f(x,y)=-(x-0.5)**2-(y-0.5)**2\")\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "où ```X``` et ```Y``` sont les vecteurs avec les coordonnées des points du plan, et ```Z``` est la matrice contenant les valeurs ```f(X,Y)```.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0043957",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6.5 64-bit",
   "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.6.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}