\n", "\\begin{align}\n", "\\left.{df(x)\\over dx} \\right| & \\approx \\frac{f(x)-f(x-\\Delta x)}{\\Delta x} \\\\\n", "\\left.{df(x)\\over dx} \\right|_{x=2} & \\approx \\frac{f(2)-f(1)}{2-1} = (5.12 - 8.38)/1.0 = -3.26\\\\\n", "\\end{align}\n", "
\n", "\n", "So, we approximate the slope at $x=2$ by the slope of the line that connects $f(1)$ and $f(2)$:\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[
\n", "\\begin{align}\n", "\\left.{df(x)\\over dx} \\right| & \\approx \\frac{f(x+\\Delta x)-f(x)}{\\Delta x} \\\\\n", "\\left.{df(x)\\over dx} \\right|_{x=2} & \\approx \\frac{f(3)-f(2)}{3-2} = (4.88 - 5.12)/1.0 = -0.24\\\\\n", "\\end{align}\n", "
\n", "\n", "Shown graphically:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[
\n", "\\begin{align}\n", "\\left.{df(x)\\over dx} \\right| & \\approx \\frac{f(x+\\Delta x)-f(x-\\Delta x)}{2\\Delta x} \\\\\n", "\\left.{df(x)\\over dx} \\right|_{x=2} & \\approx \\frac{f(3)-f(1)}{3-1} = {(4.88 - 8.38)\\over 2.0} = -1.75\\\\\n", "\\end{align}\n", "
\n", "\n", "Shown graphically:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[
\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "lines_to_next_cell": 2 }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact derivative at x = 2.00, df/dx = -0.7500\n", " FD approximation |error|\n", " forward backward central forward backward central \n", " -0.250 -3.250 -1.750 0.500 2.500 1.000\n", " -0.250 -1.750 -1.000 0.500 1.000 0.250\n", " -0.438 -1.188 -0.812 0.312 0.438 0.062\n", " -0.578 -0.953 -0.766 0.172 0.203 0.016\n", " -0.660 -0.848 -0.754 0.090 0.098 0.004\n" ] }, { "data": { "text/plain": [ "
\n", "In the application of finite differences to solve differential equations, the function that we are solving for is not known a priori, so the error in the approximation cannot be computed. For that reason, when we assess a finite-difference approximation, we look at **how quickly the error shrinks as we shrink** $\\Delta x$. To see that we have to revisit Taylor series.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Taylor Series analysis of finite-difference approximations\n", "\n", "Taylor series reference: Khan Academy https://www.khanacademy.org/math/ap-calculus-bc/bc-series-new/bc-10-11/v/maclaurin-and-taylor-series-intuition\n", "UBC Peter Wall https://www.math.ubc.ca/~pwalls/math-python/differentiation/\n", "
Video - Taylor series and finite differences https://www.coursera.org/lecture/computers-waves-simulations/w2v3-taylor-series-PEN8L\n", "
\n", "#### Approach\n", "* manipulate Taylor series approximations to develop exact expressions for derivatives.\n", "* compare these expressions to finite-difference approximations and thereby develop expressions for the **truncation error** = exact - approximate.\n", "\n", "
\n", "Recall a Taylor series gives you the value of a function at a location $x_i + dx$ using the value of the function and its derivatives at the point $x_i$. Functions that can be approximated in this way are called analytic https://en.wikipedia.org/wiki/Analytic_function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "Let's represent a function $f$ at points $x_i+\\Delta x$ and $x_i-\\Delta x$:\n", "\\begin{align}\n", "f(x_i+\\Delta x) &= f(x_i) + (\\Delta x){df(x_i)\\over dx} + {(\\Delta x)^2\\over 2!}{d^2f(x_i)\\over dx^2} + {(\\Delta x)^3\\over 3!}{d^3f(x_i)\\over dx^3} \\nonumber \\\\\n", "&+ \\cdots + {(\\Delta x)^n\\over n!}{d^nf(x_i)\\over dx^n} + \\cdots \\label{8fd884} \\\\\n", "f(x_i-\\Delta x) &= f(x_i) + (-\\Delta x){df(x_i)\\over dx} + {(-\\Delta x)^2\\over 2!}{d^2f(x_i)\\over dx^2} + {(-\\Delta x)^3\\over 3!}{d^3f(x_i)\\over dx^3} + \\nonumber\\\\\n", "&\\cdots + {(-\\Delta x)^n\\over n!}{d^nf(x_i)\\over dx^n} + \\cdots \\label{8fd885} \\\\\n", "\\end{align}\n", "
\n", "#### Forward difference Taylor analysis\n", "
\n", "Subtract $f(x_i)$ from the Taylor approximation $f(x_i+\\Delta x)$ \\ref{8fd884}\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{align}\n", "f(x_i+\\Delta x) - f(x_i) &= (\\Delta x){df(x_i)\\over dx} + {(\\Delta x)^2\\over 2!}{d^2f(x_i)\\over dx^2} + {(\\Delta x)^3\\over 3!}{d^3f(x_i)\\over dx^3} + \\cdots + {(\\Delta x)^n\\over n!}{d^nf(x_i)\\over dx^n} + \\cdots \\label{8fd886} \\\\\n", "\\end{align}\n", "
\n", "Solve for ${df(x_i)\\over dx}$, the derivative at point $x_i$:\n", "
\n", "\\begin{align}\n", "{df(x_i)\\over dx} & = \\overbrace{{f(x_i+\\Delta x) - f(x_i)\\over \\Delta x}}^\\text{forward difference approximation} \\nonumber\\\\\n", "&- \\overbrace{\\frac{1}{\\Delta x}\\left[ {(\\Delta x)^2\\over 2!}{d^2f(x_i)\\over dx^2} +\n", "{(\\Delta x)^3\\over 3!}{d^3f(x_i)\\over dx^3} + \\cdots + {(\\Delta x)^n\\over n!}{d^nf(x_i)\\over dx^n} + \\cdots \\right]}^\\text{truncation error} \\label{8fd887} \\\\\n", "\\end{align}\n", "
\n", "Let's pause for a moment. The expression \\ref{8fd887} is exact. In words this expression says that the exact derivative at point $x_i$ is given by the forward-difference approximation minus a **truncation error**, all the terms in the $[]$ braces. That is, true = approximate - error.\n", "
\n", "##### Truncation error\n", "The forward difference truncation error is\n", "
\n", "\\begin{align}\n", "\\left[ \\underbrace{\\frac{(\\Delta x)}{2!}}_\\text{leading term}{d^2f(x_i)\\over dx^2} + \\underbrace{(\\Delta x)^2\\over 3!}_\\text{2nd term}{d^3f(x_i)\\over dx^3} + \\cdots + {(\\Delta x)^{n-1}\\over n!}{d^nf(x_i)\\over dx^n} + \\cdots \\right] \\label{8fd888} \\\\\n", "\\end{align}\n", "
\n", "1. What controls the size of the truncation error?\n", "* How can we characterize this truncation error?\n", "
\n", "##### Order of approximation\n", "**Problem**: when we used finite-difference approximations to compute Euler timesteps or fluxes, we did not know the function, and therefore in general we don't know the values of the derivatives ${df(x_i)\\over dx}$, ${d^2f(x_i)\\over dx^2}$, etc, in the truncation error. So instead we focus on the effects of changing $\\Delta x$ on the truncation error.\n", "
\n", "\n", "\\begin{align*}\n", "\\begin{matrix}\n", " & & \\text{leading term} & &\\text{2nd term} \\\\\n", "\\text{Reduce }\\quad \\Delta x \\quad \\text{by}\\quad 1/2 & \\text{reduces by}\\quad &1/2& \\text{reduces by}\\quad &1/4 \\\\\n", "\\text{Reduce }\\quad \\Delta x \\quad \\text{by}\\quad 1/4 & \\text{reduces by}\\quad &1/4& \\text{reduces by}\\quad &1/16 \\\\\n", "\\text{Reduce }\\quad \\Delta x \\quad \\text{by}\\quad 1/8 & \\text{reduces by}\\quad &1/8& \\text{reduces by}\\quad &1/64 \\\\\n", "\\end{matrix}\n", "\\end{align*}\n", "\n", "
\n", "When $\\Delta x$ is reduced, the leading term $(\\Delta x){df(x_i)\\over dx}$ decreases the most slowly. When $\\Delta x$ is reduced by a factor of 8, the leading term is reduced by a factor of 8, but the second term is reduced by a factor of 64. So, the error is dominated by the leading order behavior. Higher-order terms are eventually insignificant.\n", "
\n", "So, we now write\n", "
\n", "\\begin{align}\n", "{df(x_i)\\over dx} & = \\overbrace{{f(x_i+\\Delta x) - f(x_i)\\over \\Delta x}}^\\text{forward difference approximation} \\nonumber\\\\\n", "&- \\overbrace{\\left[ {(\\Delta x)\\over 2!}{d^2f(x_i)\\over dx^2} + {(\\Delta x)^2\\over 3!}{d^3f(x_i)\\over dx^3} + \\cdots + {(\\Delta x)^{n-1}\\over n!}{d^nf(x_i)\\over dx^n} + \\cdots \\right]}^\\text{truncation error} \\label{8fd889} \\nonumber \\\\\n", "& = \\overbrace{{f(x_i+\\Delta x) - f(x_i)\\over \\Delta x}}^\\text{forward difference approximation} + \\mathcal{O}(\\Delta x) \\label{8fd8810} \\\\\n", "\\end{align}\n", "
\n", "For this reason, we say that the forward difference approximation is **First order** because as $\\Delta x \\rightarrow 0$, the truncation error shrinks by $(\\Delta x)^1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Backward difference Taylor analysis\n", "
\n", "Subtract expansion for $f(x_i-\\Delta x)$, \\ref{8fd885} from $f(x_i)$.\n", "
\n", "\\begin{align}\n", "{df(x_i)\\over dx} & = \\overbrace{{f(x_i) - f(x_i-\\Delta x)\\over \\Delta x}}^\\text{backward difference approximation} + \\mathcal{O}(\\Delta x) \\label{8fd8811} \\\\\n", "\\end{align}\n", "
\n", "The backward difference is also a **First order** approximation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Central difference Taylor analysis\n", "
\n", "Subtract expansion for $f(x_i+ \\Delta x)$, \\ref{8fd884} from $f(x_i-\\Delta x)$, \\ref{8fd885}. Some terms cancel out, which leads to:\n", "
\n", "\\begin{align}\n", "{df(x_i)\\over dx} & = \\overbrace{{f(x_i + \\Delta x) - f(x_i-\\Delta x)\\over 2\\Delta x}}^\\text{central difference approximation} + \\mathcal{O}(\\Delta x)^2 \\label{8fd88112} \\\\\n", "\\end{align}\n", "
\n", "The central difference is a **Second order** approximation because as $\\Delta x \\rightarrow 0$, the truncation error shrinks by $(\\Delta x)^2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Second derivative\n", "
\n", "Add the Taylor expansions for $f(x_i+\\Delta x)$ \\ref{8fd884} and $f(x_i-\\Delta x)$ \\ref{8fd885}\n", "
\n", "\\begin{align}\n", "f(x_i + \\Delta x) + f(x_i - \\Delta x)&= 2 f(x_i) + 2{(\\Delta x)^2\\over 2!}{d^2f(x_i)\\over dx^2} + 2\\left[{(\\Delta x)^4\\over 4!}{d^4f(x_i)\\over dx^4} + \\cdots + \\right] \\label{8fd8812} \\\\\n", "\\end{align}\n", "
\n", "Solve for ${d^2f(x_i)\\over dx^2}$:\n", "
\n", "\\begin{align}\n", "{d^2f(x_i)\\over dx^2}&= \\overbrace{{f(x_i + \\Delta x) - 2 f(x_i) + f(x_i - \\Delta x)\\over (\\Delta x)^2}}^\\text{central second derivative} + \\mathcal{O}(\\Delta x)^2 \\label{8fd8813} \\\\\n", "\\end{align}\n", "
\n", "The central second derivative is also **second order**. Does that last expression remind you of something?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Summary of finite difference approximations\n", "We derived three approximations for the first derivative: forward, backward and central. Let's express them in terms of values of, say, concentration on a grid with constant grid spacing $\\Delta x$ and\n", "\\begin{align*}\n", "c_i = c(&x_i) \\\\\n", "c_{i+1} = c(x_i &+\\Delta x) \\\\\n", "c_{i-1} = c(x_i &- \\Delta x) \\\\\n", "\\end{align*}\n", "\n", "#### First derivative\n", "
\n", "\\begin{align}\n", "\\text{Forward} \\qquad {dc(x_i)\\over dx} =& {c_{i+1} - c_{i} \\over \\Delta x} + \\mathcal{O}(\\Delta x) \\label{8fd8814} \\\\\n", "\\text{Backward}\\qquad {dc(x_i)\\over dx}= & {c_{i} - c_{i-1} \\over \\Delta x} + \\mathcal{O}(\\Delta x) \\label{8fd8815} \\\\\n", "\\text{Central} \\qquad{dc(x_i)\\over dx} =& {c_{i+1} - c_{i-1} \\over 2\\Delta x} + \\mathcal{O}(\\Delta x)^2 \\label{8fd8816}\\\\\n", "\\end{align}\n", "
\n", "#### Second derivative\n", "
\n", "\\begin{align}\n", "\\text{Central} \\qquad{d^2c(x_i)\\over dx^2}= & {c_{i+1} -2 c_i + c_{i-1} \\over (\\Delta x)^2} + \\mathcal{O}(\\Delta x)^2 \\label{8fd8817}\\\\\n", "\\end{align}\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implications for finite volumes\n", "\n", "#### Euler methods\n", "
\n", "In the forward and backward Euler (or explicit and implicit) time stepping the time derivative is approximated as:\n", "
\n", "\\begin{align}\n", "{dc\\over dt} \\approx {c(t+\\Delta t) - c(t)\\over \\Delta t} \\label{8fd8818}\\\\\n", "\\end{align}\n", "
\n", "What is the order of this approximation?\n", "
\n", "#### Fluxes\n", "
\n", "In our stencil, we approximate fluxes such as $j_{EC}$ as:\n", "
\n", "\\begin{align}\n", "j_{EC}&= D\\theta {c_E - c_C \\over \\Delta x} \\label{8fd8819}\\\\\n", "\\end{align}\n", "
\n", "Although this looks like a first order approximation, it's really a second order. Why? It's because\n", "the flux is really being approximated at the interface between the two gridblocks.\n", "
\n", "
![\"pic05\"](\"figures/flux_fd.png\")
\n",
"\n", "To emphasize this, let's use this finite-difference notation, $c_{i-1}=c_C$, $c_{i+1} = c_E$ and write:\n", "
\n", "\\begin{align}\n", "j_{EC}&= D\\theta {c_{i+1} - c_{i-1} \\over \\Delta x} \\label{8fd8820}\\\\\n", "\\end{align}\n", "
\n", "The idea is that we are approximating the flux at a point $x_i$, between $x_{i+1}$ and $x_{i-1}$. So, the term looks like a central difference, where the grid points are separated by $\\Delta x /2$. Accordingly, the order is order $\\mathcal{O}((\\Delta x)^2)$\n", "
\n", "#### Why not finite differences?\n", "Finite-difference methods are straightforward to apply to simple equations (see the assignment below) with constant coefficients. However, for cases where the material properties vary with location (heterogeneous), they can be difficult to apply, for example, when the material properties are heterogeneous (diffusion coefficient, hydraulic conductivity, thermal conductivity, etc). A naive application can lead to **non-conservative** schemes - that is, ones that do not produce discrete approximations that are conservative. In contrast, the finite-volume method is a conservative discrete approximation.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Assignment\n", "## Finite difference stencil for diffusion\n", "### 1-D Steady-state diffusion with constant coefficients\n", "#### Partial differential equation for diffusion with homogeneous diffusion coefficient\n", "To apply finite-difference approximations to solve differential equations, one simply replaces the derivatives that appear in the equations with a finite difference approximation. \n", "
\n", "The partial differential equation for diffusion in three dimensions (not in porous media) is (see `9_pdes_1.ipynb`):\n", "
\n", "\\begin{align}\n", "&{\\partial c\\over \\partial t} = {\\partial \\over \\partial x}\\left( D {\\partial c \\over \\partial x}\\right) + {\\partial \\over \\partial y}\\left( D {\\partial c \\over \\partial y}\\right)+ {\\partial \\over \\partial z}\\left( D {\\partial c \\over \\partial z}\\right) \\label{8fd8821}\\\\\n", "\\text{or}&\\nonumber\\\\\n", "&{\\partial c\\over \\partial t} = \\vec{\\nabla} \\cdot (D \\vec{\\nabla} c)\\label{8fd8822}\\\\\n", "\\end{align}\n", "
\n", "Write the partial differential equation for 1-D (in the $x$ direction) steady-state diffusion where the diffusion coefficient is spatially homogeneous.\n", "\n", "**Put your equation here**\n", "
\n", "\n", "#### Finite-difference stencil for steady-state diffusion with homogeneous diffusion coefficient\n", "
\n", "To construct a finite-difference approximation to a differential equation, we simply replace the derivatives in the equation with a finite difference approximation. For example, for an ordinary differential equation of the form\n", "
\n", "\\begin{align}\n", "{dy\\over dx} &= 2 y + 6 \\label{8fd8824}\\\\\n", "\\end{align}\n", "
\n", "we develop a finite - difference stencil by replacing the derivative by a finite difference approximation. For example a forward difference approximation would be:\n", "
\n", "\\begin{align}\n", "{y(x_{i+1}) - y(x_i)\\over \\Delta x} &= 2 y(x_i) + 6 \\label{8fd8825}\\\\\n", "\\end{align}\n", "
\n", "This is usual written using this notation: $y(x_i)\\rightarrow y_i$:\n", "
\n", "\\begin{align}\n", "{y_{i+1} - y_i\\over \\Delta x} &= 2 y_i + 6 \\label{8fd8826}\\\\\n", "\\end{align}\n", "
\n", "Write the finite-difference approximation for 1-D (in the $x$ direction) steady-state diffusion where the diffusion coefficient is spatially homogeneous (that is, for the equation you developed above). Use a centered approximation.\n", "
\n", "**Write your equation here **\n", "
\n", "\n", "#### Finite-volume stencil for steady-state diffusion with homogeneous diffusion coefficient\n", "The finite-volume approximation for one-dimensional steady diffusion (not in porous media) is:\n", "
\n", "\\begin{align}\n", "&\\left(D {c_E - c_C \\over \\Delta x} + D {c_W - c_C \\over \\Delta x} \\right) (\\Delta y) (\\Delta z) =0 \\label{8fd8823}\\\\\n", "\\end{align}\n", "
\n", "Write this equation for the case of homogeneous diffusion coefficient and constant gridblock size in the form $a_W c_W + a_C c_C + a_E c_E = rhs$ (that is, determine $a_W$, $a_C$, $a_E$ and $rhs$).\n", "
\n", "**Write your equation here.**\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Approximation error\n", "Compute the forward, backward and central difference error of `cos(x)` at $x = 0.7$, starting at $\\Delta x = 0.5$, and the halving $\\Delta x$ 10 times (see code snippet below). Write a python code to:\n", "1. compute and print the error for each value of $\\Delta x$ for the forward,centered and backwards approximations.\n", "* plot the error of each expression versus $\\Delta x$ \n", "* What should the slope of this plot be for the different approximations?\n", "\n", "
\n", "You can reuse most of the code above for this assignment." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# code fragment to compute and print the values of delta x to be used in the error computation\n", "dx_vals = 2.**-np.arange(1,11)\n", "for dx in np.nditer(dx_vals):\n", " print(f\"{dx:2.9f}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# add your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Reflection\n", "\n", "* What dictates your choice in the order of approximation used in a finite - difference approximation?\n", "* How does the finite-difference approximation differ from the finite-volume approximation for the case you considered in the assignment above?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "cell_metadata_filter": "all", "formats": "ipynb,py:percent", "notebook_metadata_filter": "all" }, "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.7.3" }, "nbsphinx": { "execute": "never" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }