{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "8ea8318a", "metadata": {}, "source": [ "# L3 Analyse Numérique – TP4\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", "Dans ce TP, nous étudierons la méthode des moindres carrés qui permet d'“approcher” une solution de $Ax= b$ lorsque $A$ n'est pas inversible puis nous étudierons les méthodes de gradient.\n", "\n", "**IMPORTANT :** L'Exercice 2 est à rendre par mail (à Pierre Le Barbenchon ET Antoine Dequay) avant le dimanche 19 mars à 23h59 en format .ipynb afin que nous puissions exécuter les cellules de votre notebook, vous pouvez nous envoyer tout le TP, nous ne regarderons que l'Exercice 2. Pour les questions dont les réponses ne contiennent pas de code, vous avez le choix entre joindre une photo/un pdf de votre réponse, ou écrire directement en $\\LaTeX$ à l'intérieur d'un bloc de texte (tutoriel disponible [ici](https://jupyter-notebook.readthedocs.io/en/latest/examples/Notebook/Working%20With%20Markdown%20Cells.html?highlight=latex#LaTeX-equations) et [ici](https://www.youtube.com/watch?v=_slr5F9dfZ4&list=PLaidQL63zBcKzwKesj3Ef_mSqZnFSDAj8)).\n", "\n", "\n", "- Exercice 1. Moindres carrés.\n", "- Exercice 2. Minimum d'une fonctionnelle quadratique strictement convexe." ] }, { "cell_type": "markdown", "id": "74bb2a34", "metadata": {}, "source": [ "Nous importons avant tout quelques librairies qui seront utiles." ] }, { "cell_type": "code", "execution_count": 1, "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. Moindres carrés\n", "\n", "On souhaite déterminer la linéarisation de l'expression trigonométrique \n", "$\n", "\\cos^5x\\,(1+\\sin^3x),% = \\dfrac{1}{128}\\Big(80\\cos(x)+40\\cos(3x)+8\\cos(5x)+6\\sin(2x)+2\\sin(4x)-2\\sin(6x)-\\sin(8x)\\Big).\n", "$\n", "c'est à dire les coefficients rationnels $(a_j)_{0\\leq j\\leq d}$ et $(b_j)_{1\\leq j \\leq d}$ tels que \n", "$$\n", "\\cos^5x\\,(1+\\sin^3x) = a_0 + \\sum_{j=1}^d a_j \\cos (jx) + b_j \\sin (jx)\n", "$$\n" ] }, { "cell_type": "markdown", "id": "5880772d", "metadata": {}, "source": [ "### Question 1)\n", "\n", "Considérer un entier $N$ suffisamment grand et constituer le vecteur colonne $F$ de taille $N$ comme suit:\n", "\n", "```python\n", "X = 2*pi/N*np.arange(N)\n", "F = np.cos(X)**5*(1+np.sin(X)**3)\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "d4f90b69", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "4ce2a7ab", "metadata": {}, "source": [ "On doit trouver le résultat suivant :\n", "$$\n", "\\begin{aligned}\n", "\\cos^5(x)(1+\\sin^3(x)) &= \\dfrac{1}{128}\\Bigl[80\\cos(x) + 40\\cos(3x) + 8\\cos(5x) + 6\\sin(2x) + 2\\sin(4x) -2\\sin(6x) -\\sin(8x)\\Bigr]\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "markdown", "id": "b578ae67", "metadata": {}, "source": [ "### Question 2)\n", "\n", "Définir une fonction ```base(d,X)``` prenant en argument un entier ```d``` supérieur ou égal à 1, et un vecteur colonne ```X``` de taille notée $N$, et renvoyant la matrice suivante\n", "$$\n", "\\texttt{A} = \\begin{pmatrix}\n", "1 & \\cos \\texttt{X}_1 & \\cdots &\\cos \\texttt{d}\\texttt{X}_1 & \\sin \\texttt{X}_1 & \\cdots &\\sin \\texttt{d}\\texttt{X}_1\\\\\n", "1 & \\cos \\texttt{X}_2 & \\cdots &\\cos \\texttt{d}\\texttt{X}_2 & \\sin \\texttt{X}_2 & \\cdots &\\sin \\texttt{d}\\texttt{X}_2\\\\\n", "\\vdots & \\vdots & & \\vdots &\\vdots & & \\vdots \\\\\n", "1 & \\cos \\texttt{X}_\\texttt{N} & \\cdots &\\cos \\texttt{d}\\texttt{X}_\\texttt{N} & \\sin \\texttt{X}_\\texttt{N} & \\cdots &\\sin \\texttt{d}\\texttt{X}_\\texttt{N}\\\\\n", "\\end{pmatrix}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "88d2c720", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "5a73912f", "metadata": {}, "source": [ "### Question 3) \n", "\n", "Pour des valeurs de ```d``` de plus en plus grandes, déterminer le vecteur ```C``` de taille ```2d+1```, solution au sens des moindres carrés de l'équation\n", "$\n", "\\texttt{A}\\, \\texttt{C} = \\texttt{F}.\n", "$ On utilisera la commande ```np.linalg.lstsq(a,b,rcond=None)[0]``` qui renvoie la solution $x$ de $ax=b$ au sens des moindres carrés." ] }, { "cell_type": "code", "execution_count": null, "id": "1238ac49", "metadata": { "scrolled": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "d55c4800", "metadata": {}, "source": [ "### Question 4)\n", "\n", "On considère désormais une valeur ```d``` pour laquelle la norme du résidu $\\|\\texttt{A}\\, \\texttt{C} - \\texttt{F}\\|$ est jugée suffisamment petite. Déterminer une expression rationnelle des coefficients de la solution retenue ```C``` et la linéarisation recherchée." ] }, { "cell_type": "code", "execution_count": null, "id": "6c2cbb96", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "ac32cd23", "metadata": {}, "source": [ "### Question 5)\n", "\n", "Déterminer numériquement la matrice symétrique définie positive intervenant dans l'équation normale." ] }, { "cell_type": "code", "execution_count": null, "id": "a0278d82", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "682807bf", "metadata": {}, "source": [ "### Question 6)\n", "\n", "Linéariser également l'expression $(1+\\cos^3(2x) + \\sin^3(2x))\\cos^{10}(x)$." ] }, { "cell_type": "code", "execution_count": null, "id": "59ef7d01", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "f21decb8", "metadata": {}, "source": [ "On doit trouver le résultat suivant :\n", "$$\n", "\\begin{aligned}\n", "(1+\\cos^3(2x) + \\sin^3(2x))\\cos^{10}(x) &= \\dfrac{1}{4096}\n", "\\Bigl[ 1683 + 2926\\cos(2x) + 1936\\cos(4x) + 1002\\cos(6x) + 428\\cos(8 \n", " x)\\\\ &+ 158\\cos(10x) + 48\\cos(12x) + 10\\cos(14x) + \\cos(16x) + 286\\sin(2x) \\\\ &+ 286\\sin(4x) + 78\\sin( \n", " 6x) -78\\sin(8x) -90\\sin(10x) -42\\sin(12x) \\\\\n", " & -10\\sin(14x) -\\sin(16x)\\Bigr]\n", " \\end{aligned}$$" ] }, { "cell_type": "markdown", "id": "9f5612fa", "metadata": {}, "source": [ "## Exercice 2. Minimum d'une fonctionnelle quadratique strictement convexe\n", "\n", "\n", "Soit $n$ un entier naturel non-nul et $J:\\mathbb R^n \\rightarrow \\mathbb R$ une fonctionnelle prenant la forme\n", "\\begin{equation}\n", " J(x) = \\frac{1}{2} \\langle Ax,x\\rangle - \\langle b,x\\rangle,\n", "\\end{equation}\n", "où $b\\in\\mathbb R^n$ et $A\\in S_n^{++}(\\mathbb R)$. Dans ce cadre, il est connu (cf. cours) que la fonctionnelle $J$ admet un minimum global sur $\\mathbb R^n$, atteint en un unique point $x\\in\\mathbb R^n$ solution de l'équation $\\nabla J(x)= 0$ où $\\nabla J(x) := Ax-b$.\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "dec49ab5", "metadata": {}, "source": [ "### Question 1)\n", "\n", "On considère pour premier exemple les données suivantes:\n", "$$\n", "b=\\begin{pmatrix}0\\\\0\\end{pmatrix},\\qquad A=\\begin{pmatrix}1 & 0\\\\0 & \\tfrac{1}{20}\\end{pmatrix}.\n", "$$\n", "\n", "Ecrire une fonction ```gradient_pas_fixe(A,b,x0,rho,itmax,eps)``` qui renvoie les itérés $(x_k)_{k\\geq 0}$ obtenus par la méthode du gradient à pas fixe, partant de l'initialisation $x_0={}^t(1,1)\\in\\mathbb R^2$ en utilisant les deux conditions :\n", "- ```itmax``` le nombre d'itérations maximal.\n", "- ```eps``` epsilon l'erreur maximal qu'on s'autorise à la fin de l'algorithme." ] }, { "cell_type": "markdown", "id": "1cd8be50", "metadata": {}, "source": [ "Méthode du gradient à pas fixe $\\rho>0$\n", "$$\n", "\\begin{cases}\n", "x_0\\in\\mathbb{R}^n\\\\\n", "x_{k+1}=x_k-\\rho \\nabla J(x_k).%,\\quad k\\geq 0.\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "1e550066", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "8ceb79c2", "metadata": {}, "source": [ "### Question 2)\n", "Représenter les itérations successives dans le plan $\\mathbb R^2$ pour ```itmax = 100``` et ```eps = 0.01```. Pour le pas $\\rho$, on choisira par exemple successivement les valeurs parmi $\\{0.10,1,1.98\\}$. Sur une autre figure on tracera également l'évolution du résidu $\\|Ax_k-b\\|_2$ en fonction de $k$ pour chacun des cas." ] }, { "cell_type": "code", "execution_count": null, "id": "1ca732e9", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "21b3cce9", "metadata": {}, "source": [ "### Question 3)\n", "\n", "Dans cette question et les suivantes, on se place en dimension $n=20$ avec les données suivantes:\n", "$$\n", "b = \\begin{pmatrix}1 \\\\ 0 \\\\ \\vdots \\\\ 0\\\\ 1\\end{pmatrix}\n", ",\\qquad \n", " A = \n", " (n+1)^2\\begin{pmatrix}\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{pmatrix}.\n", "$$\n", "Tester de nouveau la convergence de l'algorithme du gradient à pas fixe en traçant le résidu pour quelques valeurs du pas $\\rho$. Vérifier numériquement que la valeur du pas pour laquelle la vitesse de convergence semble maximale est\n", "$$\\rho_c = \\tfrac{2}{\\lambda_{\\min} + \\lambda_{\\max}}, \\textrm{ avec }\\lambda_{\\min}=\\min(\\mathrm{spec}(A)) \\textrm{ et }\\lambda_{\\max}=\\max(\\mathrm{spec}(A)).$$\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "923304c3", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "8da413c0", "metadata": {}, "source": [ "### Question 4)\n", "\n", "Retenir le choix $\\rho=\\rho_c$ et vérifier numériquement que la convergence est alors au plus géométrique, de raison $r \\overset{def}{=}\\tfrac{\\kappa -1}{\\kappa + 1}$, où $\\kappa$ est le conditionnement de $A$ en norme 2, à savoir ici $\\kappa = \\mathrm{cond}_2(A)=\\lambda_{\\max}/\\lambda_{\\min}$, autrement dit vérifier numériquement que \n", "$$\\|Ax_k-b\\| \\leqslant C r^k.$$\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9797b0d9", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "7c708769", "metadata": {}, "source": [ "### Question 5)" ] }, { "cell_type": "markdown", "id": "11159f2a", "metadata": {}, "source": [ "Méthode du gradient à pas optimal\n", "$$\n", "\\begin{cases}\n", "x_0\\in\\mathbb{R}^n\\\\\n", "\\rho_k=\\mathrm{argmin}_{\\rho\\in\\mathbb R}J(x_k-\\rho\\nabla J(x_k)),\\\\%\\quad k\\geq 0\\\\\n", "x_{k+1}=x_k-\\rho_k\\nabla J(x_k).%,\\quad k\\geq 0.\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "id": "96ca742a", "metadata": {}, "source": [ "\n", "\n", "Dans la méthode du gradient à pas optimal, démontrer que le pas optimal $\\rho_k$, solution d'un problème de minimisation monodimensionnel, est:\n", "$$\n", "\\rho_k = \\dfrac{\\langle R_k,R_k\\rangle}{\\langle A R_k,R_k\\rangle},\\quad \\textrm{ où } R_k=Ax_k-b.\n", "$$" ] }, { "cell_type": "markdown", "id": "47dc2e43", "metadata": {}, "source": [ "Réponse à rédiger ici" ] }, { "cell_type": "markdown", "id": "49d0b17d", "metadata": {}, "source": [ " ### Question 6) \n", "\n", "Programmer la méthode du gradient à pas optimal pour le problème de la Question 3) avec de nouveau les deux conditions d'arrêt.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b4588c7a", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "4009c960", "metadata": {}, "source": [ "### Question 7)\n", "\n", "Vérifier que l'algorithme converge au plus géométriquement, avec une raison égale à $r \\overset{def}{=}\\tfrac{\\kappa -1}{\\kappa + 1}$, autrement dit que\n", "$$\\|Ax_k - b\\| \\leqslant C r^k.$$" ] }, { "cell_type": "code", "execution_count": null, "id": "eba2d4fb", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "3db500de", "metadata": {}, "source": [ "### Question 8)\n", "\n", "Programmer la méthode du gradient conjugué pour le problème de la Question 3) avec de nouveau les deux conditions d'arrêt." ] }, { "cell_type": "markdown", "id": "f8f1547a", "metadata": {}, "source": [ "**Méthode du gradient conjugué pour une fonctionnelle quadratique**\n", "$$\n", "\\begin{aligned}\n", "\\textit{Initialisation :}~~~& x_0\\in\\mathbb R^n\\\\% \\textit{donn\\'e}, $\\varepsilon$\\textit{ donn\\'e}\\\\\n", "&r_0=b-A\\cdot x_0 &\\textit{(r\\'esidu initial)}\\\\\n", "&p_0=r_0 &\\textit{(direction de descente initiale)}\\\\\n", "&\\theta_0=\\langle p_0,r_0\\rangle \\\\\n", "\\\\\n", "\\textit{It\\'erations : }k\\ge 0~~~\n", "&\\alpha_k=\\theta_k/\\langle A p_k\\,,\\,p_k\\rangle \n", "&\\textit{(pas de descente)}\\\\ \n", "&x_{k+1}=x_k+\\alpha_k\\,p_k \n", "&\\textit{(mise \\`a jour de la solution)}\\\\\n", "&r_{k+1}=r_k-\\alpha_k\\,A p_k\n", "&\\textit{(résidu à l'itération } k+1)\\\\\n", "%&\\textit{Arr\\^et des it\\'erations} : $\\Vert r^{(k+1)}\\Vert\\le\\varepsilon$ ?\\\\\n", "\n", "&\\theta_{k+1}=\\langle r_{k+1},r_{k+1}\\rangle\\\\\n", "&\\beta_{k+1}=\\theta_{k+1}/\\theta_k\\\\\n", "&p_{k+1}=r_{k+1}+\\beta_{k+1}\\,p_k\n", "&\\textit{(nouvelle direction de descente)}\n", "\\end{aligned}$$\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e99ac998", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "46b98ece", "metadata": {}, "source": [ "### Question 9) \n", "\n", "Vérifier numériquement que la méthode converge en réalité en au plus $n$ itérations." ] }, { "cell_type": "code", "execution_count": null, "id": "1557fcfa", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3.6.5 64-bit", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.6.5" }, "vscode": { "interpreter": { "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49" } } }, "nbformat": 4, "nbformat_minor": 5 }