tp08-EF-VP-sujet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# TP 8: Approximation des valeurs propres du Laplacien #\n",
    "\n",
    "## 1. Formulation par éléments finis ## \n",
    "\n",
    "$$\n",
    "\\newcommand{\\sca}[2]{\\langle #1\\mid #2 \\rangle}\n",
    "\\newcommand{\\LL}{\\mathrm{L}} \n",
    "\\newcommand{\\Rsp}{\\mathbb{R}}\n",
    "\\newcommand{\\Psp}{\\mathbb{P}}\\phantom{x}\n",
    "$$\n",
    "On s'intéresse à la recherche des valeurs propres du laplacien avec conditions de Neumann d'un ouvert $\\Omega\\subseteq \\Rsp^2$ par la méthode des éléments finis en 2D. On part d'une triangulation $T_h$ de $\\Omega$, et on note:\n",
    "* $V_h$ l'espace des éléments finis $\\Psp_1$ sur $T_h$, de sorte que $V_h\\subseteq V := H^1(\\Omega)$;\n",
    "* $N_h :=\\dim(V_h)$;\n",
    "* $H := \\LL^2(\\Omega) $ ;\n",
    "\n",
    "Il s'agit alors de trouver $\\lambda_{1,h}\\leq \\dots \\lambda_{N_h,h} \\in \\Rsp$ et $w_{1,h},\\dots,w_{N_h,h}$ une base orthonormée de $V_h$ telle que \n",
    "\n",
    "$$\n",
    "\\forall m\\in\\{1,\\dots,N_h\\},~\\forall v_h\\in V_h,\\qquad a(w_{m,h}, v_h) = \\lambda_{m,h} \\sca{w_{m,h}}{v_h}_H \n",
    "$$\n",
    "\n",
    "où $a(u,v) = \\int_{\\Omega} \\sca{\\nabla u}{\\nabla v} + uv$.\n",
    "\n",
    "Pour tout sommet $x_i$ de la triangulation, on note $\\phi_i$ la fonction chapeau correspondante, i.e. la fonction $\\phi_i\\in V_h$ telle que $\\phi_i(x_j) = \\delta_{i,j}$.\n",
    "\n",
    "**Q1.** Montrer que $\\lambda\\in \\Rsp$ et $w = \\sum_{1\\leq i\\leq N_h} W_i \\phi_i\\neq 0$ est un couple de valeur/vecteur propre pour le problème variationnel  si et seulement si \n",
    "\n",
    "$$ (K+M) W = \\lambda M W$$\n",
    "\n",
    "où les matrice $K$ et $M$ sont à préciser.\n",
    "\n",
    "**Q2.** On suppose donnée une matrice $D_x$ de taille $NT_h\\times N_h$ où $NT_h$ est le nombre de triangles, et telle que si $u = \\sum_i U_i \\phi_i$, alors $(D_x U)_j = \\left.\\partial_x u\\right|_{\\sigma_j}$ où $\\sigma_j$ est le $j$ième simplexe. En d'autre terme, le vecteur $D_x U$ donne la liste des dérivées partielles par rapport à la première coordonnée $x$ de la fonction $u$ sur chaque triangle. Une seconde matrice $D_y$ est définie de même, et on note $A = (\\mathrm{aire}(\\sigma_j))_{1\\leq j\\leq NT_h}$. Donner une expression simple de $K$ en fonction de $D_x, D_y$ et $\\mathrm{diag}(A)$.\n",
    "\n",
    "**Q3.** Montrer que la matrice $M$ n'est jamais diagonale.\n",
    "\n",
    "On va approcher le problème de calcul de valeurs propres en remplaçant l'expression $M_{ij} = \\int_{\\Omega} \\phi_i\\phi_j$\n",
    "par une formule de quadrature d'ordre $1$, pour un triangle $\\sigma$ de sommets $A,B,C$\n",
    "$$ I_\\sigma(f) := \\int_\\sigma f \\simeq \\frac{\\mathrm{aire}(\\sigma)}{3}(f(A) + f(B) + f(C)).$$\n",
    "\n",
    "**Q4.** On définit la matrice de masse condensée (on parle de *mass lumping* en anglais) par\n",
    "\n",
    "$$\\tilde{M}_{ij} = \\sum_{\\sigma\\in T_h} I_\\sigma(\\psi_i\\psi_j). $$\n",
    "\n",
    "Montrer que $\\tilde{M}$ est diagonale inversible et calculer ses coefficients. On s'intéresse maintenant à la résolution du problème aux valeurs propres  $(K+\\tilde M) W = \\lambda \\tilde{M}$.\n",
    "\n",
    "**Q5.** Construire explicitement une matrice symétrique $Z$ (dépendant de $K$ et $\\tilde{M}$) telle que\n",
    "\n",
    "$$ \\exists W\\neq 0 \\hbox{ t.q. } (K+\\tilde M) W = \\lambda \\tilde{M} W \\Longleftrightarrow \\lambda \\hbox{ est valeur propre de } Z $$\n",
    "\n",
    "**Q6.** Construire les matrices $K$, $\\tilde{M}$ et $Z$ en Python, et calculer les 10 premières valeurs propres dans le cas $\\Omega = [0,1]^2$ et $h=1/30$. Pour le calcul des k plus petites valeurs propres de $Z$, on utilisera l'instruction\n",
    "*splinalg.eigsh(Z,k,sigma=0)* \n",
    "\n",
    "(NB: splinalg est le package d'algèbre linéaire sparse (creuse), eigsh est la fonction de calcul des valeurs propres pour des matrices hermitiennes (h))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.tri as tri\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "import scipy.sparse.linalg as splinalg\n",
    "\n",
    "# calcule la normale d'une famille de triangles triangle 3D, en déduit le gradient des fonctions dont \n",
    "# les triangles sont le graphe\n",
    "def tri_gradient(A,B,C):\n",
    "    N = np.cross(B - A, C-B)\n",
    "    Gx = -N[:,0] / N[:,2]\n",
    "    Gy = -N[:,1] / N[:,2]\n",
    "    return Gx, Gy\n",
    "\n",
    "# calcule les matrices Dx, Dy et A \n",
    "# NB: on utilise le package scipy.sparse pour un stockage \"creux\" des matrices,\n",
    "# c'est-à-dire qu'on ne stocke que les entrées non nulles de la matrice.\n",
    "def derivees_et_aires(X,T):\n",
    "    nt = T.shape[0];\n",
    "    npts = X.shape[0];\n",
    "  \n",
    "    A = X[T[:,0],:]; # liste des points A des triangles\n",
    "    B = X[T[:,1],:]; # liste des points B des triangles\n",
    "    C = X[T[:,2],:]; # liste des points C des triangles\n",
    "  \n",
    "    # calcul de dx et dy sur chaque triangle, en fonction de la valeur de\n",
    "    # la fonction aux sommets du triangle\n",
    "    Z = np.zeros((nt,1));\n",
    "    I = np.ones((nt,1));\n",
    "    [Gx1,Gy1] = tri_gradient(np.hstack((A,I)),np.hstack((B,Z)),np.hstack((C,Z)));\n",
    "    [Gx2,Gy2] = tri_gradient(np.hstack((A,Z)),np.hstack((B,I)),np.hstack((C,Z)));\n",
    "    [Gx3,Gy3] = tri_gradient(np.hstack((A,Z)),np.hstack((B,Z)),np.hstack((C,I)));\n",
    "  \n",
    "    # assemblage des matrices creuses\n",
    "    # pour chaque triangle 1\\leq t \\leq nt, on a trois coefficients non nuls\n",
    "    # en (t,a), (t,b), (t,c) où a = T[i,1], b=T[i,2], c=T[i,3]\n",
    "    I = np.hstack((np.arange(nt), np.arange(nt), np.arange(nt)))\n",
    "    J = np.hstack((T[:,0], T[:,1], T[:,2]))\n",
    "    Sx = np.hstack((Gx1,Gx2,Gx3))\n",
    "    Sy = np.hstack((Gy1,Gy2,Gy3))\n",
    "    \n",
    "    # les matrices creuses sont construites de la manière suivante:\n",
    "    # pour tout l dans 0, ..., len(J)-1, on pose dx[I[k],J[k]] = Sx[k]\n",
    "    dx = sp.csc_matrix((Sx,(I,J)),shape=(nt,npts))\n",
    "    dy = sp.csc_matrix((Sy,(I,J)),shape=(nt,npts))\n",
    "    \n",
    "    # calcul de l'aire des triangles:\n",
    "    AB = B-A\n",
    "    BC = C-B\n",
    "    aires = .5*(AB[:,0]*BC[:,1] - AB[:,1]*BC[:,0])\n",
    "    return dx, dy, aires\n",
    "\n",
    "# construction de la triangulation et des matrices Dx, Dy et vecteur A des aires\n",
    "def triangulation_carre(n):\n",
    "    x,y = np.meshgrid(np.linspace(0.,1.,n),\n",
    "                      np.linspace(0.,1.,n))\n",
    "    x = x.reshape(n*n,1)\n",
    "    y = y.reshape(n*n,1)\n",
    "    X = np.hstack((x,y))\n",
    "    T = tri.Triangulation(x.flatten(), y.flatten()).triangles\n",
    "    return X,T\n",
    "\n",
    "X,T = triangulation_carre(30)\n",
    "[Dx,Dy,aires] = derivees_et_aires(X,T)\n",
    "\n",
    "# assemblage de la matrice de masse condensée\n",
    "#... compléter ...\n",
    "\n",
    "# assemblage de la matrice de rigidité\n",
    "#... compléter ...\n",
    "\n",
    "# sysème équivalent\n",
    "#... compléter ...\n",
    "\n",
    "[lbd,psi] = splinalg.eigsh(Z,10,sigma=0)\n",
    "for i in xrange(1,numel(lbd)):\n",
    "    print \"valeurs propre %d: %g\" % (i, lbd[i])\n",
    "    v = Linv.dot(psi[:,i])\n",
    "    plt.tripcolor(X[:,0],X[:,1],T,v,shading='gouraud')\n",
    "    plt.gca().set_aspect('equal')\n",
    "    plt.show()\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "**Q7.** Refaire les mêmes calculs dans le cas d'une couronne triangulée."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Créer la triangulation d'une couronne\n",
    "def triangulation_couronne(nrad=8,nang=36):\n",
    "    radii = np.linspace(0.25, 1, nrad)\n",
    "    angles = np.linspace(0, 2*np.pi, nang, endpoint=False)\n",
    "    angles = np.repeat(angles[..., np.newaxis], nrad, axis=1)\n",
    "    angles[:, 1::2] += np.pi/nang\n",
    "    npts = len(radii)*len(angles)\n",
    "    x = (radii*np.cos(angles)).flatten()\n",
    "    y = (radii*np.sin(angles)).flatten()\n",
    "\n",
    "    # On construit une triangulation de Delaunay et on filtre les triangles proches de l'origine\n",
    "    triang = tri.Triangulation(x, y)\n",
    "    xmid = x[triang.triangles].mean(axis=1)\n",
    "    ymid = y[triang.triangles].mean(axis=1)\n",
    "    mask = xmid*xmid + ymid*ymid >= .25*.25\n",
    "    T = triang.triangles[mask,:]\n",
    "    X = np.hstack((x.reshape(npts,1), y.reshape(npts,1)))\n",
    "    return X,T\n",
    "\n",
    "# à compléter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 2. Recherche de valeurs propres par méthode de la puissance ## \n",
    "\n",
    "**Q8.** Soit $Z$ une matrice symétrique de taille $N$, $\\lambda_1,\\dots, \\lambda_N$ ses valeurs propres et $w_1,\\dots,w_N$ une base orthonormale de vecteurs propres de $\\Rsp^N$. On suppose que pour tout $i\\neq N$, $|\\lambda_{i}| < |\\lambda_N|$. Soit $x_0 \\in \\Rsp^N$ tel que $\\sca{x_0}{w_N} \\neq 0$ et $x_{k+1} = A x_k / \\|A x_k\\|$. Démontrer que\n",
    "\n",
    "$$ \\min_{\\varepsilon \\in \\{\\pm 1\\}} \\left\\|x_{k+1} - \\varepsilon w_N\\right\\| \\leq C(x_0) \\lambda_{N-1}. $$\n",
    "\n",
    "En déduire que $\\lim_{k\\to\\infty} \\|A x_k\\| = |\\lambda_N|$.\n",
    "\n",
    "Pourquoi la méthode fonctionne-t'elle avec probabilité $1$ si $x_0$ est tiré aléatoirement selon une loi Gaussienne (par exemple) ?\n",
    "\n",
    "**Q9.** Ici, nous souhaitons trouver la plus *petite* valeur propre. Soit encore $Z$ symétrique réelle et $\\sigma \\not\\in \\{\\lambda_1,\\dots,\\lambda_N\\}$. Pour $x_0 \\in \\Rsp^N$ on pose\n",
    "\n",
    "$$ x_{k+1} = (Z - \\sigma I_N)^{-1} x_k $$\n",
    "\n",
    "On suppose qu'il existe $i_0 \\in \\{1,\\dots,N\\}$ tel que $\\forall i\\neq i_0$, $|\\lambda_i - \\sigma| < |\\lambda_{i_0} - \\sigma|$. Démontrer que si $\\sca{x_0}{w_{i_0}} \\neq 0$ alors \n",
    "\n",
    "$$ \\lim_{N\\to\\infty} x_i = \\pm w_{i_0}. $$\n",
    "\n",
    "Cette méthode est appelée méthode de la puissance inverse.\n",
    "\n",
    "**Q10.** Programmer la méthode de la puissance inverse, et la tester pour trouver le quatrième vecteur propre du laplacien de Neumann pour le carré en choisissant $\\sigma$ une approximation numérique de $\\lambda_4$. On utilisera *splinalg.spsolve* pour résoudre les systèmes linéaires creux (et *sp.eye* pour construire la matrice identité creuse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "N = X.shape[0]\n",
    "I = sp.eye(N)\n",
    "psi = np.random.rand(N)\n",
    "sigma = 20.0\n",
    "\n",
    "# compléter \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}