{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "The famous *diffusion equation*, also known as the *heat equation*,\n", "reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial t} =\n", "\\dfc \\frac{\\partial^2 u}{\\partial x^2},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $u(x,t)$ is the unknown function to be solved for, $x$ is a coordinate\n", "in space, and $t$ is time. The coefficient $\\dfc$ is the *diffusion\n", "coefficient* and determines how fast $u$ changes in time. A quick\n", "short form for the diffusion equation is $u_t = \\dfc u_{xx}$.\n", "\n", "Compared to the wave equation, $u_{tt}=c^2u_{xx}$, which looks very similar,\n", "the diffusion equation features solutions that are very different from\n", "those of the wave equation. Also, the diffusion equation\n", "makes quite different demands to the numerical\n", "methods.\n", "\n", "\n", "Typical diffusion problems may experience rapid change in the very\n", "beginning, but then the evolution of $u$ becomes slower and slower.\n", "The solution is usually very smooth, and after some time, one cannot\n", "recognize the initial shape of $u$. This is in sharp contrast to\n", "solutions of the wave equation where the initial shape is preserved in\n", "homogeneous media - the solution is then basically a moving initial\n", "condition. The standard wave equation $u_{tt}=c^2u_{xx}$ has solutions\n", "that propagate with speed $c$ forever, without changing shape, while\n", "the diffusion equation converges to a *stationary solution* $\\bar\n", "u(x)$ as $t\\rightarrow\\infty$. In this limit, $u_t=0$, and $\\bar u$ is\n", "governed by $\\bar u''(x)=0$. This stationary limit of the diffusion\n", "equation is called the *Laplace* equation and arises in a very wide\n", "range of applications throughout the sciences.\n", "\n", "mathcal{I}_t is possible to solve for $u(x,t)$ using an explicit scheme, as we\n", "do in the section [An explicit method for the 1D diffusion equation](#diffu:pde1:FEsec), but the time step restrictions\n", "soon become much less favorable than for an explicit scheme applied to\n", "the wave equation. And of more importance, since the solution $u$ of\n", "the diffusion equation is very smooth and changes slowly, small time\n", "steps are not convenient and not required by accuracy as the diffusion\n", "process converges to a stationary state. Therefore, implicit schemes\n", "(as described in the section [Implicit methods for the 1D diffusion equation](#diffu:pde1:implicit)) are popular, but\n", "these require solutions of systems of algebraic equations. We shall\n", "use ready-made software for this purpose, but also program some simple\n", "iterative methods.\n", "The exposition is, as usual in this book, very basic and focuses on\n", "the basic ideas and how to implement. More comprehensive mathematical\n", "treatments and classical analysis\n", "of the methods are found in lots of textbooks. A favorite\n", "of ours in this respect is the one by LeVeque [[LeVeque_2007]](#LeVeque_2007).\n", "The books by Strikwerda [[Strikwerda_2007]](#Strikwerda_2007) and by\n", "Lapidus and Pinder [[Lapidus_Pinder_1982]](#Lapidus_Pinder_1982) are also highly recommended\n", "as additional material on the topic.\n", "\n", "\n", "# An explicit method for the 1D diffusion equation\n", "
\n", "\n", "Explicit finite difference methods for the wave equation $u_{tt}=c^2u_{xx}$\n", "can be used, with small modifications, for solving $u_t = \\dfc u_{xx}$\n", "as well.\n", "% if BOOK == \"book\":\n", "The exposition below assumes that the reader is familiar with the\n", "basic ideas of discretization and implementation of wave\n", "equations from the chapter [ch:wave](#ch:wave). Readers not familiar with the\n", "Forward Euler, Backward Euler, and Crank-Nicolson (or centered or\n", "midpoint) discretization methods in time should consult, e.g., Section 1.1 \n", "in [[Langtangen_decay]](#Langtangen_decay).\n", "% endif\n", "\n", "## The initial-boundary value problem for 1D diffusion\n", "\n", "To obtain a unique solution of the diffusion equation, or equivalently,\n", "to apply numerical methods, we need initial and boundary conditions.\n", "The diffusion equation goes with one initial condition $u(x,0)=I(x)$, where\n", "$I$ is a prescribed function. One boundary condition is required at\n", "each point on the boundary, which in 1D means that $u$ must be known,\n", "$u_x$ must be known, or some combination of them.\n", "\n", "\n", "We shall start with the simplest boundary condition: $u=0$. The\n", "complete initial-boundary value diffusion problem in one space\n", "dimension can then be specified as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} =\n", "\\dfc \\frac{\\partial^2 u}{\\partial x^2} + f, \\quad x\\in (0,L),\\ t\\in (0,T]\n", "\\label{diffu:pde1} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(x,0) = I(x), \\quad x\\in [0,L]\n", "\\label{diffu:pde1:ic:u} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(0,t) = 0, \\quad t>0,\n", "\\label{diffu:pde1:bc:0} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(L,t) = 0, \\quad t>0\\thinspace .\n", "\\label{diffu:pde1:bc:L} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With only a first-order derivative in time,\n", "only one *initial condition* is needed, while the second-order\n", "derivative in space leads to a demand for two *boundary conditions*.\n", "We have added a source term $f=f(x,t)$, which is\n", "convenient when testing implementations.\n", "\n", "\n", "Diffusion equations like ([1](#diffu:pde1)) have a wide range of\n", "applications throughout physical, biological, and financial sciences.\n", "One of the most common applications is propagation of heat, where\n", "$u(x,t)$ represents the temperature of some substance at point $x$ and\n", "time $t$. Other applications are listed in the section [diffu:app](#diffu:app).\n", "\n", "\n", "## Forward Euler scheme\n", "
\n", "\n", "The first step in the discretization procedure is to replace the\n", "domain $[0,L]\\times [0,T]$ by a set of mesh points. Here we apply\n", "equally spaced mesh points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_i=i\\Delta x,\\quad i=0,\\ldots,N_x,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "t_n=n\\Delta t,\\quad n=0,\\ldots,N_t \\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Moreover, $u^n_i$ denotes the mesh function that\n", "approximates $u(x_i,t_n)$ for $i=0,\\ldots,N_x$ and $n=0,\\ldots,N_t$.\n", "Requiring the PDE ([1](#diffu:pde1)) to be fulfilled at a mesh point $(x_i,t_n)$\n", "leads to the equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial}{\\partial t} u(x_i, t_n) =\n", "\\dfc\\frac{\\partial^2}{\\partial x^2} u(x_i, t_n) + f(x_i,t_n),\n", "\\label{diffu:pde1:step2} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next step is to replace the derivatives by finite difference approximations.\n", "The computationally simplest method arises from\n", "using a forward difference in time and a central difference in\n", "space:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "[D_t^+ u = \\dfc D_xD_x u + f]^n_i \\thinspace .\n", "\\label{diffu:pde1:step3a} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Written out," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{u^{n+1}_i-u^n_i}{\\Delta t} = \\dfc \\frac{u^{n}_{i+1} - 2u^n_i + u^n_{i-1}}{\\Delta x^2} + f_i^n\\thinspace .\n", "\\label{diffu:pde1:step3b} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have turned the PDE into algebraic equations, also often called\n", "discrete equations. The key property of the equations is that they\n", "are algebraic, which makes them easy to solve.\n", "As usual, we anticipate that $u^n_i$ is already computed such that\n", "$u^{n+1}_i$ is the only unknown in ([7](#diffu:pde1:step3b)).\n", "Solving with respect to this unknown is easy:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u^{n+1}_i = u^n_i + F\\left(\n", "u^{n}_{i+1} - 2u^n_i + u^n_{i-1}\\right) + \\Delta t f_i^n,\n", "\\label{diffu:pde1:step4} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where we have introduced the *mesh Fourier number*:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "F = \\dfc\\frac{\\Delta t}{\\Delta x^2}\\thinspace .\n", "\\label{_auto1} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**$F$ is the key parameter in the discrete diffusion equation.**\n", "\n", "Note that $F$ is a *dimensionless* number that lumps the key physical\n", "parameter in the problem, $\\dfc$, and the discretization parameters\n", "$\\Delta x$ and $\\Delta t$ into a single parameter. Properties\n", "of the numerical method are critically dependent upon the value of\n", "$F$ (see the section [diffu:pde1:analysis](#diffu:pde1:analysis) for details).\n", "\n", "\n", "\n", "The computational algorithm then becomes\n", "\n", "1. compute $u^0_i=I(x_i)$ for $i=0,\\ldots,N_x$\n", "\n", "2. for $n=0,1,\\ldots,N_t$:\n", "\n", "a. apply ([8](#diffu:pde1:step4)) for all the internal\n", " spatial points $i=1,\\ldots,N_x-1$\n", "\n", "b. set the boundary values\n", " $u^{n+1}_i=0$ for $i=0$ and $i=N_x$\n", "\n", "\n", "The algorithm is compactly and fully specified in Python:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'L' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mL\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNx\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# mesh points in space\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0mdx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNt\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# mesh points in time\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mdt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'L' is not defined" ] } ], "source": [ "import numpy as np\n", "\n", "I = lambda x: 1\n", "Nx = 100000\n", "a = 2.0\n", "L = 2.0\n", "dx = L/Nx\n", "dt = dx**2/(2*a)\n", "T = 100*dt\n", "Nt = int(round(T/float(dt)))\n", "\n", "\n", "x = np.linspace(0, L, Nx+1) # mesh points in space\n", "dx = x[1] - x[0]\n", "t = np.linspace(0, T, Nt+1) # mesh points in time\n", "dt = t[1] - t[0]\n", "F = a*dt/dx**2\n", "u = np.zeros(Nx+1) # unknown u at new time level\n", "u_n = np.zeros(Nx+1) # u at the previous time level\n", "\n", "# Set initial condition u(x,0) = I(x)\n", "for i in range(0, Nx+1):\n", " u_n[i] = I(x[i])\n", "\n", "for n in range(0, Nt):\n", " # Compute u at inner mesh points\n", " for i in range(1, Nx):\n", " u[i] = u_n[i] + F*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + \\\n", "\t dt*f(x[i], t[n])\n", "\n", " # Insert boundary conditions\n", " u[0] = 0; u[Nx] = 0\n", "\n", " # Update u_n before next step\n", " u_n[:]= u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we use `a` for $\\dfc$ in the code, motivated by easy visual\n", "mapping between the variable name and the mathematical symbol in formulas.\n", "\n", "We need to state already now that the shown algorithm does not\n", "produce meaningful results unless $F\\leq 1/2$. Why is explained in\n", "the section [diffu:pde1:analysis](#diffu:pde1:analysis).\n", "\n", "## Implementation\n", "
\n", "\n", "The file [`diffu1D_u0.py`](${src_diffu}/diffu1D_u0.py)\n", "contains a complete function `solver_FE_simple`\n", "for solving the 1D diffusion equation with $u=0$ on the boundary\n", "as specified in the algorithm above:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "def solver_FE_simple(I, a, f, L, dt, F, T):\n", " \"\"\"\n", " Simplest expression of the computational algorithm\n", " using the Forward Euler method and explicit Python loops.\n", " For this method F <= 0.5 for stability.\n", " \"\"\"\n", " import time; t0 = time.clock() # For measuring the CPU time\n", "\n", " Nt = int(round(T/float(dt)))\n", " t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n", " dx = np.sqrt(a*dt/F)\n", " Nx = int(round(L/dx))\n", " x = np.linspace(0, L, Nx+1) # Mesh points in space\n", " # Make sure dx and dt are compatible with x and t\n", " dx = x[1] - x[0]\n", " dt = t[1] - t[0]\n", "\n", " u = np.zeros(Nx+1)\n", " u_n = np.zeros(Nx+1)\n", "\n", " # Set initial condition u(x,0) = I(x)\n", " for i in range(0, Nx+1):\n", " u_n[i] = I(x[i])\n", "\n", " for n in range(0, Nt):\n", " # Compute u at inner mesh points\n", " for i in range(1, Nx):\n", " u[i] = u_n[i] + F*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + \\\n", " dt*f(x[i], t[n])\n", "\n", " # Insert boundary conditions\n", " u[0] = 0; u[Nx] = 0\n", "\n", " # Switch variables before next step\n", " #u_n[:] = u # safe, but slow\n", " u_n, u = u, u_n\n", "\n", " t1 = time.clock()\n", " return u_n, x, t, t1-t0 # u_n holds latest u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A faster alternative is available in the function `solver_FE`, which \n", "adds the possibility of solving the finite difference scheme by vectorization.\n", "The vectorized version replaces the explicit loop" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "for i in range(1, Nx):\n", " u[i] = u_n[i] + F*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) \\\n", " + dt*f(x[i], t[n])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "by arithmetics on displaced slices of the `u` array:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "u[1:Nx] = u_n[1:Nx] + F*(u_n[0:Nx-1] - 2*u_n[1:Nx] + u_n[2:Nx+1]) \\\n", " + dt*f(x[1:Nx], t[n])\n", "# or\n", "u[1:-1] = u_n[1:-1] + F*(u_n[0:-2] - 2*u_n[1:-1] + u_n[2:]) \\\n", " + dt*f(x[1:-1], t[n])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example,\n", "the vectorized version runs 70 times faster than the scalar version\n", "in a case with 100 time steps and a spatial mesh of $10^5$ cells.\n", "\n", "The `solver_FE` function also features a callback function such that the\n", "user can process the solution at each time level. The callback\n", "function looks like `user_action(u, x, t, n)`, where `u` is the array\n", "containing the solution at time level `n`, `x` holds all the\n", "spatial mesh points, while `t` holds all the temporal mesh points.\n", "The `solver_FE` function is very similar to `solver_FE_simple` above:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def solver_FE(I, a, f, L, dt, F, T,\n", " user_action=None, version='scalar'):\n", " \"\"\"\n", " Vectorized implementation of solver_FE_simple.\n", " \"\"\"\n", " import time; t0 = time.clock() # for measuring the CPU time\n", "\n", " Nt = int(round(T/float(dt)))\n", " t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n", " dx = np.sqrt(a*dt/F)\n", " Nx = int(round(L/dx))\n", " x = np.linspace(0, L, Nx+1) # Mesh points in space\n", " # Make sure dx and dt are compatible with x and t\n", " dx = x[1] - x[0]\n", " dt = t[1] - t[0]\n", "\n", " u = np.zeros(Nx+1) # solution array\n", " u_n = np.zeros(Nx+1) # solution at t-dt\n", "\n", " # Set initial condition\n", " for i in range(0,Nx+1):\n", " u_n[i] = I(x[i])\n", "\n", " if user_action is not None:\n", " user_action(u_n, x, t, 0)\n", "\n", " for n in range(0, Nt):\n", " # Update all inner points\n", " if version == 'scalar':\n", " for i in range(1, Nx):\n", " u[i] = u_n[i] +\\\n", " F*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) +\\\n", " dt*f(x[i], t[n])\n", "\n", " elif version == 'vectorized':\n", " u[1:Nx] = u_n[1:Nx] + \\\n", " F*(u_n[0:Nx-1] - 2*u_n[1:Nx] + u_n[2:Nx+1]) +\\\n", " dt*f(x[1:Nx], t[n])\n", " else:\n", " raise ValueError('version=%s' % version)\n", "\n", " # Insert boundary conditions\n", " u[0] = 0; u[Nx] = 0\n", " if user_action is not None:\n", " user_action(u, x, t, n+1)\n", "\n", " # Switch variables before next step\n", " u_n, u = u, u_n\n", "\n", " t1 = time.clock()\n", " return t1-t0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Verification\n", "
\n", "\n", "### Exact solution of discrete equations\n", "\n", "
\n", "\n", "Before thinking about running the functions in the previous section,\n", "we need to construct a suitable test example for verification. mathcal{I}_t\n", "appears that a manufactured solution that is linear in time and at\n", "most quadratic in space fulfills the Forward Euler scheme\n", "exactly. With the restriction that $u=0$ for $x=0,L$, we can try the\n", "solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u(x,t) = 5tx(L-x)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Inserted in the PDE, it requires a source term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x,t) = 10\\dfc t + 5x(L-x)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "% if BOOK == 'book':\n", "With the formulas from [sec:form:fdtn](#sec:form:fdtn) we can easily check\n", "% else:\n", "Let us check\n", "% endif\n", "that the manufactured `u` fulfills the scheme:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\lbrack D_t^+ u = \\dfc D_x D_x u + f\\rbrack^n_i &=\n", "\\lbrack 5x(L-x)D_t^+ t = 5 t\\dfc D_x D_x (xL-x^2) +\\\\ \n", "&\\quad\\quad 10\\dfc t + 5x(L-x)\\rbrack^n_i\\\\ \n", "&=\n", "\\lbrack 5x(L-x) = 5 t\\dfc (-2) + 10\\dfc t + 5x(L-x) \\rbrack^n_i,\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is a 0=0 expression.\n", "The computation of the source term, given any $u$,\n", "is easily automated with `sympy`:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle 10 a t + 5 x \\left(L - x\\right)$" ], "text/plain": [ "10*a*t + 5*x*(L - x)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import sympy as sym\n", "x, t, a, L = sym.symbols('x t a L')\n", "u = x*(L-x)*5*t\n", "\n", "def pde(u):\n", " return sym.diff(u, t) - a*sym.diff(u, x, x)\n", "\n", "f = sym.simplify(pde(u))\n", "f" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can choose any expression for `u` and automatically\n", "get the suitable source term `f`. However, the manufactured solution\n", "`u` will in general\n", "not be exactly reproduced by the scheme: only constant and linear\n", "functions are differentiated correctly by a forward difference, while only\n", "constant, linear, and quadratic functions are differentiated exactly by\n", "a $[D_xD_x u]^n_i$ difference.\n", "\n", "The numerical code will need to access the `u` and `f` above\n", "as Python functions. The exact solution is wanted as a Python\n", "function `u_exact(x, t)`, while the source term is wanted as\n", "`f(x, t)`. The parameters `a` and `L` in `u` and `f` above\n", "are symbols and must be replaced by `float` objects in a Python\n", "function. This can be done by redefining `a` and `L` as\n", "`float` objects and performing substitutions of symbols by\n", "numbers in `u` and `f`. The appropriate code looks like this:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "a = 0.5\n", "L = 1.5\n", "u_exact = sym.lambdify(\n", " [x, t], u.subs('L', L).subs('a', a), modules='numpy')\n", "f = sym.lambdify(\n", " [x, t], f.subs('L', L).subs('a', a), modules='numpy')\n", "I = lambda x: u_exact(x, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we also make a function `I` for the initial condition.\n", "\n", "The idea now is that our manufactured solution should be\n", "exactly reproduced by the code (to machine precision).\n", "For this purpose we make a test function for comparing\n", "the exact and numerical solutions at the end of the\n", "time interval:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def test_solver_FE():\n", " # Define u_exact, f, I as explained above\n", "\n", " dx = L/3 # 3 cells\n", " F = 0.5\n", " dt = F*dx**2\n", "\n", " u, x, t, cpu = solver_FE_simple(\n", " I=I, a=a, f=f, L=L, dt=dt, F=F, T=2)\n", " u_e = u_exact(x, t[-1])\n", " diff = abs(u_e - u).max()\n", " tol = 1E-14\n", " assert diff < tol, 'max diff solver_FE_simple: %g' % diff\n", "\n", " u, x, t, cpu = solver_FE(\n", " I=I, a=a, f=f, L=L, dt=dt, F=F, T=2,\n", " user_action=None, version='scalar')\n", " u_e = u_exact(x, t[-1])\n", " diff = abs(u_e - u).max()\n", " tol = 1E-14\n", " assert diff < tol, 'max diff solver_FE, scalar: %g' % diff\n", "\n", " u, x, t, cpu = solver_FE(\n", " I=I, a=a, f=f, L=L, dt=dt, F=F, T=2,\n", " user_action=None, version='vectorized')\n", " u_e = u_exact(x, t[-1])\n", " diff = abs(u_e - u).max()\n", " tol = 1E-14\n", " assert diff < tol, 'max diff solver_FE, vectorized: %g' % diff" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**The critical value $F=0.5$.**\n", "\n", "We emphasize that the value `F=0.5` is critical: the tests above\n", "will fail if `F` has a larger value. This is because the Forward\n", "Euler scheme is unstable for $F>1/2$.\n", "\n", "The reader may wonder if\n", "$F=1/2$ is safe or if $F<1/2$ should be required. Experiments show\n", "that $F=1/2$ works fine for $u_t=\\dfc u_{xx}$, so\n", "there is no accumulation of rounding\n", "errors in this case and hence no need to introduce any safety factor\n", "to keep $F$ away from the limiting value 0.5.\n", "\n", "\n", "\n", "\n", "### Checking convergence rates\n", "\n", "
\n", "\n", "\n", "If our chosen exact solution does not satisfy the discrete equations\n", "exactly, we are left with checking the convergence rates, just as we did\n", "previously for the wave equation. However, with the Euler scheme here,\n", "we have different accuracies in time and space, since we use a second\n", "order approximation to the spatial derivative and a first order approximation\n", "to the time derivative. Thus, we must expect different convergence rates in\n", "time and space. For the numerical error," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E = C_t\\Delta t^r + C_x\\Delta x^p,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "we should get convergence rates $r=1$ and $p=2$ ($C_t$ and $C_x$ are unknown constants).\n", "As previously,\n", "in the section [wave:pde2:fd:MMS](#wave:pde2:fd:MMS),\n", "we simplify matters by introducing a single discretization parameter $h$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "h = \\Delta t,\\quad \\Delta x = Kh^{r/p},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $K$ is any constant. This allows us to factor out only *one*\n", "discretization parameter $h$ from the formula:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E = C_t h + C_x (Kh^{r/p})^p = \\tilde C h^r,\\quad\n", "\\tilde C = C_t + C_sK^r\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The computed rate $r$ should approach 1 with increasing resolution.\n", "\n", "mathcal{I}_t is tempting, for simplicity,\n", "to choose $K=1$, which gives $\\Delta x = h^{r/p}$, expected to be\n", "$\\sqrt{\\Delta t}$. However,\n", "we have to control the stability requirement: $F\\leq\\frac{1}{2}$,\n", "which means" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\dfc\\Delta t}{\\Delta x^2}\\leq\\frac{1}{2}\\quad\\Rightarrow\n", "\\quad \\Delta x \\geq \\sqrt{2\\dfc}h^{1/2} ,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "implying that $K=\\sqrt{2\\dfc}$ is our choice in experiments where we\n", "lie on the stability limit $F=1/2$.\n", "\n", "\n", "## Numerical experiments\n", "
\n", "\n", "When a test function like the one above runs silently without errors,\n", "we have some evidence for a correct implementation of the numerical\n", "method. The next step is to do some experiments with more interesting\n", "solutions.\n", "\n", "We target a scaled diffusion problem where $x/L$ is a new spatial\n", "coordinate and $\\dfc t/L^2$ is a new time coordinate. The source term\n", "$f$ is omitted, and $u$ is scaled by $\\max_{x\\in [0,L]}|I(x)|$ (see Section 3.2 in\n", " [[Langtangen_scaling]](#Langtangen_scaling) for details).\n", "The governing PDE is then" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial t} = \\frac{\\partial^2 u}{\\partial x^2},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "in the spatial domain $[0,L]$, with boundary conditions $u(0)=u(1)=0$.\n", "Two initial conditions will be tested: a discontinuous plug," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "I(x) = \\left\\lbrace\\begin{array}{ll}\n", "0, & |x-L/2| > 0.1\\\\ \n", "1, & \\hbox{otherwise}\n", "\\end{array}\\right.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and a smooth Gaussian function," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "I(x) = e^{-\\frac{1}{2\\sigma^2}(x-L/2)^2}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The functions `plug` and `gaussian` in [`diffu1D_u0.py`](${src_diffu}/diffu1D_u0.py) run the two cases,\n", "respectively:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def plug(scheme='FE', F=0.5, Nx=50):\n", " L = 1.\n", " a = 1.\n", " T = 0.1\n", " # Compute dt from Nx and F\n", " dx = L/Nx; dt = F/a*dx**2\n", "\n", " def I(x):\n", " \"\"\"Plug profile as initial condition.\"\"\"\n", " if abs(x-L/2.0) > 0.1:\n", " return 0\n", " else:\n", " return 1\n", "\n", " cpu = viz(I, a, L, dt, F, T,\n", " umin=-0.1, umax=1.1,\n", " scheme=scheme, animate=True, framefiles=True)\n", " print('CPU time:', cpu)\n", "\n", "def gaussian(scheme='FE', F=0.5, Nx=50, sigma=0.05):\n", " L = 1.\n", " a = 1.\n", " T = 0.1\n", " # Compute dt from Nx and F\n", " dx = L/Nx; dt = F/a*dx**2\n", "\n", " def I(x):\n", " \"\"\"Gaussian profile as initial condition.\"\"\"\n", " return exp(-0.5*((x-L/2.0)**2)/sigma**2)\n", "\n", " u, cpu = viz(I, a, L, dt, F, T,\n", " umin=-0.1, umax=1.1,\n", " scheme=scheme, animate=True, framefiles=True)\n", " print('CPU time:', cpu)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These functions make use of the function `viz` for running the\n", "solver and visualizing the solution using a callback function\n", "with plotting:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def viz(I, a, L, dt, F, T, umin, umax,\n", " scheme='FE', animate=True, framefiles=True):\n", "\n", " def plot_u(u, x, t, n):\n", " plt.plot(x, u, 'r-', axis=[0, L, umin, umax],\n", " title='t=%f' % t[n])\n", " if framefiles:\n", " plt.savefig('tmp_frame%04d.png' % n)\n", " if t[n] == 0:\n", " time.sleep(2)\n", " elif not framefiles:\n", " # mathcal{I}_t takes time to write files so pause is needed\n", " # for screen only animation\n", " time.sleep(0.2)\n", "\n", " user_action = plot_u if animate else lambda u,x,t,n: None\n", "\n", " cpu = eval('solver_'+scheme)(I, a, L, dt, F, T,\n", " user_action=user_action)\n", " return cpu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that this `viz` function stores all the solutions in a\n", "list `solutions` in the callback function. Modern computers have\n", "hardly any problem with storing a lot of such solutions for moderate\n", "values of $N_x$ in 1D problems, but for 2D and 3D problems, this\n", "technique cannot be used and solutions must be stored in files.\n", "\n", "[hpl 1: Better to show the scalable file solution here?]\n", "\n", "Our experiments employ a time step $\\Delta t = 0.0002$ and\n", "simulate for $t\\in [0,0.1]$. First we try the highest value of\n", "$F$: $F=0.5$. This resolution corresponds to\n", "$N_x=50$. A possible terminal command is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Terminal> python -c 'from diffu1D_u0 import gaussian\n", " gaussian(\"solver_FE\", F=0.5, dt=0.0002)'\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $u(x,t)$ curve as a function of $x$ is shown in [Figure](#diffu:pde1:FE:fig:F=0.5) at four time levels.\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "\n", "
\n", "

\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import HTML\n", "_s = \"\"\"\n", "
\n", "\n", "
\n", "

\n", "\n", "\n", "\n", "\n", "\"\"\"\n", "HTML(_s)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "We see that the curves have saw-tooth waves in the beginning of the\n", "simulation. This non-physical noise is smoothed out with time, but\n", "solutions of the diffusion equations are known to be smooth, and\n", "this numerical solution is definitely not smooth.\n", "Lowering $F$ helps: $F\\leq 0.25$ gives a smooth solution, see\n", "% if FORMAT == \"pdflatex\":\n", "[Figure](#diffu:pde1:FE:fig:F=0.25) (and a\n", "[movie](${docraw}/mov-diffu/diffu1D_u0_FE_plug_F025/movie.ogg)).\n", "% else:\n", "[Figure](#diffu:pde1:FE:fig:F=0.25).\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "_s = \"\"\"\n", "
\n", "\n", "
\n", "

\n", "\n", "\n", "\n", "\n", "\"\"\"\n", "HTML(_s)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "% endif\n", "\n", "Increasing $F$ slightly beyond the limit 0.5, to $F=0.51$,\n", "gives growing, non-physical instabilities,\n", "as seen in [Figure](#diffu:pde1:FE:fig:F=0.51).\n", "\n", "\n", "\n", "
\n", "\n", "

Forward Euler scheme for $F=0.5$.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "

Forward Euler scheme for $F=0.25$.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "

Forward Euler scheme for $F=0.51$.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "Instead of a discontinuous initial condition we now try the smooth\n", "Gaussian function for $I(x)$. A simulation for $F=0.5$\n", "is shown in [Figure](#diffu:pde1:FE:fig:gauss:F=0.5). Now the numerical solution\n", "is smooth for all times, and this is true for any $F\\leq 0.5$.\n", "\n", "% if FORMAT != \"pdflatex\":\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "_s = \"\"\"\n", "
\n", "\n", "
\n", "

\n", "\n", "\n", "\n", "\n", "\"\"\"\n", "HTML(_s)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "% endif\n", "\n", "\n", "\n", "
\n", "\n", "

Forward Euler scheme for $F=0.5$.

\n", "\n", "\n", "\n", "\n", "\n", "Experiments with these two choices of $I(x)$ reveal some\n", "important observations:\n", "\n", " * The Forward Euler scheme leads to growing solutions if $F>\\frac{1}{2}$.\n", "\n", " * $I(x)$ as a discontinuous plug leads to a saw tooth-like noise\n", " for $F=\\frac{1}{2}$, which is absent for $F\\leq\\frac{1}{4}$.\n", "\n", " * The smooth Gaussian initial function leads to a smooth solution\n", " for all relevant $F$ values ($F\\leq \\frac{1}{2}$).\n", "\n", "# Implicit methods for the 1D diffusion equation\n", "
\n", "\n", "Simulations with the Forward Euler scheme show that the time step\n", "restriction, $F\\leq\\frac{1}{2}$, which means $\\Delta t \\leq \\Delta x^2/(2\\dfc)$,\n", "may be relevant in the beginning of the diffusion process, when the\n", "solution changes quite fast, but as time increases, the process slows\n", "down, and a small $\\Delta t$ may be inconvenient. With \n", "*implicit schemes*, which lead to coupled systems of linear equations\n", "to be solved at each time level, any size of $\\Delta t$ is possible\n", "(but the accuracy decreases with increasing $\\Delta t$).\n", "The Backward Euler scheme, derived and implemented below, is the\n", "simplest implicit scheme for the diffusion equation.\n", "\n", "## Backward Euler scheme\n", "
\n", "\n", "In ([5](#diffu:pde1:step2)), we now apply a backward difference in time,\n", "but the same central difference in space:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "[D_t^- u = D_xD_x u + f]^n_i,\n", "\\label{diffu:pde1:step3aBE} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which written out reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{u^{n}_i-u^{n-1}_i}{\\Delta t} = \\dfc\\frac{u^{n}_{i+1} - 2u^n_i + u^n_{i-1}}{\\Delta x^2} + f_i^n\\thinspace .\n", "\\label{diffu:pde1:step3bBE} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we assume $u^{n-1}_i$ is already computed, but that all quantities at the \"new\"\n", "time level $n$ are unknown. This time it is not possible to solve\n", "with respect to $u_i^{n}$ because this value couples to its neighbors\n", "in space, $u^n_{i-1}$ and $u^n_{i+1}$, which are also unknown.\n", "Let us examine this fact for the case when $N_x=3$. Equation ([11](#diffu:pde1:step3bBE)) written for $i=1,\\ldots,Nx-1= 1,2$ becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{u^{n}_1-u^{n-1}_1}{\\Delta t} = \\dfc\\frac{u^{n}_{2} - 2u^n_1 + u^n_{0}}{\\Delta x^2} + f_1^n\n", "\\label{_auto2} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{u^{n}_2-u^{n-1}_2}{\\Delta t} = \\dfc\\frac{u^{n}_{3} - 2u^n_2 + u^n_{1}}{\\Delta x^2} + f_2^n\n", "\\label{_auto3} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The boundary values $u^n_0$ and $u^n_3$ are known as zero. Collecting the\n", "unknown new values $u^n_1$ and $u^n_2$ on the left-hand side and multiplying\n", "by $\\Delta t$ gives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\left(1+ 2F\\right) u^{n}_1 - F u^{n}_{2} = u^{n-1}_1 + \\Delta t f_1^n,\n", "\\label{_auto4} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "- F u^{n}_{1} + \\left(1+ 2F\\right) u^{n}_2 = u^{n-1}_2 + \\Delta t f_2^n\\thinspace .\n", "\\label{_auto5} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a coupled $2\\times 2$ system of algebraic equations for\n", "the unknowns $u^n_1$ and $u^n_2$. The equivalent matrix form is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\left(\\begin{array}{cc}\n", "1+ 2F & - F\\\\ \n", "- F & 1+ 2F\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c}\n", "u^{n}_1\\\\ \n", "u^{n}_2\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c}\n", "u^{n-1}_1 + \\Delta t f_1^n\\\\ \n", "u^{n-1}_2 + \\Delta t f_2^n\n", "\\end{array}\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Terminology: implicit vs. explicit methods.**\n", "\n", "Discretization methods that lead to a coupled system of equations\n", "for the unknown function at a new time level are said to be\n", "*implicit methods*.\n", "The counterpart, *explicit methods*, refers to discretization\n", "methods where there is a simple explicit formula for the values of\n", "the unknown function at each of the spatial mesh points at the new\n", "time level. From an implementational point of view, implicit methods\n", "are more comprehensive to code since they require\n", "the solution of coupled equations, i.e., a matrix system, at each time level.\n", "With explicit methods we have a closed-form formula for the value of\n", "the unknown at each mesh point.\n", "\n", "Very often explicit schemes have a restriction on the size of the time\n", "step that can be relaxed by using implicit schemes. In fact,\n", "implicit schemes are frequently unconditionally stable, so the size of the\n", "time step is governed by accuracy and not by stability. This is the great\n", "advantage of implicit schemes.\n", "\n", "\n", "\n", "\n", "In the general case, ([11](#diffu:pde1:step3bBE)) gives rise to\n", "a coupled $(N_x-1)\\times (N_x-1)$ system of algebraic equations for\n", "all the unknown $u^n_i$ at the interior spatial points $i=1,\\ldots,N_x-1$.\n", "Collecting the unknowns on the left-hand side,\n", "([11](#diffu:pde1:step3bBE)) can be written" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "- F u^n_{i-1} + \\left(1+ 2F \\right) u^{n}_i - F u^n_{i+1} =\n", "u_{i-1}^{n-1},\n", "\\label{diffu:pde1:step4BE} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "for $i=1,\\ldots,N_x-1$.\n", "One can either view these equations as a system where the\n", "$u^{n}_i$ values at the internal mesh points, $i=1,\\ldots,N_x-1$, are\n", "unknown, or we may append the boundary values $u^n_0$ and $u^n_{N_x}$\n", "to the system. In the latter case, all $u^n_i$ for $i=0,\\ldots,N_x$\n", "are considered unknown, and we must add the boundary equations to\n", "the $N_x-1$ equations in ([16](#diffu:pde1:step4BE)):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_0^n = 0,\\label{diffu:pde1:step4BE:BC:0} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u_{N_x}^n = 0\\thinspace .\n", "\\label{diffu:pde1:step4BE:BC:L} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A coupled system of algebraic equations can be written on matrix form,\n", "and this is important if we want to call up ready-made software for\n", "solving the system. The equations ([16](#diffu:pde1:step4BE))\n", "and ([17](#diffu:pde1:step4BE:BC:0))--([18](#diffu:pde1:step4BE:BC:L))\n", "correspond to the matrix equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "AU = b\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $U=(u^n_0,\\ldots,u^n_{N_x})$, and\n", "the matrix $A$ has the following structure:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "A =\n", "\\left(\n", "\\begin{array}{cccccccccc}\n", "A_{0,0} & A_{0,1} & 0\n", "&\\cdots &\n", "\\cdots & \\cdots & \\cdots &\n", "\\cdots & 0 \\\\ \n", "A_{1,0} & A_{1,1} & A_{1,2} & \\ddots & & & & & \\vdots \\\\ \n", "0 & A_{2,1} & A_{2,2} & A_{2,3} &\n", "\\ddots & & & & \\vdots \\\\ \n", "\\vdots & \\ddots & & \\ddots & \\ddots & 0 & & & \\vdots \\\\ \n", "\\vdots & & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\ \n", "\\vdots & & & 0 & A_{i,i-1} & A_{i,i} & A_{i,i+1} & \\ddots & \\vdots \\\\ \n", "\\vdots & & & & \\ddots & \\ddots & \\ddots &\\ddots & 0 \\\\ \n", "\\vdots & & & & &\\ddots & \\ddots &\\ddots & A_{N_x-1,N_x} \\\\ \n", "0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots & 0 & A_{N_x,N_x-1} & A_{N_x,N_x}\n", "\\end{array}\n", "\\right)\n", "\\label{diffu:pde1:matrix:sparsity} \\tag{19}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The nonzero elements are given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "A_{i,i-1} = -F\n", "\\label{_auto6} \\tag{20}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{i,i} = 1+ 2F\n", "\\label{_auto7} \\tag{21}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{i,i+1} = -F\n", "\\label{_auto8} \\tag{22}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "in the equations for internal points, $i=1,\\ldots,N_x-1$. The first and last\n", "equation correspond to the boundary condition, where we know the solution,\n", "and therefore we must have" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "A_{0,0} = 1,\n", "\\label{_auto9} \\tag{23}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{0,1} = 0,\n", "\\label{_auto10} \\tag{24}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{N_x,N_x-1} = 0,\n", "\\label{_auto11} \\tag{25}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{N_x,N_x} = 1\\thinspace .\n", "\\label{_auto12} \\tag{26}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The right-hand side $b$ is written as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "b = \\left(\\begin{array}{c}\n", "b_0\\\\ \n", "b_1\\\\ \n", "\\vdots\\\\ \n", "b_i\\\\ \n", "\\vdots\\\\ \n", "b_{N_x}\n", "\\end{array}\\right)\n", "\\label{_auto13} \\tag{27}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "b_0 = 0,\n", "\\label{_auto14} \\tag{28}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "b_i = u^{n-1}_i,\\quad i=1,\\ldots,N_x-1,\n", "\\label{_auto15} \\tag{29}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "b_{N_x} = 0 \\thinspace . \\label{_auto16} \\tag{30}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We observe that the matrix $A$ contains quantities that do not change\n", "in time. Therefore, $A$ can be formed once and for all before we enter\n", "the recursive formulas for the time evolution.\n", "The right-hand side $b$, however, must be updated at each time step.\n", "This leads to the following computational algorithm, here sketched\n", "with Python code:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(0, L, Nx+1) # mesh points in space\n", "dx = x[1] - x[0]\n", "t = np.linspace(0, T, N+1) # mesh points in time\n", "u = np.zeros(Nx+1) # unknown u at new time level\n", "u_n = np.zeros(Nx+1) # u at the previous time level\n", "\n", "# Data structures for the linear system\n", "A = np.zeros((Nx+1, Nx+1))\n", "b = np.zeros(Nx+1)\n", "\n", "for i in range(1, Nx):\n", " A[i,i-1] = -F\n", " A[i,i+1] = -F\n", " A[i,i] = 1 + 2*F\n", "A[0,0] = A[Nx,Nx] = 1\n", "\n", "# Set initial condition u(x,0) = I(x)\n", "for i in range(0, Nx+1):\n", " u_n[i] = I(x[i])\n", "\n", "import scipy.linalg\n", "\n", "for n in range(0, Nt):\n", " # Compute b and solve linear system\n", " for i in range(1, Nx):\n", " b[i] = -u_n[i]\n", " b[0] = b[Nx] = 0\n", " u[:] = scipy.linalg.solve(A, b)\n", "\n", " # Update u_n before next step\n", " u_n[:] = u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Regarding verification, the same considerations apply as for the\n", "Forward Euler method (the section [Verification](#diffu:pde1:FE:verify)).\n", "\n", "\n", "\n", "## Sparse matrix implementation\n", "
\n", "\n", "We have seen from ([19](#diffu:pde1:matrix:sparsity)) that the matrix\n", "$A$ is tridiagonal. The code segment above used a full, dense matrix\n", "representation of $A$, which stores a lot of values we know are zero\n", "beforehand, and worse, the solution algorithm computes with all these\n", "zeros. With $N_x+1$ unknowns, the work by the solution algorithm is\n", "$\\frac{1}{3} (N_x+1)^3$ and the storage requirements $(N_x+1)^2$. By\n", "utilizing the fact that $A$ is tridiagonal and employing corresponding\n", "software tools that work with the three diagonals, the work and\n", "storage demands can be proportional to $N_x$ only. This leads to a\n", "dramatic improvement: with $N_x=200$, which is a realistic resolution,\n", "the code runs about 40,000 times faster and reduces the storage to\n", "just 1.5%! mathcal{I}_t is no doubt that we should take advantage of the fact\n", "that $A$ is tridiagonal.\n", "\n", "The key idea is to apply a data structure for a tridiagonal or sparse\n", "matrix. The `scipy.sparse` package has relevant utilities. For\n", "example, we can store only the nonzero diagonals of a matrix. The\n", "package also has linear system solvers that operate on sparse matrix\n", "data structures. The code below illustrates how we can store only the\n", "main diagonal and the upper and lower diagonals." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# Representation of sparse matrix and right-hand side\n", "main = np.zeros(Nx+1)\n", "lower = np.zeros(Nx)\n", "upper = np.zeros(Nx)\n", "b = np.zeros(Nx+1)\n", "\n", "# Precompute sparse matrix\n", "main[:] = 1 + 2*F\n", "lower[:] = -F\n", "upper[:] = -F\n", "# Insert boundary conditions\n", "main[0] = 1\n", "main[Nx] = 1\n", "\n", "A = scipy.sparse.diags(\n", " diagonals=[main, lower, upper],\n", " offsets=[0, -1, 1], shape=(Nx+1, Nx+1),\n", " format='csr')\n", "print A.todense() # Check that A is correct\n", "\n", "# Set initial condition\n", "for i in range(0,Nx+1):\n", " u_n[i] = I(x[i])\n", "\n", "for n in range(0, Nt):\n", " b = u_n\n", " b[0] = b[-1] = 0.0 # boundary conditions\n", " u[:] = scipy.sparse.linalg.spsolve(A, b)\n", " u_n[:] = u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `scipy.sparse.linalg.spsolve` function utilizes the sparse storage\n", "structure of `A` and performs, in this case, a very efficient Gaussian\n", "elimination solve.\n", "\n", "The program [`diffu1D_u0.py`](${src_diffu}/diffu1D_u0.py)\n", "contains a function `solver_BE`, which implements the Backward Euler scheme\n", "sketched above.\n", "As mentioned in the section [Forward Euler scheme](#diffu:pde1:FE),\n", "the functions `plug` and `gaussian`\n", "run the case with $I(x)$ as a discontinuous plug or a smooth\n", "Gaussian function. All experiments point to two characteristic\n", "features of the Backward Euler scheme: 1) it is always stable, and\n", "2) it always gives a smooth, decaying solution.\n", "\n", "## Crank-Nicolson scheme\n", "
\n", "\n", "The idea in the Crank-Nicolson scheme is to apply centered\n", "differences in space and time, combined with an average in time.\n", "We demand the PDE to be fulfilled at the spatial mesh points, but\n", "midway between the points in the time mesh:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial}{\\partial t} u(x_i, t_{n+\\frac{1}{2}}) =\n", "\\dfc\\frac{\\partial^2}{\\partial x^2}u(x_i, t_{n+\\frac{1}{2}}) + f(x_i,t_{n+\\frac{1}{2}}),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "for $i=1,\\ldots,N_x-1$ and $n=0,\\ldots, N_t-1$.\n", "\n", "With centered differences in space and time, we get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_t u = \\dfc D_xD_x u + f]^{n+\\frac{1}{2}}_i\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On the right-hand side we get an expression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{1}{\\Delta x^2}\\left(u^{n+\\frac{1}{2}}_{i-1} - 2u^{n+\\frac{1}{2}}_i + u^{n+\\frac{1}{2}}_{i+1}\\right) + f_i^{n+\\frac{1}{2}}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This expression is problematic since $u^{n+\\frac{1}{2}}_i$ is not one of\n", "the unknowns we compute. A possibility is to replace $u^{n+\\frac{1}{2}}_i$\n", "by an arithmetic average:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{n+\\frac{1}{2}}_i\\approx\n", "\\frac{1}{2}\\left(u^{n}_i +u^{n+1}_{i}\\right)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the compact notation, we can use the arithmetic average\n", "notation $\\overline{u}^t$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_t u = \\dfc D_xD_x \\overline{u}^t + f]^{n+\\frac{1}{2}}_i\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also use an average for $f_i^{n+\\frac{1}{2}}$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_t u = \\dfc D_xD_x \\overline{u}^t + \\overline{f}^t]^{n+\\frac{1}{2}}_i\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After writing out the differences and average, multiplying by $\\Delta t$,\n", "and collecting all unknown terms on the left-hand side, we get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{n+1}_i - \\frac{1}{2} F(u^{n+1}_{i-1} - 2u^{n+1}_i + u^{n+1}_{i+1})\n", "= u^{n}_i + \\frac{1}{2} F(u^{n}_{i-1} - 2u^{n}_i + u^{n}_{i+1})\\nonumber\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\qquad + \\frac{1}{2} f_i^{n+1} + \\frac{1}{2} f_i^n\\thinspace .\n", "\\label{_auto17} \\tag{31}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Also here, as in the Backward Euler scheme, the new unknowns\n", "$u^{n+1}_{i-1}$, $u^{n+1}_{i}$, and $u^{n+1}_{i+1}$ are coupled\n", "in a linear system $AU=b$, where $A$ has the same structure\n", "as in ([19](#diffu:pde1:matrix:sparsity)), but with slightly\n", "different entries:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "A_{i,i-1} = -\\frac{1}{2} F\n", "\\label{_auto18} \\tag{32}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{i,i} = 1 + F\n", "\\label{_auto19} \\tag{33}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{i,i+1} = -\\frac{1}{2} F\n", "\\label{_auto20} \\tag{34}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "in the equations for internal points, $i=1,\\ldots,N_x-1$. The equations\n", "for the boundary points correspond to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "A_{0,0} = 1,\n", "\\label{_auto21} \\tag{35}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{0,1} = 0,\n", "\\label{_auto22} \\tag{36}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{N_x,N_x-1} = 0,\n", "\\label{_auto23} \\tag{37}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "A_{N_x,N_x} = 1\\thinspace .\n", "\\label{_auto24} \\tag{38}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The right-hand side $b$ has entries" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "b_0 = 0,\n", "\\label{_auto25} \\tag{39}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "b_i = u^{n-1}_i + \\frac{1}{2}(f_i^n + f_i^{n+1}),\\quad i=1,\\ldots,N_x-1,\n", "\\label{_auto26} \\tag{40}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "b_{N_x} = 0 \\thinspace . \\label{_auto27} \\tag{41}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When verifying some implementation of the Crank-Nicolson scheme by convergence rate testing,\n", "one should note that the scheme is second order accurate in both space and time. The numerical\n", "error then reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E = C_t\\Delta t^r + C_x\\Delta x^r,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $r=2$ ($C_t$ and $C_x$ are unknown constants, as before).\n", "When introducing a single discretization parameter, we may now simply choose" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "h = \\Delta x = \\Delta t,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which gives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E = C_th^r + C_xh^r = (C_t + C_x)h^r,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $r$ should approach 2 as resolution is increased in the convergence rate computations.\n", "\n", "\n", "\n", "## The unifying $\\theta$ rule\n", "
\n", "\n", "For the equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial t} = G(u),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $G(u)$ is some\n", "spatial differential operator, the $\\theta$-rule\n", "looks like" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{u^{n+1}_i - u^n_i}{\\Delta t} =\n", "\\theta G(u^{n+1}_i) + (1-\\theta) G(u^{n}_i)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The important feature of this time discretization scheme is that\n", "we can implement one formula and then generate a family of\n", "well-known and widely used schemes:\n", "\n", " * $\\theta=0$ gives the Forward Euler scheme in time\n", "\n", " * $\\theta=1$ gives the Backward Euler scheme in time\n", "\n", " * $\\theta=\\frac{1}{2}$ gives the Crank-Nicolson scheme in time\n", "\n", "In the compact difference notation, we write the $\\theta$ rule\n", "as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[\\overline{D}_t u = \\dfc D_xD_x u]^{n+\\theta}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have that $t_{n+\\theta} = \\theta t_{n+1} + (1-\\theta)t_n$.\n", "\n", "Applied to the 1D diffusion problem, the $\\theta$-rule gives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\frac{u^{n+1}_i-u^n_i}{\\Delta t} &=\n", "\\dfc\\left( \\theta \\frac{u^{n+1}_{i+1} - 2u^{n+1}_i + u^{n+1}_{i-1}}{\\Delta x^2}\n", "+ (1-\\theta) \\frac{u^{n}_{i+1} - 2u^n_i + u^n_{i-1}}{\\Delta x^2}\\right)\\\\ \n", "&\\qquad + \\theta f_i^{n+1} + (1-\\theta)f_i^n\n", "\\thinspace .\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This scheme also leads to a matrix system with entries" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A_{i,i-1} = -F\\theta,\\quad A_{i,i} = 1+2F\\theta\\quad,\n", "A_{i,i+1} = -F\\theta,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "while right-hand side entry $b_i$ is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "b_i = u^n_{i} + F(1-\\theta)\n", "\\frac{u^{n}_{i+1} - 2u^n_i + u^n_{i-1}}{\\Delta x^2} +\n", "\\Delta t\\theta f_i^{n+1} + \\Delta t(1-\\theta)f_i^n\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The corresponding entries for the boundary points are as in the Backward\n", "Euler and Crank-Nicolson schemes listed earlier.\n", "\n", "Note that convergence rate testing with implementations of the theta rule must\n", "adjust the error expression according to which of the underlying schemes is actually being run.\n", "That is, if $\\theta=0$ (i.e., Forward Euler) or $\\theta=1$ (i.e., Backward Euler), there should\n", "be first order convergence, whereas with $\\theta=0.5$ (i.e., Crank-Nicolson), one should get\n", "second order convergence (as outlined in previous sections).\n", "\n", "\n", "\n", "## Experiments\n", "
\n", "\n", "\n", "We can repeat the experiments from the section [Numerical experiments](#diffu:pde1:FE:experiments)\n", "to see if the Backward Euler or Crank-Nicolson schemes have problems\n", "with sawtooth-like noise when starting with a discontinuous initial\n", "condition. We can also verify that we can have $F>\\frac{1}{2}$,\n", "which allows larger time steps than in the Forward Euler method.\n", "\n", "\n", "\n", "
\n", "\n", "

Backward Euler scheme for $F=0.5$.

\n", "\n", "\n", "\n", "\n", "\n", "The Backward Euler scheme always produces smooth solutions for any $F$.\n", "[Figure](#diffu:pde1:BE:fig:F=0.5) shows one example.\n", "Note that the mathematical discontinuity at $t=0$ leads to a linear\n", "variation on a mesh, but the approximation to a jump becomes better\n", "as $N_x$ increases. In our simulation, we specify $\\Delta t$ and $F$,\n", "and set $N_x$ to $L/\\sqrt{\\dfc\\Delta t/F}$. Since $N_x\\sim\\sqrt{F}$,\n", "the discontinuity looks sharper in the Crank-Nicolson\n", "simulations with larger $F$.\n", "\n", "The Crank-Nicolson method produces smooth solutions for small $F$,\n", "$F\\leq\\frac{1}{2}$, but small noise gets more and more evident as $F$\n", "increases. Figures [diffu:pde1:CN:fig:F=3](#diffu:pde1:CN:fig:F=3) and [diffu:pde1:CN:fig:F=10](#diffu:pde1:CN:fig:F=10)\n", "demonstrate the effect for $F=3$ and $F=10$, respectively.\n", "The section [diffu:pde1:analysis](#diffu:pde1:analysis) explains why such noise occur.\n", "\n", "\n", "\n", "
\n", "\n", "

Crank-Nicolson scheme for $F=3$.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "

Crank-Nicolson scheme for $F=10$.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## The Laplace and Poisson equation\n", "\n", "\n", "The Laplace equation, $\\nabla^2 u = 0$, and the Poisson equation,\n", "$-\\nabla^2 u = f$, occur in numerous applications throughout science and\n", "engineering. In 1D these equations read\n", "$u''(x)=0$ and $-u''(x)=f(x)$, respectively.\n", "We can solve 1D variants of the Laplace equations with the listed\n", "software, because we can interpret $u_{xx}=0$ as the limiting solution\n", "of $u_t = \\dfc u_{xx}$ when $u$ reaches a steady state limit where\n", "$u_t\\rightarrow 0$.\n", "Similarly, Poisson's equation $-u_{xx}=f$ arises from solving\n", "$u_t = u_{xx} + f$ and letting $t\\rightarrow\\infty$ so $u_t\\rightarrow 0$.\n", "\n", "Technically in a program, we can simulate $t\\rightarrow\\infty$\n", "by just taking one large time step:\n", "$\\Delta t\\rightarrow\\infty$. In the limit, the Backward Euler\n", "scheme gives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "-\\frac{u^{n+1}_{i+1} - 2u^{n+1}_i + u^{n+1}_{i-1}}{\\Delta x^2} = f^{n+1}_i,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is nothing but the discretization $[-D_xD_x u = f]^{n+1}_i=0$ of\n", "$-u_{xx}=f$.\n", "\n", "The result above means that\n", "the Backward Euler scheme can solve the limit equation directly and\n", "hence produce a solution of the 1D Laplace equation.\n", "With the Forward Euler scheme we must do the time stepping since $\\Delta t >\n", "\\Delta x^2/\\dfc$\n", "is illegal and leads to instability.\n", "We may interpret this time stepping\n", "as solving the equation system from $-u_{xx}=f$ by iterating on a\n", "pseudo time variable.\n", "\n", "[hpl 2: Better to say the last sentence when we treat iterative methods.]" ] } ], "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.8.2" } }, "nbformat": 4, "nbformat_minor": 4 }