tp03-transport-sujet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# TP 3: Équation de transport en 1D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "On s'intéresse à la résolution numérique de l'équation de transport sur $[0,T]\\times [0,1]$ avec une vitesse variable $c$ avec des conditions au bord périodiques\n",
    "\n",
    "$$\n",
    "(T)\\quad \\left\\{\\begin{aligned}\n",
    "&\\frac{\\partial u}{\\partial t}(t,x) + c(x) \\frac{\\partial u}{\\partial x}(t,x) = 0  \\hbox{ sur } [0,1]\\times ]0,T]\\\\\n",
    "&u(t,0) = u(t,1) \\hbox{ pour } t\\in [0,T]\\\\\n",
    "&u(0,\\cdot) = u_0\n",
    "\\end{aligned}\\right. \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "\n",
    "## 1. Cas d'une vitesse constante\n",
    "Soit $M$ le nombre de pas d'espace, $x_j = h j$ où $h=1/M$ et $0\\leq j\\leq M-1$. Par périodicité, on a $u(x_M) = u(x_0)$ et il est donc inutile d'encoder le dernier point. La discrétisation du temps est la même que d'habitude, c'est-à-dire\n",
    "$t_n = n \\tau$ où $\\tau = T/N$ et $0\\leq n\\leq N$. On considère deux discrétisations possibles de l'opérateur $\\partial/\\partial x$:\n",
    "\n",
    "$$ \\mathrm{D}^+_h u_j = \\frac{u_{j+1} - u_j}{h} \\hbox{ et } \\mathrm{D}^-_{h} u_j = \\frac{u_{j} - u_{j-1}}{h}, $$\n",
    "\n",
    "pour $0\\leq j\\leq M-1$, où l'on a posé implicitement $u^n_{M} = u^n_0$ et $u^n_{-1} = u^n_{M-1}$. Pour $c\\geq 0$ on considère le schéma \n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{aligned}\n",
    "& u_j^{n+1} = u_j^n - c \\tau \\mathrm{D}^-_j u & \\hbox{ pour } 1\\leq j \\leq M-1, 0\\leq n \\leq N-1 \\\\\n",
    "& u_j^0 = u_0(x_j) & \\hbox{ pour } 0\\leq j\\leq M-1\\\\\n",
    "& u_M^n = u_0^n & \\hbox{ pour } 0\\leq n\\leq N\n",
    "\\end{aligned}\\right.\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "**Q1)** Écrire une fonction DP(u,h) retournant le vecteur des dérivées discrètes à droite et DM(u,h) retournant le vecteur des dérivées discrètes à gauche (ne pas hésiter à s'aider de la fonction np.roll). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# compléter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "**Q2)** Comparer numériquement les schémas décentrés avec $\\mathrm{D}_h^-$ et $\\mathrm{D}_h^+$ (tester en particulier l'instablité inconditionnelle du schéma décentré aval)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "c = 1.0\n",
    "M = 1000\n",
    "h = 1.0/M\n",
    "x = np.linspace(0,1-h,M)\n",
    "u = np.exp(-100.0*np.power(x-.4,2))\n",
    "N = 1000\n",
    "tau = .9/N\n",
    "\n",
    "# compléter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "**Q3)** Constater la forte diffusion numérique induite par le schéma lorsque $|c|\\frac{\\tau}{h} \\ll 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "c = 1.0\n",
    "M = 100\n",
    "h = 1.0/M\n",
    "x = np.linspace(0,1-h,M)\n",
    "u = np.exp(-100.0*np.power(x-.4,2))\n",
    "N = 1000\n",
    "tau = .9/N\n",
    "#print(\"CFL=%f\" % (np.abs(c)*tau/h))\n",
    "\n",
    "# compléter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "**Q4)** Considérer le schéma décentré amont d'ordre deux, construit en posant\n",
    "\n",
    "$$\\mathrm{D}^{-,2}_{h} u_j = \\frac{3u_{j} - 4u_{j-1} + u_{j-2}}{2h} \\qquad \\mathrm{D}^{+,2}_{h} u_j = \\frac{-3u_{j} + 4u_{j+1} - u_{j+2}}{2h}$$\n",
    "\n",
    "Si $c \\geq 1$ on posera par exemple\n",
    "\n",
    "$$ u_j^{n+1} = u_j^n + \\tau \\mathrm{D}_h^{-,2} u_j $$\n",
    "\n",
    "Étudier la consistance ce schéma. Comparer l'effet du choix de ce schéma sur le phénomène de diffusion numérique rencontre dans la question précédente (écrire au préalable une fonction DM2(u,h))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# compléter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 2. Cas d'une vitesse variables\n",
    "\n",
    "Dans le cas d'une vitesse variable, le choix de la dérivée discrète (à droite ou à gauche) choisie en un point $x_j$ va dépendre de $c_j := c(x_j)$. Précisément, \n",
    "\n",
    "$$ u_j^{n+1} = \n",
    "\\begin{cases}\n",
    "u_j^n - \\tau c_j \\mathrm{D}_h^- u_j &\\hbox{ si  } c_j \\geq 0 \\\\\n",
    "u_j^n - \\tau c_j \\mathrm{D}_h^+ u_j &\\hbox{ si  } c_j < 0,\n",
    "\\end{cases}$$\n",
    "\n",
    "ce qu'on peut mettre sous la forme compacte (et propice aux calculs: on peut extraire $\\min(c_j,0)$ via np.maximum(c,0) où $c$ est le vecteur des vitesses, et vectoriser les opérations)\n",
    "\n",
    "$$ u_j^{n+1} = u_j^n + \\max(c_j,0) \\mathrm{D}_h^- u_j + \\min(c_j,0) \\mathrm{D}_h^+ u_j. $$\n",
    "\n",
    "**Q5)** Montrer que le schéma est stable en norme infinie sous la condition CFL $\\|c\\|_\\infty \\frac{\\tau}{h} \\leq 1$.\n",
    "\n",
    "**Q6)** Montrer la consistance du schéma en une solution $u \\in \\mathcal{C}^1([0,T]\\times [0,1])$ périodique ($u(t,0) = u(t,1)$ et $u'(t,0) = u'(t,1)$ pour tout $t\\in [0,T]$) de l'équation de transport.\n",
    "\n",
    "**Q7)** Tester numériquement avec les données suivantes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "M = 1000\n",
    "h = 1.0/M\n",
    "x = np.linspace(0,1-h,M)\n",
    "u = np.exp(-200.0*np.power(x-.3,2)) +  np.exp(-200.0*np.power(x-.7,2))\n",
    "c = .5 - x #0.2 + .5*x*x\n",
    "N = 1000\n",
    "tau = .9/N\n",
    "\n",
    "# compléter"
   ]
  }
 ],
 "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
}