{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "06e7bcdb-1ecc-4bd2-9aa9-7fb025fc63cd",
   "metadata": {},
   "source": [
    "## MT09 - TP4 - Automne 2024\n",
    "### Calcul approché de valeurs propres, puissances itérées, puissances inverses, déflation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "668fd69c-d772-469c-a611-76e49a3797d4",
   "metadata": {},
   "source": [
    "#### Algorithme des puissances itérées\n",
    "1. Programmer l'algorithme des puissances itérées \n",
    "\n",
    "```\n",
    "xsol, lambdasol, kout, boolcvg = puissancesIterees(A, x0, tol, kmax)\n",
    "```\n",
    "\n",
    "pour une matrice $A$ et un vecteur initial $x_0\\neq 0$, une tolérance de précision $tol$ et un nombre maximal d'itérations $k_{max}$. La sortie ```boolcvg``` retournera 1 si la tolérance est atteinte, 0 sinon. Le critère de convergence à l'itération $k$ sera\n",
    "\n",
    "$$\n",
    "\\|A x^{(k)} - \\lambda^{(k)} x^{(k)}\\|< tol.\n",
    "$$\n",
    "L'entier $k_{out}$ est l'indice d'itération de sortie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5adc3ca8-8334-440d-8639-ca3719def520",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7bee834-9014-4442-9ce8-feb93a346dc2",
   "metadata": {},
   "source": [
    "Codez une matrice générique $A$ de taille $n$ définie par\n",
    "$$\n",
    "A = \\text{tridiag}(-1, 2, -1).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2de7ec55-6104-4e7f-a078-35599fce28b1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "n = 10;\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f4a0519-5a11-4948-9885-afbc8297017e",
   "metadata": {
    "tags": []
   },
   "source": [
    "Appliquer l'algorithme de puissances itérées à la matrice $A$ en choisissant $x_0=(1,0,...,0)^T=\\mathbf{e}_1$, $kmax=1000$ et $tol=10^{-5}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f86faf0-6a32-453e-89af-e156ef7ca5b6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28a1e657-3d87-4c83-ae2d-8bfe968a724e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xsol, lambdasol, kout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1600635c-9f64-464f-9e80-8651785dbd7c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(\"lambda = \", lambdasol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15eec56f-7467-43ff-8a2f-ebf483811df4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "la.norm(A@xsol - lambdasol*xsol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab425569-5195-4efb-afdc-052c60104a1c",
   "metadata": {
    "tags": []
   },
   "source": [
    "En fait, dans ce cas particulier, on peut calculer analytiquement les valeurs propres de $A$, elles sont données par la formule\n",
    "\n",
    "$$\n",
    "\\lambda_j = 4\\sin^2\\left(\\frac{j\\pi}{2(n+1)}\\right), \\quad 1\\leq j\\leq n.\n",
    "$$\n",
    "Vérifiez que l'algorithme de puissances itérées approche bien la plus grande valeur propre de $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e41ede6c-52c5-4ad6-8d35-0b121b2cafc8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b07f1ff2-4633-4fec-b1ef-a70471440fb9",
   "metadata": {},
   "source": [
    "#### Algorithme des puissances itérées inverses\n",
    "2. Programmer l'algorithme des puissances itérées inverses\n",
    "\n",
    "```\n",
    "xsol, lambdasol, kout, boolcvg = puissancesItereesInverses(A, x0, tol, kmax)\n",
    "```\n",
    "\n",
    "pour une matrice $A$ et un vecteur initial $x_0\\neq 0$, une tolérance de précision $tol$ et un nombre maximal d'itérations $k_{max}$. La sortie ```boolcvg``` retournera 1 si la tolérance est atteinte, 0 sinon. Le critère de convergence à l'itération $k$ sera\n",
    "\n",
    "$$\n",
    "\\|A x^{(k)} - \\lambda^{(k)} x^{(k)}\\|< tol.\n",
    "$$\n",
    "\n",
    "On utilisera ```np.linsolve()``` pour résoudre les systèmes linéaires de l'algorithme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e2f392f-724e-446c-a99e-0ac9715e2ae9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "75893490-2456-49bc-984d-ed8d11353003",
   "metadata": {},
   "source": [
    "Appliquer l'algorithme de puissances itérées à la matrice tridiagonale $A$ en choisissant $x_0=(1,0,...,0)^T$, $kmax=1000$ et $tol=10^{-5}$. Vérifiez que l'algorithme de puissances itérées inverses approche bien la plus petite valeur propre de $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff1366b4-8b6f-4aeb-b9c2-a27701cf3b27",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65aaa9bf-fa4c-41c2-9b79-c65f075b4b86",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a7a34f67-4c4c-4ada-96f5-8249f749acb4",
   "metadata": {},
   "source": [
    "### Méthode de déflation\n",
    "3. La méthode de déflation permet de calculer toutes les valeurs propres d'une matrice. Dans le cas d'une matrice symétrique, l'algorithme de déflation est particulièrement simple. Une fois calculé la plus grande valeur propre $\\lambda_n$ (en valeur absolue) avec le vecteur propre associé $\\mathbf{x}_n$, on considère la matrice \"retranchée\"\n",
    "\n",
    "$$\n",
    "A \\leftarrow A - \\lambda_n\\, \\mathbf{x}_n \\mathbf{x}_n^T,\n",
    "$$\n",
    "\n",
    "et on réapplique la méthode de puissances itérées. On itère jusqu'à obtenir la matrice nulle.\n",
    "Programmer une méthode \n",
    "\n",
    "```\n",
    "D, V = deflation(A)\n",
    "```\n",
    "\n",
    "qui calcule le tableau $D$ constitué des valeurs propres de $A$ et la matrice $V$ des vecteurs propres en faisant appel à ```puissancesIterees()```. Appliquer à la matrice tridiagonale précédente. Comparer les valeurs propres calculées aux valeurs propres exactes. Pour cette matrice $A$ symétrique, vérifiez que la matrice $V$ est une matrice orthogonale aux erreurs d'approximation près (on prendra $k_{max}=1000$ et $tol=10^{-8}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a75ddc1-ff0d-4b0d-be34-fa2e41b1193b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "213e0c1b-9006-4a1b-8d0e-a6d0d6e30a2f",
   "metadata": {},
   "source": [
    "#### \"Visualisation\" des vecteurs propres\n",
    "4. Tracez les $n$ vecteurs propres  de $A$ comme s'il s'agissait de fonctions discrétisées sur un intervalle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b72f3e64-3dca-4824-9d01-0dbc99df08ac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, axs = plt.subplots(5, 2)\n",
    "axs[0,0].plot(V[:,0], '-k')\n",
    "axs[0,1].plot(V[:,1], '-k')\n",
    "axs[1,0].plot(V[:,2], '-k')\n",
    "axs[1,1].plot(V[:,3], '-k')\n",
    "axs[2,0].plot(V[:,4], '-k')\n",
    "axs[2,1].plot(V[:,5], '-k')\n",
    "axs[3,0].plot(V[:,6], '-k')\n",
    "axs[3,1].plot(V[:,7], '-k')\n",
    "axs[4,0].plot(V[:,8], '-k')\n",
    "axs[4,1].plot(V[:,9], '-k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3e0918b-f551-4dad-82b1-74de0c4b8ce9",
   "metadata": {},
   "source": [
    "#### Application du calcul de valeurs propres. Mesure d'incertitude sur un résultat\n",
    "\n",
    "Un système est gouverné par un système algébrique\n",
    "\n",
    "$$\n",
    "\\mathbf{F}(\\mathbf{u}) = \\mathbf{b}\n",
    "$$\n",
    "\n",
    "où $\\mathbf{b}$ est le vecteur représentant une force extérieure et $\\mathbf{u}$ est la variable d'état du système (par exemple le vecteur de déplacement). On suppose que $\\mathbf{F}$ est continûment différentiable et \n",
    "qu'il existe un unique $\\mathbf{u}^\\star$ tel que $\\mathbf{F}(\\mathbf{u}^\\star)=\\mathbf{b}$.\n",
    "\n",
    "Un système de capteurs permet de mesurer $b$ avec une certaine incertitude. La certification des\n",
    "capteurs de mesure garantit une erreur relative\n",
    "\n",
    "$$\n",
    "\\frac{\\|\\delta \\mathbf{b}\\|_2}{\\|\\mathbf{b}\\|_2}\n",
    "$$\n",
    "\n",
    "inférieure à $\\eta$. Si une erreur $\\delta \\mathbf{b}$ est commise sur $\\mathbf{b}$, on a une réponse $\\mathbf{u}$ telle que $\\mathbf{F}(\\mathbf{u}) = \\mathbf{b}+\\delta \\mathbf{b}$. Sous l'hypothèse de petites perturbations, on peut linéariser $\\mathbf{F}$ selon\n",
    "\n",
    "$$\n",
    "\\mathbf{F}(\\mathbf{u}) \\approx \\mathbf{F}(\\mathbf{u}^\\star) + DF(\\mathbf{u}^\\star)(\\mathbf{u}-\\mathbf{u}^\\star) = \\mathbf{b} + DF(\\mathbf{u}^\\star)(\\mathbf{u}-\\mathbf{u}^\\star) = \\mathbf{b}+\\delta \\mathbf{b}.\n",
    "$$\n",
    "\n",
    "En notant $\\delta \\mathbf{u} = \\mathbf{u}-\\mathbf{u}^\\star$. On peut donc raisonnablement considérer que\n",
    "\n",
    "$$\n",
    "DF(\\mathbf{u}^\\star)\\, \\delta \\mathbf{u} = \\delta \\mathbf{b}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7200bce-7611-41e3-9056-28b25d8646cf",
   "metadata": {},
   "source": [
    "Q : Comment évaluer numériquement une majoration optimale de $\\dfrac{\\|\\delta \\mathbf{u}\\|_2}{\\|\\mathbf{u}\\|_2}$ ('optimale' au sens où la borne de majoration peut être atteinte) ~?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8f63619-1f72-4ed8-90a7-9d9876ffd91e",
   "metadata": {},
   "source": [
    "**Application** : on considérera $\\eta=10^{-4}$, $\\mathbf{F}$ linéaire\n",
    "\n",
    "$$\n",
    "\\mathbf{F}(\\mathbf{u}) = A\\mathbf{u}\n",
    "$$\n",
    "\n",
    "avec\n",
    "\n",
    "$$\n",
    "A = \\frac{1}{\\sqrt{\\varepsilon}} \\begin{pmatrix} 1 & 1 \\\\ 1-\\varepsilon & 1 \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "avec $\\varepsilon = 10^{-7}$. Evaluer numériquement une majoration $M$ optimale de\n",
    "$\\dfrac{\\|\\delta \\mathbf{u}\\|_2}{\\|\\mathbf{u}\\|_2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5f758f6-5002-4fe1-b93c-ffca14372a0b",
   "metadata": {
    "tags": []
   },
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}