tp04-chaleur2d-sujet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# TP 3: Équation de la chaleur en 2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "On s'intéresse à la résolution numérique de l'équation de la chaleur sur $[0,T]\\times \\Omega$ où $\\Omega = ]0,1[^d$ et $d\\in\\{1,2\\}$. On considèrera l'équation avec conditions de Neumann homogènes:\n",
    "\n",
    "$$\n",
    "(1D) \\quad \\left\\{\\begin{aligned}\n",
    "&\\frac{\\partial u}{\\partial t}(t,x) = \\frac{\\partial^2 u}{\\partial x^2}(t,x)  \\hbox{ sur } ]0,1[\\\\\n",
    "&\\frac{\\partial u}{\\partial x}(t,0) = \\frac{\\partial u}{\\partial x}(t,1) = 0 \\hbox{ pour } t\\in [0,T]\\\\\n",
    "&u(0,\\cdot) = u_0\n",
    "\\end{aligned}\\right. \n",
    "$$\n",
    "\n",
    "$$\n",
    "(2D) \\quad \\left\\{\\begin{aligned}\n",
    "&\\frac{\\partial u}{\\partial t}(t,x,y) = \\frac{\\partial^2 u}{\\partial x^2}(t,x,y) + \\frac{\\partial^2 u}{\\partial y^2}(t,x,y)   \\hbox{ sur } \\Omega\\\\\n",
    "&\\frac{\\partial u}{\\partial x}(t,0,y) = \\frac{\\partial u}{\\partial x}(t,1,y) = 0 \\hbox{ pour } t\\in [0,T], y \\in[0,1] \\\\\n",
    "&\\frac{\\partial u}{\\partial y}(t,x,0) = \\frac{\\partial u}{\\partial y}(t,x,1) = 0 \\hbox{ pour } t\\in [0,T], x \\in[0,1] \\\\\n",
    "&u(0,\\cdot,\\cdot) = u_0\n",
    "\\end{aligned}\\right. \n",
    "$$\n",
    "\n",
    "La discrétisation de (1D) et de (2D) par un schéma aux différences finies conduit à la résolution d'un système linéaire $Ax=b$ dont l'inconnue est le vecteur de la solution approchée aux points de la discrétisation. Dans ce TP, nous allons résoudre ce système en utilisant une méthode itérative. On a ici fait le choix d'utliser la méthode du gradient conjugué.\n",
    "\n",
    "On commence par importer les modules qui seront utilisés et par définir les fonctions permettant de résoudre itérativement le système linéaire par la méthode du gradient conjugué."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import scipy.sparse.linalg as sparse\n",
    "def solve_iterative(Operator, b):\n",
    "    M = len(b)\n",
    "    A = sparse.LinearOperator([M,M], matvec=Operator)\n",
    "    x,info = sparse.cg(A,b)\n",
    "    return x\n",
    "\n",
    "def solve_iterative_2d(Operator, b):\n",
    "    assert(b.shape[0] == b.shape[1])\n",
    "    M = b.shape[0]\n",
    "    A = sparse.LinearOperator([M*M,M*M], matvec=lambda u: np.reshape(Operator(np.reshape(u,(M,M))), M*M));\n",
    "    x,info = sparse.cg(A,np.reshape(b,M*M))\n",
    "    return np.reshape(x,(M,M))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 1.  Cas 1D avec résolution itérative du système linéaire\n",
    "\n",
    "On considère $x = (x_0,\\dots,x_{M+1})$ une grille uniforme dans $[0,1]$ de pas $h=1/(M+1)$ et $t = (t_0,\\dots,t_N)$ une grille uniforme dans $[0,T]$ de pas $\\tau = 1/N$. On note également $u_{j}^n \\simeq u(x_j,t_n)$ l'approximation discrète de la solution du système (1D). On rappelle que le schéma implicite pour (1D) est de la forme\n",
    "$$\n",
    "\\left\\{\\begin{aligned}\n",
    "& \\frac{1}{\\tau} (u_{j}^{n+1} - u_j^n) = \\frac{1}{h^2} (u_{j-1}^{n+1} - 2 u_j^{n+1} + u_{j+1}^{n+1}) &\\hbox{ pour } 1\\leq j \\leq M, 0\\leq n \\leq N-1\\\\\n",
    "& u_0^n = u_1^n \\hbox{ et } u_{M+1}^n = u_M^n &\\hbox{ pour } 0\\leq n\\leq N \\\\\n",
    "& u_j^0 = u_0(x_j) &\\hbox{ pour } 1\\leq j\\leq M+1\n",
    "\\end{aligned}\\right.,\n",
    "$$\n",
    "En posant $v^n := (u_1^n,\\dots,u^n_M)$, vérifier que le schéma peut être réécrit de la manière suivante:\n",
    "$$ v^{n+1} = \\left(\\mathrm{Id}_M - \\frac{\\tau}{h}(D_h^+ - D_h^-)\\right)^{-1} v^n$$\n",
    "où les opérateurs/matrices \"dérivées directionnelles\" $D_h^+$ et $D_h^-$ sont définis de la manière suivante:\n",
    "$$ D_h^+ v = \\frac{1}{h}(v_2 - v_1, \\dots, v_{M} - v_{M-1}, 0) \\in \\mathbb{R}^M $$\n",
    "$$ D_h^- v = \\frac{1}{h}(0, v_2 - v_1, \\dots, v_{M} - v_{M-1}) \\in \\mathbb{R}^M, $$\n",
    "pour $v=(v_1,\\dots,v_M)\\in\\mathbb{R}^M.$ On note que les zéros correspondent aux conditions de Neumann.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "**Q1)** Soit $u = $ np.random.random$(10)$. Tester l'accès aux éléments $(u_n,\\dots,u_{m-1})$ de $u,$ pour $m>n,$ via la syntaxe $u[n:m]$. En déduire deux fonctions DP$(u,h)$ DM$(u,h)$ calculant respectivement $D_h^+ u$ et $D_h^- u$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# à compléter\n",
    "\n",
    "# def DP(u,h):\n",
    "\n",
    "# def DM(u,h):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "En utilisant DM et DP, on peut définir une fonction python $A$ calculant $u\\mapsto (\\mathrm{Id}_M - (\\tau/h)(D_h^+ - D_h^-)) u$. L'intérêt est de pouvoir résoudre le système $A(u^{n+1}) = u^n$ par une méthode itérative qui ne nécessite pas de connaître $A$ mais seulement de savoir calculer $A(u)$ (ici, on utilise la méthode du gradient conjugué). En particulier, il n'y a aucun besoin de stocker la matrice en mémoire !\n",
    "\n",
    "Ecrire un programme qui permette de calculer les itérés du schéma implicite pour (D1). Tracer sur une même figure les itérés pour M=99, N=10 et T=0.1, pour chacune des valeurs initiales pour $u_0$ définies ci-dessus. Augmenter le temps final et N (dans ce cas ne pas tracer toutes les itérées). \n",
    "\n",
    "Lorsque $u_0(x)=\\cos(\\pi x),$ la solution de (1D) est donnée par $u(t,x)=e^{-t\\pi ^2}\\cos(\\pi x).$ Afficher dans ce cas sur une même figure la solution exacte et la solution approchée au temps final."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# A completer\n",
    "#N = 10\n",
    "#M = 99\n",
    "#T = 0.1\n",
    "#h = \n",
    "#tau = \n",
    "\n",
    "#x = \n",
    "\n",
    "# Donnée initiale\n",
    "# Cas 1\n",
    "# u = np.cos(np.pi*x)\n",
    "\n",
    "# Cas 2\n",
    "# u = np.exp(-10.0*(np.power(x - 0.5,2.0)))\n",
    "# u = u/np.max(u)\n",
    "\n",
    "# Cas 3\n",
    "# u = np.random.random(M)\n",
    "\n",
    "# On définit une fonction via la syntaxe \"lambda\",\n",
    "# ici, A(u) retourne u - (tau/h)*(D^+_hu-D^-_hu)\n",
    "\n",
    "# A compléter\n",
    "# A = lambda u: ....\n",
    "\n",
    "#for i in xrange(N):\n",
    "#    plt.plot(x,u)\n",
    "#    u = solve_iterative(A, u)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 2.  Cas 2D avec résolution itérative du système linéaire\n",
    "\n",
    "Pour discrétiser le problème en 2D on considère une discrétisation $(x_i)_{0\\leq i\\leq M+1\\,}, (y_j)_{0\\leq j\\leq M+1}$ de l'axe des abcisse et des ordonnées, avec $x_i = hi$, $y_j = hj$ et $h=1/(M+1)$. La discrétisation de l'espace est donnée par la grille de points $p_{i,j} = (x_i,y_j) \\in [0,1]^2, $\n",
    "et on note $u_{i,j}^n$ l'approximation de la solution $u$ en $(x_i,y_j,t_n)$. On peut écrire un schéma implicite pour (2D) de la manière suivante:\n",
    "$$\n",
    "\\left\\{\\begin{aligned}\n",
    "& \\frac{1}{\\tau} (u_{i,j}^{n+1} - u_{i,j}^n) = \\frac{1}{h^2} (u_{i,j-1}^{n+1} + u_{i,j+1}^{n+1} + u_{i-1,j}^{n+1} + u_{i+1,j}^{n+1} +   - 4 u_{i,j}^{n+1})&\\hbox{ pour } 1\\leq i,j \\leq M, 0\\leq n \\leq N-1\\\\\n",
    "& \\frac{u_{0,j}^n - u_{1,j}^n}{h} = \\frac{u_{M+1,j}^n - u_{M,j}^n}{h} = 0 &\\hbox{ pour } 0\\leq n\\leq N \\hbox{ et } 1\\leq j\\leq M \\\\\n",
    "& \\frac{u_{i,0}^n - u_{i,1}^n}{h} = \\frac{u_{i,M+1}^n - u_{i,M}^n}{h} = 0 &\\hbox{ pour } 0\\leq n\\leq N \\hbox{ et } 1\\leq i\\leq M \\\\\n",
    "& u_{i,j}^0 = u_0(p_{ij}) &\\hbox{ pour } 0\\leq i,j\\leq M+1\n",
    "\\end{aligned}\\right.,\n",
    "$$\n",
    "On remarque qu'il manque des équations pour les valeurs de $u$ aux quatres \"coins\", i.e. $u_{0,0}^n, u_{M+1,0}^n, u_{0,M+1}^n, u_{M+1,M+1}^n$, *qui ne font donc pas partie des inconnnues*.\n",
    "\n",
    "On peut vérifier qu'en considérant $v^n := (u^n_{i,j})_{1\\leq i,j \\leq M} \\in \\mathcal{M}_{M,M}$  comme une matrice de taille $M\\times M$, le schéma peut être réécrit de la manière suivante:\n",
    "$$ A(v^{n+1}) = v^n$$\n",
    "$$ A(v) = v - \\frac{\\tau}{h}(D_{x,h}^+ v  - D_{x,h}^- v + D_{y,h}^+ v - D_{y,h}^- v)$$\n",
    "où les opérateurs $D_{x,h}^\\pm,D_{y,h}^\\pm: \\mathcal{M}_{M,M}\\to  \\mathcal{M}_{M,M}$ correspondent aux dérivées partielles à droite et à gauche. Par exemple, la $j$ième ligne de $D_{x,h}^+ v$ est donnée par:\n",
    "$$ (D_{x,h}^+ v)_{\\cdot, j} = \\frac{1}{h}(v_{2,j} - v_{1,j}, \\dots, v_{M,j} - v_{M-1,j}, 0) $$\n",
    "Là encore, le zéro final provient de la condition de Neumann.\n",
    "\n",
    "**Q2)** Écrire des fonctions DxP(u,h), DxM(u,h), DyP(u,h) et DyM(u,h) prenant en entrée une matrice $u$ de taille $M\\times M$ et retournant  les matrices des dérivées partielles directionnelles $D_{x,h}^+(u)$, $D_{x,h}^-(u), D_{y,h}^+(u)$, et  $D_{y,h}^-(u)$. Pour faire concis, on n'hésitera pas à utliser $u[m:n, r:s]$ pour extraire la sous-matrice $(u_{i,j})_{m\\leq i\\leq n-1, r\\leq j\\leq s-1}$. (Expérimenter sur de petits exemples!) \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# A compléter\n",
    "\n",
    "# def DxP(u,h):\n",
    "\n",
    "# ...."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "**Q3)** Écrire une fonction calculant $A(u)$ où l'opérateur $A$ est défini ci-dessus. On supposera toujours que $u$ est représenté par une matrice. Résoudre l'équation de la chaleur en 2D pour les paramètres ci-dessous (on utilisera v = solve_iterative_2d(A,u) pour résoudre le système linéaire $A(u) = v$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def u0(X,Y):\n",
    "    u = np.exp(-10.0*(np.power(X - 0.3,2.0) + np.power(Y - 0.5,2.0)))\n",
    "    u = u + np.exp(-30.0*(np.power(X - 0.7,2.0) + np.power(Y - 0.1,2.0)))\n",
    "    u = u/np.max(u)\n",
    "    #u = u + 0.3*np.random.random((M,M)) # -> si on veut rajouter du bruit\n",
    "    return u\n",
    "\n",
    "M = 30\n",
    "h = 1.0/M\n",
    "T = 0.05\n",
    "N = 5\n",
    "tau = T/N\n",
    "\n",
    "x = np.linspace(h,1.0-h,M)\n",
    "y = np.linspace(h,1.0-h,M)\n",
    "X, Y = np.meshgrid(x,y) # X et Y sont des matrices encodant les fonctions coordonnées (x,y) -> x  et (x,y) -> y\n",
    "\n",
    "u = u0(X,Y)\n",
    "# On peut visualiser u via la fonction plt.imshow(u, vmin=0.0, vmax=1.0, cmap='gray'). Ne pas hésiter à faire de même pour\n",
    "# les fonctions coordonnées.\n",
    "\n",
    "# A compléter\n",
    "\n",
    "# A = lambda u:\n",
    "\n",
    "#for i in xrange(4):\n",
    "#    plt.imshow(u, vmin=0.0, vmax=1.0, cmap=None)\n",
    "#    plt.colorbar()\n",
    "#    plt.show()\n",
    "#    u = solve_iterative_2d(A, u)\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'module' object has no attribute 'misc'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-f407b705a2fe>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mscipy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0msp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmisc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mascent\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mM\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m512\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m0.3\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'module' object has no attribute 'misc'"
     ]
    }
   ],
   "source": [
    "import scipy as sp\n",
    "u = sp.misc.ascent().astype(float)\n",
    "u = u / np.max(u)\n",
    "M = 512\n",
    "u = u + 0.3*np.random.random((M,M)) \n",
    "h = 1.0/M\n",
    "tau = 0.00002\n",
    "A = lambda u: u - (tau/h)*(DxP(u,h) - DxM(u,h) + DyP(u,h) - DyM(u,h))\n",
    "for i in xrange(3):\n",
    "    plt.imshow(u, cmap='gray', vmin=0.0, vmax=1.0)\n",
    "    plt.colorbar()\n",
    "    plt.show()\n",
    "    un = solve_iterative_2d(A, u)\n",
    "    print(np.max(np.abs(A(un) - u)))\n",
    "    u = un"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 3. Analyse du schéma 2D\n",
    "\n",
    "**Q1)** Montrer que les solutions du schéma (en supposant qu'elles existent) vérifient le principe du maximum discret\n",
    "$$ \\forall n< N,  \\min_{1\\leq i,j \\leq M} u_{i,j}^n \\leq \\min_{1\\leq i,j \\leq M} u_ij^{n+1} $$\n",
    "$$ \\forall n< N,  \\max_{1\\leq i,j \\leq M} u_{i,j}^{n+1} \\leq \\max_{1\\leq i,j \\leq M} u_{i,j}^{n}. $$\n",
    "En déduire la stabilité du schéma en norme infinie:\n",
    "$$ \\max_{0\\leq n\\leq N}\\max_{0\\leq i,j\\leq M+1} |u_{i,j}^n| \\leq \\|u_0\\|_\\infty $$\n",
    "\n",
    "**Q2)** Étudier la consistance du schéma.\n",
    "[comment]: %$$ \\frac{1}{\\tau}(u(x_i,x_j,t_{n+1}) - (u(x_i,x_j,t_n)) - \\frac{1}{h^2} (u(x_i,y_{j-1},t_{n+1}) + \n",
    "[comment]: %u(x_i,y_{j+1},t_{n+1})+u(x_{i+1},y_{j},t_{n+1})+u(x_{i-1},y_{j},t_{n+1}) - 4u(x_{i},y_{j},t_{n+1}))$$\n",
    "\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": 1
}