{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite volume schemes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "plt.rcParams['font.size'] = 12.0\n",
    "plt.rcParams['figure.figsize']=[8,6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Advection problem – constant or variable velocity\n",
    "\n",
    "We first solve the easy linear problem\n",
    "$$\\partial_t u + \\partial_x (a u) = 0,\\quad x\\in[0,L]$$\n",
    "with periodic boundary conditions. Firstly, the velocity $a>0$ is constant in space and time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Some tools:**  \n",
    "> The following function will be usefull to represent the discrete solutions $(U_j^n)$ as functions of space and time.  \n",
    "> General syntax: `SHOW(data=UTX,leg='scheme')`  \n",
    "> \n",
    "> `SHOW()` to draw the data stored in the 2-dimensional array `UTX`.  \n",
    ">\n",
    "> `SHOW(UTX_other)` to draw the data stored in another array `UTX_other`.\n",
    ">\n",
    "> `SHOW(UTX_other,leg='truc')` to add the title `truc` to the figures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "UTX = []\n",
    "def SHOW(data=UTX,leg='scheme'):\n",
    "    fig,(ax1, ax2) = plt.subplots(1,2,figsize=(14, 6))\n",
    "    im = ax1.imshow(data[1:,:-1],interpolation='none',\\\n",
    "                   origin=\"lower\",extent=[0,L,0,T],\\\n",
    "                   cmap=plt.get_cmap('jet'))\n",
    "    fig.colorbar(im,ax=ax1)\n",
    "    ax1.set_aspect('auto')\n",
    "    ax1.set_xlabel('$x$')\n",
    "    ax1.set_ylabel('$t$')\n",
    "    ax2.plot(data[0,:-1],u0(data[0,:-1]),'.-',data[0,:-1],data[-1,:-1],'x-')\n",
    "    ax2.set_aspect('auto')\n",
    "    ax2.set_xlabel('$x$')\n",
    "    ax2.set_ylabel('$u$')\n",
    "    ax2.legend(['Initial','Final'])\n",
    "    fig.suptitle(leg, fontsize=14, weight='bold')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "###############################\n",
    "# PDE problem\n",
    "T = 5       # final time\n",
    "L = 5       # space length\n",
    "\n",
    "def Fadv(u):\n",
    "    a = 1\n",
    "    return a*u, a*np.ones_like(u)  # returns f(u) and f'(u) (usefull for some CFL bounds)\n",
    "def u0(x):   # initial data\n",
    "    return np.array([0.5-0.5*np.sign(y-1)*np.sign(y-2) for y in x])\n",
    "Flux = Fadv # flux function\n",
    "\n",
    "###############################\n",
    "\n",
    "# Numerical parameters\n",
    "J = 20     # number of points\n",
    "mu = 1  # cfl parameter\n",
    "\n",
    "# Discretization\n",
    "dx = L/J\n",
    "x = np.linspace(0, L, J+1)\n",
    "x = x[0:-1] # removes the very last point for periodicity\n",
    "\n",
    "###############################\n",
    "\n",
    "## Initialisation\n",
    "t = 0\n",
    "U = u0(x)\n",
    "UTX = np.concatenate([[x],[U]],axis=0)\n",
    "LT = [t,t]\n",
    "\n",
    "F,dF = Flux(U)           # compute f(U_j^0)\n",
    "Vmax = max(np.abs(dF))   # get a bound for velocities (to adapt depending on the properties of f)\n",
    "dt = mu*dx/max(1,Vmax)   # max(1,Vmax) to avoid very large dt (bad for accuracy)\n",
    "\n",
    "## Time iterations ##\n",
    "while t<T:\n",
    "    dt = min([dt,T-t])\n",
    "    V = np.concatenate(([U[-1]],U,[U[0]]))    # boundary conditions\n",
    "    ##########################################\n",
    "    F,dF = Flux(V)                            # compute f(U_j)\n",
    "    # then the discrete fluxes F_{j+1/2}\n",
    "    NumFlux = F[:-1]                          # left flux\n",
    "    ##########################################\n",
    "    U = U - dt/dx*(NumFlux[1:]-NumFlux[:-1])  # Update\n",
    "    t = t+dt\n",
    "\n",
    "    LT = np.append(LT,t)\n",
    "    UTX = np.concatenate([UTX,[U]],axis=0)\n",
    "    \n",
    "LT  = np.resize(LT,(len(LT),1))\n",
    "UTX = np.concatenate([UTX,LT],axis=1)\n",
    "\n",
    "# Storage\n",
    "UTX1 = UTX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PDE problem\n",
    "T = 5       # final time\n",
    "L = 5       # space length\n",
    "def Fadv(u):\n",
    "    a = 1\n",
    "    return a*u, a*np.ones_like(u)  # returns f(u) and f'(u) (usefull for some CFL bounds)\n",
    "def u0(x):   # initial data\n",
    "    return np.array([0.5-0.5*np.sign(y-1)*np.sign(y-2) for y in x])\n",
    "Flux = Fadv # flux function\n",
    "\n",
    "\n",
    "# Numerical parameters\n",
    "J = 20     # number of points\n",
    "mu = 0.9  # cfl parameter\n",
    "\n",
    "# Discretization\n",
    "dx = L/J\n",
    "x = np.linspace(0, L, J+1)\n",
    "x = x[0:-1] # removes the very last point for periodicity\n",
    "\n",
    "## Initialize\n",
    "U = u0(x)\n",
    "t = 0\n",
    "UTX = np.concatenate([[x],[U]],axis=0)\n",
    "LT = [t,t]\n",
    "\n",
    "F,dF = Flux(U)           # compute f(U_j^0)\n",
    "Vmax = max(np.abs(dF))   # get a bound for velocities (to adapt depending on the properties of f)\n",
    "dt = mu*dx/max(1,Vmax)   # max(1,Vmax) to avoid very large dt (bad for accuracy)\n",
    "\n",
    "## Time iterations ##\n",
    "while t<T:\n",
    "    dt = min([dt,T-t])\n",
    "    V = np.concatenate(([U[-1]],U,[U[0]]))    # boundary conditions\n",
    "    ##########################################\n",
    "    F,dF = Flux(V)                            # compute f(U_j)\n",
    "    # then the discrete fluxes F_{j+1/2}\n",
    "    NumFlux = F[:-1]                          # left flux\n",
    "    ##########################################\n",
    "    U = U - dt/dx*(NumFlux[1:]-NumFlux[:-1])  # Update\n",
    "    t = t+dt\n",
    "\n",
    "    LT = np.append(LT,t)\n",
    "    UTX = np.concatenate([UTX,[U]],axis=0)\n",
    "\n",
    "LT  = np.resize(LT,(len(LT),1))\n",
    "UTX = np.concatenate([UTX,LT],axis=1)\n",
    "    \n",
    "UTX2=UTX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SHOW(UTX1,leg='Left flux for $\\mu=1$')\n",
    "SHOW(UTX2,leg='Left flux for $\\mu=0.9$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 1.** For the two above CFL parameters $\\mu\\in\\{0.9,1\\}$, compute the solution with $J=200$ cells to observe the convergence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Answer 1.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 2.** Now implement the Lax-Friedrichs flux  \n",
    "> $$F_{j+1/2}^n = \\dfrac{1}{2}(f(U_j^n)+f(U_{j+1}^n)) - \\dfrac{\\sigma}{2}(U_{j+1}^n-U_{j}^n),$$\n",
    "> with $\\sigma = \\max\\big\\{|f'(u)|,\\ u \\in[\\inf u_0,\\sup u_0]\\big\\}.$  \n",
    "> Compare the discrete solution with $\\mu=0.9$ to the one obtained with the previous left scheme."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Answer 2.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advection with time+space-depending velocities\n",
    "\n",
    "We consider now the linear problem\n",
    "$$\\partial_t u + \\partial_x (a(t,x) u) = 0$$\n",
    "with periodic boundary conditions.\n",
    "\n",
    "For the tests, we consider the following velocity field\n",
    "$$ a(t,x)= 2\\cos(2\\pi t)\\sin(2\\pi x),$$\n",
    "varying both in time and space, for $x\\in[0,4]$.  \n",
    "The initial data is\n",
    "$$u_0(x) = e^{-9(x-1)^2}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FadvTX(u,t,x):\n",
    "    a = 2*np.cos(2*np.pi*t)*np.sin(2*np.pi*x)\n",
    "    return a*u, a\n",
    "\n",
    "def u0(x):\n",
    "    return np.exp(-9*(x-2)**2)\n",
    "    #return np.array([0.5-0.5*np.sign(y-1.5)*np.sign(y-2.5) for y in x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "T = 3         # final time\n",
    "L = 4         # space length\n",
    "\n",
    "x    = np.linspace(0,L,21)\n",
    "LT   = np.linspace(0,T,17)\n",
    "U, V = np.meshgrid(x, LT)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "U = np.cos(2*np.pi*V)*np.sin(2*np.pi*U)\n",
    "V = np.ones_like(V)\n",
    "q = ax.quiver(x, LT, U, V)\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('t')\n",
    "ax.set_title('Velocity field $a(t,x)$', fontsize=14, weight='bold')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the Lax-Friedrichs scheme with $F_{j+1/2}^n \\overset{\\rm def}{=} \\dfrac{1}{2}(f(U_j^n,t^n,x_{j+1/2})+f(U_{j+1}^n,t^n,x_{j+1/2})) - \\dfrac{\\sigma}{2}(U_{j+1}^n-U_{j}^n)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 3         # final time\n",
    "L = 4         # space length\n",
    "\n",
    "# PDE problem\n",
    "Flux = FadvTX # flux function\n",
    "\n",
    "# Numerical data\n",
    "J = 36     # number of points\n",
    "mu = 1    # cfl parameter\n",
    "\n",
    "# Discretization parameters\n",
    "dx = L/J\n",
    "x = np.linspace(0, L, J+1)\n",
    "x = x[0:-1]\n",
    "\n",
    "## Initialize\n",
    "U = u0(x)\n",
    "t = 0\n",
    "UTX = np.concatenate([[x.copy()],[U]],axis=0)\n",
    "LT = [t,t]\n",
    "\n",
    "F,dF = Flux(U,t,x)       # compute f(U_j^0)\n",
    "Vmax = max(np.abs(dF))   # get a bound for velocities (to adapt depending on the properties of f)\n",
    "dt = mu*dx/max(1,Vmax)   # max(1,Vmax) to avoid very large dt (bad for accuracy)\n",
    "\n",
    "y = np.concatenate(([-dx],x,[L])) + dx/2\n",
    "\n",
    "## Time iterations ##\n",
    "while t<T:\n",
    "    dt = min([dt,T-t])\n",
    "    V = np.concatenate(([U[-1]],U,[U[0]]))\n",
    "    ##########################################\n",
    "    FL,dFL = Flux(V[:-1],t,y[:-1])                       \n",
    "    FR,dFR = Flux(V[1:] ,t,y[:-1])\n",
    "    Vmax = max(np.abs(dFL))\n",
    "    sigma = Vmax\n",
    "    NumFlux = 0.5*(FL + FR) - 0.5*sigma*(V[1:]-V[:-1])\n",
    "    ##########################################\n",
    "    U = U - dt/dx*(NumFlux[1:]-NumFlux[:-1])\n",
    "    t = t + dt\n",
    "\n",
    "    LT = np.append(LT,t)\n",
    "    UTX = np.concatenate([UTX,[U]],axis=0)\n",
    "    \n",
    "    \n",
    "LT  = np.resize(LT,(len(LT),1))\n",
    "UTX = np.concatenate([UTX,LT],axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "SHOW(UTX,leg='Lax-Friedrichs flux ; $\\mu=1$ ; '+str(J)+' cells')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 3.** From the properties of the velocity field above, explain why the choice of a grid with $J=8p+4$ cells is interesting. What properties are expected to be satisfied for the exact solution to the problem and what kind of defect is however present in the numerical solution ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We perform several computations with finer grids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for J in [8*10+4,8*20+4,8*40+4,8*80+4,8*160+4]:\n",
    "\n",
    "    mu = 1    # cfl parameter\n",
    "    \n",
    "    # Discretization parameters\n",
    "    dx = L/J\n",
    "    x = np.linspace(0, L, J+1)\n",
    "    x = x[0:-1]\n",
    "\n",
    "    ## Initialize\n",
    "    U = u0(x)\n",
    "    t = 0\n",
    "    UTX = np.concatenate([[x.copy()],[U]],axis=0)\n",
    "    LT = [t,t]\n",
    "\n",
    "    F,dF = Flux(U,t,x)       # compute f(U_j^0)\n",
    "    Vmax = max(np.abs(dF))   # get a bound for velocities (to adapt depending on the properties of f)\n",
    "    dt = mu*dx/max(1,Vmax)   # max(1,Vmax) to avoid very large dt (bad for accuracy)\n",
    "\n",
    "    y = np.concatenate(([-dx],x,[L])) + dx/2\n",
    "\n",
    "    ## Time iterations ##\n",
    "    while t<T:\n",
    "        dt = min([dt,T-t])\n",
    "        V = np.concatenate(([U[-1]],U,[U[0]]))\n",
    "        ##########################################\n",
    "        FL,dFL = Flux(V[:-1],t,y[:-1])                       \n",
    "        FR,dFR = Flux(V[1:] ,t,y[:-1])\n",
    "        Vmax = max(np.abs(dFL))\n",
    "        sigma = Vmax\n",
    "        NumFlux = 0.5*(FL + FR) - 0.5*sigma*(V[1:]-V[:-1])\n",
    "        ##########################################\n",
    "        U = U - dt/dx*(NumFlux[1:]-NumFlux[:-1])\n",
    "        t = t + dt\n",
    "\n",
    "        LT = np.append(LT,t)\n",
    "        UTX = np.concatenate([UTX,[U]],axis=0)\n",
    "    \n",
    "    LT  = np.resize(LT,(len(LT),1))\n",
    "    UTX = np.concatenate([UTX,LT],axis=1)\n",
    "    \n",
    "    plt.plot(x,U,label='J='+str(J))\n",
    "    \n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "SHOW(UTX,leg='Lax-Friedrichs flux ; $\\mu=1$ ; '+str(J)+' cells')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 4.** Did the defect disappear for large values of $J$ ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider now another scheme for the resolution of the same problem, using an upwind strategy.\n",
    "The numerical flux is\n",
    "$$F_{j+1/2}^n = \\begin{cases}f(U_j^n,t_n,x_{j+1/2}), & \\textrm{ if } a(t_n,x_{j+1/2})>0 ; \\\\ \n",
    "f(U_{j+1}^n,t_n,x_{j+1/2}), & \\textrm{ if } a(t_n,x_{j+1/2})<0.\\end{cases}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 3         # final time\n",
    "L = 4         # space length\n",
    "\n",
    "# PDE problem\n",
    "Flux = FadvTX # flux function\n",
    "\n",
    "# Numerical data\n",
    "J = 36     # number of points\n",
    "mu = 1    # cfl parameter\n",
    "\n",
    "# Discretization parameters\n",
    "dx = L/J\n",
    "x = np.linspace(0, L, J+1)\n",
    "x = x[0:-1]\n",
    "\n",
    "## Initialize\n",
    "U = u0(x)\n",
    "t = 0\n",
    "UTX = np.concatenate([[x.copy()],[U]],axis=0)\n",
    "LT = [t,t]\n",
    "\n",
    "F,dF = Flux(U,t,x)       # compute f(U_j^0)\n",
    "Vmax = max(np.abs(dF))   # get a bound for velocities (to adapt depending on the properties of f)\n",
    "dt = mu*dx/max(1,Vmax)   # max(1,Vmax) to avoid very large dt (bad for accuracy)\n",
    "\n",
    "y = np.concatenate(([-dx],x,[L])) + dx/2\n",
    "\n",
    "## Time iterations ##\n",
    "while t<T:\n",
    "    dt = min([dt,T-t])\n",
    "    V = np.concatenate(([U[-1]],U,[U[0]]))\n",
    "    ##########################################\n",
    "    FL,dFL = Flux(V[:-1],t,y[:-1])                       \n",
    "    FR,dFR = Flux(V[1:],t,y[:-1])\n",
    "    NumFlux = FL*(dFL>0) + FR*(dFR<0)\n",
    "    ##########################################\n",
    "    U = U - dt/dx*(NumFlux[1:]-NumFlux[:-1])\n",
    "    t = t + dt\n",
    "\n",
    "    LT = np.append(LT,t)\n",
    "    UTX = np.concatenate([UTX,[U]],axis=0)\n",
    "    \n",
    "    \n",
    "LT  = np.resize(LT,(len(LT),1))\n",
    "UTX = np.concatenate([UTX,LT],axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "SHOW(UTX,leg='Upwind flux ; $\\mu=1$ ; '+str(J)+' cells')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for J in [8*10+4,8*20+4,8*40+4,8*80+4,8*160+4]:\n",
    "\n",
    "    mu = 1    # cfl parameter\n",
    "    \n",
    "    # Discretization parameters\n",
    "    dx = L/J\n",
    "    x = np.linspace(0, L, J+1)\n",
    "    x = x[0:-1]\n",
    "\n",
    "    ## Initialize\n",
    "    U = u0(x)\n",
    "    t = 0\n",
    "    UTX = np.concatenate([[x.copy()],[U]],axis=0)\n",
    "    LT = [t,t]\n",
    "\n",
    "    F,dF = Flux(U,t,x)       # compute f(U_j^0)\n",
    "    Vmax = max(np.abs(dF))   # get a bound for velocities (to adapt depending on the properties of f)\n",
    "    dt = mu*dx/max(1,Vmax)   # max(1,Vmax) to avoid very large dt (bad for accuracy)\n",
    "\n",
    "    y = np.concatenate(([-dx],x,[L])) + dx/2\n",
    "\n",
    "    ## Time iterations ##\n",
    "    while t<T:\n",
    "        dt = min([dt,T-t])\n",
    "        V = np.concatenate(([U[-1]],U,[U[0]]))\n",
    "        ##########################################\n",
    "        FL,dFL = Flux(V[:-1],t,y[:-1])                       \n",
    "        FR,dFR = Flux(V[1:],t,y[:-1])\n",
    "        NumFlux = FL*(dFL>0) + FR*(dFR<0)\n",
    "        ##########################################\n",
    "        U = U - dt/dx*(NumFlux[1:]-NumFlux[:-1])\n",
    "        t = t + dt\n",
    "\n",
    "        LT = np.append(LT,t)\n",
    "        UTX = np.concatenate([UTX,[U]],axis=0)\n",
    "    \n",
    "    LT  = np.resize(LT,(len(LT),1))\n",
    "    UTX = np.concatenate([UTX,LT],axis=1)\n",
    "    \n",
    "    plt.plot(x,U,label='J='+str(J))\n",
    "    \n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "SHOW(UTX,leg='Upwind flux ; $\\mu=1$ ; '+str(J)+' cells')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 5.** Does the upwind scheme provide good results ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Nonlinear Burgers problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Fburgers(u):\n",
    "    return u**2/2, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def u0(x):\n",
    "    return np.array([0.5-np.sign(y-1)*np.sign(y-2) for y in x])\n",
    "plt.plot(x,u0(x),label='u0')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 0.5         # final time\n",
    "L = 4           # space length\n",
    "Flux = Fburgers # flux function\n",
    "\n",
    "mu = 0.9        # cfl parameter\n",
    "J = 100         # number of points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 6.** Program the Lax-Friedrichs scheme for the Burgers equation:\n",
    "> $$F_{j+1/2}^n = \\dfrac{1}{2}(f(U_j^n)+f(U_{j+1}^n)) - \\dfrac{\\sigma}{2}(U_{j+1}^n-U_{j}^n),$$\n",
    "> with $\\sigma = \\max\\big\\{|f'(u)|,\\ u \\in[\\inf u_0,\\sup u_0]\\big\\}.$  \n",
    "> Compare the discrete solution with $\\mu=0.9$ for several values of the grid size $J$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **Godunov scheme** is defined from the exact solution of local Riemann problems.  \n",
    "The discrete flux is then obtained as $$F_{j+1/2}^n = f(w(0,U_j^n,U_{j+1}^n))$$ where $w(\\xi,U,V)$ denotes the self-similar (and entropy weak) solution to the Riemann problem with left state $U$ and right state $V$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 7.** Solve the Riemann problem for the Burgers flux and then program the Godunov scheme for the Burgers problem. Compare the solution to the Lax-Friedrichs scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 8.** Program the local Lax-Friedrichs scheme for the Burgers equation:\n",
    "> $$F_{j+1/2}^n = \\dfrac{1}{2}(f(U_j^n)+f(U_{j+1}^n)) - \\dfrac{\\sigma_j^n}{2}(U_{j+1}^n-U_{j}^n),$$\n",
    "> with $\\sigma_j^n = \\max\\big\\{|f'(u)|,\\ u \\in[\\min(U_j^n,U_{j+1}^n),\\max(U_j^n,U_{j+1}^n)]\\big\\}.$  \n",
    "> Compare the solution to the previous schemes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question 9.** Solve numerically the Cauchy problem for the nonlinear equations with flux $f(u)= 2 u^3 -u$ and initial data $u_0$ given below over $x\\in[0,4]$ and $t\\in[0,1.25]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FLUX(u):\n",
    "    return 2*u**3-u, 6*u**2-1\n",
    "up = np.linspace(-1,1)\n",
    "plt.plot(up,FLUX(up)[0])\n",
    "plt.show()\n",
    "\n",
    "\n",
    "def u0(x):\n",
    "    return -1*(x<1)-0.75*(x>=1)*(x<2)+0.25*(x>=2)*(x<3)-0.5*(x>=3)\n",
    "x = np.linspace(0,4)\n",
    "plt.plot(x,u0(x),label='u0')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}