{ "cells": [ { "cell_type": "markdown", "id": "d8c29c40", "metadata": {}, "source": [ "# L3 Analyse Numérique – TP2\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 la résolution de systèmes linéaires.\n", "\n", "- Exercice 1. *Conditionnement* d'une matrice de taille fixée dépendant d'un paramètre.\n", "- Exercice 2. Résolution d'un problème d'élasticité par une méthode de différences finies.\n", "- Exercice 3. Approche optimisée par stockage creux 'sparse'" ] }, { "cell_type": "markdown", "id": "1b118759", "metadata": {}, "source": [ "## Exercice 1. Conditionnement d'une matrice" ] }, { "cell_type": "markdown", "id": "74bb2a34", "metadata": {}, "source": [ "Nous importons avant tout quelques librairies qui seront utiles." ] }, { "cell_type": "code", "execution_count": 1, "id": "185b8caf", "metadata": {}, "outputs": [], "source": [ "from math import *\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "be1602bd", "metadata": {}, "outputs": [], "source": [ "# Au besoin (?)\n", "# import sys\n", "# !{sys.executable} -m pip install scipy" ] }, { "cell_type": "markdown", "id": "c1b67926", "metadata": {}, "source": [ "Étant donnée une matrice inversible $M\\in\\mathsf{GL}_n(\\mathbb{R})$ et une donnée $b\\in\\mathbb{R}^n$, la résolution du système linéaire $Mx=b$ peut s'avérer particulièrement sensible aux imprécisions sur la donnée $b$, *ceci même pour des calculs menés en arithmétique exacte*. Un résultat du cours concerne l'erreur relative commise sur la solution $x$ qui peut être contrôlée à l'aide du conditionnement de $M$ défini comme étant la quantité: \n", "\n", "$$\\mathsf{cond}\\,M:= \\|M\\|\\,\\|M^{-1}\\|.$$\n", "\n", "En l'occurence, si $Mx=b$ et $My=b+\\delta b$ avec $x\\neq 0$, alors l'inégalité suivante est satisfaite: \n", "\n", "$$\n", "\\dfrac{\\|y-x\\|}{\\|x\\|}\\leq \\mathsf{cond}\\, M \\dfrac{\\|\\delta b\\|}{\\|b\\|}.\n", "\\tag{Err relat.}\n", "$$\n", "\n", "\n", "En Python, le conditionnement s'obtient à l'aide de la commande [`numpy.linalg.cond`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.cond.html).\n", "\n", "Dans cet exercice, on s'intéressera à la matrice $M_\\alpha$ suivante, dépendant d'un paramètre réel $\\alpha$:\n", "$$\n", "M_\\alpha=\n", "\\begin{pmatrix}\n", " 1+\\alpha&1&2\\\\\n", " 3&1&\\alpha\\\\\n", " 1&2\\alpha&1\n", "\\end{pmatrix}.\n", "$$\n", "**Attention.** Cette matrice n'est pas inversible pour toutes les valeurs du paramètre $\\alpha$. \n", "L'espace vectoriel $\\mathbb{R}^3$ est muni de la norme euclidienne usuelle notée $\\|\\cdot\\|$, qui s'obtient en Python par la fonction [`numpy.linalg.norm`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html#numpy.linalg.norm). Les erreurs seront calculés dans cette norme ainsi que le conditionnement de la matrice." ] }, { "cell_type": "markdown", "id": "67e411c8", "metadata": {}, "source": [ "### Question 1)\n", "\n", "Programmer une fonction Python `def Malpha(alpha):` qui renvoie un `numpy.array` de `float` représentant la matrice $M_\\alpha$. \n", "Calculer ensuite $M_2$ et vérifier qu'elle n'est pas inversible." ] }, { "cell_type": "code", "execution_count": null, "id": "5db160f6", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "ff987ec7", "metadata": {}, "source": [ "### Question 2)\n", "\n", "Soit $\\alpha=2.01$ pour lequel la matrice $M_\\alpha$ est inversible. \n", "Étant donné $x=(1,1,1)^T$, on considère le second membre $b=M_\\alpha x$.\n", "\n", "**a)** Calculer $y\\in\\mathbb{R}^3$ la solution du système linéaire $M_\\alpha y=b+\\delta b$ où $\\delta b\\in\\mathbb{R}^3$ est une perturbation vectorielle obtenue par une réalisation d'un aléa de petite amplitude `0.01*np.random.rand(3,1)`.\n", "\n", "**b)** Calculer alors numériquement l'erreur relative sur la solution\n", "$$\n", "\\dfrac{\\|y-x\\|}{\\|x\\|},\n", "$$\n", "puis vérifier numériquement que l'inégalité (Err relat.) est bien réalisée." ] }, { "cell_type": "code", "execution_count": null, "id": "218d273f", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "1dc4f375", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "8100c222", "metadata": {}, "source": [ "### Question 3) \n", "\n", "Reprendre la question précédente de façon à traiter un grand nombre de valeurs de $\\alpha\\in[-5,5]$ (choisies aléatoirement) et des données $x$ et $\\delta b$ aléatoires pour chaque $\\alpha$. L'objectif est de visualiser graphiquement l'inégalité (Err relat.). On pourra considérer la quantité: \n", "$$\n", "\\tau = \\left(\\dfrac{\\|y-x\\|}{\\|x\\|}\\right)\\left(\\dfrac{\\|\\delta b\\|}{\\|b\\|}\\right)^{-1},\n", "$$ \n", "et représenter séparément en fonction de $\\alpha$ les quantités $\\tau$, $\\det(M_\\alpha)$ et $\\dfrac{\\tau}{\\mathsf{cond}\\, M_\\alpha}$.\n", "\n", "Le code suivant peut vous servir de base.\n", "```python\n", "n = 50\n", "\n", "fig, (ax1, ax2, ax3) = plt.subplots(3,figsize=(10,10))\n", "Lalpha = np.zeros(n)\n", "Lx, Ly, Lz = np.zeros(n), np.zeros(n), np.zeros(n)\n", "\n", "for i in range(n):\n", " alpha = np.random.rand() * 10 - 5\n", " \n", " x = alpha\n", " y = alpha**2\n", " z = alpha**3\n", "\n", " Lalpha[i] = alpha\n", " Lx[i] = x\n", " Ly[i] = y\n", " Lz[i] = z\n", " \n", "ax1.plot(Lalpha,Lx,'.')\n", "ax2.plot(Lalpha,Ly,'.')\n", "ax3.plot(Lalpha,Lz,'.')\n", " \n", "plt.show()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "9ce3fa1d", "metadata": { "scrolled": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "b9565898", "metadata": {}, "source": [ "## Exercice 2. Problème d'élasticité\n", "\n", "On souhaite déterminer une approximation de la solution $u\\in\\mathcal{C}^2(]0,1[)\\cap\\mathcal{C}^0([0,1])$ du **problème aux limites** suivant: \n", "\n", "$$\n", "\\begin{cases}\n", "\\dfrac{d}{dx}\\left(a(x)\\dfrac{du}{dx}\\right) = f(x),& 00$. La corde est fixée aux deux points extrémaux de coordonnées $(0,0)$ et $(1,0)$ et est soumise à la charge linéique $f$ qui en modifie la géométrie par l'effet de la gravité." ] }, { "cell_type": "markdown", "id": "7a2421b0", "metadata": {}, "source": [ "Nous décrivons une stratégie de résolution numérique approchée qui s'appuie sur **une approximation en dimension finie du problème**. \n", "\n", "Pour ce faire, considérons un entier naturel non nul $n\\in\\mathbb{N}$ à partir duquel est définie une discrétisation de l'intervalle d'espace $x\\in[0,1]$: notons $h=(n+1)^{-1}$ le pas de discrétisation et pour tout entier $i$, posons $x_i=ih$ et $x_{i+1/2}=(i+1/2)h$. \n", "\n", "L'inconnue du problème à résoudre est alors le vecteur $U=(u_1,\\ldots,u_n)^T\\in\\mathbb{R}^n$, solution des équations aux différences finies suivantes: \n", "\n", "$$\n", "\\begin{aligned}\n", "&\\dfrac{1}{h}\\left(a(x_{i+1/2})\\dfrac{u_{i+1}-u_{i}}{h}-a(x_{i-1/2})\\dfrac{u_{i}-u_{i-1}}{h}\\right) = f(x_i),\\quad 1\\leq i \\leq n,\\\\\n", "&u_0 = u_{n+1} = 0.\n", "\\end{aligned}\n", "\\tag{Pb}\n", "$$\n", "\n", "On admettra que la solution $U$ de ce problème discret approche bien celle du problème continu $u$ au sens par exemple où la quantité \n", "\n", "$$\\max_{1\\leq i \\leq n} \\vert u_i-u(x_i)\\vert$$\n", "\n", "converge vers $0$ lorsque la dimension $n$ tend vers l'infini." ] }, { "cell_type": "markdown", "id": "f7e74f13", "metadata": {}, "source": [ "### Question 1)\n", "\n", "Identifier *sur le papier* avec soin une matrice $A\\in\\mathsf{M}_n(\\mathbb{R})$ et un vecteur $F\\in\\mathbb{R}^n$ permettant de réécrire le problème (Pb) d'inconnue $U\\in\\mathbb{R}^n$ sous la forme d'un système linéaire $AU=F$. " ] }, { "cell_type": "markdown", "id": "a0390bad", "metadata": {}, "source": [ "### Question 2)\n", "\n", "Programmer la construction de la matrice $A$ et du second membre $F$, à travers une fonction dont la syntaxe sera la suivante, à compléter:\n", "```python\n", "def construction(n,a,f):\n", "\t# n entier, a et f le nom des fonctions concernees\n", "\th = 1./(n+1)\n", "\tx = np.arange(0,n+2)*h\n", " vec = a(x[:-1]+0.5*h)\n", " ... # Astuce: utiliser l'instruction np.diag vue au TP1\n", " return A,F\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "c52dff3f", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "c89586e3", "metadata": {}, "source": [ "### Question 3)\n", "\n", "Vérifier que, dans le cas des fonctions $a$ et $f$ constantes égales à 1:\n", "$$\n", "\\begin{aligned}\n", "a(x) &= 1,\\\\\n", "f(x) &= 1, \n", "\\end{aligned}\n", "$$\n", "et pour une valeur de $n$ à votre convenance, le vecteur suivant\n", "$$\n", "U=\\left(\\dfrac{x_i(x_i-1)}{2}\\right)_{1\\leq i \\leq n}\n", "$$\n", "est bien solution du problème (Pb) à l'erreur machine près. \n", "On pourra faire usage de la commande `np.linalg.solve` documentée dans l'aide." ] }, { "cell_type": "code", "execution_count": null, "id": "64cc3107", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "bac4cada", "metadata": {}, "source": [ "### Question 4)\n", "\n", "Désormais, les fonctions $a$ et $f$ sont choisies selon les paramètres suivants: \n", "\n", "$$\n", "\\begin{aligned}\n", "a(x) &= 1 -0.95\\mathrm{e}^{-500(x-0.25)^2},\\\\\n", "f(x) &= 0.1 + 0.75\\mathrm{e}^{-1000(x-\\xi_1)^2} + 2\\mathrm{e}^{-1000(x-\\xi_2)^2}.\n", "\\end{aligned}\n", "$$\n", "\n", "Dans la fonction $f$ ainsi définie, le terme constant peut s'interpréter comme la contribution de la masse propre de la corde, tandis que les deux termes exponentiellement localisés modélisent deux masses suspendues, placées au voisinage des points d'abscisses respectives $\\xi_1=0.2$ et $\\xi_2=0.8$.\n", "\n", "```python\n", "def a(x):\n", " return 1 - 0.95*np.exp(-500*(x-0.25)**2)\n", "\n", "xi1 = 0.2\n", "xi2 = 0.8\n", "def f(x):\n", " return 0.1 + 0.75*np.exp(-1000*(x-xi1)**2) + 2*np.exp(-1000*(x-xi2)**2)\n", "\n", "fig, ax1 = plt.subplots()\n", "ax2 = ax1.twinx()\n", "x = np.linspace(0, 1, 200)\n", "ax1.plot(x, a(x), 'g-')\n", "ax2.plot(x, f(x), 'b-')\n", "\n", "ax1.set_xlabel('position')\n", "ax1.set_ylabel('elasticité', color='g')\n", "ax2.set_ylabel('chargement', color='b')\n", "plt.show()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "28015694", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "386f4268", "metadata": {}, "source": [ "### Question 4a)\n", "\n", "Calculer la matrice $A$ qui intervient pour une discrétisation de $n=5$ points intérieurs et vérifier numériquement que $A$ est alors bien symétrique, définie et négative. Si ce n'est pas le cas, c'est qu'il y a une erreur dans votre fonction `construction` à la question 2 !" ] }, { "cell_type": "code", "execution_count": null, "id": "0b3f1137", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "11d26cd6", "metadata": { "scrolled": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "239ce5de", "metadata": {}, "source": [ "### Question 4b)\n", "\n", "Résoudre le problème (Pb) avec les nouvelles données $a$ et $f$ et représenter la corde pour différentes valeurs de $n$. Si possible, on ajoutera dans la représentation graphique les conditions de bord $u(0)=0=u(1)$ qui ont été initialement écartées du vecteur $U$ (utiliser pour cela par exemple `np.concatenate`)." ] }, { "cell_type": "code", "execution_count": null, "id": "5c0313f7", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "e3636306", "metadata": { "scrolled": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "e470465e", "metadata": {}, "source": [ "## Exercice 3\n", "\n", "Cet exercice fait suite au précédent.\n", "\n", "Pour de grandes valeurs de $n$, le temps de calcul nécessaire pour construire la matrice $A$ puis pour résoudre le système linéaire $AU=F$ devient important. Nous allons dans la suite utiliser une stratégie de résolution adaptée à la structure particulière de la matrice $A$. De nombreuses stratégies sont possibles mais un gain d'efficacité conséquent s'obtient en mettant à profit la **structure creuse** de $A$, c'est à dire en adaptant le stockage des variables et les algorithmes de résolution de façon à ne pas tenir compte des très nombreux éléments nuls de $A$, ici localisés \"loin\" de la diagonale. La librairie `scipy` fournit différents stockages creux et des algorithmes d'algèbre linéaires adaptés. " ] }, { "cell_type": "markdown", "id": "75fab052", "metadata": {}, "source": [ "### Question 1)\n", "\n", "Tester les trois cellules de code suivantes et commentez." ] }, { "cell_type": "code", "execution_count": null, "id": "9158e3da", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0. 0. 0. 0. 0.]\n", " [0. 0. 1. 0. 0. 0.]\n", " [0. 0. 0. 2. 0. 0.]\n", " [0. 0. 0. 0. 3. 0.]\n", " [0. 0. 0. 0. 0. 4.]\n", " [0. 0. 0. 0. 0. 0.]]\n", " (1, 2)\t1.0\n", " (2, 3)\t2.0\n", " (3, 4)\t3.0\n", " (4, 5)\t4.0\n" ] } ], "source": [ "import numpy as np\n", "import scipy.sparse as sp\n", "from scipy.sparse.linalg import spsolve\n", "\n", "print(np.diag(np.arange(5.),1))\n", "print(sp.diags(np.arange(5),1))" ] }, { "cell_type": "code", "execution_count": null, "id": "8541b606", "metadata": { "scrolled": true }, "outputs": [], "source": [ "%%time\n", "n = 2000\n", "A = np.eye(n,n)\n", "x = np.ones(n)\n", "y = np.linalg.solve(A,x)" ] }, { "cell_type": "code", "execution_count": null, "id": "b46c2052", "metadata": {}, "outputs": [], "source": [ "%%time\n", "n = 200000\n", "A = sp.eye(n,n,format='csr')\n", "x = np.ones((n,1))\n", "y = spsolve(A,x)" ] }, { "cell_type": "markdown", "id": "b76ff0ae", "metadata": {}, "source": [ "Dans la suite, nous utiliserons systématiquement ce stockage creux. La fonction `construction` a été réécrite de façon adaptée en une nouvelle fonction `constructionSP` définie par:\n", "```python\n", "def constructionSP(n,a,f):\n", " h = 1./(n+1)\n", " x = np.arange(0,n+2)*h\n", " vec = a((0.5+np.arange(n+1))*h)\n", " A = (sp.diags(vec[1:-1],1,format='csr')+sp.diags(vec[1:-1],-1,format='csr')-sp.diags(vec[:-1]+vec[1:],format='csr'))/h**2\n", " F = f(x[1:-1])\n", " return A,F\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "aa7a4eb3", "metadata": { "scrolled": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "b617651c", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "2df25567", "metadata": {}, "source": [ "### Question 2)\n", "\n", "Revenons à notre problème d'élasticité (Pb). La seconde charge est supposée fixée à la position $\\xi_2=0.8$. On souhaite maintenant optimiser la position $\\xi_1$ de la première masse de façon à maximiser l'amplitude de la corde, autrement dit, à minimiser la quantité négative $\\min_{x\\in[0,1]} u(x)$. \n", "Proposer et mettre en œuvre un procédé numérique de votre choix qui aborde ce problème. \n", "On pourra tester et décrypter au préalable la commande suivante:\n", "```python\n", "test = (np.linspace(0,10,11)-5.)**2\n", "print(test)\n", "np.where(test == np.min(test))[0][0]\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "70698dac", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "dae91bb5", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "7c847d86", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "3c051ac5", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "ee61aa52", "metadata": {}, "source": [ "### Question 3)\n", "\n", "Reprendre la question précédente, avec cette fois-ci les deux positions $\\xi_1\\in[0,1]$ et $\\xi_2\\in[0,1]$ à optimiser. \n", "\n", "On pourra jeter un œil aux exemples d'utilisation de la fonction [`matplotlib.contour`](https://matplotlib.org/stable/gallery/images_contours_and_fields/contour_demo.html) et utiliser la commande `numpy.meshgrid`." ] }, { "cell_type": "code", "execution_count": null, "id": "dfcf4238", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "9cc49d61", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "7c4c6c57", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "78ca9923", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "349e6e06", "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.9.2 (v3.9.2:1a79785e3e, Feb 19 2021, 09:09:00) \n[Clang 12.0.0 (clang-1200.0.32.29)]" }, "vscode": { "interpreter": { "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49" } } }, "nbformat": 4, "nbformat_minor": 5 }